home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 3 / PDCD_3.iso / utilities / utilsp / pgplot / Examples / f77 / PGDemo1 next >
Encoding:
Text File  |  1994-02-24  |  29.5 KB  |  1,032 lines
  1.       PROGRAM PGDEM1
  2. C-----------------------------------------------------------------------
  3. C Demonstration program for PGPLOT. The main program opens the output
  4. C device and calls a series of subroutines, one for each sample plot.
  5. C-----------------------------------------------------------------------
  6.       INTEGER PGBEG
  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. Always
  10. C check the return code from PGBEG.
  11. C
  12.       IF (PGBEG(0,'?',1,1) .NE. 1) STOP
  13. C
  14. C Print information about device.
  15. C
  16.       CALL PGEX0
  17. C
  18. C Call the demonstration subroutines.
  19. C
  20.       CALL PGEX1
  21.       CALL PGEX2
  22.       CALL PGEX3
  23.       CALL PGEX4
  24.       CALL PGEX5
  25.       CALL PGEX6
  26.       CALL PGEX7
  27.       CALL PGEX8
  28.       CALL PGEX9
  29.       CALL PGEX10
  30.       CALL PGEX11
  31.       CALL PGEX12
  32.       CALL PGEX13
  33.       CALL PGEX14
  34.       CALL PGEX15
  35. C
  36. C Finally, call PGEND to terminate things properly.
  37. C
  38.       CALL PGEND
  39. C-----------------------------------------------------------------------
  40.       END
  41.  
  42.       SUBROUTINE PGEX0
  43. C-----------------------------------------------------------------------
  44. C This subroutine tests PGQINF and displays the information returned on
  45. C the standard output.
  46. C-----------------------------------------------------------------------
  47.       CHARACTER*64 VALUE
  48.       INTEGER LENGTH
  49.       REAL X, Y, X1, X2, Y1, Y2
  50. C
  51. C Information available from PGQINF:
  52. C
  53.       CALL PGQINF('version',  VALUE, LENGTH)
  54.       WRITE (*,*) 'version=', VALUE(:LENGTH)
  55.       CALL PGQINF('state',    VALUE, LENGTH)
  56.       WRITE (*,*) 'state=',   VALUE(:LENGTH)
  57.       CALL PGQINF('user',     VALUE, LENGTH)
  58.       WRITE (*,*) 'user=',    VALUE(:LENGTH)
  59.       CALL PGQINF('now',      VALUE, LENGTH)
  60.       WRITE (*,*) 'now=',     VALUE(:LENGTH)
  61.       CALL PGQINF('device',   VALUE, LENGTH)
  62.       WRITE (*,*) 'device=',  VALUE(:LENGTH)
  63.       CALL PGQINF('file',     VALUE, LENGTH)
  64.       WRITE (*,*) 'file=',    VALUE(:LENGTH)
  65.       CALL PGQINF('type',     VALUE, LENGTH)
  66.       WRITE (*,*) 'type=',    VALUE(:LENGTH)
  67.       CALL PGQINF('dev/type', VALUE, LENGTH)
  68.       WRITE (*,*) 'dev/type=',VALUE(:LENGTH)
  69.       CALL PGQINF('hardcopy', VALUE, LENGTH)
  70.       WRITE (*,*) 'hardcopy=',VALUE(:LENGTH)
  71.       CALL PGQINF('terminal', VALUE, LENGTH)
  72.       WRITE (*,*) 'terminal=',VALUE(:LENGTH)
  73.       CALL PGQINF('cursor',   VALUE, LENGTH)
  74.       WRITE (*,*) 'cursor=',  VALUE(:LENGTH)
  75. C
  76. C Get view surface dimensions:
  77. C
  78.       CALL PGQVSZ(1, X1, X2, Y1, Y2)
  79.       X = X2-X1
  80.       Y = Y2-Y1
  81.       WRITE (*,100) X, Y, X*25.4, Y*25.4
  82.   100 FORMAT (' Plot dimensions (x,y; inches): ',F9.2,', ',F9.2/
  83.      1        '                          (mm): ',F9.2,', ',F9.2)
  84. C-----------------------------------------------------------------------
  85.       END
  86.  
  87.       SUBROUTINE PGEX1
  88. C-----------------------------------------------------------------------
  89. C This example illustrates the use of PGENV, PGLAB, PGPT, PGLINE.
  90. C-----------------------------------------------------------------------
  91.       INTEGER I
  92.       REAL XS(5),YS(5), XR(100), YR(100)
  93.       DATA XS/1.,2.,3.,4.,5./
  94.       DATA YS/1.,4.,9.,16.,25./
  95. C
  96. C Call PGENV to specify the range of the axes and to draw a box, and
  97. C PGLAB to label it. The x-axis runs from 0 to 10, and y from 0 to 20.
  98. C
  99.       CALL PGENV(0.,10.,0.,20.,0,1)
  100.       CALL PGLAB('(x)', '(y)', 'PGPLOT Example 1:  y = x\u2')
  101. C
  102. C Mark five points (coordinates in arrays XS and YS), using symbol
  103. C number 9.
  104. C
  105.       CALL PGPT(5,XS,YS,9)
  106. C
  107. C Compute the function at 60 points, and use PGLINE to draw it.
  108. C
  109.       DO 10 I=1,60
  110.           XR(I) = 0.1*I
  111.           YR(I) = XR(I)**2
  112.    10 CONTINUE
  113.       CALL PGLINE(60,XR,YR)
  114. C-----------------------------------------------------------------------
  115.       END
  116.  
  117.       SUBROUTINE PGEX2
  118. C-----------------------------------------------------------------------
  119. C Repeat the process for another graph. This one is a graph of the
  120. C sinc (sin x over x) function.
  121. C-----------------------------------------------------------------------
  122.       INTEGER I
  123.       REAL XR(100), YR(100)
  124. C
  125.       CALL PGENV(-2.,10.,-0.4,1.2,0,1)
  126.       CALL PGLAB('(x)', 'sin(x)/x', 
  127.      $             'PGPLOT Example 2:  Sinc Function')
  128.       DO 20 I=1,100
  129.           XR(I) = (I-20)/6.
  130.           YR(I) = 1.0
  131.           IF (XR(I).NE.0.0) YR(I) = SIN(XR(I))/XR(I)
  132.    20 CONTINUE
  133.       CALL PGLINE(100,XR,YR)
  134. C-----------------------------------------------------------------------
  135.       END
  136.  
  137.       SUBROUTINE PGEX3
  138. C----------------------------------------------------------------------
  139. C This example illustrates the use of PGBOX and attribute routines to
  140. C mix colors and line-styles.
  141. C----------------------------------------------------------------------
  142.       REAL PI
  143.       PARAMETER (PI=3.14159265)
  144.       INTEGER I
  145.       REAL XR(360), YR(360)
  146.       REAL ARG
  147. C
  148. C Call PGENV to initialize the viewport and window; the
  149. C AXIS argument is -2, so no frame or labels will be drawn.
  150. C
  151.       CALL PGENV(0.,720.,-2.0,2.0,0,-2)
  152.       CALL PGSAVE
  153. C
  154. C Set the color index for the axes and grid (index 5 = cyan).
  155. C Call PGBOX to draw first a grid at low brightness, and then a
  156. C frame and axes at full brightness. Note that as the x-axis is
  157. C to represent an angle in degrees, we request an explicit tick 
  158. C interval of 90 deg with subdivisions at 30 deg, as multiples of
  159. C 3 are a more natural division than the default.
  160. C
  161.       CALL PGSCI(14)
  162.       CALL PGBOX('G',30.0,0,'G',0.2,0)
  163.       CALL PGSCI(5)
  164.       CALL PGBOX('ABCTSN',90.0,3,'ABCTSNV',0.0,0)
  165. C
  166. C Call PGLAB to label the graph in a different color (3=green).
  167. C
  168.       CALL PGSCI(3)
  169.       CALL PGLAB('x (degrees)','f(x)','PGPLOT Example 3')
  170. C
  171. C Compute the function to be plotted: a trig function of an
  172. C angle in degrees, computed every 2 degrees.
  173. C
  174.       DO 20 I=1,360
  175.           XR(I) = 2.0*I
  176.           ARG = XR(I)/180.0*PI
  177.           YR(I) = SIN(ARG) + 0.5*COS(2.0*ARG) + 
  178.      1                0.5*SIN(1.5*ARG+PI/3.0)
  179.    20 CONTINUE
  180. C
  181. C Change the color (6=magenta), line-style (2=dashed), and line
  182. C width and draw the function.
  183. C
  184.       CALL PGSCI(6)
  185.       CALL PGSLS(2)
  186.       CALL PGSLW(3)
  187.       CALL PGLINE(360,XR,YR)
  188. C
  189. C Restore attributes to defaults.
  190. C
  191.       CALL PGUNSA
  192. C-----------------------------------------------------------------------
  193.       END
  194.  
  195.       SUBROUTINE PGEX4
  196. C-----------------------------------------------------------------------
  197. C Demonstration program for PGPLOT: draw histograms.
  198. C-----------------------------------------------------------------------
  199.       INTEGER  I, ISEED
  200.       REAL     DATA(1000), X(620), Y(620)
  201.       REAL     RNGAUS, RAN5, JUNK
  202. C
  203. C Call RNGAUS to obtain 1000 samples from a normal distribution.
  204. C
  205.       ISEED = -5678921
  206.       JUNK = RAN5(ISEED)
  207.       DO 10 I=1,1000
  208.           DATA(I) = RNGAUS(ISEED)
  209.    10 CONTINUE
  210. C
  211. C Draw a histogram of these values.
  212. C
  213.       CALL PGSAVE
  214.       CALL PGHIST(1000,DATA,-3.1,3.1,31,0)
  215. C
  216. C Samples from another normal distribution.
  217. C
  218.       DO 15 I=1,200
  219.           DATA(I) = 1.0+0.5*RNGAUS(ISEED)
  220.    15 CONTINUE
  221. C
  222. C Draw another histogram (filled) on same axes.
  223. C
  224.       CALL PGSCI(15)
  225.       CALL PGHIST(200,DATA,-3.1,3.1,31,3)
  226.       CALL PGSCI(0)
  227.       CALL PGHIST(200,DATA,-3.1,3.1,31,1)
  228.       CALL PGSCI(1)
  229. C
  230. C Redraw the box which may have been clobbered by the histogram.
  231. C
  232.       CALL PGBOX('BST', 0.0, 0, ' ', 0.0, 0)
  233. C
  234. C Label the plot.
  235. C
  236.       CALL PGLAB('Variate', ' ',
  237.      $             'PGPLOT Example 4:  Histograms (Gaussian)')
  238. C
  239. C Superimpose the theoretical distribution.
  240. C
  241.       DO 20 I=1,620
  242.           X(I) = -3.1 + 0.01*(I-1)
  243.           Y(I) = 0.2*1000./SQRT(2.*3.14159265)*EXP(-0.5*X(I)*X(I))
  244.    20 CONTINUE
  245.       CALL PGLINE(620,X,Y)
  246.       CALL PGUNSA
  247. C-----------------------------------------------------------------------
  248.       END
  249.  
  250.       SUBROUTINE PGEX5
  251. C----------------------------------------------------------------------
  252. C Demonstration program for the PGPLOT plotting package.  This example
  253. C illustrates how to draw a log-log plot.
  254. C PGPLOT subroutines demonstrated:
  255. C    PGENV, PGERRY, PGLAB, PGLINE, PGPT, PGSCI.
  256. C----------------------------------------------------------------------
  257.       INTEGER   RED, GREEN, CYAN
  258.       PARAMETER (RED=2)
  259.       PARAMETER (GREEN=3)
  260.       PARAMETER (CYAN=5)
  261.       INTEGER   I
  262.       REAL      X, YLO, YHI
  263.       REAL      FREQ(15), FLUX(15), XP(100), YP(100), ERR(15)
  264.       DATA FREQ / 26., 38., 80., 160., 178., 318.,
  265.      1            365., 408., 750., 1400., 2695., 2700.,
  266.      2            5000., 10695., 14900. /
  267.       DATA FLUX / 38.0, 66.4, 89.0, 69.8, 55.9, 37.4,
  268.      1            46.8, 42.4, 27.0, 15.8, 9.09, 9.17,
  269.      2            5.35, 2.56, 1.73 /
  270.       DATA ERR  / 6.0, 6.0, 13.0, 9.1, 2.9, 1.4,
  271.      1            2.7, 3.0, 0.34, 0.8, 0.2, 0.46,
  272.      2            0.15, 0.08, 0.01 /
  273. C
  274. C Call PGENV to initialize the viewport and window; the AXIS argument 
  275. C is 30 so both axes will be logarithmic. The X-axis (frequency) runs 
  276. C from 0.01 to 100 GHz, the Y-axis (flux density) runs from 0.3 to 300
  277. C Jy. Note that it is necessary to specify the logarithms of these
  278. C quantities in the call to PGENV. We request equal scales in x and y
  279. C so that slopes will be correct.  Use PGLAB to label the graph.
  280. C
  281.       CALL PGSAVE
  282.       CALL PGSCI(CYAN)
  283.       CALL PGENV(-2.0,2.0,-0.5,2.5,1,30)
  284.       CALL PGLAB('Frequency, \gn (GHz)',
  285.      1             'Flux Density, S\d\gn\u (Jy)',
  286.      2             'PGPLOT Example 5:  Log-Log plot')
  287. C
  288. C Draw a fit to the spectrum (don't ask how this was chosen). This 
  289. C curve is drawn before the data points, so that the data will write 
  290. C over the curve, rather than vice versa.
  291. C
  292.       DO 10 I=1,100
  293.           X = 1.3 + I*0.03
  294.           XP(I) = X-3.0
  295.           YP(I) = 5.18 - 1.15*X -7.72*EXP(-X)
  296.    10 CONTINUE
  297.       CALL PGSCI(RED)
  298.       CALL PGLINE(100,XP,YP)
  299. C
  300. C Plot the measured flux densities: here the data are installed with a
  301. C DATA statement; in a more general program, they might be read from a
  302. C file. We first have to take logarithms (the -3.0 converts MHz to GHz).
  303. C
  304.       DO 20 I=1,15
  305.           XP(I) = ALOG10(FREQ(I))-3.0
  306.           YP(I) = ALOG10(FLUX(I))
  307.    20 CONTINUE
  308.       CALL PGSCI(GREEN)
  309.       CALL PGPT(15, XP, YP, 17)
  310. C
  311. C Draw +/- 2 sigma error bars: take logs of both limits.
  312. C
  313.       DO 30 I=1,15
  314.           YHI = ALOG10(FLUX(I)+2.*ERR(I))
  315.           YLO = ALOG10(FLUX(I)-2.*ERR(I))
  316.           CALL PGERRY(1,XP(I),YLO,YHI,1.0)
  317.    30 CONTINUE
  318.       CALL PGUNSA
  319. C-----------------------------------------------------------------------
  320.       END
  321.  
  322.       SUBROUTINE PGEX6
  323. C----------------------------------------------------------------------
  324. C Demonstration program for the PGPLOT plotting package.  This example
  325. C illustrates the use of PGPOLY using SOLID and HOLLOW fill-area
  326. C attributes.
  327. C----------------------------------------------------------------------
  328.       REAL TWOPI
  329.       PARAMETER (TWOPI=2.0*3.14159265)
  330.       INTEGER I, J
  331.       REAL X(10), Y(10)
  332. C
  333. C Call PGENV to initialize the viewport and window; the
  334. C AXIS argument is -2, so no frame or labels will be drawn.
  335. C
  336.       CALL PGSAVE
  337.       CALL PGENV(0.,8.,0.,8.0,1,-2)
  338. C
  339. C Call PGLAB to label the graph.
  340. C
  341.       CALL PGSCI(3)
  342.       CALL PGLAB(' ',' ','PGPLOT Example 6:  routine PGPOLY')
  343. C
  344. C Draw assorted regular convex polygons (solid).
  345. C
  346.       CALL PGSFS(1)
  347.       DO 20 I=3,9
  348.           CALL PGSCI(I-2)
  349.           DO 10 J=1,I
  350.             X(J) = I-2 + 0.5*COS(TWOPI*(J-1)/I)
  351.             Y(J) = 6 + 0.5*SIN(TWOPI*(J-1)/I)
  352.    10     CONTINUE
  353.           CALL PGPOLY (I,X,Y)
  354.    20 CONTINUE
  355. C
  356. C Draw assorted regular convex polygons (hollow).
  357. C
  358.       CALL PGSFS(2)
  359.       DO 40 I=3,9
  360.           CALL PGSCI(I-2)
  361.           DO 30 J=1,I
  362.             X(J) = I-2 + 0.5*COS(TWOPI*(J-1)/I)
  363.             Y(J) = 3 + 0.5*SIN(TWOPI*(J-1)/I)
  364.    30     CONTINUE
  365.           CALL PGPOLY (I,X,Y)
  366.    40 CONTINUE
  367. C
  368.       CALL PGUNSA
  369. C-----------------------------------------------------------------------
  370.       END
  371.  
  372.       SUBROUTINE PGEX7
  373. C-----------------------------------------------------------------------
  374. C A plot with a large number of symbols; plus test of PGERRB.
  375. C-----------------------------------------------------------------------
  376.       INTEGER I, ISEED
  377.       REAL XS(300),YS(300), XR(101), YR(101), XP, YP, XSIG, YSIG
  378.       REAL RAN5, RNGAUS
  379. C
  380. C Window and axes.
  381. C
  382.       CALL PGBBUF
  383.       CALL PGSAVE
  384.       CALL PGSCI(1)
  385.       CALL PGENV(0.,5.,-1.0,6.0,0,1)
  386.       CALL PGLAB('\fix', '\fiy', 'PGPLOT Example 7: scatter plot')
  387. C
  388. C Random data points.
  389. C
  390.       ISEED = -45678921
  391.       DO 10 I=1,300
  392.           XS(I) = 5.0*RAN5(ISEED)
  393.           YS(I) = SQRT(XS(I)) + 0.5*RNGAUS(ISEED)
  394.    10 CONTINUE
  395.       CALL PGSCI(3)
  396.       CALL PGPT(100,XS,YS,3)
  397.       CALL PGPT(100,XS(101),YS(101),17)
  398.       CALL PGPT(100,XS(201),YS(101),21)
  399. C
  400. C Curve defining parent distribution.
  401. C
  402.       DO 20 I=1,101
  403.           XR(I) = 0.05*(I-1)
  404.           YR(I) = SQRT(XR(I))
  405.    20 CONTINUE
  406.       CALL PGSCI(2)
  407.       CALL PGLINE(101,XR,YR)
  408. C
  409. C Test of PGERRB.
  410. C
  411.       XP = 3.0
  412.       YP = 3.0
  413.       XSIG = 0.75
  414.       YSIG = 1.2
  415.       CALL PGSCI(5)
  416.       CALL PGSCH(3.0)
  417.       CALL PGERRB(1, 1, XP, YP, XSIG, 1.0)
  418.       CALL PGERRB(2, 1, XP, YP, YSIG, 1.0)
  419.       CALL PGERRB(3, 1, XP, YP, XSIG, 1.0)
  420.       CALL PGERRB(4, 1, XP, YP, YSIG, 1.0)
  421.       CALL PGPT(1,XP,YP,21)
  422. C
  423.       CALL PGUNSA
  424.       CALL PGEBUF
  425. C-----------------------------------------------------------------------
  426.       END
  427.  
  428.       SUBROUTINE PGEX8
  429. C-----------------------------------------------------------------------
  430. C Demonstration program for PGPLOT. This program shows some of the
  431. C possibilities for overlapping windows and viewports.
  432. C T. J. Pearson  1986 Nov 28
  433. C-----------------------------------------------------------------------
  434.       INTEGER I
  435.       REAL XR(720), YR(720)
  436. C-----------------------------------------------------------------------
  437. C Color index:
  438.       INTEGER BLACK, WHITE, RED, GREEN, BLUE, CYAN, MAGENT, YELLOW
  439.       PARAMETER (BLACK=0)
  440.       PARAMETER (WHITE=1)
  441.       PARAMETER (RED=2)
  442.       PARAMETER (GREEN=3)
  443.       PARAMETER (BLUE=4)
  444.       PARAMETER (CYAN=5)
  445.       PARAMETER (MAGENT=6)
  446.       PARAMETER (YELLOW=7)
  447. C Line style:
  448.       INTEGER FULL, DASHED, DOTDSH, DOTTED, FANCY
  449.       PARAMETER (FULL=1)
  450.       PARAMETER (DASHED=2)
  451.       PARAMETER (DOTDSH=3)
  452.       PARAMETER (DOTTED=4)
  453.       PARAMETER (FANCY=5)
  454. C Character font:
  455.       INTEGER NORMAL, ROMAN, ITALIC, SCRIPT
  456.       PARAMETER (NORMAL=1)
  457.       PARAMETER (ROMAN=2)
  458.       PARAMETER (ITALIC=3)
  459.       PARAMETER (SCRIPT=4)
  460. C Fill-area style:
  461.       INTEGER SOLID, HOLLOW
  462.       PARAMETER (SOLID=1)
  463.       PARAMETER (HOLLOW=2)
  464. C-----------------------------------------------------------------------
  465. C
  466.       CALL PGPAGE
  467.       CALL PGBBUF
  468.       CALL PGSAVE
  469. C
  470. C Define the Viewport
  471. C
  472.       CALL PGSVP(0.1,0.6,0.1,0.6)
  473. C
  474. C Define the Window
  475. C
  476.       CALL PGSWIN(0.0, 630.0, -2.0, 2.0)
  477. C
  478. C Draw a box
  479. C
  480.       CALL PGSCI(CYAN)
  481.       CALL PGBOX ('ABCTS', 90.0, 3, 'ABCTSV', 0.0, 0)
  482. C
  483. C Draw labels
  484. C
  485.       CALL PGSCI (RED)
  486.       CALL PGBOX ('N',90.0, 3, 'VN', 0.0, 0)
  487. C
  488. C Draw SIN line
  489. C
  490.       DO 10 I=1,360
  491.           XR(I) = 2.0*I
  492.           YR(I) = SIN(XR(I)/57.29577951)
  493.    10 CONTINUE
  494.       CALL PGSCI (MAGENT)
  495.       CALL PGSLS (DASHED)
  496.       CALL PGLINE (360,XR,YR)
  497. C
  498. C Draw COS line by redefining the window
  499. C
  500.       CALL PGSWIN (90.0, 720.0, -2.0, 2.0)
  501.       CALL PGSCI (YELLOW)
  502.       CALL PGSLS (DOTTED)
  503.       CALL PGLINE (360,XR,YR)
  504.       CALL PGSLS (FULL)
  505. C
  506. C Re-Define the Viewport
  507. C
  508.       CALL PGSVP(0.45,0.85,0.45,0.85)
  509. C
  510. C Define the Window, and erase it
  511. C
  512.       CALL PGSWIN(0.0, 180.0, -2.0, 2.0)
  513.       CALL PGSCI(0)
  514.       CALL PGRECT(0.0, 180., -2.0, 2.0)
  515. C
  516. C Draw a box
  517. C
  518.       CALL PGSCI(BLUE)
  519.       CALL PGBOX ('ABCTSM', 60.0, 3, 'VABCTSM', 1.0, 2)
  520. C
  521. C Draw SIN line
  522. C
  523.       CALL PGSCI (WHITE)
  524.       CALL PGSLS (DASHED)
  525.       CALL PGLINE (360,XR,YR)
  526. C
  527.       CALL PGUNSA
  528.       CALL PGEBUF
  529. C-----------------------------------------------------------------------
  530.       END
  531.  
  532.       SUBROUTINE PGEX9
  533. C----------------------------------------------------------------------
  534. C Demonstration program for the PGPLOT plotting package.  This example
  535. C illustrates curve drawing with PGFUNT; the parametric curve drawn is
  536. C a simple Lissajous figure.
  537. C                              T. J. Pearson  1983 Oct 5
  538. C----------------------------------------------------------------------
  539.       REAL     FX, FY
  540.       EXTERNAL FX, FY
  541. C
  542. C Call PGFUNT to draw the function (autoscaling).
  543. C
  544.       CALL PGBBUF
  545.       CALL PGSAVE
  546.       CALL PGSCI(5)
  547.       CALL PGFUNT(FX,FY,360,0.0,2.0*3.14159265,0)
  548. C
  549. C Call PGLAB to label the graph in a different color.
  550. C
  551.       CALL PGSCI(3)
  552.       CALL PGLAB('x','y','PGPLOT Example 9:  routine PGFUNT')
  553.       CALL PGUNSA
  554.       CALL PGEBUF
  555. C
  556.       END
  557.  
  558.       REAL FUNCTION FX(T)
  559.       REAL T
  560.       FX = SIN(T*5.0)
  561.       RETURN
  562.       END
  563.  
  564.       REAL FUNCTION FY(T)
  565.       REAL T
  566.       FY = SIN(T*4.0)
  567.       RETURN
  568.       END
  569.  
  570.       SUBROUTINE PGEX10
  571. C----------------------------------------------------------------------
  572. C Demonstration program for the PGPLOT plotting package.  This example
  573. C illustrates curve drawing with PGFUNX.
  574. C                              T. J. Pearson  1983 Oct 5
  575. C----------------------------------------------------------------------
  576. C The following define mnemonic names for the color indices and
  577. C linestyle codes.
  578.       INTEGER   BLACK, WHITE, RED, GREEN, BLUE, CYAN, MAGENT, YELLOW
  579.       PARAMETER (BLACK=0)
  580.       PARAMETER (WHITE=1)
  581.       PARAMETER (RED=2)
  582.       PARAMETER (GREEN=3)
  583.       PARAMETER (BLUE=4)
  584.       PARAMETER (CYAN=5)
  585.       PARAMETER (MAGENT=6)
  586.       PARAMETER (YELLOW=7)
  587.       INTEGER   FULL, DASH, DOTD
  588.       PARAMETER (FULL=1)
  589.       PARAMETER (DASH=2)
  590.       PARAMETER (DOTD=3)
  591. C
  592. C The Fortran functions to be plotted must be declared EXTERNAL.
  593. C
  594.       REAL     BESJ0, BESJ1
  595.       EXTERNAL BESJ0, BESJ1
  596. C
  597. C Call PGFUNX twice to draw two functions (autoscaling the first time).
  598. C
  599.       CALL PGBBUF
  600.       CALL PGSAVE
  601.       CALL PGSCI(YELLOW)
  602.       CALL PGFUNX(BESJ0,500,0.0,10.0*3.14159265,0)
  603.       CALL PGSCI(RED)
  604.       CALL PGSLS(DASH)
  605.       CALL PGFUNX(BESJ1,500,0.0,10.0*3.14159265,1)
  606. C
  607. C Call PGLAB to label the graph in a different color. Note the
  608. C use of "\f" to change font.  Use PGMTXT to write an additional
  609. C legend inside the viewport.
  610. C
  611.       CALL PGSCI(GREEN)
  612.       CALL PGSLS(FULL)
  613.       CALL PGLAB('\fix', '\fiy',
  614.      2           '\frPGPLOT Example 10: routine PGFUNX')
  615.       CALL PGMTXT('T', -4.0, 0.5, 0.5,
  616.      1     '\frBessel Functions')
  617. C
  618. C Call PGARRO to label the curves.
  619. C
  620.       CALL PGARRO(8.0, 0.7, 1.0, BESJ0(1.0))
  621.       CALL PGARRO(12.0, 0.5, 9.0, BESJ1(9.0))
  622.       CALL PGSTBG(GREEN)
  623.       CALL PGSCI(0)
  624.       CALL PGPTXT(8.0, 0.7, 0.0, 0.0, ' \fiy = J\d0\u(x)')
  625.       CALL PGPTXT(12.0, 0.5, 0.0, 0.0, ' \fiy = J\d1\u(x)')
  626.       CALL PGUNSA
  627.       CALL PGEBUF
  628. C-----------------------------------------------------------------------
  629.       END
  630.  
  631.       SUBROUTINE PGEX11
  632. C-----------------------------------------------------------------------
  633. C Test routine for PGPLOT: draws a skeletal dodecahedron.
  634. C-----------------------------------------------------------------------
  635.       INTEGER NVERT
  636.       REAL T, T1, T2, T3
  637.       PARAMETER (NVERT=20)
  638.       PARAMETER (T=1.618)
  639.       PARAMETER (T1=1.0+T)
  640.       PARAMETER (T2=-1.0*T)
  641.       PARAMETER (T3=-1.0*T1)
  642.       INTEGER I, J, K
  643.       REAL VERT(3,NVERT), R, ZZ
  644.       REAL X(2),Y(2)
  645. C
  646. C Cartesian coordinates of the 20 vertices.
  647. C
  648.       DATA VERT/ T, T, T,       T, T,T2,
  649.      3           T,T2, T,       T,T2,T2,
  650.      5          T2, T, T,      T2, T,T2,
  651.      7          T2,T2, T,      T2,T2,T2,
  652.      9          T1,1.0,0.0,    T1,-1.0,0.0,
  653.      B          T3,1.0,0.0,    T3,-1.0,0.0,
  654.      D          0.0,T1,1.0,    0.0,T1,-1.0,
  655.      F          0.0,T3,1.0,    0.0,T3,-1.0,
  656.      H          1.0,0.0,T1,    -1.0,0.0,T1,
  657.      J          1.0,0.0,T3,   -1.0,0.0,T3 /
  658. C
  659. C Initialize the plot (no labels).
  660. C
  661.       CALL PGBBUF
  662.       CALL PGSAVE
  663.       CALL PGENV(-4.,4.,-4.,4.,1,-2)
  664.       CALL PGSCI(2)
  665.       CALL PGSLS(1)
  666.       CALL PGSLW(1)
  667. C
  668. C Write a heading.
  669. C
  670.       CALL PGLAB(' ',' ','PGPLOT Example 11:  Dodecahedron')
  671. C
  672. C Mark the vertices.
  673. C
  674.       DO 2 I=1,NVERT
  675.           ZZ = VERT(3,I)
  676.           CALL PGPT(1,VERT(1,I)+0.2*ZZ,VERT(2,I)+0.3*ZZ,9)
  677.     2 CONTINUE
  678. C
  679. C Draw the edges - test all vertex pairs to find the edges of the 
  680. C correct length.
  681. C
  682.       CALL PGSLW(3)
  683.       DO 20 I=2,NVERT
  684.           DO 10 J=1,I-1
  685.               R = 0.
  686.               DO 5 K=1,3
  687.                   R = R + (VERT(K,I)-VERT(K,J))**2
  688.     5         CONTINUE
  689.               R = SQRT(R)
  690.               IF(ABS(R-2.0).GT.0.1) GOTO 10
  691.               ZZ = VERT(3,I)
  692.               X(1) = VERT(1,I)+0.2*ZZ
  693.               Y(1) = VERT(2,I)+0.3*ZZ
  694.               ZZ = VERT(3,J)
  695.               X(2) = VERT(1,J)+0.2*ZZ
  696.               Y(2) = VERT(2,J)+0.3*ZZ
  697.               CALL PGLINE(2,X,Y)
  698.    10     CONTINUE
  699.    20 CONTINUE
  700.       CALL PGUNSA
  701.       CALL PGEBUF
  702. C-----------------------------------------------------------------------
  703.       END
  704.  
  705.       SUBROUTINE PGEX12
  706. C-----------------------------------------------------------------------
  707. C Test routine for PGPLOT: draw arrows with PGARRO.
  708. C-----------------------------------------------------------------------
  709.       INTEGER NV, I, K
  710.       REAL A, D, X, Y, XT, YT
  711. C
  712. C Number of arrows.
  713. C
  714.       NV =16
  715. C
  716. C Select a square viewport.
  717. C
  718.       CALL PGBBUF
  719.       CALL PGSAVE
  720.       CALL PGSCH(0.7)
  721.       CALL PGSCI(2)
  722.       CALL PGENV(-1.05,1.05,-1.05,1.05,1,-1)
  723.       CALL PGLAB(' ', ' ', 'PGPLOT Example 12: PGARRO')
  724.       CALL PGSCI(1)
  725. C
  726. C Draw the arrows
  727. C
  728.       K = 1
  729.       D = 360.0/57.29577951/NV
  730.       A = -D
  731.       DO 20 I=1,NV
  732.           A = A+D
  733.           X = COS(A)
  734.           Y = SIN(A)
  735.           XT = 0.2*COS(A-D)
  736.           YT = 0.2*SIN(A-D)
  737.           CALL PGSAH(K, 80.0-3.0*I, 0.5*REAL(I)/REAL(NV))
  738.           CALL PGSCH(0.25*I)
  739.           CALL PGARRO(XT, YT, X, Y)
  740.           K = K+1
  741.           IF (K.GT.2) K=1
  742.    20 CONTINUE
  743. C
  744.       CALL PGUNSA
  745.       CALL PGEBUF
  746. C-----------------------------------------------------------------------
  747.       END
  748.  
  749.       SUBROUTINE PGEX13
  750. C----------------------------------------------------------------------
  751. C This example illustrates the use of PGTBOX.
  752. C----------------------------------------------------------------------
  753.       INTEGER N
  754.       PARAMETER (N=10)
  755.       INTEGER I
  756.       REAL X1(N), X2(N)
  757.       CHARACTER*20 XOPT(N), BSL*1
  758.       DATA X1 /   4*0.0, -8000.0, 100.3, 205.3, -45000.0, 2*0.0/
  759.       DATA X2 /4*8000.0,  8000.0, 101.3, 201.1, 3*-100000.0/
  760.       DATA XOPT / 'BSTN', 'BSTNZ', 'BSTNZH', 'BSTNZD', 'BSNTZHFO', 
  761.      :      'BSTNZD', 'BSTNZHI', 'BSTNZHP', 'BSTNZDY', 'BSNTZHFOY'/
  762. C
  763.       BSL = CHAR(92)
  764.       CALL PGPAGE
  765.       CALL PGSAVE
  766.       CALL PGBBUF
  767.       CALL PGSCH(0.7)
  768.       DO 100 I=1,N
  769.         CALL PGSVP(0.15, 0.85, (0.7+REAL(N-I))/REAL(N), 
  770.      :                         (0.7+REAL(N-I+1))/REAL(N)) 
  771.         CALL PGSWIN(X1(I), X2(I), 0.0, 1.0)
  772.         CALL PGTBOX(XOPT(I),0.0,0,' ',0.0,0)
  773.         CALL PGLAB('Option = '//XOPT(I), ' ', ' ')
  774.         IF (I.EQ.1) CALL PGMTXT('B', -1.0, 0.5, 0.5,
  775.      :             BSL//'fiAxes drawn with PGTBOX')
  776.   100 CONTINUE
  777.       CALL PGEBUF
  778.       CALL PGUNSA
  779. C-----------------------------------------------------------------------
  780.       END
  781.  
  782.       SUBROUTINE PGEX14
  783. C-----------------------------------------------------------------------
  784. C Test routine for PGPLOT: polygon fill and color representation.
  785. C-----------------------------------------------------------------------
  786.       INTEGER I, J, N, M
  787.       REAL PI, THINC, R, G, B, THETA
  788.       REAL XI(100),YI(100),XO(100),YO(100),XT(3),YT(3)
  789.       PARAMETER(PI = 3.14159265)
  790. C
  791.       N = 33
  792.       M = 8
  793.       THINC=2.0*PI/N
  794.       DO 10 I=1,N
  795.         XI(I) = 0.0
  796.         YI(I) = 0.0
  797.    10 CONTINUE
  798.       CALL PGBBUF
  799.       CALL PGSAVE
  800.       CALL PGENV(-1.,1.,-1.,1.,1,-2)
  801.       CALL PGLAB(' ', ' ', 'PGPLOT Example 14: PGPOLY and PGSCR')
  802.       DO 50 J=1,M
  803.         R = 1.0
  804.         G = 1.0 - REAL(J)/REAL(M)
  805.         B = G
  806.         CALL PGSCR(J, R, G, B)
  807.         THETA = -REAL(J)*PI/REAL(N)
  808.         R = REAL(J)/REAL(M)
  809.         DO 20 I=1,N
  810.           THETA = THETA+THINC
  811.           XO(I) = R*COS(THETA)
  812.           YO(I) = R*SIN(THETA)
  813.    20   CONTINUE
  814.         DO 30 I=1,N
  815.           XT(1) = XO(I)
  816.           YT(1) = YO(I)
  817.           XT(2) = XO(MOD(I,N)+1)
  818.           YT(2) = YO(MOD(I,N)+1)
  819.           XT(3) = XI(I)
  820.           YT(3) = YI(I)
  821.           CALL PGSCI(J)
  822.           CALL PGSFS(1)
  823.           CALL PGPOLY(3,XT,YT)
  824.           CALL PGSFS(2)
  825.           CALL PGSCI(1)
  826.           CALL PGPOLY(3,XT,YT)
  827.    30   CONTINUE
  828.         DO 40 I=1,N
  829.           XI(I) = XO(I)
  830.           YI(I) = YO(I)
  831.    40   CONTINUE
  832.    50 CONTINUE
  833.       CALL PGUNSA
  834.       CALL PGEBUF
  835. C-----------------------------------------------------------------------
  836.       END
  837.  
  838.       SUBROUTINE PGEX15
  839. C----------------------------------------------------------------------
  840. C This is a line-drawing test; it draws a regular n-gon joining
  841. C each vertex to every other vertex. It is not optimized for pen
  842. C plotters.
  843. C----------------------------------------------------------------------
  844.       INTEGER I, J, NV
  845.       REAL A, D, X(100), Y(100)
  846. C
  847. C Set the number of vertices, and compute the 
  848. C coordinates for unit circumradius.
  849. C
  850.       NV = 17
  851.       D = 360.0/NV
  852.       A = -D
  853.       DO 20 I=1,NV
  854.           A = A+D
  855.           X(I) = COS(A/57.29577951)
  856.           Y(I) = SIN(A/57.29577951)
  857.    20 CONTINUE
  858. C
  859. C Select a square viewport.
  860. C
  861.       CALL PGBBUF
  862.       CALL PGSAVE
  863.       CALL PGSCH(0.5)
  864.       CALL PGSCI(2)
  865.       CALL PGENV(-1.05,1.05,-1.05,1.05,1,-1)
  866.       CALL PGLAB(' ', ' ', 'PGPLOT Example 15: PGMOVE and PGDRAW')
  867.       CALL PGSCR(0,0.2,0.3,0.3)
  868.       CALL PGSCR(1,1.0,0.5,0.2)
  869.       CALL PGSCR(2,0.2,0.5,1.0)
  870.       CALL PGSCI(1)
  871. C
  872. C Draw the polygon.
  873. C
  874.       DO 40 I=1,NV-1
  875.           DO 30 J=I+1,NV
  876.             CALL PGMOVE(X(I),Y(I))
  877.             CALL PGDRAW(X(J),Y(J))
  878.    30     CONTINUE
  879.    40 CONTINUE
  880. C
  881. C Flush the buffer.
  882. C
  883.       CALL PGUNSA
  884.       CALL PGEBUF
  885. C-----------------------------------------------------------------------
  886.       END
  887.  
  888.       REAL FUNCTION BESJ0(XX)
  889.       REAL XX
  890. C     BESSEL FUNCTION J0(XX)
  891. C     REVISED SEPT, 1971  -  TRANSFERED TO VAX JULY 1979.
  892. C     J0(-XX) = J0(XX)
  893. C-----------------------------------------------------------------------
  894.       REAL X, XO3, T, F0, THETA0
  895. C
  896.       X = ABS(XX)
  897.       IF (X .LE. 3.0) THEN
  898.           XO3 = X/3.
  899.           T   = XO3*XO3
  900.           BESJ0 = 1.+T*(-2.2499997+T*(1.2656208+T*(-.3163866+T*(.0444479
  901.      1            +T*(-.003944+T*.0002100)))))
  902.       ELSE
  903.           T = 3./X
  904.           F0 = .79788456+T*(-.00000077+T*(-.00552740+T*(-.00009512
  905.      1         +T*(.00137237+T*(-.00072805+T*.00014476)))))
  906.           THETA0 = X-.78539816+T*(-.04166397+T*(-.00003954
  907.      1      +T*(.00262573+T*(-.00054125+T*(-.00029333+T*.00013558)))))
  908.           BESJ0 = F0*COS(THETA0)/SQRT(X)
  909.       END IF
  910.       RETURN
  911.       END
  912.  
  913.       REAL FUNCTION BESJ1(XX)
  914.       REAL XX
  915. C     BESSEL FUNCTION J1(XX)
  916. C     REVISED SEPT,1971
  917. C     J1(-XX)=-J1(XX)
  918. C     TRANSFERED TO VAX    JULY 1979
  919. C-----------------------------------------------------------------------
  920.       REAL X, XO3, T, F1, THETA1
  921. C
  922.       X = ABS(XX)
  923.       IF (X .LE. 3.0) THEN
  924.           XO3 = X/3.
  925.           T = XO3*XO3
  926.           BESJ1 = .5+T*(-.56249985+T*(.21093573+T*
  927.      1                 (-.03954289+T*(.00443319+T*
  928.      2                 (-.00031761+T*.00001109)))))
  929.           BESJ1 = BESJ1*XX
  930.       ELSE
  931.           T = 3./X
  932.           F1 = .79788456+T*(.00000156+T*(.01659667+T*(.00017105
  933.      1         +T*(-.00249511+T*(.00113653-T*.00020033)))))
  934.           THETA1 = X-2.35619449+T*(.12499612+T*(.00005650+T*(-.00637879
  935.      1             +T*(.00074348+T*(.00079824-T*.00029166)))))
  936.           BESJ1 = F1*COS(THETA1)/SQRT(X)
  937.       END IF
  938.       IF (XX .LT. 0.0) BESJ1 = -BESJ1
  939.       RETURN
  940.       END
  941.  
  942. C*RNGAUS -- random number from Gaussian (normal) distribution
  943. C+
  944.       REAL FUNCTION RNGAUS (ISEED)
  945.       INTEGER ISEED
  946. C
  947. C Returns a normally distributed deviate with zero mean and unit 
  948. C variance. The routine uses the Box-Muller transformation of uniform
  949. C deviates. Reference: Press et al., Numerical Recipes, Sec. 7.2.
  950. C
  951. C Arguments:
  952. C  ISEED  (in/out) : seed used for RAN5 random-number generator.
  953. C
  954. C Subroutines required:
  955. C  RAN5 -- return a uniform random deviate between 0 and 1.
  956. C
  957. C History:
  958. C  1987 Nov 13 - TJP.
  959. C-----------------------------------------------------------------------
  960.       INTEGER ISET
  961.       REAL R, V1, V2, FAC, GSET
  962.       REAL RAN5
  963.       SAVE ISET, GSET
  964.       DATA ISET/0/
  965. C
  966.       IF (ISET.EQ.0) THEN
  967.    10     V1 = 2.*RAN5(ISEED)-1.
  968.           V2 = 2.*RAN5(ISEED)-1.
  969.           R = V1**2+V2**2
  970.           IF (R.GE.1.) GOTO 10
  971.           FAC = SQRT(-2.*LOG(R)/R)
  972.           GSET = V1*FAC
  973.           RNGAUS = V2*FAC
  974.           ISET = 1
  975.       ELSE
  976.           RNGAUS = GSET
  977.           ISET = 0
  978.       END IF
  979. C
  980.       END
  981.  
  982. C*RAN5 -- random number from uniform distribution
  983. C+
  984.       REAL FUNCTION RAN5(IDUM)
  985.       INTEGER IDUM
  986. C
  987. C Returns a uniform random deviate between 0.0 and 1.0.
  988. C
  989. C Park and Miller's "Minimal Standard" random number generator (Comm.
  990. C ACM, 31, 1192, 1988) with Bays-Durham shuffle. Call with IDUM negative
  991. C to initialize. Be sure to preserve IDUM between calls.
  992. C Reference: Press and Farrar, Computers in Physics, Vol.4, No. 2, 
  993. C p.190, 1990.
  994. C
  995. C Arguments:
  996. C  IDUM  (in/out) : seed.
  997. C
  998. C History:
  999. C  1990 Apr 6 - TJP.
  1000. C-----------------------------------------------------------------------
  1001.       INTEGER IA, IM, IQ, IR, NTAB
  1002.       REAL    AM, ATAB
  1003.       PARAMETER (IA=16807, IM=2147483647, AM=1.0/IM)
  1004.       PARAMETER (IQ=127773, IR=2836, NTAB=32, ATAB=NTAB-1)
  1005. C-
  1006.       INTEGER J, K
  1007.       REAL V(NTAB), Y
  1008.       SAVE V, Y
  1009.       DATA V/NTAB*0.0/, Y/0.5/
  1010. C-
  1011.       IF (IDUM.LE.0) THEN
  1012.           IDUM = MAX(-IDUM,1)
  1013.           DO 12 J=NTAB,1,-1
  1014.               K = IDUM/IQ
  1015.               IDUM = IA*(IDUM-K*IQ) - IR*K
  1016.               IF (IDUM.LT.0) IDUM = IDUM+IM
  1017.               V(J) = AM*IDUM
  1018.    12     CONTINUE
  1019.           Y = V(1)
  1020.       END IF
  1021.     1 CONTINUE
  1022.           K = IDUM/IQ
  1023.           IDUM = IA*(IDUM-K*IQ) - IR*K
  1024.           IF (IDUM.LT.0) IDUM = IDUM+IM
  1025.           J = 1 + INT(ATAB*Y)
  1026.           Y = V(J)
  1027.           RAN5 = Y
  1028.           V(J) = AM*IDUM
  1029.       IF (RAN5.EQ.0.0 .OR. RAN5.EQ.1.0) GOTO 1
  1030.       RETURN
  1031.       END
  1032.