home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d2xx / d267 / diglib.lha / Diglib / diglib.zoo / diglib / PURJOY.FOR < prev    next >
Encoding:
Text File  |  1989-06-20  |  30.8 KB  |  923 lines
  1.  
  2.         SUBROUTINE PURJOY(Z,IZDIM1,IZ,KX,KY,CAMLOC,XYLIM,
  3.      1   XLAB,YLAB,ZLAB,CSIZE,MARPLT)
  4. C
  5. CPurpose:   This subroutine will plot a function Z=F(X,Y) as a lined surface.
  6. C       The function must be defined on a regular grid.   This routine will
  7. C       optionally remove hidden lines.
  8. C
  9. CArguments:
  10. C
  11. C  Input
  12. C
  13. C       Z               * Type: real array.
  14. C                       * The function values: Z(I,J)=F(Xi,Yj), where
  15. C                               Xi = XMIN + (i-1)*(XMAX-XMIN)/(KX-1)
  16. C                               Yj = YMIN + (j-1)*(YMAX-YMIN)/(KY-1)
  17. C
  18. C       IZDIM1          * Type: integer constant or variable.
  19. C                       * The first dimension of the Z array - not
  20. C                               necessarily the number of X values.
  21. C
  22. C       IZ              * Type: byte array.
  23. C                       * A working array of bytes dimensioned atleast
  24. C                               KX*KY long.
  25. C
  26. C       KX              * Type: integer constant or variable.
  27. C                       * The number of X values in the Z array.
  28. C                               KX <= IZDIM1 ofcourse.
  29. C
  30. C       KY              * Type: integer constant or variable.
  31. C                       * The number of Y values in the Z array.
  32. C
  33. C       CAMLOC          * Type: real array.
  34. C                       * The relative location of the viewer in space.
  35. C                               The viewer always faces toward the center
  36. C                               of the surface.
  37. C                               CAMLOC(1) = distance from surface in units
  38. C                                the same as those of Z.
  39. C                               CAMLOC(2) = angle between the viewer and the
  40. C                                X axis in degrees.   Usually, multiples of
  41. C                                30 or 45 degrees are best.
  42. C                               CAMLOC(3) = angle between the viewer and the
  43. C                                XY plane located at Z=(ZMIN+ZMAX)/2 in
  44. C                                degrees.   Thus 90 degrees is directly above
  45. C                                the surface - an unexciting picture!   Usual
  46. C                                the angle is selected near 45 degrees.
  47. C
  48. C       XYLIM           * Type: real two dimensional array dimensioned (2,6).
  49. C                       * General parameters:
  50. C                               XYLIM(1,1) = XMIN ==> the minimum value of X.
  51. C                               XYLIM(2,1) = XMAX ==> the maximum value of X.
  52. C                               XYLIM(1,2) = YMIN ==> the minimum value of Y.
  53. C                               XYLIM(2,2) = YMAX ==> the maximum value of Y.
  54. C                                Note: Z(I,J) = F(Xi,Yj) where:
  55. C                                  Xi = XMIN + (i-1)*(XMAX-XMIN)/(KX-1)
  56. C                                  Yj = YMIN + (j-1)*(YMAX-YMIN)/(KY-1)
  57. C                               XYLIM(1,3) = ZMIN ==> the minimum value of Z.
  58. C                               XYLIM(2,3) = ZMAX ==> the maximum value of Z.
  59. C                                These Z values define the range of Z values
  60. C                                to fit on the screen.   It is strongly
  61. C                                advised that ZMIN and ZMAX bound Z(I,J).
  62. C                               XYLIM(1,4) = X/Z axis length ratio.   If this
  63. C                                parameter is 0, then X and Z are assumed to
  64. C                                have the same units, so their relative
  65. C                                lengths will be in proportion to their
  66. C                                ranges.   If this parameter is nonzero, then
  67. C                                the X axis will be XYLIM(1,4) times as long
  68. C                                as the Z axis.
  69. C                               XYLIM(2,4) = Y/Z axis length ratio.   Same as
  70. C                                XYLIM(1,4), but for Y axis.
  71. C                               XYLIM(1,5) = plot width in virtual coordinate
  72. C                               XYLIM(2,5) = plot height in virtual coord.
  73. C                                Note: The plot is expanded/contracted until
  74. C                                it all fits within the box defined by
  75. C                                XYLIM(1,5) and XYLIM(2,5).
  76. C                               XYLIM(1,6) = virtual X coord. of the lower
  77. C                                left corner of the plot box.
  78. C                               XYLIM(2,6) = virtual Y coord. of the lower
  79. C                                left corner of the box.
  80. C
  81. C       XLAB            * Type: string constant or variable.
  82. C                       * The X axis lable.
  83. C
  84. C       YLAB            * Type: string constant or variable.
  85. C                       * The Y axis lable.
  86. C
  87. C       ZLAB            * Type: string constant or variable.
  88. C                       * The Z axis lable.
  89. C
  90. C       CSIZE           * Type: real constant or variable.
  91. C                       * The character size in virtual coord. for the tick
  92. C                               mark lables and the axis lables.
  93. C
  94. C       MARPLT          * Type: integer constant or variable.
  95. C                       * Hidden line flag:
  96. C                         0 ==> draw all lines, hidden or not.
  97. C                         1 ==> suppress all lines hidden by the surface, but
  98. C                               display both the top and bottom of the surfac
  99. C                         3 ==> suppress all lines hidden by the surface, and
  100. C                               all lines showing the bottom of the surface.
  101. C                         Add 4 to MARPLT if you do not want the axes nor the
  102. C                          ticks labled.   This is useful on small plots.
  103. C
  104.         REAL*4 DEVID,XLENCM,YLENCM,XRES,YRES,NDCLRS,IDVBTS,NFLINE
  105.      1   XCLIPD,YCLIPD
  106.         COMMON /GCDCHR/ DEVID, XLENCM, YLENCM, XRES, YRES,
  107.      1   NDCLRS, IDVBTS, NFLINE, XCLIPD, YCLIPD
  108.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  109.         DIMENSION Z(IZDIM1,KY), CAMLOC(3), XYLIM(2,6)
  110.         CHARACTER*1 IZ(KX,KY), XLAB(2), YLAB(2), ZLAB(2)
  111.         EXTERNAL LEN
  112. C
  113. C       COMMON STORAGE DESCRIPTOR
  114.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  115.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  116.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  117.      3   AMTX(3,3),FOCALL
  118.         DIMENSION LIMIT(2),FLIM(2)
  119.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  120.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  121. C       END CDE
  122. C
  123.         LOGICAL LSOLID
  124.         COMMON /COMDP1/ LSOLID
  125. C
  126. C       LOCAL CDE
  127.         DIMENSION XMINA(2,2)
  128.         LOGICAL LLABLE
  129.         EQUIVALENCE(XMIN,XMINA(1,1))
  130.     COMMON/PLTPRM/CXSIZE,CYSIZE,TICKLN,YVINI
  131.         COMMON /DBASE/VX,VY,VOLDX,VOLDY
  132. C
  133. C       TYPE 800, KX,KY,XYLIM
  134. C800    FORMAT(' DEBUG-- KX,KY = ',2I5/' XYLIM='6(/1X,2E14.6))
  135. C       PICK UP XY LIMITS, BOX SIZES, ETC.
  136.         DO 9 J=1,2
  137.         XMINA(1,J) = XYLIM(1,J)
  138. 9       XMINA(2,J) = XYLIM(2,J)
  139.         ZMIN = XYLIM(1,3)
  140.         ZMAX = XYLIM(2,3)
  141.         PLOTX = XYLIM(1,5)
  142.         PLOTY = XYLIM(2,5)
  143.         AXISR(1) = XYLIM(1,4)
  144.         AXISR(2) = XYLIM(2,4)
  145.         PLTORG(1) = XYLIM(1,6)
  146.         PLTORG(2) = XYLIM(2,6)
  147. C
  148. C       NOW SET UP LIMITS IF AXIS RATIOS ARE REQUESTED
  149. C
  150.         IF (AXISR(1) .EQ. 0.0) GO TO 260
  151.         DO 255 I=1,2
  152. 255     XMINA(1,I)=AXISR(I)*ZMIN
  153. 260     IF (AXISR(2) .EQ. 0.0) GO TO 266
  154.         DO 265 I=1,2
  155. 265     XMINA(2,I)=AXISR(I)*ZMAX
  156. C       SET TOLERANCE FOR VISIBLE TESTS = HALF PLOTTER STEP SIZE
  157. 266     PQLMT = AMIN1(0.5/XRES,0.5/YRES)
  158. C       TYPE 504, XMINA
  159. C
  160. C       CONVERT R, PHI, THETA TO DX, DY, DZ
  161. C
  162.         RAD = 3.14159/180.0
  163.         PHI = CAMLOC(2)*RAD
  164.         THETA = CAMLOC(3)*RAD
  165.         CAMWKG(1)=CAMLOC(1)*COS(PHI)*COS(THETA)
  166.         CAMWKG(2)=CAMLOC(1)*SIN(PHI)*COS(THETA)
  167.         CAMWKG(3)=CAMLOC(1)*SIN(THETA)
  168. C       PICK UP CAMERA DATA
  169.         DO 3 J=1,2
  170.         CAMWKG(J+3)=(XMINA(1,J)+XMINA(2,J))/2.0
  171. 3       CAMWKG(J)=CAMWKG(J+3)+CAMWKG(J)
  172.         CAMWKG(6)=ZMIN+ZMAX/2.0
  173.         CAMWKG(3)=CAMWKG(6)+CAMWKG(3)
  174.         CALL CAMROT
  175.         MX=KX
  176.         FMX=FLOAT(KX)
  177.         NY=KY
  178.         FNY=FLOAT(NY)
  179. C       OPTION FOR SCALING Z
  180. C       SCALE FACTORS TO CONVERT USER VALUES TO INDICES
  181.         GX(1) = (XMAX-XMIN)/(FMX-1.0)
  182.         GX(2) = (YMAX-YMIN)/(FNY-1.0)
  183. C
  184. C       FIND Z SCALE FACTOR
  185. C
  186.         GX(3)=1.0
  187.         ZORG=0.0
  188. C
  189. C       FIND SCALE FACTORS FOR PLOT
  190. C
  191.         CYSIZE = CSIZE
  192.         CXSIZE = 9.0*CYSIZE/8.0
  193.         XA=1.0E30
  194.         XB=-1.0E30
  195.         YA=1.0E30
  196.         YB=-1.0E30
  197.         IF (CAMWKG(3) .LT. CAMWKG(6)) GO TO 16
  198.         DX=FLOAT(MX-1)/20.0
  199.         DY=FLOAT(NY-1)/20.0
  200.         IF=MX
  201.         XZ = XMAX
  202.         IB=1
  203.         JF=NY
  204.         YZ = YMIN
  205.         JB=1
  206.         IF (CAMWKG(1) .GE. CAMWKG(4)) GO TO 120
  207.         IF=1
  208.         XZ = XMIN
  209.         IB=MX
  210.         DX=-DX
  211. 120     IF (CAMWKG(2) .GE. CAMWKG(5)) GO TO 130
  212.         JF=1
  213.         YZ = YMAX
  214.         JB=NY
  215.         DY=-DY
  216. 130     FRX=IF
  217.         BKX=IB
  218.         FRY=JF
  219.         BKY=JB
  220.         VX = XMIN + (FRX-1.0)*GX(1) - CAMWKG(1)
  221.         VY = YMIN + (BKY-1.0-DY)*GX(2) - CAMWKG(2)
  222.         CALL EXTRMA(VX,VY,ZMAX-CAMWKG(3),XA,XB,YA,YB,IERR)
  223.         IF (IERR .NE. 0) GO TO 50
  224.         TEMP = ZMIN - CAMWKG(3)
  225.         CALL EXTRMA(VX,VY,TEMP,XA,XB,YA,YB,IERR)
  226.         IF (IERR .NE. 0) GO TO 50
  227.         VY = YMIN + (FRY-1.0+DY)*GX(2) - CAMWKG(2)
  228.         CALL EXTRMA(VX,VY,TEMP,XA,XB,YA,YB,IERR)
  229.         IF (IERR .NE. 0) GO TO 50
  230.         CALL EXTRMA(XMIN+(BKX-1.0)*GX(1)-CAMWKG(1),VY,TEMP,
  231.      1   XA,XB,YA,YB,IERR)
  232.         IF (IERR .NE. 0) GO TO 50
  233.         VX = VX + DX*GX(1)
  234.         CALL EXTRMA(VX,YMIN+(BKY-1.0)*GX(2)-CAMWKG(2),TEMP,
  235.      1   XA,XB,YA,YB,IERR)
  236.         IF (IERR .NE. 0) GO TO 50
  237.         CALL EXTRMA(VX,VY-DY*GX(2),TEMP,XA,XB,YA,YB,IERR)
  238.         IF (IERR .NE. 0) GO TO 50
  239. 16      DO 20 J=1,NY
  240.         VY = YMIN + (J-1)*GX(2) - CAMWKG(2)
  241.         DO 20 I=1,MX
  242.         VX = XMIN + (I-1)*GX(1) - CAMWKG(1)
  243.         CALL EXTRMA(VX,VY,Z(I,J)-CAMWKG(3),XA,XB,YA,YB,IERR)
  244.         IF (IERR .NE. 0) GO TO 50
  245. 20      CONTINUE
  246. C
  247. C       SCALE X AND Y RANGES TO FIT ON PLOT
  248. C
  249.         TEMP = 12.5*CXSIZE
  250.         LLABLE = .TRUE.
  251.         IF ((MARPLT.AND.4) .NE. 0) LLABLE = .FALSE.
  252.         IF (.NOT. LLABLE) TEMP = 0.0
  253.         FX(1) = (PLOTX-TEMP)/(XB-XA)
  254.         TEMP = 2.0*CYSIZE
  255.         IF (.NOT. LLABLE) TEMP = 0.0
  256.         FX(2) = (PLOTY-TEMP)/(YB-YA)
  257. C       CHOOSE MINIMUM FOCAL LENGTH OF THE TWO
  258.         FOCALL = AMIN1(FX(1),FX(2))
  259. C       SET X,Y ORIGINS (BEFORE SCALING TO FOCAL LGTH)
  260.         XORG(1) = XA
  261.         XORG(2) = YA
  262. C       SIZES IN X,Y (NOT INCLUDING OUT-OF-BOX POIINTS THAT GET IN PIC)
  263.         XB = (XB-XA)*FOCALL
  264.         YB = (YB-YA)*FOCALL
  265. C       CENTER FOR NOW, BUT LATER MAKE OPTIONAL
  266.         CENTER(1) = (PLOTX-XB)/2.0
  267.         CENTER(2) = (PLOTY-YB)/2.0
  268. C       TYPE 602, FX,XB,YB,PLOTX,PLOTY,CENTER
  269. C
  270. C       CAMERA LOCATION EXPRESSED AS XY INDICES
  271.         U = 1.0+(FMX-1.0)*(CAMWKG(1)-XMIN)/(XMAX-XMIN)
  272.         V = 1.0+(FNY-1.0)*(CAMWKG(2)-YMIN)/(YMAX-YMIN)
  273. C       FOR VISIBILITY CHECKING, SCALE CAMERA Z COORDINATE OPPOSITE TO THE
  274. C       WAY Z WILL BE SCALED FOR PLOTTING - RATHER THAN SCALING ALL THE
  275. C       Z-S ON THE SURFACE WHEN CHECKING.
  276.         W = (CAMWKG(3)-ZORG)/GX(3)
  277. C       CALCULATE VISIBILITIES
  278. C
  279. C       IF LSB OF MARPLT IS SET, SUPRESS ALL HIDDEN LINES
  280.         IF ((MARPLT.AND.1) .NE. 0) GO TO 7
  281.         DO 8 K = 1,NY
  282.         DO 8 J = 1,MX
  283. 8       IZ(J,K)=0
  284.         GO TO 40
  285. 7       LSOLID = .FALSE.
  286.         IF ((MARPLT.AND.2) .NE. 0) LSOLID = .TRUE.
  287.         DO 1 K = 1,NY
  288.         ETA = FLOAT(K)
  289.         DO 1 J =1,MX
  290.          L = IVIS(FLOAT(J),ETA,Z(J,K),Z,IZDIM1)+1
  291. 1       IZ(J,K)=L
  292. C
  293. C       NOW PLOT
  294. 40      CALL DRAW3D(Z,IZDIM1,IZ,KX)
  295.         IF (CAMWKG(3) .LT. CAMWKG(6)) GO TO 45
  296.         CALL GSSETC(CYSIZE,0.0)
  297.         CALL XYPRM(FRX,BKY,ZMAX,0)
  298.         VOLDX=VX
  299.         VOLDY=VY
  300.         VXT=VX
  301.         VYT=VY
  302.         CALL XYPRM(FRX,BKY-DY,ZMAX,1)
  303.         IF (LLABLE) CALL TICKL(ZMAX,-0.5)
  304.         CALL GSMOVE(VXT,VYT)
  305.         CALL XYPRM(FRX,BKY,ZMIN,1)
  306.         VOLDX=VX
  307.         VOLDY=VY
  308.         CALL XYPRM(FRX,BKY-DY,ZMIN,1)
  309.         IF (.NOT. LLABLE) GO TO 140
  310.         CALL TICKL(ZMIN,0.25)
  311.         TEMP = AMAX1(VOLDX,VXT)+1.5*CYSIZE
  312.         IF (VX .LT. VOLDX) TEMP = AMIN1(VOLDX,VXT)-0.5*CYSIZE
  313.         CALL GSMOVE(TEMP,(VOLDY+VYT-CXSIZE*LEN(ZLAB))/2.0)
  314.         CALL GSSETC(CYSIZE,90.0)
  315.         CALL GSPSTR(ZLAB)
  316.         CALL GSSETC(CYSIZE,0.0)
  317. 140     CALL GSMOVE(VOLDX,VOLDY)
  318.         CALL XYPRM(FRX+DX,BKY,ZMIN,1)
  319.         IF (LLABLE) CALL TICKL(XYLIM(1+JB/NY,2),-0.5)
  320.         CALL GSMOVE(VOLDX,VOLDY)
  321.         CALL XYPRM(FRX,FRY+DY,ZMIN,1)
  322.         IF (.NOT. LLABLE) GO TO 150
  323.         CALL TICKL(XYLIM(1+IF/MX,1),-0.5)
  324.         TEMP = CXSIZE*(LEN(YLAB)+0.25)
  325.         IF (VX .LT. VOLDX) TEMP = -0.25*CXSIZE
  326.         CALL GSMOVE((VX+VOLDX)/2.0-TEMP,(VY+VOLDY)/2.0-CYSIZE)
  327.         CALL GSPSTR(YLAB)
  328. 150     CALL XYPRM(FRX,FRY,Z(IF,JF),-1)
  329.         CALL GSMOVE(VX,VY)
  330.         CALL XYPRM(FRX,FRY,ZMIN,1)
  331.         VOLDX=VX
  332.         VOLDY=VY
  333.         CALL XYPRM(FRX+DX,FRY,ZMIN,1)
  334.         IF (LLABLE) CALL TICKL(XYLIM(1+JF/NY,2),-0.5)
  335.         CALL GSMOVE(VOLDX,VOLDY)
  336.         CALL XYPRM(BKX,FRY,ZMIN,1)
  337.         IF (.NOT. LLABLE) GO TO 160
  338.         TEMP = CXSIZE*(LEN(XLAB)+0.25)
  339.         IF (VX .GT. VOLDX) TEMP = -0.25*CXSIZE
  340.         CALL GSMOVE((VX+VOLDX)/2.0-TEMP,(VY+VOLDY)/2.0-CYSIZE)
  341.         CALL GSPSTR(XLAB)
  342. 160     VOLDX=VX
  343.         VOLDY=VY
  344.         CALL GSMOVE(VX,VY)
  345.         CALL XYPRM(BKX,FRY+DY,ZMIN,1)
  346.         IF (LLABLE) CALL TICKL(XYLIM(1+IB/MX,1),-0.5)
  347.         CALL GSMOVE(VOLDX,VOLDY)
  348.         CALL XYPRM(BKX,FRY,Z(IB,JF),1)
  349. 45      RETURN
  350. C
  351. C       POINT ON THE SURFACE IS BEHIND THE CAMERA. QUIT.
  352. C
  353.   50    WRITE(9,603)
  354.         RETURN
  355. C
  356. C       Z IS A FLAT PLANE, DO NOT DRAW (FOR NOW)
  357. C
  358.   60    WRITE(9,604)
  359.         RETURN
  360. C
  361. C 503 FORMAT(' Z MULTIPLIER',E15.6,', Z ORIGIN SHIFT',E15.6)
  362. C504    FORMAT('0X LIMITS',2F10.3/' Y LIMITS',2F10.3/' Z LIMITS',2F10.3/
  363. C       1' Z CUTOFF',2E15.6/
  364. C       2                 ' PLOT SIZE',2F10.3/' PLOT ORIGIN',2F10.3)
  365. C 602 FORMAT('0FOCAL LENGTHS TO FILL X,Y PLOTTER SPACE',2E15.6,
  366. C    1 ', LESSER VALUE CHOSEN'/'0PICTURE SIZE IN X,Y =',2F9.3,
  367. C    2 ', REQUESTED SIZES',2F9.3/' CENTERS = ',2G14.7)
  368.   603 FORMAT('PART OF SURFACE IS BEHIND THE CAMERA, UNABLE TO PLOT. SOR
  369.      1RY.')
  370.   604 FORMAT('FUNCTION IS LEVEL PLANE, NO USE PLOTTING IT')
  371.         END
  372.         SUBROUTINE EXTRMA(XV,YV,ZV,XA,XB,YA,YB,IERR)
  373. C
  374. C       COMMON STORAGE DESCRIPTOR
  375.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  376.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  377.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  378.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  379.      3   AMTX(3,3),FOCALL
  380.         DIMENSION LIMIT(2),FLIM(2)
  381.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  382.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  383. C       END CDE
  384. C
  385.         DIMENSION XS(3), XC(3)
  386. C
  387. C
  388.         XS(1) = XV
  389.         XS(2) = YV
  390.         XS(3) = ZV
  391.         CALL ROTATE(XS,AMTX,XC)
  392. C       QUIT IF POINT IS BEHIND CAMERA
  393.         IF(XC(3).LE.0.0) GO TO 50
  394.         XC(1) = XC(1)/XC(3)
  395.         XC(2) = XC(2)/XC(3)
  396.         XA = AMIN1(XA,XC(1))
  397.         XB = AMAX1(XB,XC(1))
  398.         YA = AMIN1(YA,XC(2))
  399.         YB = AMAX1(YB,XC(2))
  400.         IERR = 0
  401.         RETURN
  402. 50      IERR = -1
  403.         RETURN
  404.         END
  405.         SUBROUTINE XYPRM(X,Y,ZETA,ILINE)
  406. C
  407. C       COMMON STORAGE DESCRIPTOR
  408.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  409.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  410.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  411.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  412.      3   AMTX(3,3),FOCALL
  413.         DIMENSION LIMIT(2),FLIM(2)
  414.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  415.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  416. C       END CDE
  417. C
  418.     COMMON/PLTPRM/CXSIZE,CYSIZE,TICKLN,YVINI
  419.         COMMON /DBASE/VX,VY,VOLDX,VOLDY
  420.         DIMENSION XS(3),XC(3)
  421.         XS(1)=XMIN+(X-1.0)*GX(1)-CAMWKG(1)
  422.         XS(2)=YMIN+(Y-1.0)*GX(2)-CAMWKG(2)
  423.         XS(3)=ZORG + ZETA*GX(3)-CAMWKG(3)
  424.         CALL ROTATE(XS,AMTX,XC)
  425.         VX=(XC(1)/XC(3)-XORG(1))*FOCALL+PLTORG(1)+CENTER(1)
  426.         VY=(XC(2)/XC(3)-XORG(2))*FOCALL+PLTORG(2)+CENTER(2)
  427.         IF (ILINE) 30, 20, 10
  428. 10      CALL GSDRAW(VX,VY)
  429.         GO TO 30
  430. 20      CALL GSMOVE(VX,VY)
  431. 30      RETURN
  432.         END
  433.         SUBROUTINE TICKL(ANUM,UP)
  434.     COMMON/PLTPRM/CXSIZE,CYSIZE,TICKLN,YVINI
  435.         COMMON /DBASE/VX,VY,VOLDX,VOLDY
  436.         CHARACTER*8 TEMP1,TEMP2
  437.         CHARACTER*1 NUMBR(8)
  438.         EQUIVALENCE (NUMBR(1),TEMP2)
  439. C       ENCODE(6,10,NUMBR) ANUM
  440. C10     FORMAT(F6.2)
  441.         WRITE(TEMP2,451) ANUM
  442. C       READ(TEMP1,452) TEMP2
  443. 451     FORMAT(F6.2)
  444. 452     FORMAT(A8)
  445.         DO 20 IS=1,6
  446.         IF(NUMBR(IS) .NE. 32) GO TO 30
  447. 20      CONTINUE
  448. 30      NUMBR(7)=0
  449.         TEMP = CXSIZE*((7-IS)+0.25)
  450.         IF (VX .GT. VOLDX) TEMP = -0.25*CXSIZE
  451.         CALL GSMOVE(VX-TEMP,VY+UP*CYSIZE)
  452.         CALL GSPSTR(NUMBR(IS))
  453.         RETURN
  454.         END
  455.         SUBROUTINE CAMROT
  456. C       MAKE UP CAMERA ROTATION MATRIX
  457. C
  458. C       ROTATION IS DONE SO THAT Z PRIME AXIS IS DIRECTED FROM THE
  459. C       CAMERA TO THE AIMING POINT.   NOTE ALSO THAT THE PRIMED
  460. C       COORDINATE SYSTEM IS LEFT-HANDED IF EPSLON=-1.
  461. C       THIS IS SO THAT THE PICTURE COMES OUT RIGHT WHEN PROJECTED
  462. C       ON THE PRIMED COORDINATE SYSTEM.
  463. C
  464. C       COMMON STORAGE DESCRIPTOR
  465.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  466.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  467.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  468.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  469.      3   AMTX(3,3),FOCALL
  470.         DIMENSION LIMIT(2),FLIM(2)
  471.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  472.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  473. C       END CDE
  474. C
  475. C       LOCAL CDE
  476.         DIMENSION AU(3),AV(3),AW(3)
  477. C       HANDEDNESS PARAMETER, -1 FOR LEFT-HANDED USUALLY
  478.         DATA EPSLON/-1.0/
  479. C
  480.         S = 0.0
  481.         DO 1 J = 1,3
  482.         AV(J) = 0.0
  483.         AW(J) = 0.0
  484.         AU(J) = CAMWKG(J+3)-CAMWKG(J)
  485. 1       S = S + AU(J)**2
  486.         S = SQRT(S)
  487.         DO 2 J = 1,3
  488. 2       AU(J) = AU(J)/S
  489.         SIGMA = SQRT(AU(1)**2 + AU(2)**2)
  490. C       PREPARE LOOKING STRAIGHT UP OR DOWN
  491.         AV(1) = 1.0
  492.         AW(2) = -EPSLON
  493.         IF(AU(3) .GT. 0.0) AW(2) = -AW(2)
  494.         IF(SIGMA .LT. 1.0E-3) GO TO 4
  495. C       X AXIS
  496.         AV(1) = AU(2)/SIGMA
  497.         AV(2) = -AU(1)/SIGMA
  498.         AV(3) = 0.0
  499. C       Y AXIS
  500.         AW(1) = EPSLON*AU(1)*AU(3)/SIGMA
  501.         AW(2) = EPSLON*AU(2)*AU(3)/SIGMA
  502.         AW(3) = -EPSLON*SIGMA
  503. C       TRANSFER AXIS DIRECTION COSINES TO ROTATION MATRIX ROWS
  504. 4       DO 3 J = 1,3
  505.         AMTX(1,J) = AV(J)
  506.         AMTX(2,J) = AW(J)
  507. 3       AMTX(3,J) = AU(J)
  508.         RETURN
  509.         END
  510.         SUBROUTINE DRAWPQ(Z,IZDIM1)
  511. C       DRAW VISIBLE PART OF SEGMENT PC-QC
  512. C
  513.         DIMENSION Z(IZDIM1,2)
  514. C
  515. C       COMMON STORAGE DESCRIPTOR
  516.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  517.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  518.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  519.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  520.      3   AMTX(3,3),FOCALL
  521.         DIMENSION LIMIT(2),FLIM(2)
  522.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  523.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  524. C       END CDE
  525. C
  526. C
  527. C       CDE PACKAGE FOR DRAW3D,DRAWPQ
  528.         COMMON/COMDPA/PC(3),QC(3),P(3),Q(3),ENDA(6),ENDB(6),OLDQ(3),
  529.      1   PW(3),QW(3),T(6),PK(3),QK(3),PHIP,PHIQ,PHIA,IBEAM,ICOLOR
  530.         INTEGER PHIP,PHIQ,PHIA
  531. C       END OF CDE PACKAGE
  532. C
  533.         P(1) = PC(1)
  534.         P(2) = PC(2)
  535.         P(3) = PC(3)
  536.         Q(1) = QC(1)
  537.         Q(2) = QC(2)
  538.         Q(3) = QC(3)
  539. C       TEST IF P VISIBLE
  540. 2       IF(PHIP .EQ. 0) GO TO 30
  541. C       YES, TEST Q
  542. 7       IF(PHIP*PHIQ)10,4,3
  543. C       BOTH VISIBLE   SEGMENT DRAWABLE, PLOT   EXIT
  544. 3       KGOTO = 0
  545.         GO TO 300
  546. C       Q IS INVISIBLE, FIND LAST VISIBLE POINT ON SEGMENT PQ
  547. 4       JGOTO = 1
  548.         GO TO 200
  549. C       GIVE UP IF NOT FOUND IN MAXCUT1 BISECTIONS
  550. 5       IF(KFLAG .NE. 0) GO TO 6
  551. C       NEXT POINT
  552.         IBEAM = 0
  553.         RETURN
  554. C       POINT FOUND
  555. 6       Q(1) = ENDA(1)
  556.         Q(2) = ENDA(2)
  557.         Q(3) = ENDA(3)
  558.         GO TO 3
  559. C
  560. C       GAP IN SEGMENT, FIND LAST POINT TO CONNECT P.
  561. 10      JGOTO = 2
  562.         GO TO 200
  563. C       IF NOT FOUND (CANNOT FIND POINT WITH SAME VISIBILITY FN). TRY 2ND
  564. 11      IF(KFLAG .EQ. 0) GO TO 15
  565. C       SAVE OLD Q, RESET POINT   PLOT THIS PIECE.
  566.         OLDQ(1) = Q(1)
  567.         OLDQ(2) = Q(2)
  568.         OLDQ(3) = Q(3)
  569.         Q(1) = ENDA(1)
  570.         Q(2) = ENDA(2)
  571.         Q(3) = ENDA(3)
  572. C       DRAW FIRST PART OF SEGMENT AND COME BACK HERE
  573.         KGOTO = 2
  574.         GO TO 300
  575. C       RESTORE Q   FIND LOWER LIMIT OF UPPER SEGMENT.
  576. C       LIMITS FOR SEARCH
  577. 12      P(1) = Q(1)
  578.         P(2) = Q(2)
  579.         P(3) = Q(3)
  580.         Q(1) = OLDQ(1)
  581.         Q(2) = OLDQ(2)
  582.         Q(3) = OLDQ(3)
  583. C       BEAM OFF FIRST
  584. 15      IBEAM = 0
  585.         JGOTO = 3
  586.         GO TO 201
  587. C       IF SEGMENT TOO SHORT, GIVE UP.
  588. 13      IF(KFLAG .EQ. 0) GO TO 50
  589. C       LOWER END NOW NEWLY FOUND POINT.
  590. 14      P(1) = ENDA(1)
  591.         P(2) = ENDA(2)
  592.         P(3) = ENDA(3)
  593.         GO TO 3
  594. C       P INVISIBLE, CHECK Q. IF INVISIBLE, ADVANCE.
  595. 30      IBEAM = 0
  596.         IF(PHIQ .EQ. 0) GO TO 50
  597. C       FIND P
  598.         JGOTO = 4
  599.         GO TO 201
  600. C       IF NO POINT, GIVE UP.
  601. 31      IF(KFLAG) 14,50,14
  602. C
  603. C
  604. C       P VISIBLE, Q INVISIBLE, FIND Q.
  605. C       ENDB = INVISIBLE END OF INTERVAL, ENDA = VISIBLE
  606. 200     ENDB(1) = Q(1)
  607.         ENDB(2) = Q(2)
  608.         ENDB(3) = Q(3)
  609.         ENDA(1) = P(1)
  610.         ENDA(2) = P(2)
  611.         ENDA(3) = P(3)
  612. C       REQUIRED IVIS FUNCTION
  613. C       IN CASE OF GAP IN SEGMENT, CONSIDER POINT VISIBLE IF ITS VISIB.
  614. C       FUNCTION MATCHES THIS ONE AND UPDATE ENDA, ELSE ENDB.
  615.         PHIA = PHIP
  616.         GO TO 205
  617. C       P INVISIBLE, Q VISIBLE. FIND P.
  618. 201     ENDB(1) = P(1)
  619.         ENDB(2) = P(2)
  620.         ENDB(3) = P(3)
  621.         ENDA(1) = Q(1)
  622.         ENDA(2) = Q(2)
  623.         ENDA(3) = Q(3)
  624.         PHIA = PHIQ
  625. 205     KFLAG = 0
  626. C       GET PROJECTED LENGTH OF SEGMENT
  627.         PK(1) = XMIN + (ENDA(1)-1.0)*GX(1) - CAMWKG(1)
  628.         PK(2) = YMIN + (ENDA(2)-1.0)*GX(2) - CAMWKG(2)
  629.         PK(3) = ENDA(3)*GX(3) + ZORG - CAMWKG(3)
  630.         CALL ROTATE(PK,AMTX,ENDA(4))
  631.         PK(1) = XMIN + (ENDB(1)-1.0)*GX(1) - CAMWKG(1)
  632.         PK(2) = YMIN + (ENDB(2)-1.0)*GX(2) - CAMWKG(2)
  633.         PK(3) = ENDB(3)*GX(3) + ZORG - CAMWKG(3)
  634.         CALL ROTATE(PK,AMTX,ENDB(4))
  635. C       NEXT STEP
  636. 210     T(1) = (ENDA(1)+ENDB(1))/2.0
  637.         T(2) = (ENDA(2)+ENDB(2))/2.0
  638.         T(3) = (ENDA(3)+ENDB(3))/2.0
  639.         T(4) = (ENDA(4)+ENDB(4))/2.0
  640.         T(5) = (ENDA(5)+ENDB(5))/2.0
  641.         T(6) = (ENDA(6)+ENDB(6))/2.0
  642.         MFLAG = IVIS(T(1),T(2),T(3),Z,IZDIM1)
  643.         IF(MFLAG .EQ. PHIA) GO TO 220
  644. C       NOT VISIBLE, RESET INVISIBLE END.
  645.         ENDB(1) = T(1)
  646.         ENDB(2) = T(2)
  647.         ENDB(3) = T(3)
  648.         ENDB(4) = T(4)
  649.         ENDB(5) = T(5)
  650.         ENDB(6) = T(6)
  651. C       CHECK SEGMENT LENGTH (USE MAX OF X, Y DIFFERENCES)
  652. 216     SL = FOCALL*AMAX1(ABS(ENDA(4)/ENDA(6)-ENDB(4)/ENDB(6)),
  653.      1    ABS(ENDA(5)/ENDA(6)-ENDB(5)/ENDB(6)))
  654.         IF(SL .GE. PQLMT) GO TO 210
  655.         GO TO (5,11,13,31), JGOTO
  656. C       RECORD VISIBLE, UPDATE ENDA
  657. 220     KFLAG = MFLAG
  658.         ENDA(1) = T(1)
  659.         ENDA(2) = T(2)
  660.         ENDA(3) = T(3)
  661.         ENDA(4) = T(4)
  662.         ENDA(5) = T(5)
  663.         ENDA(6) = T(6)
  664.         GO TO 216
  665. C
  666. C
  667. C       DRAW P TO Q
  668. C
  669. C       IF BEAM IS ON, JUST MOVE IT TO Q.
  670. 300     IF(IBEAM .GT. 0) GO TO 310
  671. C       MOVE TO P, BEAM OFF.
  672.         PK(1) = XMIN + (P(1)-1.0)*GX(1) - CAMWKG(1)
  673.         PK(2) = YMIN + (P(2)-1.0)*GX(2) - CAMWKG(2)
  674.         PK(3) = P(3)*GX(3) + ZORG - CAMWKG(3)
  675.         CALL ROTATE(PK,AMTX,PW)
  676.         PW(1) = (PW(1)/PW(3)-XORG(1))*FOCALL + PLTORG(1) + CENTER(1)
  677.         PW(2) = (PW(2)/PW(3)-XORG(2))*FOCALL + PLTORG(2) + CENTER(2)
  678.         CALL GSMOVE(PW(1),PW(2))
  679. C       MOVE TO Q, BEAM ON. BEAM IS LEFT AND AT POINT Q.
  680. 310     QK(1) = XMIN + (Q(1)-1.0)*GX(1) - CAMWKG(1)
  681.         QK(2) = YMIN + (Q(2)-1.0)*GX(2) - CAMWKG(2)
  682.         QK(3) = Q(3)*GX(3) + ZORG - CAMWKG(3)
  683.         CALL ROTATE(QK,AMTX,QW)
  684.         QW(1) = (QW(1)/QW(3)-XORG(1))*FOCALL + PLTORG(1) + CENTER(1)
  685.         QW(2) = (QW(2)/QW(3)-XORG(2))*FOCALL + PLTORG(2) + CENTER(2)
  686.         CALL GSDRAW(QW(1),QW(2))
  687.         IBEAM = 1
  688.         IF(KGOTO .NE. 0) GO TO 12
  689. C
  690. 50      RETURN
  691.         END
  692.         SUBROUTINE DRAW3D(Z,IZDIM1,IZ,KX)
  693. C       DRAW PLOT
  694. C
  695.         DIMENSION Z(IZDIM1,2)
  696.         CHARACTER*1 IZ(KX,2)
  697. C
  698. C       COMMON STORAGE DESCRIPTOR
  699.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  700.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  701.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  702.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  703.      3   AMTX(3,3),FOCALL
  704.         DIMENSION LIMIT(2),FLIM(2)
  705.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  706.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  707. C       END CDE
  708. C
  709. C
  710. C       CDE PACKAGE FOR DRAW3D,DRAWPQ
  711.         INTEGER PHIP,PHIQ,PHIA
  712.         COMMON/COMDPA/PC(3),QC(3),P(3),Q(3),ENDA(6),ENDB(6),OLDQ(3),
  713.      1   PW(3),QW(3),T(6),PK(3),QK(3),PHIP,PHIQ,PHIA,IBEAM,ICOLOR
  714. C       END OF CDE PACKAGE
  715. C       END OF CDE PACKAGE
  716. C
  717. C       SAVE Z DIMENSION IN COMMON TO PASS ALONG THROUGH DRAWPQ TO IVIS
  718. C       SCAN ALONG X FIRST AT CONSTANT Y
  719. C
  720. C       INDEX OF COORDINATE BEING STEPPED ALONG A LINE
  721.         KSCAN = 1
  722. C       INDEX OF COORDINATE BEING HELD FIXED
  723.         KFIX = 2
  724. C       SET FIXED COORDINATE   INCREMENT
  725.         PC(KFIX) = 1.0
  726.         DELFIX = 1.0
  727. C       SET ROVING COORDINATE   INCREMENT INITIALLY
  728.         DELSCN = 1.0
  729.         QC(KSCAN) = 1.0
  730. C       BEGIN SCANNING A LINE
  731. 101     QC(KFIX) = PC(KFIX)
  732.         IBEAM = 0
  733. C       NEXT POINT IN LINE SCAN
  734. 102     PC(KSCAN) = QC(KSCAN)
  735.         QC(KSCAN) = PC(KSCAN) + DELSCN
  736. C       WORKING INDICES
  737.         JPC = IFIX(PC(1))
  738.         KPC = IFIX(PC(2))
  739.         JQC = IFIX(QC(1))
  740.         KQC = IFIX(QC(2))
  741. C       PHI FUNCTIONS
  742.         PC(3)=Z(JPC,KPC)
  743.         QC(3)=Z(JQC,KQC)
  744.         PHIP=ICHAR(IZ(JPC,KPC))-1
  745.         PHIQ=ICHAR(IZ(JQC,KQC))-1
  746. 200     CALL DRAWPQ(Z,IZDIM1)
  747. C       TEST IF LINE IS DONE
  748.         IF((QC(KSCAN)-1.0)*(QC(KSCAN)-FLIM(KSCAN)) .LT. 0.0) GO TO 102
  749. C       LINE DONE. ADVANCE FIXED COORDINATE.
  750.         PC(KFIX) = PC(KFIX) + DELFIX
  751. C       TEST IF FIXED COORDINATE NOW OFF LIMITS
  752.         IF((PC(KFIX)-1.0)*(PC(KFIX)-FLIM(KFIX)) .GT. 0.0) GO TO 55
  753. C       FLIP INCREMENT. SCAN BEGINS AT QC OF PREVIOUS LINE.
  754.         DELSCN = -DELSCN
  755.         GO TO 101
  756. C       TEST IF WE HAVE DONE Y SCAN YET.
  757. 55      IF(KSCAN .EQ. 2) RETURN
  758. C       NO, SCAN Y DIRECTION AT FIXED X.
  759.         KSCAN = 2
  760.         KFIX = 1
  761. C       START FIXED X AT X OF LAST TRAVERSE
  762.         PC(1) = QC(1)
  763. C       THEN STEP X IN OPPOSITE DIRECTION
  764.         DELFIX = -DELSCN
  765. C       WE ENDED UP AT MAX. Y, SO FIRST Y SCAN GOES BACKWARDS
  766.         DELSCN = -1.0
  767. C       INITIAL Y FOR FIRST LINE
  768.         QC(2) = FNY
  769.         GO TO 101
  770.         END
  771.         FUNCTION IVIS(XI,ETA,ZETA,Z,IZDIM1)
  772. C                  CORRECTED VERSION, 24FEB69
  773. C       DETERMINE IF POINT XI, ETA IS VISIBLE
  774. C       POINT IS GIVEN BY XI, ETA, ZIN
  775. C       AND VISIBILITY IS TESTED WITH RESPECT TO SURFACE Z(X,Y)
  776. C        XI, ETA COORDINATES EXPRESSED AS INDICES OF ARRAY Z, BUT NEED NOT
  777. C       BE INTEGERS IN GENERAL. FOR ENTRY IVIS, THEY MUST BE.
  778. C
  779.         DIMENSION Z(IZDIM1,2)
  780. C
  781. C       COMMON STORAGE DESCRIPTOR
  782.     COMMON/PLTCLP/XMIN,XMAX,YMIN,YMAX
  783.         COMMON/COMDP/ZMIN,ZMAX,AXISR(2),PLOTX,
  784.      1   PLOTY,PLTORG(2),CAMXYZ(3),MX,NY,FMX,FNY,CAMWKG(6),XORG(3),
  785.      2   GX(3),FX(2),KSCALE,ZORG,CENTER(2),PQLMT,
  786.      3   AMTX(3,3),FOCALL
  787.         DIMENSION LIMIT(2),FLIM(2)
  788.         EQUIVALENCE(U,CAMXYZ(1)),(V,CAMXYZ(2)),(W,CAMXYZ(3)),
  789.      1   (MX,LIMIT(1)),(FMX,FLIM(1))
  790. C       END CDE
  791. C
  792.         LOGICAL LSOLID
  793.         COMMON /COMDP1/ LSOLID
  794.         EQUIVALENCE (CX,CY), (DXI,DETA), (XIW,ETAW),
  795.      1  (XIEND,ETAEND), (KDXI,KDETA), (KXIEND,KETEND),
  796.      2  (DX,DY)
  797. C
  798. C
  799. C
  800. C       INITIAL P FUNCTION
  801. 5       IVIS = 0
  802.         R = U-XI
  803.         S = V-ETA
  804.         T = W-ZETA
  805. C       TEST IF WE CHECK ALONG X
  806.         IF(ABS(R) .LT. 1.0) GO TO 20
  807. C       CONSTANTS FOR Y(X),Z(X)
  808.         CY = S/R
  809.         CZ = T/R
  810.         DXI = SIGN(1.0,R)
  811. C       INITIAL POINT. TAKE AINT(XI) IF .NE. XI AND STEPS IN RIGHT DIRECTI
  812.         XIW = AINT(XI)
  813.         IF((XIW-XI)*DXI .LE. 0.0) XIW = XIW+DXI
  814. C       SKIP IF OFF LIMITS (WE ARE ON EDGE OF PLOT REGION)
  815.         IF((XIW-1.0)*(XIW-FMX) .GT. 0.0) GO TO 20
  816. C       FINAL POINT. TAKE AINT(U) IF IT MOVES OPPOSITE DXI, ELSE ROUND
  817.         XIEND = AINT(U)
  818.         IF((XIEND-U)*DXI .GE. 0.0) XIEND = XIEND-DXI
  819. C       BUT DO NOT GO BEYOND EDGES
  820.         XIEND = AMAX1(1.0,AMIN1(XIEND,FMX))
  821. C
  822. C                AFTER TESTING, RE-PRDER THESE STATEMENTS
  823.         J = IFIX(XIW)
  824.         KDXI = IFIX(DXI)
  825.         KXIEND = IFIX(XIEND)
  826.         XW = XIW-U
  827. C
  828. C       IF LIMITS CROSS, NO TEST
  829.         IF((XIEND-XIW)*DXI .LE. 0.0) GO TO 20
  830. C       GET Y(X)
  831. 3       YW = V + XW*CY
  832. C       IF Y IS OFF LIMITS, DONE
  833.         IF((YW-1.0)*(YW-FNY)) 21,25,20
  834. C       ON EDGE EXACTLY, NO INTERPOLATION
  835. 25      K = IFIX(YW)
  836.         IF(W + XW*CZ - Z(J,K)) 4,10,7
  837. C       INDEX FOR LOWER Y OF INTERVAL
  838. 21      K = IFIX(YW)
  839.         DY = YW-FLOAT(K)
  840. C       TEST Z OF LINE - Z OF SURFACE. ACCEPT ZERO DIFFERENCE.
  841.         IF((W + XW*CZ)-(Z(J,K) + DY*(Z(J,K+1)-Z(J,K)))) 4,10,7
  842. C       NEGATIVE. OK IF IVIS NEG. OR ZERO, ELSE REJECT
  843. 4       IF(IVIS) 10,6,40
  844. C       IVIS WAS ZERO, SET NEG.
  845. 6       IVIS = -1
  846.         GO TO 10
  847. C       PLUS. OK IF IVIS + OR ZERO, ELSE, REJECT
  848. 7       IF(IVIS) 40,8,10
  849. C       SET PLUS
  850. 8       IVIS = 1
  851. C       TEST IF DONE. ADVANCE IF NOT
  852. 10      IF(J .EQ. KXIEND) GO TO 20
  853.         J = J+KDXI
  854.         XW = XW+DXI
  855.         GO TO 3
  856. C
  857. C       CHECK IF WE TEST IN Y DIRECTION
  858. 20      IF(ABS(S) .LT. 1.0) GO TO 45
  859. C       CONSTANTS FOR X(Y),Z(Y)
  860.         CX = R/S
  861.         CZ = T/S
  862.         DETA = SIGN(1.0,S)
  863.         ETAW = AINT(ETA)
  864.         IF((ETAW-ETA)*DETA .LE. 0.0) ETAW = ETAW+DETA
  865. C       CHECK WHETHER ON LIMITS
  866.         IF((ETAW-1.0)*(ETAW-FNY) .GT. 0.0) GO TO 45
  867.         ETAEND = AINT(V)
  868.         IF((ETAEND-V)*DETA .GE. 0.0) ETAEND = ETAEND-DETA
  869.         ETAEND = AMAX1(1.0,AMIN1(FNY,ETAEND))
  870.         K = IFIX(ETAW)
  871.         KDETA = IFIX(DETA)
  872.         YW = ETAW-V
  873.         KETEND = IFIX(ETAEND)
  874. C       IF LIMITS CROSS, NO TEST, BUT TEST SINGLE POINT IF WE HAVE ALREADY
  875. C       TESTED X
  876.         A = ETAEND-ETAW
  877.         IF(A*DETA .LT. 0.0) GO TO 45
  878.         IF(A .EQ. 0.0 .AND. IVIS .EQ. 0) GO TO 45
  879. C       GET X(Y)
  880. 23      XW = U + YW*CX
  881. C       IF X OFF LIMITS, DONE
  882.         IF((XW-1.0)*(XW-FMX)) 44,46,45
  883. 46      J = IFIX(XW)
  884.         IF(W + YW*CZ - Z(J,K)) 24,30,27
  885. 44      J = IFIX(XW)
  886.         DX = XW-FLOAT(J)
  887.         IF((W + YW*CZ) - (Z(J,K)+DX*(Z(J+1,K)-Z(J,K)))) 24,30,27
  888. C       NEG., IVIS MUST BE NEG OR ZERO ELSE REJCT
  889. 24      IF(IVIS) 30,26,40
  890. C       SET IVIS NEG
  891. 26      IVIS = -1
  892.         GO TO 30
  893. C       POS, IVIS MUST BE ZERO OR + ELSE REJECT
  894. 27      IF(IVIS) 40,28,30
  895. 28      IVIS = 1
  896. C       TEST IF DONE, ADVANCE IF NOT.
  897. 30      IF(K .EQ. KETEND) GO TO 45
  898.         K = K+KDETA
  899.         YW = YW+DETA
  900.         GO TO 23
  901. C
  902. C       REJECT THIS POINT, RETURN ZERO.
  903. 40      IVIS = 0
  904.         RETURN
  905. C
  906. C       ACCEPT. RETURN +/- 1
  907. C       IF IVIS ZERO, CAMERA WAS RIGHT OVER XI, ETA.
  908. 45      IF(IVIS .EQ. 0) IVIS = ISIGN(1,IFIX(T))
  909.         IF (LSOLID .AND. (IVIS .EQ. -1)) GO TO 40
  910.         RETURN
  911.         END
  912.         SUBROUTINE ROTATE(XIN,A,XOUT)
  913. C       ROTATE VECTOR XIN BY MATRIX A TO GET XOUT
  914. C
  915. C
  916.         DIMENSION XIN(3),A(9),XOUT(3)
  917.         XOUT(1) = A(1)*XIN(1) + A(4)*XIN(2) + A(7)*XIN(3)
  918.         XOUT(2) = A(2)*XIN(1) + A(5)*XIN(2) + A(8)*XIN(3)
  919.         XOUT(3) = A(3)*XIN(1) + A(6)*XIN(2) + A(9)*XIN(3)
  920.         RETURN
  921.         END
  922.  
  923.