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