home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 3 / PDCD_3.iso / utilities / utilsp / pgplot / Examples / f77 / PGDemo7 < prev    next >
Encoding:
Text File  |  1994-02-24  |  7.7 KB  |  292 lines
  1.       PROGRAM PGDEM7
  2. C
  3. C Demonstration program for 3D surface plotting routine FREDDY.
  4. C
  5.       INTEGER PGBEG
  6.       REAL A(51,51), R, SIZE
  7.       INTEGER I, J
  8. C
  9.       DO 20 I=1,51
  10.       DO 10 J=1,51
  11.           R = (I-26)**2 + (J-26)**2
  12.           R = 0.5*SQRT(R)
  13.           A(I,J) = 1.0
  14.           IF (R.GT.0.0) A(I,J) = SIN(R)/R
  15.    10 CONTINUE
  16.    20 CONTINUE
  17.       IF (PGBEG(0, '?', 1, 1) .NE. 1) STOP
  18. C        FREDDY assumes the window is square of size SIZE.
  19.       SIZE = 1.0
  20.       CALL PGENV(0., SIZE, 0., SIZE, 1, -2)
  21.       CALL FREDDY(A,51,51, SIZE, 25.0)
  22.       CALL PGEND
  23.       END
  24. C-----------------------------------------------------------------------
  25.       SUBROUTINE FREDDY(ARRAY,KX,NY,SIZE,ANGLE)
  26.       INTEGER KX, NY
  27.       REAL ARRAY(KX,NY), SIZE, ANGLE
  28. C
  29. C Draws isometric plot of array
  30. C
  31.       REAL FMAX,FMIN,DELTAX,DELTAY,DELTAV,SINE,PEAK,X,DX,HEIGHT
  32.       INTEGER I,J,KI,KJ,NX,MX,MY,STEP,LEFT,RIGHT,IT,MN,INCX
  33.       LOGICAL VISBLE
  34.       COMMON /FREDCM/ DELTAX,X,STEP,LEFT,RIGHT,IT,NX,VISBLE
  35. C
  36.       MN = KX*NY
  37.       NX = KX
  38. C     Check array size:
  39.       IF(NX.LT.2 .OR. NY.LT.2) RETURN
  40.       FMAX = ARRAY(1,1)
  41.       FMIN = FMAX
  42.       DO 20 J=1,NY
  43.           DO 10 I=1,NX
  44.               FMIN = AMIN1(ARRAY(I,J),FMIN)
  45.               FMAX = AMAX1(ARRAY(I,J),FMAX)
  46.    10     CONTINUE
  47.    20 CONTINUE
  48.       DELTAX = SIZE/(NX+NY)
  49.       SINE = SIN(ANGLE/58.)
  50.       DELTAY = DELTAX*SINE
  51.       HEIGHT = SIZE*(1.-ABS(SINE))
  52.       DELTAV = HEIGHT
  53.       FMAX = FMAX-FMIN
  54.       IF(FMAX.LT.0.0001) FMAX = DELTAV
  55.       DELTAV = DELTAV/FMAX
  56.       MX = NX+1
  57.       MY = NY+1
  58.       STEP = MX
  59. C
  60. C Start PGPLOT buffering.
  61. C
  62.       CALL PGBBUF
  63. C
  64. C Work our way down the Y axis, then up the X axis,
  65. C calculating the Y plotter coordinates for each
  66. C column of the plot, doing the hidden-line suppression
  67. C at the same time.
  68. C
  69.       DO 50 J=1,NY
  70.           KJ = MY-J
  71.           KI = 1
  72. C               ( KI,KJ are coordinates of bottom of column)
  73.           ARRAY(KI,KJ) = DELTAY*(KI+KJ) + DELTAV*(ARRAY(KI,KJ)-FMIN)
  74.    30     PEAK = ARRAY(KI,KJ)
  75.    40     KI = KI+1
  76.           KJ = KJ+1
  77.           IF(KI.GT.NX .OR. KJ.GT.NY) GOTO 50
  78.           ARRAY(KI,KJ) = DELTAY*(KI+KJ) + DELTAV*(ARRAY(KI,KJ)-FMIN)
  79.           IF(ARRAY(KI,KJ).GT.PEAK) GOTO 30
  80.           IF(ARRAY(KI,KJ).LE.PEAK) ARRAY(KI,KJ) = -ABS(ARRAY(KI,KJ))
  81.           GOTO 40
  82.    50 CONTINUE
  83. C
  84. C Now to work our way up the X axis
  85. C
  86.       DO 80 I=2,NX
  87.           KI = I
  88.           KJ = 1
  89.           ARRAY(KI,KJ) = DELTAY*(KI+KJ)+DELTAV*(ARRAY(KI,KJ)-FMIN)
  90.    60     PEAK = ARRAY(KI,KJ)
  91.    70     KI = KI+1
  92.           KJ = KJ+1
  93.           IF(KI.GT.NX .OR. KJ.GT.NY) GOTO 80
  94.           ARRAY(KI,KJ) = DELTAY*(KI+KJ)+DELTAV*(ARRAY(KI,KJ)-FMIN)
  95.           IF(ARRAY(KI,KJ).GT.PEAK) GOTO 60
  96.           IF(ARRAY(KI,KJ).LE.PEAK) ARRAY(KI,KJ) = -ABS(ARRAY(KI,KJ))
  97.           GOTO 70
  98.    80 CONTINUE
  99. C
  100. C Draw a line along the bottom of the vertical faces
  101. C
  102.       CALL PGMOVE(DELTAX*(NX+NY-2), DELTAY*(MX))
  103.       CALL PGDRAW(DELTAX*(NY-1),    DELTAY*2)
  104.       CALL PGDRAW(0.0,              DELTAY*MY)
  105. C
  106. C Array is now ready for plotting.  If a point is
  107. C positive, then it is to be plotted at that Y
  108. C coordinate; if it is negative, then it is
  109. C invisible, but at minus that Y coordinate (the point
  110. C where the line heading towards it disappears has to
  111. C be determined by finding the intersection of it and
  112. C the cresting line).
  113. C
  114. C Plot rows:
  115. C
  116.       DO 110 J=1,NY,2
  117.           KJ = MY-J
  118.           DX = DELTAX*(J-2)
  119.           X = DX+DELTAX
  120.           CALL PGMOVE(X,DELTAY*(KJ+1))
  121.           CALL PGDRAW(X,ARRAY(1,KJ))
  122.           VISBLE = .TRUE.
  123.           DO 90 I=2,NX
  124.               RIGHT = I+NX*(KJ-1)
  125.               LEFT = RIGHT-1
  126.               IT = RIGHT
  127.               X = DX+DELTAX*I
  128.               CALL FREDGO(ARRAY,MN)
  129.    90     CONTINUE
  130. C
  131. C Now at far end of row so come back
  132. C
  133.           KJ = KJ-1
  134.           IF(KJ.LE.0) GOTO 170
  135.           VISBLE = ARRAY(NX,KJ).GE.0.0
  136.           DX = DELTAX*(NX+J)
  137.           IF(VISBLE) CALL PGMOVE(DX-DELTAX,ARRAY(NX,KJ))
  138.           DELTAX = -DELTAX
  139.           DO 100 I=2,NX
  140.               KI = MX-I
  141.               LEFT = KI+NX*(KJ-1)
  142.               RIGHT = LEFT+1
  143.               IT = LEFT
  144.               X = DX+DELTAX*I
  145.               CALL FREDGO(ARRAY,MN)
  146.   100     CONTINUE
  147. C
  148.           X = DX+DELTAX*NX
  149.           IF(.NOT.VISBLE) CALL PGMOVE(X,ARRAY(1,KJ))
  150.           CALL PGDRAW(X,DELTAY*(KJ+1))
  151. C               (set DELTAX positive for return trip)
  152.           DELTAX = -DELTAX
  153.   110 CONTINUE
  154. C
  155. C Now do the columns:
  156. C as we fell out of the last DO-loop we do the
  157. C columns in ascending-X order
  158. C
  159.       INCX = 1
  160.       KI = 1
  161. C               (set DELTAX -ve since scanning R to L)
  162.   120 DX = DELTAX*(KI+NY-1)
  163.       DELTAX = -DELTAX
  164.       X = DX+DELTAX
  165.       CALL PGMOVE(X,ARRAY(1,1))
  166.   130 VISBLE = .TRUE.
  167.       DO 140 J=2,NY
  168.           LEFT = KI+NX*(J-1)
  169.           RIGHT = LEFT-NX
  170.           IT = LEFT
  171.           X = DX+DELTAX*J
  172.           CALL FREDGO(ARRAY,MN)
  173.   140 CONTINUE
  174. C
  175. C At far end, increment X and check still inside array
  176. C
  177.       KI = KI+INCX
  178.       IF(KI.LE.0 .OR. KI.GT.NX) GOTO 180
  179.       VISBLE = ARRAY(KI,NY).GE.0.0
  180.       DELTAX = -DELTAX
  181.       DX = DELTAX*(KI-2)
  182.       X = DX+DELTAX
  183.       IF(VISBLE) CALL PGMOVE(X,ARRAY(KI,NY))
  184.       DO 150 J=2,NY
  185.           KJ = MY-J
  186.           RIGHT = KI+NX*(KJ-1)
  187.           LEFT = RIGHT+NX
  188.           IT = RIGHT
  189.           X = DX+DELTAX*J
  190.           CALL FREDGO(ARRAY,MN)
  191.   150 CONTINUE
  192. C
  193.       X = DX+DELTAX*NY
  194.       IF(.NOT.VISBLE) CALL PGMOVE(X,ARRAY(KI,1))
  195.       IF(KI.EQ.1) GOTO 180
  196.       CALL PGDRAW(X,DELTAY*(KI+1))
  197.       KI = KI+INCX
  198.       IF(KI.GT.NX) GOTO 180
  199.       IF(KI.EQ.1) GOTO 120
  200.   160 DELTAX = -DELTAX
  201.       DX = DELTAX*(1-KI-NY)
  202.       X = DX+DELTAX
  203.       CALL PGMOVE(X,DELTAY*(KI+1))
  204.       CALL PGDRAW(X,ARRAY(KI,1))
  205.       GOTO 130
  206. C
  207. C Do columns backwards because ended rows at far end of X
  208. C
  209.   170 KI = NX
  210.       INCX = -1
  211.       DX = DELTAX*(KI+NY)
  212.       GOTO 160
  213. C
  214. C
  215.   180 CALL PGEBUF
  216.       END
  217. C-----------------------------------------------------------------------
  218.       SUBROUTINE FREDGO(ARRAY,MN)
  219.       INTEGER MN
  220.       REAL ARRAY(MN)
  221. C
  222.       INTEGER STEP,LEFT,RIGHT,IT,NX
  223.       LOGICAL VISBLE
  224.       REAL AL,AR,BL,EM,XX,X,Y,DELTAX
  225.       COMMON /FREDCM/ DELTAX,X,STEP,LEFT,RIGHT,IT,NX,VISBLE
  226. C
  227. C Test visibility
  228. C
  229.       IF(ARRAY(IT).LT.0.0) GOTO 80
  230. C
  231. C This point is visible - was last?
  232. C
  233.       IF(VISBLE) GOTO 50
  234. C
  235. C No: calculate point where this line vanishes
  236. C
  237.    10 IF(LEFT.LE.NX .OR. MOD(LEFT-1,NX).EQ.0 .OR.
  238.      1     RIGHT.LE.NX .OR. MOD(RIGHT-1,NX).EQ.0) GOTO 100
  239.       AL = ABS(ARRAY(LEFT))
  240.       AR = ABS(ARRAY(RIGHT))
  241.       IF(ARRAY(LEFT).LT.0.0) GOTO 70
  242. C               Right-hand point is crested
  243.    20 RIGHT = RIGHT-STEP
  244.       IF(ARRAY(RIGHT).LT.0.0) GOTO 20
  245. C               Left-hand end of cresting line is either
  246. C               RIGHT+NX or RIGHT-1
  247.       LEFT = RIGHT+NX
  248.       IF(ARRAY(LEFT).LT.0.0) LEFT = RIGHT-1
  249. C
  250. C               RIGHT and LEFT index into the endpoints of the
  251. C               cresting line
  252.    30 BL = ABS(ARRAY(LEFT))
  253.       EM = ABS(ARRAY(RIGHT))-BL
  254.       XX = EM-AR+AL
  255.       IF(ABS(XX).LT.0.0001) GOTO 60
  256.       XX = (AL-BL)/XX
  257.    40 Y = EM*XX+BL
  258.       IF(DELTAX.GT.0.0) XX = 1.0-XX
  259.       XX = X-XX*DELTAX
  260.       IF(VISBLE) GOTO 90
  261. C               Drawing a line from an invisible point
  262. C               to a visible one
  263.       CALL PGMOVE(XX,Y)
  264.       VISBLE = .TRUE.
  265.    50 CALL PGDRAW(X,ARRAY(IT))
  266.       RETURN
  267. C
  268.    60 XX = 0.5
  269.       GOTO 40
  270. C
  271. C Left-hand point crested
  272. C
  273.    70 LEFT = LEFT-STEP
  274.       IF(ARRAY(LEFT).LT.0.0) GOTO 70
  275. C
  276. C Right-hand end of cresting line is either LEFT+1 or LEFT-NX
  277. C
  278.       RIGHT = LEFT+1
  279.       IF(ARRAY(RIGHT).LT.0.0) RIGHT = LEFT-NX
  280.       GOTO 30
  281. C
  282. C This point is invisible; if last one was too, then forget it;
  283. C else draw a line towards it
  284. C
  285.    80 IF(.NOT.VISBLE) RETURN
  286.       GOTO 10
  287. C
  288.    90 CALL PGDRAW(XX,Y)
  289.   100 VISBLE = .FALSE.
  290.       RETURN
  291.       END
  292.