home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / fortran / pgplot5.1 / pgplot5 / pgplot5.1.0 / examples-src / pgdemo4.f < prev    next >
Encoding:
Text File  |  1994-11-29  |  8.6 KB  |  275 lines
  1.       PROGRAM PGDEM4
  2. C-----------------------------------------------------------------------
  3. C Test program for PGPLOT: test of imaging routine PGGRAY (with 
  4. C PGCONT). This program serves to both demonstrate and test PGGRAY.
  5. C It computes an (arbitrary) function on a 2D array, and uses both
  6. C PGGRAY and PGCONT to display it. An irregular transformation (TR 
  7. C matrix) is used, to test (a) that the routine works when the array
  8. C pixels are not aligned with the device pixels, and (b) that the
  9. C image is clipped correctly at the edge of the viewport. The program
  10. C also draws the bounding quadrilateral of the contour map. The contours
  11. C should end on this quadrilateral, but note that the grayscale image 
  12. C extends one half pixel beyond this quadrilateral (if the subarray to 
  13. C be displayed contains N pixels in one dimension, the width of the
  14. C image is N units, while the width of the contour map is N-1 pixels).
  15. C-----------------------------------------------------------------------
  16.       INTEGER PGBEG
  17.       INTEGER   MXI, MXJ
  18.       PARAMETER (MXI=40, MXJ=40)
  19.       INTEGER I,J
  20.       REAL F(MXI,MXJ)
  21.       REAL ANGLE,FMIN,FMAX,ALEV,RADIUS,TR(6)
  22. C
  23. C Open device for graphics.
  24. C
  25.       IF (PGBEG(0,'?',1,1) .NE. 1) STOP
  26. C
  27. C Compute a suitable function.
  28. C
  29.       FMIN = F(1,1)
  30.       FMAX = F(1,1)
  31.       DO 20 I=1,MXI
  32.           DO 10 J=1,MXJ
  33.               F(I,J) = COS(0.6*SQRT(I*80./MXI)-16.0*J/(3.*MXJ))*
  34.      :                 COS(16.0*I/(3.*MXI))+(I/FLOAT(MXI)-J/FLOAT(MXJ))
  35.               FMIN = MIN(F(I,J),FMIN)
  36.               FMAX = MAX(F(I,J),FMAX)
  37.    10     CONTINUE
  38.    20 CONTINUE
  39. C
  40. C Compute the coordinate transformation matrix.
  41. C ANGLE is the image rotation angle in radians.
  42. C
  43.       ANGLE=-120./57.29578
  44.       RADIUS= 40.*SQRT(2.)
  45.       TR(1) =-(40./MXI)*COS(ANGLE)+(40./MXJ)*SIN(ANGLE)-
  46.      :         RADIUS*COS(ANGLE+45/57.29578)
  47.       TR(2) = 80.*COS(ANGLE)/MXI
  48.       TR(3) =-80.*SIN(ANGLE)/MXJ
  49.       TR(4) =-(40./MXI)*SIN(ANGLE)-(40./MXJ)*COS(ANGLE)-
  50.      :         RADIUS*SIN(ANGLE+45/57.29578)
  51.       TR(5) = 80.*SIN(ANGLE)/MXI
  52.       TR(6) = 80.*COS(ANGLE)/MXJ
  53. C-----------------------------------------------------------------------
  54. C PGGRAY
  55. C-----------------------------------------------------------------------
  56. C Clear the screen. Set up window and viewport.
  57. C
  58.       CALL PGPAGE
  59.       CALL SETVP
  60.       CALL PGWNAD(-40., 40.,-40., 40.)
  61.       CALL PGSCI(1)
  62. C
  63. C Draw the map with PGGRAY.  
  64. C
  65.       CALL PGGRAY(F,MXI,MXJ,1,MXI,1,MXJ,FMAX,FMIN,TR)
  66. C
  67. C Overlay contours in red.
  68. C
  69.       CALL PGSCI(2)
  70.       DO 30 I=1,21
  71.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  72.           IF (MOD(I,5).EQ.0) THEN
  73.               CALL PGSLW(3)
  74.           ELSE
  75.               CALL PGSLW(1)
  76.           END IF
  77.           IF (I.LT.10) THEN
  78.               CALL PGSLS(2)
  79.           ELSE
  80.               CALL PGSLS(1)
  81.           END IF
  82.           CALL PGCONT(F,MXI,MXJ,1,MXI,1,MXJ,ALEV,-1,TR)
  83.    30 CONTINUE
  84.       CALL PGSLS(1)
  85.       CALL PGSLW(1)
  86. C
  87. C Annotate the plot.
  88. C
  89.       CALL PGSCI(2)
  90.       CALL OUTLIN(1,MXI,1,MXJ,TR)
  91.       CALL PGSCI(5)
  92.       CALL PGMTXT('t',1.0,0.0,0.0,'Routines PGGRAY, PGCONT, PGWEDG')
  93.       CALL PGBOX('bcnts',0.0,0,'bcnts',0.0,0)
  94. C
  95. C Draw a wedge.
  96. C
  97.       CALL PGSCH(0.8)
  98.       CALL PGWEDG('BG', 3.0, 4.0, FMAX, FMIN, 'Elevation')
  99.       CALL PGSCH(1.0)
  100. C-----------------------------------------------------------------------
  101. C PGIMAG
  102. C-----------------------------------------------------------------------
  103. C Clear the screen. Set up window and viewport.
  104. C
  105.       CALL PGPAGE
  106.       CALL SETVP
  107.       CALL PGWNAD(-40., 40.,-40., 40.)
  108.       CALL PGSCI(1)
  109. C
  110. C Set up the color map.
  111. C
  112.       CALL PALETT(2)
  113. C
  114. C Draw the map with PGIMAG.  
  115. C
  116.       CALL PGIMAG(F,MXI,MXJ,1,MXI,1,MXJ,FMIN,FMAX,TR)
  117. C
  118. C Overlay contours in white.
  119. C
  120.       CALL PGSCI(1)
  121.       DO 40 I=1,21
  122.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  123.           IF (MOD(I,5).EQ.0) THEN
  124.               CALL PGSLW(3)
  125.           ELSE
  126.               CALL PGSLW(1)
  127.           END IF
  128.           IF (I.LT.10) THEN
  129.               CALL PGSLS(2)
  130.           ELSE
  131.               CALL PGSLS(1)
  132.           END IF
  133.           CALL PGCONT(F,MXI,MXJ,1,MXI,1,MXJ,ALEV,-1,TR)
  134.    40 CONTINUE
  135.       CALL PGSLS(1)
  136.       CALL PGSLW(1)
  137. C
  138. C Annotate the plot.
  139. C
  140.       CALL PGSCI(2)
  141.       CALL OUTLIN(1,MXI,1,MXJ,TR)
  142.       CALL PGSCI(5)
  143.       CALL PGMTXT('t',1.0,0.0,0.0,'Routines PGIMAG, PGCONT, PGWEDG')
  144.       CALL PGBOX('bcnts',0.0,0,'bcnts',0.0,0)
  145. C
  146. C Draw a wedge.
  147. C
  148.       CALL PGSCH(0.8)
  149.       CALL PGWEDG('BI', 3.0, 4.0, FMIN, FMAX, 'Elevation')
  150.       CALL PGSCH(1.0)
  151. C-----------------------------------------------------------------------
  152. C Close the device and exit.
  153. C
  154.       CALL PGEND
  155.       END
  156.  
  157.       SUBROUTINE OUTLIN(I1,I2,J1,J2,TR)
  158.       INTEGER I1,I2,J1,J2
  159.       REAL TR(6)
  160. C-----------------------------------------------------------------------
  161. C Draw the enclosing rectangle of the subarray to be contoured,
  162. C applying the transformation TR.
  163. C
  164. C For a contour map, the corners are (I1,J1) and (I2,J2); for
  165. C a gray-scale map, they are (I1-0.5,J1-0.5), (I2+0.5, J2+0.5).
  166. C-----------------------------------------------------------------------
  167.       INTEGER K
  168.       REAL XW(5), YW(5), T
  169. C
  170.       XW(1) = I1
  171.       YW(1) = J1
  172.       XW(2) = I1
  173.       YW(2) = J2
  174.       XW(3) = I2
  175.       YW(3) = J2
  176.       XW(4) = I2
  177.       YW(4) = J1
  178.       XW(5) = I1
  179.       YW(5) = J1
  180.       DO 10 K=1,5
  181.           T = XW(K)
  182.           XW(K) = TR(1) + TR(2)*T + TR(3)*YW(K)
  183.           YW(K) = TR(4) + TR(5)*T + TR(6)*YW(K)
  184.    10 CONTINUE
  185.       CALL PGLINE(5,XW,YW)
  186.       END
  187.  
  188.       SUBROUTINE PALETT(TYPE)
  189. C-----------------------------------------------------------------------
  190. C Set a "palette" of colors in the range of color indices used by
  191. C PGIMAG.
  192. C-----------------------------------------------------------------------
  193.       INTEGER TYPE
  194. C
  195.       REAL GL(2), GR(2), GG(2), GB(2)
  196.       REAL RL(9), RR(9), RG(9), RB(9)
  197.       REAL HL(5), HR(5), HG(5), HB(5)
  198.       REAL WL(10), WR(10), WG(10), WB(10)
  199.       REAL AL(20), AR(20), AG(20), AB(20)
  200.       REAL TL(4), TR(4), TG(4), TB(4)
  201. C
  202.       DATA GL /0.0, 1.0/
  203.       DATA GR /0.0, 1.0/
  204.       DATA GG /0.0, 1.0/
  205.       DATA GB /0.0, 1.0/
  206. C
  207.       DATA RL /-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7/
  208.       DATA RR / 0.0, 0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0, 1.0/
  209.       DATA RG / 0.0, 0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0, 1.0/
  210.       DATA RB / 0.0, 0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0, 1.0/
  211. C
  212.       DATA HL /0.0, 0.2, 0.4, 0.6, 1.0/
  213.       DATA HR /0.0, 0.5, 1.0, 1.0, 1.0/
  214.       DATA HG /0.0, 0.0, 0.5, 1.0, 1.0/
  215.       DATA HB /0.0, 0.0, 0.0, 0.3, 1.0/
  216. C
  217.       DATA WL /0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0/
  218.       DATA WR /0.0, 1.0, 0.0, 0.0, 0.3,  0.8,  0.3,  1.0,  1.0, 1.0/
  219.       DATA WG /0.0, 0.5, 0.4, 1.0, 0.0,  0.0,  0.2,  0.7,  1.0, 1.0/
  220.       DATA WB /0.0, 0.0, 0.0, 0.0, 0.4,  1.0,  0.0,  0.0, 0.95, 1.0/
  221. C
  222.       DATA AL /0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5,
  223.      :         0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0/
  224.       DATA AR /0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0,
  225.      :         0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/
  226.       DATA AG /0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8,
  227.      :         0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0/
  228.       DATA AB /0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9,
  229.      :         0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
  230. C
  231.       DATA TL /0.0, 0.5, 0.5, 1.0/
  232.       DATA TR /0.2, 0.6, 0.6, 1.0/
  233.       DATA TG /0.0, 0.0, 0.5, 1.0/
  234.       DATA TB /1.0, 0.0, 0.0, 0.0/
  235. C
  236.       IF (TYPE.EQ.1) THEN
  237. C        -- gray scale
  238.          CALL PGCTAB(GL, GR, GG, GB, 2, 1.0, 0.5)
  239.       ELSE IF (TYPE.EQ.2) THEN
  240. C        -- rainbow
  241.          CALL PGCTAB(RL, RR, RG, RB, 9, 1.0, 0.5)
  242.       ELSE IF (TYPE.EQ.3) THEN
  243. C        -- heat
  244.          CALL PGCTAB(HL, HR, HG, HB, 5, 1.0, 0.5)
  245.       ELSE IF (TYPE.EQ.4) THEN
  246. C        -- weird IRAF
  247.          CALL PGCTAB(WL, WR, WG, WB, 10, 1.0, 0.5)
  248.       ELSE IF (TYPE.EQ.5) THEN
  249. C        -- AIPS
  250.          CALL PGCTAB(AL, AR, AG, AB, 20, 1.0, 0.5)
  251.       ELSE IF (TYPE.EQ.6) THEN
  252. C        -- TJP
  253.          CALL PGCTAB(TL, TR, TG, TB, 4, 1.0, 0.5)
  254.       END IF
  255.       END
  256.  
  257.       SUBROUTINE SETVP
  258. C-----------------------------------------------------------------------
  259. C Set the viewport, allowing margins around the edge for annotation.
  260. C (This is similar in effect to PGVSTD, but has different margins.)
  261. C The routine determines the view-surface size and allocates margins
  262. C as fractions of the minimum of width and height.
  263. C-----------------------------------------------------------------------
  264.       REAL D, VPX1, VPX2, VPY1, VPY2
  265. C
  266.       CALL PGSVP(0.0, 1.0, 0.0, 1.0)
  267.       CALL PGQVP(1, VPX1, VPX2, VPY1, VPY2)
  268.       D = MIN(VPX2-VPX1, VPY2-VPY1)/40.0
  269.       VPX1 = VPX1 + 5.0*D
  270.       VPX2 = VPX2 - 2.0*D
  271.       VPY1 = VPY1 + 8.0*D
  272.       VPY2 = VPY2 - 2.0*D
  273.       CALL PGVSIZ(VPX1, VPX2, VPY1, VPY2)
  274.       END
  275.