home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / maths / pgplot_1 / Examples_f77_PGDEMO3 < prev    next >
Encoding:
Text File  |  1997-06-10  |  18.0 KB  |  602 lines
  1.       PROGRAM PGDEM3
  2. C-----------------------------------------------------------------------
  3. C Demonstration program for PGPLOT contouring routines.
  4. C-----------------------------------------------------------------------
  5.       INTEGER PGBEG
  6.       WRITE (*,'(A)') ' Demonstration of PGPLOT contouring routines'
  7. C
  8. C Call PGBEG to initiate PGPLOT and open the output device; PGBEG
  9. C will prompt the user to supply the device name and type.
  10. C
  11.       IF (PGBEG(0,'?',1,1) .NE. 1) STOP
  12. C
  13. C Call the demonstration subroutines.
  14. C
  15.       WRITE (*,'(A)') ' Routine PGCONT'
  16.       CALL PGEX31
  17.       WRITE (*,'(A)') ' Routine PGCONS'
  18.       CALL PGEX32
  19.       WRITE (*,'(A)') ' Routine PGCONT with PGCONL labels'
  20.       CALL PGEX36
  21.       WRITE (*,'(A)') ' Routine PGCONB'
  22.       CALL PGEX33
  23.       WRITE (*,'(A)') ' Routine PGCONX with arrow labels'
  24.       CALL PGEX37
  25.       WRITE (*,'(A)') ' Routine PGCONX'
  26.       CALL PGEX34
  27.       WRITE (*,'(A)') ' Routine PGCONF'
  28.       CALL PGEXX1
  29. C
  30. C Finally, call PGEND to terminate things properly.
  31. C
  32.       CALL PGEND
  33. C-----------------------------------------------------------------------
  34.       END
  35.  
  36.       SUBROUTINE PGEX31
  37. C-----------------------------------------------------------------------
  38. C Demonstration of contouring routine PGCONT.
  39. C-----------------------------------------------------------------------
  40.       INTEGER I,J
  41.       REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6)
  42.       DATA TR/0.,1.,0.,0.,0.,1./
  43. C
  44. C Compute a suitable function.
  45. C
  46.       FMIN = 0.0
  47.       FMAX = 0.0
  48.       DO 20 I=1,40
  49.           DO 10 J=1,40
  50.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  51.      1                     (I-J)/40.0
  52.               FMIN = MIN(F(I,J),FMIN)
  53.               FMAX = MAX(F(I,J),FMAX)
  54.    10     CONTINUE
  55.    20 CONTINUE
  56. C
  57. C Clear the screen. Set up window and viewport.
  58. C
  59.       CALL PGPAGE
  60.       CALL PGSVP(0.05,0.95,0.05,0.95)
  61.       CALL PGSWIN(1.0,40.0,1.0,40.0)
  62.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  63.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONT')
  64. C
  65. C Draw the map.  PGCONT is called once for each contour, using
  66. C different line attributes to distinguish contour levels.
  67. C
  68.       CALL PGBBUF
  69.       DO 30 I=1,21
  70.           ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
  71.           IF (MOD(I,5).EQ.0) THEN
  72.               CALL PGSLW(3)
  73.           ELSE
  74.               CALL PGSLW(1)
  75.           END IF
  76.           IF (I.LT.10) THEN
  77.               CALL PGSCI(2)
  78.               CALL PGSLS(2)
  79.           ELSE
  80.               CALL PGSCI(3)
  81.               CALL PGSLS(1)
  82.           END IF
  83.           CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
  84.    30 CONTINUE
  85.       CALL PGSLW(1)
  86.       CALL PGSLS(1)
  87.       CALL PGSCI(1)
  88.       CALL PGEBUF
  89.       END
  90.  
  91.       SUBROUTINE PGEX32
  92. C-----------------------------------------------------------------------
  93. C Demonstration of contouring routine PGCONS.
  94. C-----------------------------------------------------------------------
  95.       INTEGER I,J
  96.       REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6)
  97.       DATA TR/0.,1.,0.,0.,0.,1./
  98. C
  99. C Compute a suitable function.
  100. C
  101.       FMIN = 0.0
  102.       FMAX = 0.0
  103.       DO 20 I=1,40
  104.           DO 10 J=1,40
  105.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  106.      1                     (I-J)/40.0
  107.               FMIN = MIN(F(I,J),FMIN)
  108.               FMAX = MAX(F(I,J),FMAX)
  109.    10     CONTINUE
  110.    20 CONTINUE
  111. C
  112. C Clear the screen. Set up window and viewport.
  113. C
  114.       CALL PGPAGE
  115.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  116.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONS')
  117. C
  118. C Draw the map.  PGCONS is called once for each contour, using
  119. C different line attributes to distinguish contour levels.
  120. C
  121.       CALL PGBBUF
  122.       DO 40 I=1,21
  123.           ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
  124.           IF (MOD(I,5).EQ.0) THEN
  125.               CALL PGSLW(3)
  126.           ELSE
  127.               CALL PGSLW(1)
  128.           END IF
  129.           IF (I.LT.10) THEN
  130.               CALL PGSCI(2)
  131.               CALL PGSLS(2)
  132.           ELSE
  133.               CALL PGSCI(3)
  134.               CALL PGSLS(1)
  135.           END IF
  136.           CALL PGCONS(F,40,40,1,40,1,40,ALEV,-1,TR)
  137.    40 CONTINUE
  138.       CALL PGSLW(1)
  139.       CALL PGSLS(1)
  140.       CALL PGSCI(1)
  141.       CALL PGEBUF
  142.       END
  143.  
  144.       SUBROUTINE PGEX33
  145. C-----------------------------------------------------------------------
  146. C Demonstration of contouring routine PGCONB.
  147. C-----------------------------------------------------------------------
  148.       REAL BLANK
  149.       PARAMETER (BLANK=-1.2E20)
  150.       INTEGER I,J
  151.       REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6),X,Y,R
  152.       DATA TR/0.,1.,0.,0.,0.,1./
  153. C
  154. C Compute a suitable function.
  155. C
  156.       FMIN = 0.0
  157.       FMAX = 0.0
  158.       DO 20 I=1,40
  159.           DO 10 J=1,40
  160.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  161.      1                     (I-J)/40.0
  162.               FMIN = MIN(F(I,J),FMIN)
  163.               FMAX = MAX(F(I,J),FMAX)
  164.    10     CONTINUE
  165.    20 CONTINUE
  166. C
  167. C "Blank" the data outside an annulus.
  168. C
  169.       DO 60 I=1,40
  170.           DO 50 J=1,40
  171.               R = SQRT((I-20.5)**2 + (J-20.5)**2)
  172.               IF (R.GT.20.0 .OR. R.LT.3.0) F(I,J) = BLANK
  173.    50     CONTINUE
  174.    60 CONTINUE
  175. C
  176.       CALL PGPAGE
  177.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  178.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONB')
  179.       CALL PGBBUF
  180.       DO 80 I=1,21
  181.           ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
  182.           IF (MOD(I,5).EQ.0) THEN
  183.               CALL PGSLW(3)
  184.           ELSE
  185.               CALL PGSLW(1)
  186.           END IF
  187.           IF (I.LT.10) THEN
  188.               CALL PGSCI(2)
  189.               CALL PGSLS(2)
  190.           ELSE
  191.               CALL PGSCI(3)
  192.               CALL PGSLS(1)
  193.           END IF
  194.           CALL PGCONB(F,40,40,1,40,1,40,ALEV,-1,TR,BLANK)
  195.    80 CONTINUE
  196.       CALL PGEBUF
  197. C
  198. C Mark the blanked points for easy identification.
  199. C
  200.       CALL PGBBUF
  201.       CALL PGSCI(1)
  202.       DO 100 I=1,40
  203.           DO 90 J=1,40
  204.               IF (F(I,J).EQ.BLANK) THEN
  205.                   X = TR(1) + REAL(I)*TR(2) + REAL(J)*TR(3)
  206.                   Y = TR(4) + REAL(I)*TR(5) + REAL(J)*TR(6)
  207.                   CALL PGPT1(X, Y, -1)
  208.               END IF
  209.    90     CONTINUE
  210.   100 CONTINUE
  211.       CALL PGEBUF
  212.       END
  213.  
  214.       SUBROUTINE PGEX34
  215. C-----------------------------------------------------------------------
  216. C This program is intended to demonstrate the use of the PGPLOT routine
  217. C PGCONX. As an example, we take data defined on a sphere. We want to
  218. C draw a contour map of the data on an equal-area projection of the
  219. C surface of the sphere; we choose the Hammer-Aitoff equal-area
  220. C projection centered on Declination (latitude) 0, Right Ascension
  221. C (longitude) 0. The data are defined at 2-degree intervals in both
  222. C coordinates. We thus need a data array dimensioned 181 by 91; the
  223. C array index runs from -90 to +90 in declination (91 elements) and
  224. C from -180 to +180 in right ascension (181 elements). The data at -180 
  225. C and +180 must be identical, of course, but they need to be duplicated 
  226. C in the array as these two longitudes appear on opposite sides of the
  227. C map. 
  228. C-----------------------------------------------------------------------
  229.       REAL  PI, RPDEG
  230.       PARAMETER (PI=3.14159265359)
  231.       PARAMETER (RPDEG=PI/180.0)
  232.       INTEGER I, J
  233.       REAL RA, DEC, B, L, XC(181), YC(181)
  234.       REAL Q(181,91), C(9)
  235.       EXTERNAL PLCAIT
  236. C
  237. C Call PGENV to create a rectangular window of 4 x 2 units. This is 
  238. C the bounding rectangle of the plot. The JUST argument is 1
  239. C to get equal scales in x and y.  
  240. C
  241.       CALL PGBBUF
  242.       CALL PGENV(-2.0, 2.0, -1.0, 1.0, 1, -2)
  243.       CALL PGLAB('Right Ascension', 'Declination', ' ')
  244.       CALL PGMTXT('t',2.0,0.0,0.0,
  245.      1           'Contouring on a non-Cartesian grid using PGCONX')
  246.       CALL PGSCH(0.6)
  247.       CALL PGMTXT('b',8.0,0.0,0.0,
  248.      1            'Hammer-Aitoff Equal-Area Projection of the Sphere')
  249.       CALL PGSCH(1.0)
  250. C
  251. C Draw 7 lines of constant longitude at longitude 0, 60, 120, ..., 
  252. C 360 degrees. Each line is made up of 90 straight-line segments.
  253. C
  254.       DO 20 J=1,7
  255.           RA = (-180.+(J-1)*60.)*RPDEG
  256.           DO 10 I=1,91
  257.               DEC = 2*(I-46)*RPDEG
  258.               CALL AITOFF(DEC,RA,XC(I),YC(I))
  259.    10     CONTINUE
  260.           CALL PGLINE(91,XC,YC)
  261.    20 CONTINUE
  262. C
  263. C Draw 5 lines of constant latitude at latitudes -60, -30, 0, 30, 
  264. C 60 degrees. Each line is made up of 360 straight-line segments.
  265. C
  266.       DO 40 J=1,5
  267.           DEC = (-60.+(J-1)*30.)*RPDEG
  268.           DO 30 I=1,181
  269.               RA = 2*(I-91)*RPDEG
  270.               CALL AITOFF(DEC,RA,XC(I),YC(I))
  271.    30     CONTINUE
  272.           CALL PGLINE(181,XC,YC)
  273.    40 CONTINUE
  274.       CALL PGEBUF
  275. C
  276. C Compute the data to be contoured. In practice the data might be read
  277. C in from an external file. In this example the data are computed: they
  278. C are the galactic latitudes of the points on the sphere. Thus the 
  279. C contours will be lines of constant galactic latitude.
  280. C
  281.       DO 60 J=1,91
  282.           DEC = 2*(J-46)*RPDEG
  283.           DO 50 I=1,181
  284.               RA = 2*(I-91)*RPDEG
  285.               CALL GALACT(RA, DEC, B,L)
  286.               Q(I,J) = B
  287.    50     CONTINUE
  288.    60 CONTINUE
  289. C
  290. C Draw the contour map using PGCONX. Contours at 0, 20, 40, 60, 80.
  291. C
  292.       DO 70 I=1,9
  293.           C(I) = -100.0 +I*20.0
  294.    70 CONTINUE
  295.       CALL PGBBUF
  296.       CALL PGSCI(2)
  297.       CALL PGCONX(Q, 181, 91, 1, 181, 1, 91, C, 9, PLCAIT)
  298.       CALL PGSCI(1)
  299.       CALL PGEBUF
  300.       END
  301.  
  302.       SUBROUTINE PLCAIT(VISBLE, X, Y, Z)
  303.       INTEGER VISBLE
  304.       REAL X,Y,Z
  305. C-----------------------------------------------------------------------
  306. C Plotting subroutine for PGCONX. This routine converts the input
  307. C coordinates (latitude and longitude) into the projected coordinates
  308. C (x and y), and moves or draws as requested by VISBLE.
  309. C-----------------------------------------------------------------------
  310.       REAL  PI, RPDEG
  311.       PARAMETER (PI=3.14159265359)
  312.       PARAMETER (RPDEG=PI/180.0)
  313.       REAL B, L, XWORLD, YWORLD
  314.       B = 2.0*(Y-46.0)*RPDEG
  315.       L = 2.0*(X-91.0)*RPDEG
  316.       CALL AITOFF(B, L, XWORLD, YWORLD)
  317.       IF (VISBLE.EQ.0) THEN
  318.           CALL PGMOVE(XWORLD, YWORLD)
  319.       ELSE
  320.           CALL PGDRAW(XWORLD, YWORLD)
  321.       END IF
  322.       END
  323.  
  324.       SUBROUTINE AITOFF(B,L,X,Y)
  325. C-----------------------------------------------------------------------
  326. C Hammer-Aitoff projection.
  327. C
  328. C       Input: latitude and longitude (B,L) in radians
  329. C       Output: cartesian (X,Y) in range +/-2, +/-1
  330. C-----------------------------------------------------------------------
  331.       REAL L,B,X,Y,L2,DEN
  332. C
  333.       L2 = L/2.0
  334.       DEN = SQRT(1.0+COS(B)*COS(L2))
  335.       X = 2.0*COS(B)*SIN(L2)/DEN
  336.       Y = SIN(B)/DEN
  337.       END
  338.  
  339.       SUBROUTINE GALACT(RA,DEC,GLAT,GLONG)
  340. C-----------------------------------------------------------------------
  341. C Convert 1950.0 equatorial coordinates (RA, DEC) to galactic
  342. C latitude and longitude (GLAT, GLONG).
  343. C
  344. C Arguments:
  345. C  RA, DEC (input): 1950.0 RA and Dec (radians).
  346. C  GLAT, GLONG (output): galactic latitude and longitude 
  347. C      (degrees).
  348. C
  349. C Reference: e.g., D. R. H. Johnson and D. R. Soderblom, A. J. v93, 
  350. C  p864 (1987).
  351. C-----------------------------------------------------------------------
  352.       REAL RA, RRA, DEC, RDEC, CDEC, R(3,3), E(3), G(3)
  353.       REAL RADDEG, GLAT, GLONG
  354.       INTEGER I, J
  355.       DATA R/-.066988740D0, .492728466D0,-.867600811D0,-.872755766D0,
  356.      $       -.450346958D0,-.188374601D0,-.483538915D0, .744584633D0,
  357.      $        .460199785D0/
  358.       DATA RADDEG/57.29577951D0/
  359. C-----------------------------------------------------------------------
  360.       RRA = RA
  361.       RDEC = DEC
  362.       CDEC = COS(RDEC)
  363.       E(1) = CDEC*COS(RRA)
  364.       E(2) = CDEC*SIN(RRA)
  365.       E(3) = SIN(RDEC)
  366.       DO 20 I=1,3
  367.           G(I) = 0.0
  368.           DO 10 J=1,3
  369.               G(I) = G(I) + E(J)*R(I,J)
  370.    10     CONTINUE
  371.    20 CONTINUE
  372.       GLAT  = ASIN(G(3))*RADDEG
  373.       GLONG = ATAN2(G(2),G(1))*RADDEG
  374.       IF (GLONG.LT.0.0) GLONG = GLONG+360.0
  375.       RETURN
  376. C-----------------------------------------------------------------------
  377.       END
  378.  
  379.       SUBROUTINE PGEX37
  380. C-----------------------------------------------------------------------
  381. C Demonstration of contouring routine PGCONX.
  382. C-----------------------------------------------------------------------
  383.       INTEGER I,J
  384.       REAL F(40,40),FMIN,FMAX,ALEV(1)
  385.       EXTERNAL PLCARO
  386. C
  387. C Compute a suitable function.
  388. C
  389.       FMIN = 0.0
  390.       FMAX = 0.0
  391.       DO 20 I=1,40
  392.           DO 10 J=1,40
  393.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  394.      1                     (I-J)/40.0
  395.               FMIN = MIN(F(I,J),FMIN)
  396.               FMAX = MAX(F(I,J),FMAX)
  397.    10     CONTINUE
  398.    20 CONTINUE
  399. C
  400. C Clear the screen. Set up window and viewport.
  401. C
  402.       CALL PGPAGE
  403.       CALL PGSVP(0.05,0.95,0.05,0.95)
  404.       CALL PGSWIN(1.0,40.0,1.0,40.0)
  405.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  406.       CALL PGMTXT('t',1.0,0.0,0.0,
  407.      :            'Contouring using PGCONX with arrows')
  408. C
  409. C Draw the map.  PGCONX is called once for each contour, using
  410. C different line attributes to distinguish contour levels.
  411. C
  412.       CALL PGBBUF
  413.       DO 30 I=1,21
  414.           ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
  415.           IF (MOD(I,5).EQ.0) THEN
  416.               CALL PGSLW(3)
  417.           ELSE
  418.               CALL PGSLW(1)
  419.           END IF
  420.           IF (I.LT.10) THEN
  421.               CALL PGSCI(2)
  422.               CALL PGSLS(2)
  423.           ELSE
  424.               CALL PGSCI(3)
  425.               CALL PGSLS(1)
  426.           END IF
  427.           CALL PGCONX(F,40,40,1,40,1,40,ALEV,-1,PLCARO)
  428.    30 CONTINUE
  429.       CALL PGSLW(1)
  430.       CALL PGSLS(1)
  431.       CALL PGSCI(1)
  432.       CALL PGEBUF
  433.       END
  434.  
  435.       SUBROUTINE PGEX36
  436. C-----------------------------------------------------------------------
  437. C Demonstration of contouring routine PGCONT and PGCONL.
  438. C-----------------------------------------------------------------------
  439.       INTEGER I,J
  440.       REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6)
  441.       CHARACTER*32 LABEL
  442.       DATA TR /0.0, 1.0, 0.0, 0.0, 0.0, 1.0/
  443. C
  444. C Compute a suitable function.
  445. C
  446.       FMIN = 0.0
  447.       FMAX = 0.0
  448.       DO 20 I=1,40
  449.           DO 10 J=1,40
  450.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  451.      1                     (I-J)/40.0
  452.               FMIN = MIN(F(I,J),FMIN)
  453.               FMAX = MAX(F(I,J),FMAX)
  454.    10     CONTINUE
  455.    20 CONTINUE
  456. C
  457. C Clear the screen. Set up window and viewport.
  458. C
  459.       CALL PGPAGE
  460.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  461.       CALL PGMTXT('t',1.0,0.0,0.0,
  462.      1            'Contouring using PGCONT and PGCONL labels')
  463. C
  464. C Draw the map.  PGCONT is called once for each contour, using
  465. C different line attributes to distinguish contour levels.
  466. C
  467.       CALL PGBBUF
  468.       DO 40 I=1,21
  469.           ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
  470.           IF (MOD(I,5).EQ.0) THEN
  471.               CALL PGSLW(3)
  472.           ELSE
  473.               CALL PGSLW(1)
  474.           END IF
  475.           IF (I.LT.10) THEN
  476.               CALL PGSCI(2)
  477.               CALL PGSLS(2)
  478.           ELSE
  479.               CALL PGSCI(3)
  480.               CALL PGSLS(1)
  481.           END IF
  482.           CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
  483.    40 CONTINUE
  484.       CALL PGSLW(1)
  485.       CALL PGSLS(1)
  486.       CALL PGEBUF
  487. C
  488. C Label the contours with PGCONL. Only even-numbered contours
  489. C are labelled.
  490. C
  491.       CALL PGBBUF
  492.       DO 50 I=2,21,2
  493.           ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
  494.           WRITE (LABEL,'(I2)') I
  495. C         WRITE (LABEL,'(F8.2)') ALEV
  496.           IF (I.LT.10) THEN
  497.               CALL PGSCI(2)
  498.           ELSE
  499.               CALL PGSCI(3)
  500.           END IF
  501.           CALL PGCONL(F,40,40,1,40,1,40,ALEV,TR,LABEL,16,8)
  502.  50   CONTINUE
  503.       CALL PGSCI(1)
  504.       CALL PGEBUF
  505.       END
  506.  
  507.       SUBROUTINE PLCARO(VISBLE, X, Y, Z)
  508.       INTEGER VISBLE
  509.       REAL X,Y,Z
  510. C-----------------------------------------------------------------------
  511. C Plotting subroutine for PGCONX. This routine labels contours with
  512. C arrows. Arrows point clockwise around minima, anticlockwise around
  513. C maxima. Arrows are drawn on 1/16 of contour line segments.
  514. C-----------------------------------------------------------------------
  515.       REAL XP, YP
  516.       INTEGER I
  517.       SAVE I
  518.       DATA I /0/
  519. C
  520.       I = MOD(I+1,16)
  521.       IF (VISBLE.EQ.0) THEN
  522.           I = 0
  523.           CALL PGMOVE(X, Y)
  524.       ELSE IF (I.EQ.8) THEN
  525. C         -- Draw line segment with arrow at midpoint
  526.           CALL PGQPOS(XP,YP)
  527.           CALL PGARRO(XP, YP, (X+XP)*0.5, (Y+YP)*0.5)
  528.           CALL PGDRAW(X, Y)
  529.       ELSE
  530. C         -- Draw plain line segment
  531.           CALL PGDRAW(X, Y)
  532.       END IF
  533.       END
  534.  
  535.       SUBROUTINE PGEXX1
  536. C-----------------------------------------------------------------------
  537. C Demonstration of contouring routine PGCONF.
  538. C-----------------------------------------------------------------------
  539.       INTEGER NX, NY, NC
  540.       PARAMETER (NX=51, NY=51, NC=9)
  541.       INTEGER I, J
  542.       REAL Z(NX,NY),TR(6), R
  543.       REAL X, Y, XMIN, XMAX, YMIN, YMAX, DX, DY, MU, C(NC)
  544.       DATA C /3.0, 3.2, 3.5, 3.6, 3.766413, 4.0 ,5.0, 10.0, 100.0/     
  545. C
  546. C Compute a suitable function. This is the test function used by
  547. C W. V. Snyder, Algorithm 531, Contour Plotting, ACM Trans. Math.
  548. C Softw. v.4, pp.290-294 (1978).
  549. C
  550.       XMIN = -2.0
  551.       XMAX = 2.0
  552.       YMIN =-2.0
  553.       YMAX = 2.0
  554.       MU = 0.3
  555.       DX = (XMAX-XMIN)/FLOAT(NX-1)                                      
  556.       DY = (YMAX-YMIN)/FLOAT(NY-1)
  557.       TR(1) = XMIN - DX
  558.       TR(2) = DX
  559.       TR(3) = 0.0
  560.       TR(4) = YMIN - DY
  561.       TR(5) = 0.0
  562.       TR(6) = DY
  563.       DO 20 I=1,NX
  564.          X = TR(1) + I*TR(2)
  565.          DO 10 J=1,NY     
  566.             Y = TR(4) + J*TR(6)
  567.             Z(I,J) = (1.0-MU)*(2.0/SQRT((X-MU)**2+Y**2)+(X-MU)**2+Y**2)   
  568.      *           + MU*(2.0/SQRT((X+1.0-MU)**2+Y**2)+(X+1.0-MU)**2+Y**2)      
  569.  10      CONTINUE                                   
  570.    20 CONTINUE                                                          
  571. C
  572. C Clear the screen. Set up window and viewport.
  573. C
  574.       CALL PGPAGE
  575.       CALL PGVSTD(0.05,0.95,0.05,0.95)
  576.       CALL PGWNAD(XMIN, XMAX, YMIN, YMAX)
  577. C
  578. C Fill contours with PGCONF.
  579. C
  580.       CALL PGSFS(1)
  581.       DO 30 I=1, NC-1
  582.          R = 0.5+0.5*REAL(I-1)/REAL(NC-1)
  583.          CALL PGSCR(I+10, R, R, R)
  584.          CALL PGSCI(I+10)
  585.          CALL PGCONF(Z,NX,NY,1,NX,1,NY,C(I),C(I+1),TR)
  586.  30   CONTINUE
  587. C
  588. C Draw the contour lines with PGCONT.
  589. C
  590.       CALL PGSCI(3)
  591.       CALL PGCONT(Z,NX,NY,1,NX,1,NY,C,NC,TR)
  592. C
  593. C Labels and box.
  594. C
  595.       CALL PGSCI(1)
  596.       CALL PGSCH(0.6)
  597.       CALL PGBOX('bctsin',1.0,10,'bctsinv',1.0,10)
  598.       CALL PGSCH(1.0)
  599.       CALL PGMTXT('t',1.0,0.0,0.0,'Contour filling using PGCONF')
  600. C
  601.       END
  602.