home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / fortran77_210 / acmtoms / f77 / contour531 < prev    next >
Encoding:
Text File  |  1992-01-17  |  17.7 KB  |  623 lines
  1. ***********************  COPYRIGHT NOTICE  ********************************
  2. *  The material in this library is copyrighted by the ACM, which grants   *
  3. *  general permission to distribute provided the copies are not made for  *
  4. *  direct commercial advantage.  For details of the copyright and         *
  5. *  dissemination agreement, consult the current issue of TOMS.            *
  6. ***************************************************************************
  7.       PROGRAM CONTOUR
  8. C              from C.Johnson  but original routines from  ACM531
  9.       DIMENSION Z(51,51),C(10),WORK(1680)
  10. C     DIMENSION OF WORK IS LARGE ENOUGH TO CONTAIN
  11. C     2*(DIMENSION OF C)*(TOTAL DIMENSION OF Z) USEFUL BITS.  SEE THE
  12. C     BITMAP ROUTINES ACCESSED BY GCONTR.
  13.       REAL MU
  14.       EXTERNAL DRAW
  15.       COMMON /CUR/ XCUR,YCUR
  16.       COMMON /GRAF/ XPMIN,XPMAX,YPMIN,YPMAX
  17.       DATA C(1),C(2),C(3),C(4),C(5) /3.05,3.2,3.5,3.50135,3.6/
  18.       DATA C(6),C(7),C(8),C(9),C(10) /3.766413,4.0,4.130149,5.0,10.0/
  19.       DATA NX /51/, NY /51/, NF /10/
  20. C
  21.       XMIN=-2.0
  22.       XMAX=2.0
  23.       YMIN=-2.0
  24.       YMAX=2.0
  25.       MU=0.3
  26. C
  27.       XPMIN=0.0
  28.       XPMAX=FLOAT(NX)
  29.       YPMIN=0.0
  30.       YPMAX=FLOAT(NY)
  31. C
  32. C      CALL MODE(0)
  33.       CALL CLG
  34.       DX=(XMAX-XMIN)/FLOAT(NX-1)
  35.       DY=(YMAX-YMIN)/FLOAT(NY-1)
  36.       XCUR=1.0
  37.       YCUR=1.0
  38.       IF (MOD(NX,2).NE.0) YCUR=FLOAT(NY)
  39.       IF (MOD(NY,2).NE.0) XCUR=FLOAT(NX)
  40.       X=XMIN-DX
  41.       DO 20 I=1,NX
  42.        Y=YMIN-DY
  43.        X=X+DX
  44.        DO 10 J=1,NY
  45.         Y=Y+DY
  46.         Z(I,J)=(1.0-MU)*(2.0/SQRT((X-MU)**2+Y**2)+(X-MU)**2+Y**2)
  47.      *         +MU*(2.0/SQRT((X+1.0-MU)**2+Y**2)+(X+1.0-MU)**2+Y**2)
  48.   10   CONTINUE
  49.   20  CONTINUE
  50. C
  51.       NOPT=0
  52.       SX=1.0
  53.       SY=1.0
  54.       CALL CGRID(NOPT,NX,SX,XMIN,XMAX,NY,SY,YMIN,YMAX)
  55.       CALL GCONTR(Z,51,NX,NY,C,NF,2.E6,WORK,DRAW)
  56.       STOP
  57.       END
  58.       SUBROUTINE GCONTR(Z,NRZ,NX,NY,CV,NCV,ZMAX,BITMAP,DRAW)
  59. C
  60. C     THIS SUBROUTINE DRAWS A CONTOUR THROUGH EQUAL VALUES OF AN ARRAY.
  61. C
  62. C     *****     FORMAL ARGUMENTS     ***********************************
  63. C
  64. C     Z IS THE ARRAY FOR WHICH CONTOURS ARE TO BE DRAWN.  THE ELEMENTS
  65. C     OF Z ARE ASSUMED TO LIE UPON THE NODES OF A TOPOLOGICALLY
  66. C     RECTANGULAR COORDINATE SYSTEM - E.G. CARTESIAN, POLAR (EXCEPT
  67. C     THE ORIGIN), ETC.
  68. C
  69. C     NRZ IS THE NUMBER OF ROWS DECLARED FOR Z IN THE CALLING PROGRAM.
  70. C
  71. C     NX IS THE LIMIT FOR THE FIRST SUBSCRIPT OF Z.
  72. C
  73. C     NY IS THE LIMIT FOR THE SECOND SUBSCRIPT OF Z.
  74. C
  75. C     CV ARE THE VALUES OF THE CONTOURS TO BE DRAWN.
  76. C
  77. C     NCV IS THE NUMBER OF CONTOUR VALUES IN CV.
  78. C
  79. C     ZMAX IS THE MAXIMUM VALUE OF Z FOR CONSIDERATION.  A VALUE OF
  80. C     Z(I,J) GREATER THAN ZMAX IS A SIGNAL THAT THAT POINT AND THE
  81. C     GRID LINE SEGMENTS RADIATING FROM THAT POINT TO IT^S NEIGHBORS
  82. C     ARE TO BE EXCLUDED FROM CONTOURING.
  83. C
  84. C     BITMAP IS A WORK AREA LARGE ENOUGH TO HOLD 2*NX*NY*NCV BITS.  IT
  85. C     IS ACCESSED BY LOW-LEVEL ROUTINES, WHICH ARE DESCRIBED BELOW.
  86. C     LET J BE THE NUMBER OF USEFUL BITS IN EACH WORD OF BITMAP,
  87. C     AS DETERMINED BY THE USER MACHINE AND IMPLEMENTATION OF
  88. C     THE BITMAP MANIPULATION SUBPROGRAMS DESCRIBED BELOW.  THEN
  89. C     THE NUMBER OF WORDS REQUIRED FOR THE BITMAP IS THE FLOOR OF
  90. C         (2*NX*NY*NCV+J-1)/J.
  91. C
  92. C     DRAW IS A USER-PROVIDED SUBROUTINE USED TO DRAW CONTOURS.
  93. C     THE CALLING SEQUENCE FOR DRAW IS
  94. C
  95. C         CALL DRAW (X,Y,IFLAG)
  96. C         LET NX = INTEGER PART OF X, FX = FRACTIONAL PART OF X.
  97. C         THEN X SHOULD BE INTERPRETED SUCH THAT INCREASES IN NX
  98. C         CORRESPOND TO INCREASES IN THE FIRST SUBSCRIPT OF Z, AND
  99. C         FX IS THE FRACTIONAL DISTANCE FROM THE ABSCISSA CORRESPONDING
  100. C         TO NX TO THE ABSCISSA CORRESPONDING TO NX+1,
  101. C         AND Y SHOULD BE INTERPRETED SIMILARLY FOR THE SECOND
  102. C         SUBSCRIPT OF Z.
  103. C         THE LOW-ORDER DIGIT OF IFLAG WILL HAVE ONE OF THE VALUES
  104. C             1 - CONTINUE A CONTOUR,
  105. C             2 - START A CONTOUR AT A BOUNDARY,
  106. C             3 - START A CONTOUR NOT AT A BOUNDARY,
  107. C             4 - FINISH A CONTOUR AT A BOUNDARY,
  108. C             5 - FINISH A CLOSED CONTOUR (NOT AT A BOUNDARY).
  109. C                 NOTE THAT REQUESTS 1, 4 AND 5 ARE FOR PEN-DOWN
  110. C                 MOVES, AND THAT REQUESTS 2 AND 3 ARE FOR PEN-UP
  111. C                 MOVES.
  112. C             6 - SET X AND Y TO THE APPROXIMATE ^PEN^ POSITION, USING
  113. C                 THE NOTATION DISCUSSED ABOVE.  THIS CALL MAY BE
  114. C                 IGNORED, THE RESULT BEING THAT THE ^PEN^ POSITION
  115. C                 IS TAKEN TO CORRESPOND TO Z(1,1).
  116. C         IFLAG/10 IS THE CONTOUR NUMBER.
  117. C
  118. C     *****     EXTERNAL SUBPROGRAMS     *******************************
  119. C
  120. C     DRAW IS THE USER-SUPPLIED LINE DRAWING SUBPROGRAM DESCRIBED ABOVE.
  121. C     DRAW MAY BE SENSITIVE TO THE HOST COMPUTER AND TO THE PLOT DEVICE.
  122. C     FILL0 IS USED TO FILL A BITMAP WITH ZEROES.  CALL FILL0 (BITMAP,N)
  123. C     FILLS THE FIRST N BITS OF BITMAP WITH ZEROES.
  124. C     MARK1 IS USED TO PLACE A 1 IN A SPECIFIC BIT OF THE BITMAP.
  125. C     CALL MARK1 (BITMAP,N) PUTS A 1 IN THE NTH BIT OF THE BITMAP.
  126. C     IGET IS USED TO DETERMINE THE SETTING OF A PARTICULAR BIT IN THE
  127. C     BITMAP.  I=IGET(BITMAP,N) SETS I TO ZERO IF THE NTH BIT OF THE
  128. C     BITMAP IS ZERO, AND SETS I TO ONE IF THE NTH BIT IS ONE.
  129. C     FILL0, MARK1 AND IGET ARE MACHINE SENSITIVE.
  130. C
  131. C     ******************************************************************
  132. C
  133.       REAL Z(NRZ,*),CV(*)
  134.       INTEGER BITMAP(*)
  135.       INTEGER L1(4),L2(4),IJ(2)
  136. C
  137. C     L1 AND L2 CONTAIN LIMITS USED DURING THE SPIRAL SEARCH FOR THE
  138. C     BEGINNING OF A CONTOUR.
  139. C     IJ STORES SUBCRIPTS USED DURING THE SPIRAL SEARCH.
  140. C
  141.       INTEGER I1(2),I2(2),I3(6)
  142. C
  143. C     I1, I2 AND I3 ARE USED FOR SUBSCRIPT COMPUTATIONS DURING THE
  144. C     EXAMINATION OF LINES FROM Z(I,J) TO IT^S NEIGHBORS.
  145. C
  146.       REAL XINT(4)
  147. C
  148. C     XINT IS USED TO MARK INTERSECTIONS OF THE CONTOUR UNDER
  149. C     CONSIDERATION WITH THE EDGES OF THE CELL BEING EXAMINED.
  150. C
  151.       REAL XY(2)
  152. C
  153. C     XY IS USED TO COMPUTE COORDINATES FOR THE DRAW SUBROUTINE.
  154. C
  155.       EQUIVALENCE (L2(1),IMAX),(L2(2),JMAX),(L2(3),IMIN),(L2(4),JMIN)
  156.       EQUIVALENCE (IJ(1),I),(IJ(2),J)
  157.       EQUIVALENCE (XY(1),X),(XY(2),Y)
  158. C
  159.       DATA L1(3) /-1/, L1(4) /-1/
  160.       DATA I1 /1,0/, I2 /1,-1/, I3 /1,0,0,1,1,0/
  161. C
  162.       L1(1)=NX
  163.       L1(2)=NY
  164.       DMAX=ZMAX
  165. C
  166. C     SET THE CURRENT PEN POSITION.  THE DEFAULT POSITION CORRESPONDS
  167. C     TO Z(1,1).
  168. C
  169.       X=1.0
  170.       Y=1.0
  171.       CALL DRAW(X,Y,6)
  172.       ICUR=MAX0(1,MIN0(INT(X),NX))
  173.       JCUR=MAX0(1,MIN0(INT(Y),NY))
  174. C
  175. C     CLEAR THE BITMAP
  176. C
  177.       CALL FILL0(BITMAP,2*NX*NY*NCV)
  178. C
  179. C     SEARCH ALONG A RECTANGULAR SPIRAL PATH FOR A LINE SEGMENT HAVING
  180. C     THE FOLLOWING PROPERTIES
  181. C          1.  THE END POINTS ARE NOT EXCLUDED,
  182. C          2.  NO MARK HAS BEEN RECORDED FOR THE SEGMENT,
  183. C          3.  THE VALUES OF Z AT THE ENDS OF THE SEGMENT ARE SUCH THAT
  184. C              ONE Z IS LESS THAN THE CURRENT CONTOUR VALUE, AND THE
  185. C              OTHER IS GREATER THAN OR EQUAL TO THE CURRENT CONTOUR
  186. C              VALUE.
  187. C
  188. C     SEARCH ALL BOUNDARIES FIRST, THEN SEARCH INTERIOR LINE SEGMENTS.
  189. C     NOTE THAT THE INTERIOR LINE SEGMENTS NEAR EXCLUDED POINTS MAY BE
  190. C     BOUNDARIES.
  191. C
  192.       IBKEY=0
  193.   10  I=ICUR
  194.       J=JCUR
  195.   20  IMAX=I
  196.       IMIN=-I
  197.       JMAX=J
  198.       JMIN=-J
  199.       IDIR=0
  200. C     DIRECTION ZERO IS +I, 1 IS +J, 2 IS -I, 3 IS -J.
  201.   30  NXIDIR=IDIR+1
  202.       K=NXIDIR
  203.       IF (NXIDIR.GT.3) NXIDIR=0
  204.   40  I=IABS(I)
  205.       J=IABS(J)
  206.       IF (Z(I,J).GT.DMAX) GOTO 140
  207.       L=1
  208. C     L=1 MEANS HORIZONTAL LINE, L=2 MEANS VERTICAL LINE.
  209.   50  IF (IJ(L).GE.L1(L)) GOTO 130
  210.       II=I+I1(L)
  211.       JJ=J+I1(3-L)
  212.       IF (Z(II,JJ).GT.DMAX) GOTO 130
  213.       ASSIGN 100 TO JUMP
  214. C     THE NEXT 15 STATEMENTS (OR SO) DETECT BOUNDARIES.
  215.   60  IX=1
  216.       IF (IJ(3-L).EQ.1) GOTO 80
  217.       II=I-I1(3-L)
  218.       JJ=J-I1(L)
  219.       IF (Z(II,JJ).GT.DMAX) GOTO 70
  220.       II=I+I2(L)
  221.       JJ=J+I2(3-L)
  222.       IF (Z(II,JJ).LT.DMAX) IX=0
  223.   70  IF (IJ(3-L).GE.L1(3-L)) GOTO 90
  224.   80  II=I+I1(3-L)
  225.       JJ=J+I1(L)
  226.       IF (Z(II,JJ).GT.DMAX) GOTO 90
  227.       IF (Z(I+1,J+1).LT.DMAX) GOTO JUMP, (100, 280)
  228.   90  IX=IX+2
  229.       GOTO JUMP,(100, 280)
  230.  100  IF (IX.EQ.3) GOTO 130
  231.       IF (IX+IBKEY.EQ.0) GOTO 130
  232. C     NOW DETERMINE WHETHER THE LINE SEGMENT IS CROSSED BY THE CONTOUR.
  233.       II=I+I1(L)
  234.       JJ=J+I1(3-L)
  235.       Z1=Z(I,J)
  236.       Z2=Z(II,JJ)
  237.       DO 120 ICV=1,NCV
  238.       IF (IGET(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L).NE.0) GOTO 120
  239.       IF (CV(ICV).LE.AMIN1(Z1,Z2)) GOTO 110
  240.       IF (CV(ICV).LE.AMAX1(Z1,Z2)) GOTO 190
  241.  110  CALL MARK1(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L)
  242.  120  CONTINUE
  243.  130  L=L+1
  244.       IF (L.LE.2) GOTO 50
  245.  140  L=MOD(IDIR,2)+1
  246.       IJ(L)=ISIGN(IJ(L),L1(K))
  247. C
  248. C     LINES FROM Z(I,J) TO Z(I+1,J) AND Z(I,J+1) ARE NOT SATISFACTORY.
  249. C     CONTINUE THE SPIRAL.
  250. C
  251.  150  IF (IJ(L).GE.L1(K)) GOTO 170
  252.       IJ(L)=IJ(L)+1
  253.       IF (IJ(L).GT.L2(K)) GOTO 160
  254.       GOTO 40
  255.  160  L2(K)=IJ(L)
  256.       IDIR=NXIDIR
  257.       GOTO 30
  258.  170  IF (IDIR.EQ.NXIDIR) GOTO 180
  259.       NXIDIR=NXIDIR+1
  260.       IJ(L)=L1(K)
  261.       K=NXIDIR
  262.       L=3-L
  263.       IJ(L)=L2(K)
  264.       IF (NXIDIR.GT.3) NXIDIR=0
  265.       GOTO 150
  266.  180  IF (IBKEY.NE.0) RETURN
  267.       IBKEY=1
  268.       GOTO 10
  269. C
  270. C     AN ACCEPTABLE LINE SEGMENT HAS BEEN FOUND.
  271. C     FOLLOW THE CONTOUR UNTIL IT EITHER HITS A BOUNDARY OR CLOSES.
  272. C
  273.  190  IEDGE=L
  274.       CVAL=CV(ICV)
  275.       IF (IX.NE.1) IEDGE=IEDGE+2
  276.       IFLAG=2+IBKEY
  277.       XINT(IEDGE)=(CVAL-Z1)/(Z2-Z1)
  278.  200  XY(L)=FLOAT(IJ(L))+XINT(IEDGE)
  279.       XY(3-L)=FLOAT(IJ(3-L))
  280.       CALL MARK1(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L)
  281.       CALL DRAW(X,Y,IFLAG+10*ICV)
  282.       IF (IFLAG.LT.4) GOTO 210
  283.       ICUR=I
  284.       JCUR=J
  285.       GOTO 20
  286. C
  287. C     CONTINUE A CONTOUR.  THE EDGES ARE NUMBERED CLOCKWISE WITH
  288. C     THE BOTTOM EDGE BEING EDGE NUMBER ONE.
  289. C
  290.  210  NI=1
  291.       IF (IEDGE.LT.3) GOTO 220
  292.       I=I-I3(IEDGE)
  293.       J=J-I3(IEDGE+2)
  294.  220  DO 250 K=1,4
  295.       IF (K.EQ.IEDGE) GOTO 250
  296.       II=I+I3(K)
  297.       JJ=J+I3(K+1)
  298.       Z1=Z(II,JJ)
  299.       II=I+I3(K+1)
  300.       JJ=J+I3(K+2)
  301.       Z2=Z(II,JJ)
  302.       IF (CVAL.LE.AMIN1(Z1,Z2)) GOTO 250
  303.       IF (CVAL.GT.AMAX1(Z1,Z2)) GOTO 250
  304.       IF (K.EQ.1) GOTO 230
  305.       IF (K.NE.4) GOTO 240
  306.  230  ZZ=Z1
  307.       Z1=Z2
  308.       Z2=ZZ
  309.  240  XINT(K)=(CVAL-Z1)/(Z2-Z1)
  310.       NI=NI+1
  311.       KS=K
  312.  250  CONTINUE
  313.       IF (NI.EQ.2) GOTO 260
  314. C
  315. C     THE CONTOUR CROSSES ALL FOUR EDGES OF THE CELL BEING EXAMINED.
  316. C     CHOOSE THE LINES TOP-TO-LEFT AND BOTTOM-TO-RIGHT IF THE
  317. C     INTERPOLATION POINT ON THE TOP EDGE IS LESS THAN THE INTERPOLATION
  318. C     POINT ON THE BOTTOM EDGE.  OTHERWISE, CHOOSE THE OTHER PAIR.  THIS
  319. C     METHOD PRODUCES THE SAME RESULTS IF THE AXES ARE REVERSED.  THE
  320. C     CONTOUR MAY CLOSE AT ANY EDGE, BUT MUST NOT CROSS ITSELF INSIDE
  321. C     ANY CELL.
  322. C
  323.       KS=5-IEDGE
  324.       IF (XINT(3).LT.XINT(1)) GOTO 260
  325.       KS=3-IEDGE
  326.       IF (KS.LE.0) KS=KS+4
  327. C
  328. C     DETERMINE WHETHER THE CONTOUR WILL CLOSE OR RUN INTO A BOUNDARY
  329. C     AT EDGE KS OF THE CURRENT CELL.
  330. C
  331.  260  L=KS
  332.       IFLAG=1
  333.       ASSIGN 280 TO JUMP
  334.       IF (KS.LT.3) GOTO 270
  335.       I=I+I3(KS)
  336.       J=J+I3(KS+2)
  337.       L=KS-2
  338.  270  IF (IGET(BITMAP,2*(NX*(NY*(ICV-1)+J-1)+I-1)+L).EQ.0) GOTO 60
  339.       IFLAG=5
  340.       GOTO 290
  341.  280  IF (IX.NE.0) IFLAG=4
  342.  290  IEDGE=KS+2
  343.       IF (IEDGE.GT.4) IEDGE=IEDGE-4
  344.       XINT(IEDGE)=XINT(KS)
  345.       GOTO 200
  346. C
  347.       END
  348.       SUBROUTINE FILL0(BITMAP,N)
  349. C
  350. C     FILL THE FIRST N BITS OF BITMAP WITH ZEROES.
  351. C
  352.       INTEGER BITMAP(*),N
  353. C
  354. C     DATA NBPW /35/
  355.       DATA NBPW /31/
  356. C     NBPW IS THE MINIMUM NUMBER OF SIGNIFICANT BITS PER WORD USED
  357. C     BY INTEGER ARITHMETIC.  THIS IS USUALLY ONE LESS THAN THE
  358. C     ACTUAL NUMBER OF BITS PER WORD, BUT AN IMPORTANT EXCEPTION IS
  359. C     THE CDC-6000 SERIES OF MACHINES, WHERE NBPW SHOULD BE 48.
  360. C
  361.       LOOP=N/NBPW
  362.       NBLW=MOD(N,NBPW)
  363.       IF (LOOP.EQ.0) GOTO 20
  364.       DO 10 I=1,LOOP
  365.       BITMAP(I)=0
  366.   10  CONTINUE
  367.   20  IF (NBLW.NE.0) BITMAP(LOOP+1)=MOD(BITMAP(LOOP+1),2**(NBPW-NBLW))
  368.       RETURN
  369.       END
  370.       SUBROUTINE MARK1(BITMAP,N)
  371. C
  372. C     PUT A ONE IN THE NTH BIT OF BITMAP.
  373. C
  374.       INTEGER BITMAP(*),N
  375. C
  376. C     DATA NBPW /35/
  377.       DATA NBPW /31/
  378. C     NBPW IS THE MINIMUM NUMBER OF SIGNIFICANT BITS PER WORD USED
  379. C     BY INTEGER ARITHMETIC.  THIS IS USUALLY ONE LESS THAN THE
  380. C     ACTUAL NUMBER OF BITS PER WORD, BUT AN IMPORTANT EXCEPTION IS
  381. C     THE CDC-6000 SERIES OF MACHINES, WHERE NBPW SHOULD BE 48.
  382. C
  383.       NWORD=(N-1)/NBPW
  384.       NBIT=MOD(N-1,NBPW)
  385.       I=2**(NBPW-NBIT-1)
  386.       BITMAP(NWORD+1)=BITMAP(NWORD+1)+I*(1-MOD(BITMAP(NWORD+1)/I,2))
  387.       RETURN
  388.       END
  389.       FUNCTION IGET(BITMAP,N)
  390. C
  391. C     IGET=0 IF THE NTH BIT OF BITMAP IS ZERO, ELSE IGET IS ONE.
  392. C
  393.       INTEGER BITMAP(*),N
  394. C     DATA NBPW /35/
  395. C
  396.       DATA NBPW /31/
  397. C     NBPW IS THE MINIMUM NUMBER OF SIGNIFICANT BITS PER WORD USED
  398. C     BY INTEGER ARITHMETIC.  THIS IS USUALLY ONE LESS THAN THE
  399. C     ACTUAL NUMBER OF BITS PER WORD, BUT AN IMPORTANT EXCEPTION IS
  400. C     THE CDC-6000 SERIES OF MACHINES, WHERE NBPW SHOULD BE 48.
  401. C
  402.       NWORD=(N-1)/NBPW
  403.       NBIT=MOD(N-1,NBPW)
  404.       IGET=MOD(BITMAP(NWORD+1)/2**(NBPW-NBIT-1),2)
  405.       RETURN
  406.       END
  407.       SUBROUTINE DRAW(X,Y,IFLAG)
  408. C     THIS SUBROUTINE USES CALCOMP PLOT ROUTINES TO DRAW LINES FOR THE
  409. C     CONTOUR PLOTTING ROUTINE GCONTR.
  410.  
  411.       DATA IBLANK /32/
  412.       IH=IFLAG/10
  413.       IL=IFLAG-10*IH
  414.       IF (IL.EQ.6) GOTO 40
  415.       IPEN=2
  416.       IF (IL.EQ.2) IPEN=3
  417.       IF (IL.EQ.3) IPEN=3
  418.       XCUR=X
  419.       YCUR=Y
  420. C      XX=(X-1.0)*XL
  421. C      YY=(Y-1.0)*YL
  422.       XX=X
  423.       YY=Y
  424.       CALL CPLOT(XX,YY,IPEN)
  425.       IF (IL.LT.2) GOTO 30
  426.       IF (IL.GT.4) GOTO 30
  427.       IF (IH.NE.0) GOTO 10
  428.       CALL NUMBER(XX,YY-0.03,0.07,IH,0.0,-1)
  429.       GOTO 20
  430.   10  CALL SYMBOL(XX,YY-0.03,0.07,IH,0.0,NCH)
  431.   20  CALL CPLOT(XX,YY,3)
  432.   30  RETURN
  433.   40  X=XCUR
  434.       Y=YCUR
  435.       RETURN
  436.       END
  437.       SUBROUTINE CPLOT(X,Y,IPEN)
  438.       IF (IPEN.EQ.2) CALL LINETO(X,Y)
  439.       IF (IPEN.EQ.3) CALL MOVETO(X,Y)
  440.       RETURN
  441.       END
  442.       SUBROUTINE NUMBER(X,Y,PENSIZE,ICH,P,NCH)
  443.       CALL VDU(5)
  444.       CALL MOVETO(X,Y)
  445.       CALL VDU(48+ICH)
  446.       CALL VDU(4)
  447.       RETURN
  448.       END
  449.       SUBROUTINE SYMBOL(X,Y,PENSIZE,ICH,P,NCH)
  450. C
  451. C First version = same as SUBROUTINE NUMBER above.
  452. C
  453.       CALL VDU(5)
  454.       CALL MOVETO(X,Y)
  455.       CALL VDU(48+ICH)
  456.       CALL VDU(4)
  457.       RETURN
  458.       END
  459.       SUBROUTINE CLG
  460.       CALL VDU(16)
  461.       RETURN
  462.       END
  463.       SUBROUTINE LINETO(XPT,YPT)
  464.       INTEGER FNX,FNY
  465.       CALL PLOT(5,FNX(XPT),FNY(YPT))
  466.       RETURN
  467.       END
  468.       SUBROUTINE MOVETO(XPT,YPT)
  469.       INTEGER FNX,FNY
  470.       CALL PLOT(4,FNX(XPT),FNY(YPT))
  471.       RETURN
  472.       END
  473.       INTEGER FUNCTION FNX(Z)
  474.       COMMON /GRAF/ XPMIN,XPMAX,YPMIN,YPMAX
  475.       REAL Z
  476.       FNX=INT(1023.0*(Z-XPMIN)/(XPMAX-XPMIN)+0.5)
  477.       RETURN
  478.       END
  479.       INTEGER FUNCTION FNY(Z)
  480.       COMMON /GRAF/ XPMIN,XPMAX,YPMIN,YPMAX
  481.       REAL Z
  482.       FNY=INT(1023.0*(Z-YPMIN)/(YPMAX-YPMIN)+0.5)
  483.       RETURN
  484.       END
  485.       SUBROUTINE CGRID(NOPT,NX,SX,XS,XF,NY,SY,YS,YF)
  486. C
  487. C SUBROUTINE WHICH DRAWS A FRAME AROUND THE PLOT AND DRAWS
  488. C EITHER TICK MARKS OR GRID LINES.
  489. C
  490. C PARAMETERS   NOPT -- =0, DRAW TICKS ONLY
  491. C                      =1, DRAW GRID LINES
  492. C                      =2, DRAW GRID LINES TO EDGE OF FRAME.
  493. C              NX -- NUMBER OF INTERVALS IN X DIRECTION
  494. C              SX -- SPACING IN INCHES BETWEEN TICK MARKS OR GRID LINES
  495. C                    ALONG THE X AXIS
  496. C              XS -- LOCATION OF FIRST TICK OR GRID LINE ON X AXIS
  497. C              XF -- LOCATION OF RIGHT EDGE OF FRAME
  498. C              NY -- NUMBER OF INTERVALS IN Y DIRECTION
  499. C              SY -- SPACING IN INCHES BETWEEN TICK MARKS OR GRID LINES
  500. C                    ALONG THE Y AXIS
  501. C              YS -- LOCATION OF FIRST TICK OR GRID LINE ON Y AXIS
  502. C              YF -- LOCATION OF TOP EDGE OF FRAME
  503. C ASSUMPTIONS  NX, SX, NY, SY ALL POSITIVE.
  504. C              THE LOWER LEFT-HAND CORNER OF THE FRAME IS DRAWN AT (0,0)
  505. C              IF XS<0, USE 0; IF YS<0, USE 0
  506. C              IF XF<=0, USE NX*SX; IF YF<=0, USE NY*SY.
  507. C
  508.       XINC=SX
  509.       YINC=SY
  510.       XLGTH=FLOAT(NX)*SX
  511.       YLGTH=FLOAT(NY)*SY
  512.       XMIN=AMAX1(XS,0.0)
  513.       YMIN=AMAX1(YS,0.0)
  514.       XMAX=AMAX1(XF,XLGTH+XMIN)
  515.       YMAX=AMAX1(YF,YLGTH+YMIN)
  516. C
  517. C     DRAW FRAME.
  518. C
  519.       CALL CPLOT(0.0,0.0,3)
  520.       CALL CPLOT(XMAX,0.0,2)
  521.       CALL CPLOT(XMAX,YMAX,2)
  522.       CALL CPLOT(0.0,YMAX,2)
  523.       CALL CPLOT(0.0,0.0,2)
  524.       IF (NOPT.NE.0) GOTO 130
  525. C
  526. C     DRAW TICK MARKS.
  527. C
  528.       DO 120 J=1,4
  529.       GOTO (10,50,20,40),J
  530.   10  X2=0.0
  531.       IF (XMIN.NE.0.0) X2=XMIN-SX
  532.       Y2=0.0
  533.       GOTO 30
  534.   20  XINC=-SX
  535.       X2=XMIN+XLGTH+SX
  536.       IF (XMAX.EQ.XMIN+XLGTH) X2=XMAX
  537.       Y2=YMAX
  538.   30  Y1=Y2
  539.       Y2=Y2+SIGN(0.125,XINC)
  540.       N=NX
  541.       IF (ABS(XMAX-XMIN-XLGTH)+ABS(XMIN)) 70,80,70
  542.   40  YINC=-SY
  543.       Y2=YMIN+YLGTH+SY
  544.       IF (YMAX.EQ.YMIN+YLGTH) Y2=YMAX
  545.       X2=0.0
  546.       GOTO 60
  547.   50  Y2=0.0
  548.       IF (YMIN.NE.0.0) Y2=YMIN-SY
  549.       X2=XMAX
  550.   60  X1=X2
  551.       N=NY
  552.       X2=X2-SIGN(0.125,YINC)
  553.       IF (ABS(YMAX-YMIN-YLGTH)+ABS(YMIN)) 70,80,70
  554.   70  N=N+1
  555.   80  DO 110 I=1,N
  556.       IF (MOD(J,2).EQ.0) GOTO 90
  557.       X2=X2+XINC
  558.       X1=X2
  559.       GOTO 100
  560.   90  Y2=Y2+YINC
  561.       Y1=Y2
  562.  100  CALL CPLOT(X1,Y1,3)
  563.       CALL CPLOT(X2,Y2,2)
  564.  110  CONTINUE
  565.  120  CONTINUE
  566.       GOTO 240
  567. C
  568. C     DRAW GRID LINES
  569. C
  570.  130  X1=XMIN
  571.       X2=XMIN+XLGTH
  572.       IF (NOPT.NE.2) GOTO 140
  573.       X1=0.0
  574.       X2=XMAX
  575.  140  Y1=YMIN-SY
  576.       N=NY+1
  577.       IF (YMAX.EQ.YMIN+YLGTH) N=N-1
  578.       IF (YMIN.NE.0.0) GOTO 150
  579.       Y1=0.0
  580.       N=N-1
  581.  150  IF (N.LE.0) GOTO 170
  582.       J=1
  583.       DO 160 I=1,N
  584.       J=-J
  585.       Y1=Y1+SY
  586.       CALL CPLOT(X1,Y1,3)
  587.       CALL CPLOT(X2,Y1,2)
  588.       XX=X1
  589.       X1=X2
  590.       X2=XX
  591.  160  CONTINUE
  592.  170  Y1=YMIN+YLGTH
  593.       Y2=YMIN
  594.       IF (NOPT.NE.2) GOTO 180
  595.       Y1=YMAX
  596.       Y2=0.0
  597.  180  N=NX+1
  598.       IF (J.LT.0) GOTO 200
  599.       X1=XMIN-SX
  600.       IF (XMAX.EQ.XMIN+XLGTH) N=N-1
  601.       IF (XMIN.NE.0.0) GOTO 190
  602.       X1=0.0
  603.       N=N-1
  604.  190  IF (N.LE.0) GOTO 240
  605.       XINC=SX
  606.       GOTO 220
  607.  200  X1=XMIN+XLGTH+SX
  608.       IF (XMIN.EQ.0.0) N=N-1
  609.       IF (XMAX.NE.XLGTH+XMIN) GOTO 210
  610.       N=N-1
  611.       X1=XMAX
  612.  210  XINC=-SX
  613.  220  DO 230 I=1,N
  614.       X1=X1+XINC
  615.       CALL CPLOT(X1,Y1,3)
  616.       CALL CPLOT(X1,Y2,2)
  617.       XX=Y1
  618.       Y1=Y2
  619.       Y2=XX
  620.  230  CONTINUE
  621.  240  RETURN
  622.       END
  623.