SUBROUTINE CONTOR(Z,NZ,IZ,MX,MY,X1,XMX,Y1,YMY,NL,CL) IMPLICIT NONE C C THIS SUBROUTINE WILL PRODUCE A CONTOUR PLOT OF THE FUNCTION C DEFINED BY Z(I,J) = F(X(I),Y(J)). IT IS ASSUMED THAT C A CALL TO "MAPIT" HAS ALREADY BEEN MADE TO ESTABLISH THE C COORDINATE AXIS (X,Y), WITH X LIMITS COVERING THE RANGE C X1 TO XMX, AND Y LIMITS COVERING THE RANGE Y1 TO YMY. C C Bad bug fixed: 23-APR-1985 by Hal R. Brand C CArguments: C C Input C C Z * Type: real array. C * The values of the function to contour: C Z(I,J) = F(Xi,Yj) where: C Xi = X1 + (i-1)*(XMX-X1)/(MX-1) C Yj = Y1 + (j-1)*(YMX-Y1)/(MY-1) C C NZ * Type: integer constant or variable. C * The first dimension of the array Z - not necessaril C equal to MX, but MX <= NZ. C C IZ * Type: byte array. C * Used internally for working storage. C C MX * Type: integer constant or variable. C * The number of X grid points. C C MY * Type: integer constant or variable. C * The number of Y grid points. C C X1 * Type: real constant or variable. C * The minimum X value. C C XMX * Type: real constant or variable. C * The maximum X value. C C Y1 * Type: real constant or variable. C * The minimum Y value. C C YMY * Type: real constant or variable. C * The maximum Y value. C C NL * Type: integer constant or variable. C * The number of contour levels. C C CL * Type: real array. C * The coutour levels to draw. (Same units as C F() or Z().) C C Output C C None. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REAL*4 Z(NZ,MY),X1,Y1,YMY REAL*4 CL(NL) INTEGER*1 IZ(MX,MY) REAL*4 CLEVEL,XL,DX,YL,DY INTEGER IOLD,JOLD,IN,JN,NX,NY,NLOOP COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN, 1 NX,NY,XL,DX,YL,DY INTEGER I,J,NL,ICIR,IU,JU,NC C C INITIALIZE ROUTINE C XL = X1 YL = Y1 DX = XMX-X1 DY = YMY-Y1 NX=MX NY=MY NLOOP=MIN1(FLOAT(NX)/2.0+.5,FLOAT(NY)/2.0+.5) C START SEARCHING FOR PLUS-MINUS TRANSITIONS C TO START A CONTOR ON. DO 50 NC=1,NL C C ZERO ARRAY SHOWING WHERE WE HAVE BEEN C DO 2 J=1,NY DO 1 I=1,NX IZ(I,J)=0 1 CONTINUE 2 CONTINUE CLEVEL=CL(NC) DO 50 ICIR=1,NLOOP IU=NX+1-ICIR JU=NY+1-ICIR DO 10 J=ICIR,JU-1 CALL LOOK(Z,ICIR,J,1,IZ,NZ,NX) 10 CONTINUE DO 20 I=ICIR,IU-1 CALL LOOK(Z,I,JU,2,IZ,NZ,NX) 20 CONTINUE DO 30 J=JU,ICIR+1,-1 CALL LOOK(Z,IU,J,3,IZ,NZ,NX) 30 CONTINUE DO 40 I=IU,ICIR+1,-1 CALL LOOK(Z,I,ICIR,4,IZ,NZ,NX) 40 CONTINUE 50 CONTINUE RETURN END C C C SUBROUTINE LOOK(Z,II,JJ,M,IZ,NZ,IZDIM) INTEGER*1 IZ(IZDIM,2) DIMENSION Z(NZ,2) C C This subroutine looks for a contour starting at the point (II,JJ) C with the contour being oriented such that the point (II,JJ) is C greater than the current contouring level, and its neighbor (specifie C by M) is less than the current contouring level. C COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN, 1 NX,NY,XL,DX,YL,DY DIMENSION IDMODE(3,4) DATA IDMODE/4,1,2, 1,2,3, 2,3,4, 3,4,1/ C IOLD=II JOLD=JJ MODE=M CALL NEWP(1,MODE) C C LOOK FOR CONTOR STARTING HERE. THE "OLD" POINT IS ALWAYS THE C POSITIVE ONE, SO THE TEST IS EASY. C IF (Z(IOLD,JOLD) .GE. CLEVEL .AND. Z(IN,JN) .LT. CLEVEL) GOTO 20 100 RETURN C C CHECK FOR CONTOR PREVIOUSLY THRU HERE - "SEGMNT" RETURNS THE POINT C WE MARK WHEN EVER A CONTOUR PASSES THRU THE POINTS "(IOLD,JOLD)" C AND "(IN,JN)". "SEGMNT" ALSO RETURNS THE MARK THAT SHOULD BE C PLACED GIVEN THE ORIENTATION OF THE CONTOUR. C 20 CALL SEGMNT(ICI,ICJ,ISEG) IF ((IZ(ICI,ICJ).AND.ISEG) .NE. 0) RETURN C C NEW CONTOUR. TRACE IT TILL IT ENDS BY LOOPING BACK ON ITSELF, OR C RUNNING OFF THE GRID. C CALL ZPNT(XX,YY,Z,NZ) CALL SCALE(XX,YY,VX,VY) CALL GSMOVE(VX,VY) IOLD=IN JOLD=JN 30 CONTINUE DO 50 N=2,4 CALL NEWP(N,MODE) IF (IN .LT. 1 .OR. IN .GT. NX) GO TO 100 IF (JN .LT. 1 .OR. JN .GT. NY) GO TO 100 IF (SIGN(1.0,Z(IOLD,JOLD)-CLEVEL) .NE. 1 SIGN(1.0,Z(IN,JN)-CLEVEL)) GO TO 60 IOLD=IN JOLD=JN 50 CONTINUE C C IT IS IMPOSSIBLE TO JUST FALL THRU DUE TO THE ALGORITHM C C STOP 'SERIOUS ERROR' 60 CONTINUE C C FOUND THE NEXT INTERSECTION. SEE IF IT HAS ALREADY BEEN MARKED. C IF SO, THEN WE ARE DONE, ELSE, MARK IT, DRAW TO IT, AND CONTINUE ON. C CALL SEGMNT(ICI,ICJ,ISEG) IF ((IZ(ICI,ICJ).AND.ISEG) .NE. 0) RETURN IZ(ICI,ICJ)=(IZ(ICI,ICJ).OR.ISEG) CALL ZPNT(XX,YY,Z,NZ) CALL SCALE(XX,YY,VX,VY) CALL GSDRAW(VX,VY) MODE=IDMODE(N-1,MODE) GO TO 30 END C C C SUBROUTINE SEGMNT(ICI,ICJ,ISEG) COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN, 1 NX,NY,XL,DX,YL,DY ICI=MIN0(IOLD,IN) ICJ=MIN0(JOLD,JN) ISEG=1 IF (IOLD .EQ. IN) ISEG=2 RETURN END C C C SUBROUTINE NEWP(I,M) COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN, 1 NX,NY,XL,DX,YL,DY DIMENSION IDELI(4),JDELJ(4) DATA IDELI,JDELJ / 0,1,0,-1, 1,0,-1,0/ INDEX=MOD(2+I+M,4)+1 IN=IOLD+IDELI(INDEX) JN=JOLD+JDELJ(INDEX) RETURN END C C C SUBROUTINE ZPNT(X,Y,Z,NZ) DIMENSION Z(NZ,2) COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN, 1 NX,NY,XL,DX,YL,DY A=Z(IN,JN)-Z(IOLD,JOLD) C IF NO CHANGE IN Z'S, PICK OLD POINT SO AS TO STAY TO RIGHT IF (A .EQ. 0.0) GO TO 10 A=(CLEVEL-Z(IOLD,JOLD))/A 10 X=A*(IN-IOLD)+IOLD Y=A*(JN-JOLD)+JOLD C NOW CONVERT INDEXS TO X,Y VALUES X=(X-1.0)*DX/(NX-1)+XL Y=(Y-1.0)*DY/(NY-1)+YL RETURN END