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