C C *********************************************************** C ******************** PSI/88 - PART 2 ******************** C *********************************************************** C C Version 1.0 Any questions to the author should specify C the version being used. C PROGRAM PSICON C C THIS ROUTINE "HACKED OUT OF" PSI77, PART 1, WRITTEN BY: C C WILLIAM L. JORGENSEN C Yale University, DEPARTMENT OF CHEMISTRY C New Haven, CT 06511, USA C PHONE 203-432-6278 C C HACKING OUT AND MODIFICATION TO ALLOW HIDDEN LINE ELIMINATION C OF INDENTATIONS AND DOUGHNUTS. The modifications are rather C minor, so the code is mostly unchanged from PSI/77. C DAN SEVERANCE, PURDUE C C Redistribution and use in source and binary forms are permitted C provided that the above paragraphs and this one are duplicated in C all such forms and that any documentation, advertising materials, C and other materials related to such distribution and use acknowledge C that the software was developed by William Jorgensen at Purdue University C The name of the University or William Jorgensen may not be used to endorse C or promote products derived from this software without specific prior C written permission. The author is now at Yale University. C THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR C IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED C WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. C COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH, * SCALE,PERZ,KA,KO,ZO COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 C C THIS ROUTINE READS THE 3-D GRID OF ORBITAL VALUES GENERATED C BY THE PSI1 PROGRAM AND STORED IN DISK FILE 17. (UNFORMATTED) C ALL DATA NEEDED BY THE PROGRAM IS ALSO STORED IN FILE 17. C C THE SAME INPUT DECK IS USED AS WAS SPECIFIED FOR PSI1. C DIFFERENT CONTOUR LEVELS MAY BE SELECTED AND THIS PROGRAM RE-RUN C WITHOUT REXECUTION OF PSI1. C C CLOSED CONTOURS AT THE SPECIFIED LEVEL(S) C ARE DRAWN THROUGH THE DENSITY OR ORBITAL VALUE C ARRAYS AND THE POINTS(IN ANGSTROMS) FOR EACH CONTOUR ARE STORED C AS OUTPUT IN DISK FILE 22. THE TOTAL C NUMBER OF POINTS AND THE NUMBER OF CURVES ALONG C WITH THE NUMBER OF POINTS IN EACH CURVE AND AN C INDICATOR FOR WHETHER THE CONTOUR LEVEL IS POSITIVE C OR NEGATIVE ARE STORED IN THE OUTPUT DISK FILE C 23. FOR PROPER HIDDEN LINE ELIMINATION C IN THE 3-D DRAWING PROGRAM, IT IS ESSENTIAL THAT ALL CONTOURS C WHICH ARE CONSTRUCTED BE CLOSED. IF THIS CONDITION IS NOT MET, C AN ERROR MESSAGE IS PRINTED. C CARD INPUT IS DESCRIBED IN SUBROUTINE THREED. C C NOT ALL OF THE VARIABLES IN COMMON/RPLOT ARE C ESSENTIAL TO THIS PROGRAM, BUT HAVE BEEN RETAINED C FOR COMPATIBILITY WITH THE HIDDEN LINE PLOTTING C PROGRAM IN CASE THE TWO ARE MERGED C IRD = 5 ILST = 6 IDSK1 = 22 IDSK2 = 23 IDSK3 = 24 THE = 0.0E+0 GAM = 0.0E+0 PHI = 0.0E+0 NC = 0 NCURV = 0 C C MAKE 3-D CONTOUR MAP AND STORE ON DISK C CALL THREED STOP END C C SUBROUTINE THREED PARAMETER (MXPTS=51,MXPTSQ=MXPTS*MXPTS) C C THE CONTOURING PORTIONS SEPARATED FROM THE ORBITAL COMPUTATION C SECTION AND PLACED IN THIS PROGRAM. OTHER THAN NECCESSARY C CHANGES IN THE INPUT, THIS ROUTINE IS THAT WRITTEN BY C WILLIAM JORGENSEN. C DAN SEVERANCE 12-87 C C MODIFIED TO OUTPUT ADDITIONAL INFORMATION TO BE USED IN HIDDEN C LINE PORTION OF THE PROGRAM. DAN SEVERANCE 2-88 C COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,XX,YY,ZZ,IDASH,SCALE,PERZ, * KA,KO,ZZZ COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 COMMON // DENS(MXPTS,MXPTS,MXPTS),SAV(MXPTSQ),CTR(15) C C PUT DENS IN A COMMON BLOCK FOR THOSE MACHINES WHICH CAN MAP C IT TO A LARGER MEMORY AREA. THIS MAY BE CONVERTED TO A DIMENSION C STATEMENT IF PREFERRED. C REAL C(50,3),XX(850),YY(850),ZZ(850),DMIN,DMAX INTEGER IAN(49) CHARACTER CALC*20 DATA AU,PI,RT3 / 0.5291670E+0,3.14159265360E+0,1.73205080E+0 / C C DENS IS THE CHARGE DENSITY OR ORBITAL VALUE MATRIX. IAN IS THE C ARRAY OF ATOMIC NUMBERS FOR THE ATOMS. ATOMIC COORDINATES ARE C READ IN ANGSTROMS AND STORED IN ARRAY C. C CTR IS THE ARRAY OF CONTOUR LEVELS IN AU. C C CARD INPUT IS AS FOLLOWS - THIS IS NOT ALL USED IN THIS PROGRAM C BUT THE SAME INPUT DECK IS USED BY THE ORBITAL VALUE C PACKAGE C C CARD 1 - TYPE OF WAVEFUNCTION = STO-3G , 3-21G, 3-21+G, C 6-31G TO 6-31++G**. (A20) C (THIS IS FOR INFORMATIONAL USE ONLY, THE PROGRAM C DOESN'T USE IT) C C CARD 2 - NO. OF CONTOURS, ICONN, ICONZ, NORB -4I2 C ICONN = 1 FOR NEGATIVE CONTOURS, TOO. C ICONZ=1 FOR ZERO CONTOUR AND NORB=2 FOR CHARGE DENSITY PLOT. C CARD 3 - CONTOUR LEVELS (ONLY POSITIVE NEED BE SPECIFIED WHEN C ICONN AND ICONZ ARE USED - 8F10.6). C C INITIALIZE FOR AUTOMATIC PROCESSING C NCT = 1 CTR(1) = 0.0750E+0 NORB = 1 ICONN = 1 ICONZ = 0 READ (IRD,10) CALC 10 FORMAT (A) C C READ THE NO. OF CONTOURS AND THEN THE LEVELS. C READ (IRD,20) NCT,ICONN,ICONZ,NORB 20 FORMAT (4I2) READ (IRD,30) (CTR(I),I=1,NCT) 30 FORMAT (8F10.6) C C DEFAULT VALUES C IF (SCALE.EQ.0.0E+0) SCALE = 1.00E+0 IF (NCT.EQ.0) CTR(1) = 0.0750E+0 IF (NCT.EQ.0) NCT = 1 IF (NORB.EQ.0) NORB = 1 IF (NORB.EQ.1) ICONN = 1 C C READ THE ORBITAL VALUE MATRIX C WRITE (ILST,*) 'READING 17 ' OPEN (17,FILE='FOR017',FORM='UNFORMATTED',STATUS='OLD') READ (17) MAXPTS,NAT READ (17) (IAN(I),I=1,NAT) READ (17) ((C(I,J),J=1,3),I=1,NAT) READ (17) (((DENS(I,J,K),I=1,MAXPTS),J=1,MAXPTS),K=1,MAXPTS) READ (17) XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX C C DETERMINE DEFAULT RANGES, INCREASING SCALE C INCREASES THE RANGE OF THE PLOT. C SPACES = FLOAT(MAXPTS-1)*AU XMI = XMIN/AU YMI = YMIN/AU ZMI = ZMIN/AU XINC = (XMAX-XMIN)/SPACES YINC = (YMAX-YMIN)/SPACES ZINC = (ZMAX-ZMIN)/SPACES DMIN = 1000.00E+0 DMAX = -1000.00E+0 DO 40 K = 1, MAXPTS DO 40 J = 1, MAXPTS DO 40 I = 1, MAXPTS DMIN = MIN(DMIN,DENS(I,J,K)) DMAX = MAX(DMAX,DENS(I,J,K)) 40 CONTINUE WRITE (ILST,*) 'MIN, MAX DENSITY(VALUE) COMPUTED IS ',DMIN,DMAX OPEN (IDSK3,FILE='FOR024',STATUS='UNKNOWN') OPEN (IDSK2,FILE='FOR023',STATUS='UNKNOWN') OPEN (IDSK1,FILE='FOR022',STATUS='UNKNOWN') KA = 1 KO = 2 ZZZ = ZMIN ZINC = ZINC*AU DO 70 I = 1, MAXPTS KNT = 0 DO 60 J = 1, MAXPTS DO 50 K = 1, MAXPTS KNT = KNT+1 SAV(KNT) = DENS(K,J,I) 50 CONTINUE 60 CONTINUE CALL TRDPLT (SAV,MAXPTS,MAXPTS,XMIN,XMAX,YMIN,YMAX,NCT,CTR, * ICONN,ICONZ) ZZZ = ZZZ+ZINC 70 CONTINUE NC1 = NC NCURV1 = NCURV KA = 1 KO = 3 ZZZ = YMIN YINC = YINC*AU DO 100 I = 1, MAXPTS KNT = 0 DO 90 J = 1, MAXPTS DO 80 K = 1, MAXPTS KNT = KNT+1 SAV(KNT) = DENS(K,I,J) 80 CONTINUE 90 CONTINUE CALL TRDPLT (SAV,MAXPTS,MAXPTS,XMIN,XMAX,ZMIN,ZMAX,NCT,CTR, * ICONN,ICONZ) ZZZ = ZZZ+YINC 100 CONTINUE NC2 = NC NCURV2 = NCURV KA = 2 KO = 3 ZZZ = XMIN XINC = XINC*AU DO 130 I = 1, MAXPTS KNT = 0 DO 120 J = 1, MAXPTS DO 110 K = 1, MAXPTS KNT = KNT+1 SAV(KNT) = DENS(I,K,J) 110 CONTINUE 120 CONTINUE CALL TRDPLT (SAV,MAXPTS,MAXPTS,YMIN,YMAX,ZMIN,ZMAX,NCT,CTR, * ICONN,ICONZ) ZZZ = ZZZ+XINC 130 CONTINUE NC3 = NC NCURV3 = NCURV WRITE (IDSK2,140) NC,NCURV,NC1,NCURV1,NC2,NCURV2,NC3,NCURV3 NCMAX = 0 140 FORMAT (10I6) 150 FORMAT (20I4) WRITE (IDSK2,150) (NAC(I),NCDASH(I),I=1,NCURV) NCMAX = 0 DO 160 I = 1, NCURV NCMAX = MAX(NCMAX,NAC(I)) 160 CONTINUE WRITE (ILST,170) NC,NCURV,NCMAX 170 FORMAT (6X,'NC = ',I5,' NCURV = ',I4,' NCMAX =',I4/) RETURN END C C SUBROUTINE TRDPLT (ARRAY,IXDIM,IYDIM,XXMIN,XXMAX,YYMIN,YYMAX,NCTR, * CTR,ICONN,ICONZ) COMMON /SCLDAT/ XMAX,XMIN,YMAX,YMIN,XINCR,YINCR,CL,NPLT COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH, * SCALE,PERZ,KA,KO,ZZ COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 DIMENSION ARRAY(IXDIM,IYDIM),CTR(15) DATA NCSAVE / 0 / C C THIS PROGRAM PRODUCES CONTOUR PLOTS FROM THE MATRIX, C ARRAY, WHICH HAS DIMENSIONS IXDIM AND IYDIM USING THE C CONTOURS,CTR. NCTR = THE NUMBER OF CONTOURS. KO = 1,2,3 C IF THE ORDINATE IS THE X,Y OR Z AXIS% KA IS SIMILAR C FOR THE ABSCISSA. ICONN = 1 IF NEGATIVE CONTOURS ARE C DESIRED THAT ARE THE SAME AS THE POSITIVE CONTOURS (ONLY C THE POSITIVE ONES NEED BE SPECIFIED IN CTR AND NCTR). C ICONZ = 1 IF A DOTTED ZERO CONTOUR IS DESIRED. C W. L. JORGENSEN - DECEMBER 6, 1971 C XMAX = XXMAX XMIN = XXMIN YMAX = YYMAX YMIN = YYMIN XINCR = (XMAX-XMIN)/FLOAT(IXDIM-1) YINCR = (YMAX-YMIN)/FLOAT(IYDIM-1) C C CALL CONTUR FOR EACH CONTOUR. C NPLT = 0 IDASH = 0 DO 10 NCL = 1, NCTR CL = CTR(NCL) NPLT = NPLT+1 CALL CONTUR (ARRAY,IXDIM,IYDIM,CL,IERROR) IF (IERROR.EQ.1) THEN WRITE (ILST,40) STOP ENDIF IDASH = IDASH + 2*NCL 10 CONTINUE C C FOR ICONN = 1, DO AGAIN WITH - CONTOURS. C IF (ICONN.EQ.1) THEN NPLT = 0 IDASH = 1 DO 20 NCL = 1, NCTR CL = -CTR(NCL) NPLT = NPLT+1 CALL CONTUR (ARRAY,IXDIM,IYDIM,CL,IERROR) IF (IERROR.EQ.1) THEN WRITE (ILST,40) STOP ENDIF IDASH = IDASH + 2*NCL 20 CONTINUE ENDIF C C FOR ICONZ=1, PLOT 0 CONTOUR. C IF (ICONZ.EQ.1) THEN CL = 0.00E+0 NPLT = 0 CALL CONTUR (ARRAY,IXDIM,IYDIM,0.00E+0,IERROR) IF (IERROR.EQ.1) THEN WRITE (ILST,40) STOP ENDIF ENDIF IF (NCURV.NE.NCSAVE) THEN WRITE (IDSK3,30) NCURV NCSAVE = NCURV ENDIF 30 FORMAT (I6) RETURN 40 FORMAT ('0** CONTOUR PASSING OUTSIDE THE PLOTTING BOX **'/ * '0** INCREASE THE CONTOUR LEVEL AND RERUN PSICON OR **'/ * '0** INCREASE THE SCALE FACTOR FOR PSI1 AND RERUN IT **') END C C SUBROUTINE ROTPLT COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH, * SCALE,PERZ,KA,KO,ZO COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 DIMENSION SAV(3) C C THIS ROUTINE STORES THE CURVES WHICH HAVE C BEEN GENERATED BY THE CONTUR PLOTTING C ROUTINES IN THE DISK FILE NAMED CURVE C C TRANSFORM CURVE TO ORIGINAL SYSTEM C KZ = 6-KA-KO DO 10 I = 1, N SAV(KA) = X(I) SAV(KO) = Y(I) SAV(KZ) = ZO X(I) = SAV(1) Y(I) = SAV(2) Z(I) = SAV(3) 10 CONTINUE C C STORE CURVE ON DISK C NCURV = NCURV+1 NCDASH(NCURV) = IDASH NC = NC+N NAC(NCURV) = N IF (N.GT.850) THEN WRITE (ILST,*) '>>>>> WARNING:XYZ REC NEED TO BE REDIMENSIONED' WRITE (ILST,*) '>>>>>>THEY NEED ',N ENDIF WRITE (IDSK1,20) (X(I),Y(I),Z(I),I=1,N) 20 FORMAT (8F10.6) RETURN END C C SUBROUTINE CONTUR (A,NX,NY,CLP,JERROR) COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,X,Y,Z,IDASH,SCALE,PERZ,KA, * KO,ZO COMMON /CRPTPM/ NRMAX,NPMAX,LX,UX,LY,UY,CL,NR,NP,IXO,IYO,ISO,PC, * LOPC COMMON /TEMP/ REC DIMENSION INX(8),INY(8),A(NX,NY),X(850),Y(850),REC(850),Z(850) INTEGER UX,UY,UX1,UY1,COM,REC LOGICAL PC,LOPC DATA INX / -1,-1,0,3*1,0,-1 / DATA INY / 0,3*1,0,3*-1 / NF = 1 LX = NF LY = NF UX = NX UY = NY CL = CLP NRMAX = 850 NPMAX = 850 C C THE 3 ROUTINES, CONTUR,CURVE AND INTERP FORM A C GENERAL CONTOUR MAPPING PACKAGE THAT SIMPLY REQUIRES C A PLOTTING ROUTINE PLOT(N,X,Y) TO PLOT N POINTS C WITH COORDINATES IN THE ARRAYS X AND Y. C IF CONTOURS WITH MORE THAN 400 POINTS ARE TO BE C DRAWN (THIS IS VERY LARGE) , NRMAX AND NPMAX ABOVE C MUST BE ADJUSTED ALONG WITH THE DIMENSION STATEMENTS C THIS PROGRAM IS BASED ON M.O. DAYHOFF S PAPER C A CONTOUR MAP PROGRAM FOR X-RAY CRYSTALLOGRAPHY C OCT. 1963 COMMUNICATIONS OF THE ACM. FIRST C PROGRAMMED BY BRUCE LANGDON, PLASMA PHYSICS C LAB,PRINCETON UNIVERSITY, NOV. 1966. C NR = 0 JERROR = 0 C C SCAN RIGHT EDGE AND BOTTOM C LY1 = LY+1 DO 10 J = LY1, UY IF (A(UX,J-1).LT.CL.AND.A(UX,J).GE.CL) THEN CALL CURVE (A,NX,UX,J,7,IER) JERROR = 1 C C WRITE (ILST,*) ' RIGHT EDGE ' C IF (IER.EQ.1) GO TO 70 ENDIF 10 CONTINUE C C SCAN TOP EDGE, LEFT TO RIGHT C UX1 = UX-1 DO 20 I1 = LX, UX1 I = UX1+LX-I1 IF (A(I+1,UY).LT.CL.AND.A(I,UY).GE.CL) THEN CALL CURVE (A,NX,I,UY,5,IER) JERROR = 1 C C WRITE (ILST,*) ' TOP EDGE ' C IF (IER.EQ.1) GO TO 70 ENDIF 20 CONTINUE C C SCAN LEFT EDGE TOP TO BOTTOM C UY1 = UY-1 DO 30 J1 = LY, UY1 J = UY1+LY-J1 IF (A(LX,J+1).LT.CL.AND.A(LX,J).GE.CL) THEN CALL CURVE (A,NX,LX,J,3,IER) C C WRITE (ILST,*) ' LEFT EDGE EDGE ' C JERROR = 1 IF (IER.EQ.1) GO TO 70 ENDIF 30 CONTINUE C C SCAN BOTTOM EDGE AND INTERIOR POINTS C LX1 = LX+1 DO 60 J = LY, UY1 DO 50 I = LX1, UX IF (A(I-1,J).LT.CL.AND.A(I,J).GE.CL) THEN IF (NR.NE.0) THEN COM = (UY-LY+1)*I+J DO 40 ID = 1, NR IF (REC(ID).EQ.COM) GO TO 50 40 CONTINUE ENDIF CALL CURVE (A,NX,I,J,1,IER) IF (J.EQ.LY) JERROR = 1 IF (IER.EQ.1) GO TO 70 ENDIF 50 CONTINUE 60 CONTINUE 70 RETURN END C C SUBROUTINE CURVE (A,IDIM,IXP,IYP,ISP,IER) INTEGER UX,UY,REC,DX,DY LOGICAL PC,LOPC,LPC DIMENSION A(IDIM,IDIM),X(850),Y(850),REC(850),INX(8),INY(8) DIMENSION Z(850) COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,X,Y,Z,IDASH,SCALE,PERZ,KA, * KO,ZO COMMON /CRPTPM/ NRMAX,NPMAX,LX,UX,LY,UY,CL,NR,NP,IXO,IYO,ISO,PC, * LOPC COMMON /TEMP/ REC COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 DATA INX / -1,-1,0,3*1,0,-1 / DATA INY / 0,3*1,0,3*-1 / C C IER = 0 IS = ISP IX = IXP IY = IYP IXO = IX IYO = IY ISO = IS NP = 1 LOPC = .FALSE. CALL INTERP (IX,IY,IS,A,IDIM,JER) GO TO (10,100,100,40), JER C C ROTATE C 10 IS = IS+1 IF (IS.EQ.9) IS = 1 DX = INX(IS) DY = INY(IS) C C FIND PLOT POINT C 20 CALL INTERP (IX,IY,IS,A,IDIM,JER) GO TO (30,60,80,110), JER 30 IF (IS.NE.1) GO TO 10 C C RECORD A CENTER C 40 IF (NR.GE.NRMAX) THEN C C REC ARRAY OVERFLOW C WRITE (ILST,50) 50 FORMAT ('0TOO MANY POINTS IN CONTOUR. IT IS TERMINATED') NNN = NP-1 CALL ROTPLT IER = 1 RETURN ENDIF NR = NR+1 REC(NR) = UY*IX+IY GO TO 10 C C DIAG FAIL C 60 IS = IS+1 IF (IS.EQ.9) IS = 1 LOPC = .TRUE. CALL INTERP (IX+INX(IS),IY+INY(IS),IS-3,A,IDIM,JER) GO TO (70,100,100,110), JER 70 IX = IX+DX IY = IY+DY IS = IS+4 IF (IS.GT.8) IS = IS-8 GO TO 20 C C NON-DIAG FAIL C 80 IX = IX+DX IY = IY+DY IS = IS+5 IF (IS.GT.8) IS = IS-8 LPC = PC CALL INTERP (IX,IY,IS,A,IDIM,JER) GO TO (90,90,100,110), JER 90 IF (.NOT.(LPC.AND.PC)) GO TO 10 C C PLOT POINT SWITCH C TEM = X(NP-2) X(NP-2) = X(NP-1) X(NP-1) = TEM TEM = Y(NP-2) Y(NP-2) = Y(NP-1) Y(NP-1) = TEM GO TO 10 100 STOP 110 NNN = NP-1 CALL ROTPLT RETURN END C C SUBROUTINE INTERP (IX,IY,ISP,A,IDIM,JER) INTEGER UX,UY,REC,DX,DY LOGICAL PC,LOPC DIMENSION A(IDIM,IDIM),X(850),Y(850),REC(850),INX(8),INY(8) DIMENSION Z(850) COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,X,Y,Z,IDASH,SCALE,PERZ,KA, * KO,ZO COMMON /CRPTPM/ NRMAX,NPMAX,LX,UX,LY,UY,CL,NR,NP,IXO,IYO,ISO,PC, * LOPC COMMON /SCLDAT/ XMAX,XMIN,YMAX,YMIN,XINCR,YINCR,CLX,NPLT COMMON /TEMP/ REC DATA INX / -1,-1,0,3*1,0,-1 / DATA INY / 0,3*1,0,3*-1 / C JER = 1 IS = ISP IF (IS.LT.1) IS = IS+8 DX = INX(IS) DY = INY(IS) FDX = FLOAT(DX)*XINCR FDY = FLOAT(DY)*YINCR FIX = FLOAT(IX-1)*XINCR+XMIN FIY = FLOAT(IY-1)*YINCR+YMIN IX1 = IX+DX IY1 = IY+DY IF (IX1.GE.LX.AND.IX1.LE.UX.AND.IY1.GE.LY.AND.IY1.LE.UY) THEN IF (DX*DY.EQ.0) THEN C C NON-DIAGONAL CASE C CHECK FOR FAIL C IF (A(IX1,IY1).GE.CL) GO TO 10 IF (DX.EQ.0) THEN X(NP) = FIX Y(NP) = (A(IX,IY)-CL)/(A(IX,IY)-A(IX,IY1))*FDY+FIY ELSE Y(NP) = FIY X(NP) = (A(IX,IY)-CL)/(A(IX,IY)-A(IX1,IY))*FDX+FIX ENDIF PC = .FALSE. ELSE C C DIAGONAL CASE C CP = (A(IX,IY)+A(IX1,IY)+A(IX,IY1)+A(IX1,IY1))*0.250E+0 IF ((CP.GE.CL.OR.LOPC)) THEN C C CONTOUR PASSES ON FAR SIDE OF CENTER POINT C IF (A(IX1,IY1).GE.CL.AND.CP.GE.CL) GO TO 20 V = (A(IX1,IY1)-CL)/(A(IX1,IY1)-CP)*0.50E+0 X(NP) = (1.0E+0-V)*FDX+FIX Y(NP) = (1.0E+0-V)*FDY+FIY PC = .TRUE. LOPC = .FALSE. ELSE C C CONTOUR PASSES ON NEAR SIDE OF CENTER POINT C V = (A(IX,IY)-CL)/(A(IX,IY)-CP)*0.50E+0 X(NP) = V*FDX+FIX Y(NP) = V*FDY+FIY PC = .FALSE. ENDIF ENDIF NP = NP+1 IF (IX.NE.IXO.OR.IY.NE.IYO.OR.IS.NE.ISO) THEN IF (NP.LE.NPMAX) RETURN C C PLOT PART OF CURVE AND CONTINUE C NNN = NP-1 CALL ROTPLT X(1) = X(NP-1) Y(1) = Y(NP-1) NP = 2 RETURN ENDIF ENDIF JER = JER+1 10 JER = JER+1 20 JER = JER+1 RETURN END