home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_disks / 200-299 / ff267.lzh / Diglib / diglib.zoo / diglib / CONTOR.FOR < prev    next >
Text File  |  1989-06-20  |  7KB  |  219 lines
  1.         SUBROUTINE CONTOR(Z,NZ,IZ,MX,MY,X1,XMX,Y1,YMY,NL,CL)
  2. C
  3. C       THIS SUBROUTINE WILL PRODUCE A CONTOUR PLOT OF THE FUNCTION
  4. C       DEFINED BY Z(I,J) = F(X(I),Y(J)).   IT IS ASSUMED THAT
  5. C       A CALL TO "MAPIT" HAS ALREADY BEEN MADE TO ESTABLISH THE
  6. C       COORDINATE AXIS (X,Y), WITH X LIMITS COVERING THE RANGE
  7. C       X1 TO XMX, AND Y LIMITS COVERING THE RANGE Y1 TO YMY.
  8. C
  9. C       Bad bug fixed: 23-APR-1985 by Hal R. Brand
  10. C
  11. CArguments:
  12. C
  13. C  Input
  14. C
  15. C       Z               * Type: real array.
  16. C                       * The values of the function to contour:
  17. C                          Z(I,J) = F(Xi,Yj) where:
  18. C                            Xi = X1 + (i-1)*(XMX-X1)/(MX-1)
  19. C                            Yj = Y1 + (j-1)*(YMX-Y1)/(MY-1)
  20. C
  21. C       NZ              * Type: integer constant or variable.
  22. C                       * The first dimension of the array Z - not necessaril
  23. C                               equal to MX, but MX <= NZ.
  24. C
  25. C       IZ              * Type: byte array.
  26. C                       * Used internally for working storage.
  27. C
  28. C       MX              * Type: integer constant or variable.
  29. C                       * The number of X grid points.
  30. C
  31. C       MY              * Type: integer constant or variable.
  32. C                       * The number of Y grid points.
  33. C
  34. C       X1              * Type: real constant or variable.
  35. C                       * The minimum X value.
  36. C
  37. C       XMX             * Type: real constant or variable.
  38. C                       * The maximum X value.
  39. C
  40. C       Y1              * Type: real constant or variable.
  41. C                       * The minimum Y value.
  42. C
  43. C       YMY             * Type: real constant or variable.
  44. C                       * The maximum Y value.
  45. C
  46. C       NL              * Type: integer constant or variable.
  47. C                       * The number of contour levels.
  48. C
  49. C       CL              * Type: real array.
  50. C                       * The coutour levels to draw.   (Same units as
  51. C                               F() or Z().)
  52. C
  53. C  Output
  54. C
  55. C    None.
  56. C
  57. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  58. C
  59.         DIMENSION Z(NZ,MY)
  60.         DIMENSION CL(NL)
  61.         INTEGER*1 IZ(MX,MY)
  62.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  63.      1   NX,NY,XL,DX,YL,DY
  64. C
  65. C       INITIALIZE ROUTINE
  66. C
  67.         XL = X1
  68.         YL = Y1
  69.         DX = XMX-X1
  70.         DY = YMY-Y1
  71.         NX=MX
  72.         NY=MY
  73.         NLOOP=MIN1(FLOAT(NX)/2.0+.5,FLOAT(NY)/2.0+.5)
  74. C       START SEARCHING FOR PLUS-MINUS TRANSITIONS
  75. C       TO START A CONTOR ON.
  76.         DO 50 NC=1,NL
  77. C
  78. C       ZERO ARRAY SHOWING WHERE WE HAVE BEEN
  79. C
  80.         DO 2 J=1,NY
  81.                 DO 1 I=1,NX
  82.                         IZ(I,J)=0
  83. 1                       CONTINUE
  84. 2               CONTINUE
  85.         CLEVEL=CL(NC)
  86.         DO 50 ICIR=1,NLOOP
  87.                 IU=NX+1-ICIR
  88.                 JU=NY+1-ICIR
  89.                 DO 10 J=ICIR,JU-1
  90.                         CALL LOOK(Z,ICIR,J,1,IZ,NZ,NX)
  91. 10                      CONTINUE
  92.                 DO 20 I=ICIR,IU-1
  93.                         CALL LOOK(Z,I,JU,2,IZ,NZ,NX)
  94. 20                      CONTINUE
  95.                 DO 30 J=JU,ICIR+1,-1
  96.                         CALL LOOK(Z,IU,J,3,IZ,NZ,NX)
  97. 30                      CONTINUE
  98.                 DO 40 I=IU,ICIR+1,-1
  99.                         CALL LOOK(Z,I,ICIR,4,IZ,NZ,NX)
  100. 40                      CONTINUE
  101. 50              CONTINUE
  102.         RETURN
  103.         END
  104. C
  105. C
  106. C
  107.         SUBROUTINE LOOK(Z,II,JJ,M,IZ,NZ,IZDIM)
  108.         INTEGER*1 IZ(IZDIM,2)
  109.         DIMENSION Z(NZ,2)
  110. C
  111. C       This subroutine looks for a contour starting at the point (II,JJ)
  112. C       with the contour being oriented such that the point (II,JJ) is
  113. C       greater than the current contouring level, and its neighbor (specifie
  114. C       by M) is less than the current contouring level.
  115. C
  116.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  117.      1   NX,NY,XL,DX,YL,DY
  118.         DIMENSION IDMODE(3,4)
  119.         DATA IDMODE/4,1,2,  1,2,3,  2,3,4,  3,4,1/
  120. C
  121.         IOLD=II
  122.         JOLD=JJ
  123.         MODE=M
  124.         CALL NEWP(1,MODE)
  125. C
  126. C       LOOK FOR CONTOR STARTING HERE.   THE "OLD" POINT IS ALWAYS THE
  127. C        POSITIVE ONE, SO THE TEST IS EASY.
  128. C
  129.         IF (Z(IOLD,JOLD) .GE. CLEVEL .AND. Z(IN,JN) .LT. CLEVEL) GOTO 20
  130. 100     RETURN
  131. C
  132. C       CHECK FOR CONTOR PREVIOUSLY THRU HERE - "SEGMNT" RETURNS THE POINT
  133. C        WE MARK WHEN EVER A CONTOUR PASSES THRU THE POINTS "(IOLD,JOLD)"
  134. C        AND "(IN,JN)".   "SEGMNT" ALSO RETURNS THE MARK THAT SHOULD BE
  135. C        PLACED GIVEN THE ORIENTATION OF THE CONTOUR.
  136. C
  137. 20      CALL SEGMNT(ICI,ICJ,ISEG)
  138.         IF ((IZ(ICI,ICJ).AND.ISEG) .NE. 0) RETURN
  139. C
  140. C       NEW CONTOUR.   TRACE IT TILL IT ENDS BY LOOPING BACK ON ITSELF, OR
  141. C        RUNNING OFF THE GRID.
  142. C
  143.         CALL ZPNT(XX,YY,Z,NZ)
  144.         CALL SCALE(XX,YY,VX,VY)
  145.         CALL GSMOVE(VX,VY)
  146.         IOLD=IN
  147.         JOLD=JN
  148. 30              CONTINUE
  149.                 DO 50 N=2,4
  150.                         CALL NEWP(N,MODE)
  151.                         IF (IN .LT. 1 .OR. IN .GT. NX) GO TO 100
  152.                         IF (JN .LT. 1 .OR. JN .GT. NY) GO TO 100
  153.                         IF (SIGN(1.0,Z(IOLD,JOLD)-CLEVEL) .NE.
  154.      1            SIGN(1.0,Z(IN,JN)-CLEVEL)) GO TO 60
  155.                         IOLD=IN
  156.                         JOLD=JN
  157. 50                      CONTINUE
  158. C
  159. C       IT IS IMPOSSIBLE TO JUST FALL THRU DUE TO THE ALGORITHM
  160. C
  161. C               STOP 'SERIOUS ERROR'
  162. 60              CONTINUE
  163. C
  164. C       FOUND THE NEXT INTERSECTION.   SEE IF IT HAS ALREADY BEEN MARKED.
  165. C        IF SO, THEN WE ARE DONE, ELSE, MARK IT, DRAW TO IT, AND CONTINUE ON.
  166. C
  167.                 CALL SEGMNT(ICI,ICJ,ISEG)
  168.                 IF ((IZ(ICI,ICJ).AND.ISEG) .NE. 0) RETURN
  169.                 IZ(ICI,ICJ)=(IZ(ICI,ICJ).OR.ISEG)
  170.                 CALL ZPNT(XX,YY,Z,NZ)
  171.                 CALL SCALE(XX,YY,VX,VY)
  172.                 CALL GSDRAW(VX,VY)
  173.                 MODE=IDMODE(N-1,MODE)
  174.                 GO TO 30
  175.         END
  176. C
  177. C
  178. C
  179.         SUBROUTINE SEGMNT(ICI,ICJ,ISEG)
  180.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  181.      1   NX,NY,XL,DX,YL,DY
  182.         ICI=MIN0(IOLD,IN)
  183.         ICJ=MIN0(JOLD,JN)
  184.         ISEG=1
  185.         IF (IOLD .EQ. IN) ISEG=2
  186.         RETURN
  187.         END
  188. C
  189. C
  190. C
  191.         SUBROUTINE NEWP(I,M)
  192.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  193.      1   NX,NY,XL,DX,YL,DY
  194.         DIMENSION IDELI(4),JDELJ(4)
  195.         DATA IDELI,JDELJ / 0,1,0,-1,   1,0,-1,0/
  196.         INDEX=MOD(2+I+M,4)+1
  197.         IN=IOLD+IDELI(INDEX)
  198.         JN=JOLD+JDELJ(INDEX)
  199.         RETURN
  200.         END
  201. C
  202. C
  203. C
  204.         SUBROUTINE ZPNT(X,Y,Z,NZ)
  205.         DIMENSION Z(NZ,2)
  206.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  207.      1   NX,NY,XL,DX,YL,DY
  208.         A=Z(IN,JN)-Z(IOLD,JOLD)
  209. C       IF NO CHANGE IN Z'S, PICK OLD POINT SO AS TO STAY TO RIGHT
  210.         IF (A .EQ. 0.0) GO TO 10
  211.         A=(CLEVEL-Z(IOLD,JOLD))/A
  212. 10      X=A*(IN-IOLD)+IOLD
  213.         Y=A*(JN-JOLD)+JOLD
  214. C       NOW CONVERT INDEXS TO X,Y VALUES
  215.         X=(X-1.0)*DX/(NX-1)+XL
  216.         Y=(Y-1.0)*DY/(NY-1)+YL
  217.         RETURN
  218.         END
  219.