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