PROGRAM PGDE13
C-----------------------------------------------------------------------
C Demonstration program for PGPLOT with multiple devices.
C It requires an interactive device which presents a menu of graphs
C to be displayed on the second device, which may be interactive or
C hardcopy.
C-----------------------------------------------------------------------
      INTEGER PGOPEN, ID, ID1, ID2, NP
C
C Call PGOPEN to initiate PGPLOT and open the output device; PGBEG
C will prompt the user to supply the device name and type. Always
C check the return code from PGBEG.

      WRITE (*,*) 'This program demonstrates the use of two devices'
      WRITE (*,*) 'in PGPLOT.  An interactive device is used to'
      WRITE (*,*) 'present a menu of graphs that may be displayed on'
      WRITE (*,*) 'a second device. Use the cursor or mouse to select'
      WRITE (*,*) 'the graph to be displayed. It is also possible'
      WRITE (*,*) 'to display either 1 graph per page or 4 graphs'
      WRITE (*,*) 'per page.'
      WRITE (*,*) 'If you have an X-Window display, try specifying'
      WRITE (*,*) '/XWIN for both devices.'
      WRITE (*,*)
C
      ID1 = PGOPEN('?Graphics device for menu (eg, /XWIN): ')    
      IF (ID1.LE.0) STOP
      CALL INIT
      CALL PGASK(.FALSE.)
      ID2 = PGOPEN('?Graphics device for graphs (eg, file/PS): ')
      IF (ID2.LE.0) STOP
      CALL PGASK(.FALSE.)
C
C Select a plot.
C
      NP = 1
 100  CALL PGSLCT(ID1)
      CALL MENU(NP, ID)
      CALL PGSLCT(ID2)
      CALL PGSAVE
      CALL PGBBUF
      IF (ID.EQ.1) THEN
         CALL PGEX1
      ELSE IF (ID.EQ.2) THEN
         CALL PGEX2
      ELSE IF (ID.EQ.3) THEN
         CALL PGEX3
      ELSE IF (ID.EQ.4) THEN
         CALL PGEX4
      ELSE IF (ID.EQ.5) THEN
         CALL PGEX5
      ELSE IF (ID.EQ.6) THEN
         CALL PGEX6
      ELSE IF (ID.EQ.7) THEN
         CALL PGEX7
      ELSE IF (ID.EQ.8) THEN
         CALL PGEX8
      ELSE IF (ID.EQ.9) THEN
         CALL PGEX9
      ELSE IF (ID.EQ.10) THEN
         CALL PGEX10
      ELSE IF (ID.EQ.11) THEN
         CALL PGEX11
      ELSE IF (ID.EQ.12) THEN
         CALL PGEX12
      ELSE IF (ID.EQ.13) THEN
         CALL PGEX13
      ELSE IF (ID.EQ.14) THEN
         CALL PGSUBP(1,1)
         NP = 1
      ELSE IF (ID.EQ.15) THEN
         CALL PGSUBP(2,2)
         NP = 4
      ELSE
         GOTO 200
      END IF
      CALL PGEBUF
      CALL PGUNSA
      GOTO 100
C
C Done: close devices.
C
 200  CALL PGEND
C-----------------------------------------------------------------------
      END

      SUBROUTINE INIT
C
C Set up graphics device to display menu.
C-----------------------------------------------------------------------
      CALL PGPAP(2.5, 2.0)
      CALL PGPAGE
      CALL PGSVP(0.0,1.0,0.0,1.0)
      CALL PGSWIN(0.0,0.5,0.0,1.0)
      CALL PGSCR(0, 0.4, 0.4, 0.4)
      RETURN
C-----------------------------------------------------------------------
      END

      SUBROUTINE MENU(NP, ID)
      INTEGER NP, ID
C
C Display menu of plots.
C-----------------------------------------------------------------------
      INTEGER NBOX
      PARAMETER (NBOX=16)
      CHARACTER*12 VALUE(NBOX)
      INTEGER I, JUNK, K
      REAL X1, X2, Y(NBOX), XX, YY, R
      CHARACTER CH
      INTEGER PGCURS
C
      DATA VALUE / '1', '2', '3', '4', '5', '6', '7', '8', '9',
     :             '10', '11', '12', '13',
     :             'One panel', 'Four panels', 'EXIT' /

      DATA XX/0.5/, YY/0.5/
C
      X1 = 0.1
      X2 = 0.2
      DO 5 I=1,NBOX
         Y(I) = 1.0 - REAL(I+1)/REAL(NBOX+2)
 5    CONTINUE
C
C Display buttons.
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGERAS
      CALL PGSCI(1)
      CALL PGSCH(2.5)
      CALL PGPTXT(X1, 1.0-1.0/REAL(NBOX+2), 0.0, 0.0, '\fiMENU')
      CALL PGSLW(1)
      CALL PGSCH(2.0)
      DO 10 I=1,NBOX
         CALL PGSCI(1)
         CALL PGSFS(1)
         CALL PGCIRC(X1, Y(I), 0.02)
         CALL PGSCI(2)
         CALL PGSFS(2)
         CALL PGCIRC(X1, Y(I), 0.02)
         CALL PGSCI(1)
         CALL PGPTXT(X2, Y(I), 0.0, 0.0, VALUE(I))
 10   CONTINUE
      K = 14
      IF (NP.EQ.4) K = 15
      CALL PGSCI(2)
      CALL PGSFS(1)
      CALL PGCIRC(X1, Y(K), 0.02)
      CALL PGUNSA
      CALL PGEBUF
C
C Cursor input.
C
 20   JUNK = PGCURS(XX, YY, CH)
      IF (ICHAR(CH).EQ.0) GOTO 50
C
C Find which box and highlight it
C
      DO 30 I=1,NBOX
         R = (XX-X1)**2 +(YY-Y(I))**2
         IF (R.LT.(0.03**2)) THEN
            ID = I
            CALL PGSAVE
            CALL PGSCI(2)
            CALL PGSFS(1)
            CALL PGCIRC(X1, Y(I), 0.02)
            CALL PGUNSA
            RETURN
         END IF
 30   CONTINUE
      GOTO 20
 50   ID = 0
      RETURN
      END


      SUBROUTINE PGEX1
C-----------------------------------------------------------------------
C This example illustrates the use of PGENV, PGLAB, PGPT, PGLINE.
C-----------------------------------------------------------------------
      INTEGER I
      REAL XS(5),YS(5), XR(100), YR(100)
      DATA XS/1.,2.,3.,4.,5./
      DATA YS/1.,4.,9.,16.,25./
C
C Call PGENV to specify the range of the axes and to draw a box, and
C PGLAB to label it. The x-axis runs from 0 to 10, and y from 0 to 20.
C
      CALL PGENV(0.,10.,0.,20.,0,1)
      CALL PGLAB('(x)', '(y)', 'PGPLOT Example 1:  y = x\u2')
C
C Mark five points (coordinates in arrays XS and YS), using symbol
C number 9.
C
      CALL PGPT(5,XS,YS,9)
C
C Compute the function at 60 points, and use PGLINE to draw it.
C
      DO 10 I=1,60
          XR(I) = 0.1*I
          YR(I) = XR(I)**2
   10 CONTINUE
      CALL PGLINE(60,XR,YR)
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX2
C-----------------------------------------------------------------------
C Repeat the process for another graph. This one is a graph of the
C sinc (sin x over x) function.
C-----------------------------------------------------------------------
      INTEGER I
      REAL XR(100), YR(100)
C
      CALL PGENV(-2.,10.,-0.4,1.2,0,1)
      CALL PGLAB('(x)', 'sin(x)/x', 
     $             'PGPLOT Example 2:  Sinc Function')
      DO 20 I=1,100
          XR(I) = (I-20)/6.
          YR(I) = 1.0
          IF (XR(I).NE.0.0) YR(I) = SIN(XR(I))/XR(I)
   20 CONTINUE
      CALL PGLINE(100,XR,YR)
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX3
C----------------------------------------------------------------------
C This example illustrates the use of PGBOX and attribute routines to
C mix colors and line-styles.
C----------------------------------------------------------------------
      REAL PI
      PARAMETER (PI=3.14159265359)
      INTEGER I
      REAL XR(360), YR(360)
      REAL ARG
C
C Call PGENV to initialize the viewport and window; the
C AXIS argument is -2, so no frame or labels will be drawn.
C
      CALL PGENV(0.,720.,-2.0,2.0,0,-2)
      CALL PGSAVE
C
C Set the color index for the axes and grid (index 5 = cyan).
C Call PGBOX to draw first a grid at low brightness, and then a
C frame and axes at full brightness. Note that as the x-axis is
C to represent an angle in degrees, we request an explicit tick 
C interval of 90 deg with subdivisions at 30 deg, as multiples of
C 3 are a more natural division than the default.
C
      CALL PGSCI(14)
      CALL PGBOX('G',30.0,0,'G',0.2,0)
      CALL PGSCI(5)
      CALL PGBOX('ABCTSN',90.0,3,'ABCTSNV',0.0,0)
C
C Call PGLAB to label the graph in a different color (3=green).
C
      CALL PGSCI(3)
      CALL PGLAB('x (degrees)','f(x)','PGPLOT Example 3')
C
C Compute the function to be plotted: a trig function of an
C angle in degrees, computed every 2 degrees.
C
      DO 20 I=1,360
          XR(I) = 2.0*I
          ARG = XR(I)/180.0*PI
          YR(I) = SIN(ARG) + 0.5*COS(2.0*ARG) + 
     1                0.5*SIN(1.5*ARG+PI/3.0)
   20 CONTINUE
C
C Change the color (6=magenta), line-style (2=dashed), and line
C width and draw the function.
C
      CALL PGSCI(6)
      CALL PGSLS(2)
      CALL PGSLW(3)
      CALL PGLINE(360,XR,YR)
C
C Restore attributes to defaults.
C
      CALL PGUNSA
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX4
C-----------------------------------------------------------------------
C Demonstration program for PGPLOT: draw histograms.
C-----------------------------------------------------------------------
      REAL     PI
      PARAMETER (PI=3.14159265359)
      INTEGER  I, ISEED
      REAL     DATA(1000), X(620), Y(620)
      REAL     PGRNRM
C
C Call PGRNRM to obtain 1000 samples from a normal distribution.
C
      ISEED = -5678921
      DO 10 I=1,1000
          DATA(I) = PGRNRM(ISEED)
   10 CONTINUE
C
C Draw a histogram of these values.
C
      CALL PGSAVE
      CALL PGHIST(1000,DATA,-3.1,3.1,31,0)
C
C Samples from another normal distribution.
C
      DO 15 I=1,200
          DATA(I) = 1.0+0.5*PGRNRM(ISEED)
   15 CONTINUE
C
C Draw another histogram (filled) on same axes.
C
      CALL PGSCI(15)
      CALL PGHIST(200,DATA,-3.1,3.1,31,3)
      CALL PGSCI(0)
      CALL PGHIST(200,DATA,-3.1,3.1,31,1)
      CALL PGSCI(1)
C
C Redraw the box which may have been clobbered by the histogram.
C
      CALL PGBOX('BST', 0.0, 0, ' ', 0.0, 0)
C
C Label the plot.
C
      CALL PGLAB('Variate', ' ',
     $             'PGPLOT Example 4:  Histograms (Gaussian)')
C
C Superimpose the theoretical distribution.
C
      DO 20 I=1,620
          X(I) = -3.1 + 0.01*(I-1)
          Y(I) = 0.2*1000./SQRT(2.0*PI)*EXP(-0.5*X(I)*X(I))
   20 CONTINUE
      CALL PGLINE(620,X,Y)
      CALL PGUNSA
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX5
C----------------------------------------------------------------------
C Demonstration program for the PGPLOT plotting package.  This example
C illustrates how to draw a log-log plot.
C PGPLOT subroutines demonstrated:
C    PGENV, PGERRY, PGLAB, PGLINE, PGPT, PGSCI.
C----------------------------------------------------------------------
      INTEGER   RED, GREEN, CYAN
      PARAMETER (RED=2)
      PARAMETER (GREEN=3)
      PARAMETER (CYAN=5)
      INTEGER   NP
      PARAMETER (NP=15)
      INTEGER   I
      REAL      X, YLO(NP), YHI(NP)
      REAL      FREQ(NP), FLUX(NP), XP(100), YP(100), ERR(NP)
      DATA FREQ / 26., 38., 80., 160., 178., 318.,
     1            365., 408., 750., 1400., 2695., 2700.,
     2            5000., 10695., 14900. /
      DATA FLUX / 38.0, 66.4, 89.0, 69.8, 55.9, 37.4,
     1            46.8, 42.4, 27.0, 15.8, 9.09, 9.17,
     2            5.35, 2.56, 1.73 /
      DATA ERR  / 6.0, 6.0, 13.0, 9.1, 2.9, 1.4,
     1            2.7, 3.0, 0.34, 0.8, 0.2, 0.46,
     2            0.15, 0.08, 0.01 /
C
C Call PGENV to initialize the viewport and window; the AXIS argument 
C is 30 so both axes will be logarithmic. The X-axis (frequency) runs 
C from 0.01 to 100 GHz, the Y-axis (flux density) runs from 0.3 to 300
C Jy. Note that it is necessary to specify the logarithms of these
C quantities in the call to PGENV. We request equal scales in x and y
C so that slopes will be correct.  Use PGLAB to label the graph.
C
      CALL PGSAVE
      CALL PGSCI(CYAN)
      CALL PGENV(-2.0,2.0,-0.5,2.5,1,30)
      CALL PGLAB('Frequency, \gn (GHz)',
     1             'Flux Density, S\d\gn\u (Jy)',
     2             'PGPLOT Example 5:  Log-Log plot')
C
C Draw a fit to the spectrum (don't ask how this was chosen). This 
C curve is drawn before the data points, so that the data will write 
C over the curve, rather than vice versa.
C
      DO 10 I=1,100
          X = 1.3 + I*0.03
          XP(I) = X-3.0
          YP(I) = 5.18 - 1.15*X -7.72*EXP(-X)
   10 CONTINUE
      CALL PGSCI(RED)
      CALL PGLINE(100,XP,YP)
C
C Plot the measured flux densities: here the data are installed with a
C DATA statement; in a more general program, they might be read from a
C file. We first have to take logarithms (the -3.0 converts MHz to GHz).
C
      DO 20 I=1,NP
          XP(I) = ALOG10(FREQ(I))-3.0
          YP(I) = ALOG10(FLUX(I))
   20 CONTINUE
      CALL PGSCI(GREEN)
      CALL PGPT(NP, XP, YP, 17)
C
C Draw +/- 2 sigma error bars: take logs of both limits.
C
      DO 30 I=1,NP
          YHI(I) = ALOG10(FLUX(I)+2.*ERR(I))
          YLO(I) = ALOG10(FLUX(I)-2.*ERR(I))
   30 CONTINUE
      CALL PGERRY(NP,XP,YLO,YHI,1.0)
      CALL PGUNSA
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX6
C----------------------------------------------------------------------
C Demonstration program for the PGPLOT plotting package.  This example
C illustrates the use of PGPOLY, PGCIRC, and PGRECT using SOLID, 
C OUTLINE, HATCHED, and CROSS-HATCHED fill-area attributes.
C----------------------------------------------------------------------
      REAL PI, TWOPI
      PARAMETER (PI=3.14159265359)
      PARAMETER (TWOPI=2.0*PI)
      INTEGER NPOL
      PARAMETER (NPOL=6)
      INTEGER I, J, N1(NPOL), N2(NPOL), K
      REAL X(10), Y(10), Y0, ANGLE
      CHARACTER*32 LAB(4)
      DATA N1 / 3, 4, 5, 5, 6, 8 /
      DATA N2 / 1, 1, 1, 2, 1, 3 /
      DATA LAB(1) /'Fill style 1 (solid)'/
      DATA LAB(2) /'Fill style 2 (outline)'/
      DATA LAB(3) /'Fill style 3 (hatched)'/
      DATA LAB(4) /'Fill style 4 (cross-hatched)'/
C
C Initialize the viewport and window.
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGPAGE
      CALL PGSVP(0.0, 1.0, 0.0, 1.0)
      CALL PGWNAD(0.0, 10.0, 0.0, 10.0)
C
C Label the graph.
C
      CALL PGSCI(1)
      CALL PGMTXT('T', -2.0, 0.5, 0.5, 
     :     'PGPLOT fill area: routines PGPOLY, PGCIRC, PGRECT')
C
C Draw assorted polygons.
C
      DO 30 K=1,4
         CALL PGSCI(1)
         Y0 = 10.0 - 2.0*K
         CALL PGTEXT(0.2, Y0+0.6, LAB(K))
         CALL PGSFS(K)
         DO 20 I=1,NPOL
            CALL PGSCI(I)
            DO 10 J=1,N1(I)
               ANGLE = REAL(N2(I))*TWOPI*REAL(J-1)/REAL(N1(I))
               X(J) = I + 0.5*COS(ANGLE)
               Y(J) = Y0 + 0.5*SIN(ANGLE)
 10         CONTINUE
            CALL PGPOLY (N1(I),X,Y)
 20      CONTINUE
         CALL PGSCI(7)
         CALL PGCIRC(7.0, Y0, 0.5)
         CALL PGSCI(8)
         CALL PGRECT(7.8, 9.5, Y0-0.5, Y0+0.5)
 30   CONTINUE
C
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX7
C-----------------------------------------------------------------------
C A plot with a large number of symbols; plus test of PGERR1.
C-----------------------------------------------------------------------
      INTEGER I, ISEED
      REAL XS(300),YS(300), XR(101), YR(101), XP, YP, XSIG, YSIG
      REAL PGRAND, PGRNRM
C
C Window and axes.
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGSCI(1)
      CALL PGENV(0.,5.,-0.3,0.6,0,1)
      CALL PGLAB('\fix', '\fiy', 'PGPLOT Example 7: scatter plot')
C
C Random data points.
C
      ISEED = -45678921
      DO 10 I=1,300
          XS(I) = 5.0*PGRAND(ISEED)
          YS(I) = XS(I)*EXP(-XS(I)) + 0.05*PGRNRM(ISEED)
   10 CONTINUE
      CALL PGSCI(3)
      CALL PGPT(100,XS,YS,3)
      CALL PGPT(100,XS(101),YS(101),17)
      CALL PGPT(100,XS(201),YS(201),21)
C
C Curve defining parent distribution.
C
      DO 20 I=1,101
          XR(I) = 0.05*(I-1)
          YR(I) = XR(I)*EXP(-XR(I))
   20 CONTINUE
      CALL PGSCI(2)
      CALL PGLINE(101,XR,YR)
C
C Test of PGERR1/PGPT1.
C
      XP = XS(101)
      YP = YS(101)
      XSIG = 0.2
      YSIG = 0.1
      CALL PGSCI(5)
      CALL PGSCH(3.0)
      CALL PGERR1(5, XP, YP, XSIG, 1.0)
      CALL PGERR1(6, XP, YP, YSIG, 1.0)
      CALL PGPT1(XP,YP,21)
C
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX8
C-----------------------------------------------------------------------
C Demonstration program for PGPLOT. This program shows some of the
C possibilities for overlapping windows and viewports.
C T. J. Pearson  1986 Nov 28
C-----------------------------------------------------------------------
      INTEGER I
      REAL XR(720), YR(720)
C-----------------------------------------------------------------------
C Color index:
      INTEGER BLACK, WHITE, RED, GREEN, BLUE, CYAN, MAGENT, YELLOW
      PARAMETER (BLACK=0)
      PARAMETER (WHITE=1)
      PARAMETER (RED=2)
      PARAMETER (GREEN=3)
      PARAMETER (BLUE=4)
      PARAMETER (CYAN=5)
      PARAMETER (MAGENT=6)
      PARAMETER (YELLOW=7)
C Line style:
      INTEGER FULL, DASHED, DOTDSH, DOTTED, FANCY
      PARAMETER (FULL=1)
      PARAMETER (DASHED=2)
      PARAMETER (DOTDSH=3)
      PARAMETER (DOTTED=4)
      PARAMETER (FANCY=5)
C Character font:
      INTEGER NORMAL, ROMAN, ITALIC, SCRIPT
      PARAMETER (NORMAL=1)
      PARAMETER (ROMAN=2)
      PARAMETER (ITALIC=3)
      PARAMETER (SCRIPT=4)
C Fill-area style:
      INTEGER SOLID, HOLLOW
      PARAMETER (SOLID=1)
      PARAMETER (HOLLOW=2)
C-----------------------------------------------------------------------
C
      CALL PGPAGE
      CALL PGBBUF
      CALL PGSAVE
C
C Define the Viewport
C
      CALL PGSVP(0.1,0.6,0.1,0.6)
C
C Define the Window
C
      CALL PGSWIN(0.0, 630.0, -2.0, 2.0)
C
C Draw a box
C
      CALL PGSCI(CYAN)
      CALL PGBOX ('ABCTS', 90.0, 3, 'ABCTSV', 0.0, 0)
C
C Draw labels
C
      CALL PGSCI (RED)
      CALL PGBOX ('N',90.0, 3, 'VN', 0.0, 0)
C
C Draw SIN line
C
      DO 10 I=1,360
          XR(I) = 2.0*I
          YR(I) = SIN(XR(I)/57.29577951)
   10 CONTINUE
      CALL PGSCI (MAGENT)
      CALL PGSLS (DASHED)
      CALL PGLINE (360,XR,YR)
C
C Draw COS line by redefining the window
C
      CALL PGSWIN (90.0, 720.0, -2.0, 2.0)
      CALL PGSCI (YELLOW)
      CALL PGSLS (DOTTED)
      CALL PGLINE (360,XR,YR)
      CALL PGSLS (FULL)
C
C Re-Define the Viewport
C
      CALL PGSVP(0.45,0.85,0.45,0.85)
C
C Define the Window, and erase it
C
      CALL PGSWIN(0.0, 180.0, -2.0, 2.0)
      CALL PGSCI(0)
      CALL PGRECT(0.0, 180., -2.0, 2.0)
C
C Draw a box
C
      CALL PGSCI(BLUE)
      CALL PGBOX ('ABCTSM', 60.0, 3, 'VABCTSM', 1.0, 2)
C
C Draw SIN line
C
      CALL PGSCI (WHITE)
      CALL PGSLS (DASHED)
      CALL PGLINE (360,XR,YR)
C
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX9
C----------------------------------------------------------------------
C Demonstration program for the PGPLOT plotting package.  This example
C illustrates curve drawing with PGFUNT; the parametric curve drawn is
C a simple Lissajous figure.
C                              T. J. Pearson  1983 Oct 5
C----------------------------------------------------------------------
      REAL PI
      PARAMETER (PI=3.14159265359)
      REAL     FX, FY
      EXTERNAL FX, FY
C
C Call PGFUNT to draw the function (autoscaling).
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGSCI(5)
      CALL PGFUNT(FX,FY,360,0.0,2.0*PI,0)
C
C Call PGLAB to label the graph in a different color.
C
      CALL PGSCI(3)
      CALL PGLAB('x','y','PGPLOT Example 9:  routine PGFUNT')
      CALL PGUNSA
      CALL PGEBUF
C
      END

      REAL FUNCTION FX(T)
      REAL T
      FX = SIN(T*5.0)
      RETURN
      END

      REAL FUNCTION FY(T)
      REAL T
      FY = SIN(T*4.0)
      RETURN
      END

      SUBROUTINE PGEX10
C----------------------------------------------------------------------
C Demonstration program for the PGPLOT plotting package.  This example
C illustrates curve drawing with PGFUNX.
C                              T. J. Pearson  1983 Oct 5
C----------------------------------------------------------------------
      REAL PI
      PARAMETER (PI=3.14159265359)
C The following define mnemonic names for the color indices and
C linestyle codes.
      INTEGER   BLACK, WHITE, RED, GREEN, BLUE, CYAN, MAGENT, YELLOW
      PARAMETER (BLACK=0)
      PARAMETER (WHITE=1)
      PARAMETER (RED=2)
      PARAMETER (GREEN=3)
      PARAMETER (BLUE=4)
      PARAMETER (CYAN=5)
      PARAMETER (MAGENT=6)
      PARAMETER (YELLOW=7)
      INTEGER   FULL, DASH, DOTD
      PARAMETER (FULL=1)
      PARAMETER (DASH=2)
      PARAMETER (DOTD=3)
C
C The Fortran functions to be plotted must be declared EXTERNAL.
C
      REAL     PGBSJ0, PGBSJ1
      EXTERNAL PGBSJ0, PGBSJ1
C
C Call PGFUNX twice to draw two functions (autoscaling the first time).
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGSCI(YELLOW)
      CALL PGFUNX(PGBSJ0,500,0.0,10.0*PI,0)
      CALL PGSCI(RED)
      CALL PGSLS(DASH)
      CALL PGFUNX(PGBSJ1,500,0.0,10.0*PI,1)
C
C Call PGLAB to label the graph in a different color. Note the
C use of "\f" to change font.  Use PGMTXT to write an additional
C legend inside the viewport.
C
      CALL PGSCI(GREEN)
      CALL PGSLS(FULL)
      CALL PGLAB('\fix', '\fiy',
     2           '\frPGPLOT Example 10: routine PGFUNX')
      CALL PGMTXT('T', -4.0, 0.5, 0.5,
     1     '\frBessel Functions')
C
C Call PGARRO to label the curves.
C
      CALL PGARRO(8.0, 0.7, 1.0, PGBSJ0(1.0))
      CALL PGARRO(12.0, 0.5, 9.0, PGBSJ1(9.0))
      CALL PGSTBG(GREEN)
      CALL PGSCI(0)
      CALL PGPTXT(8.0, 0.7, 0.0, 0.0, ' \fiy = J\d0\u(x)')
      CALL PGPTXT(12.0, 0.5, 0.0, 0.0, ' \fiy = J\d1\u(x)')
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX11
C-----------------------------------------------------------------------
C Test routine for PGPLOT: draws a skeletal dodecahedron.
C-----------------------------------------------------------------------
      INTEGER NVERT
      REAL T, T1, T2, T3
      PARAMETER (NVERT=20)
      PARAMETER (T=1.618)
      PARAMETER (T1=1.0+T)
      PARAMETER (T2=-1.0*T)
      PARAMETER (T3=-1.0*T1)
      INTEGER I, J, K
      REAL VERT(3,NVERT), R, ZZ
      REAL X(2),Y(2)
C
C Cartesian coordinates of the 20 vertices.
C
      DATA VERT/ T, T, T,       T, T,T2,
     3           T,T2, T,       T,T2,T2,
     5          T2, T, T,      T2, T,T2,
     7          T2,T2, T,      T2,T2,T2,
     9          T1,1.0,0.0,    T1,-1.0,0.0,
     B          T3,1.0,0.0,    T3,-1.0,0.0,
     D          0.0,T1,1.0,    0.0,T1,-1.0,
     F          0.0,T3,1.0,    0.0,T3,-1.0,
     H          1.0,0.0,T1,    -1.0,0.0,T1,
     J          1.0,0.0,T3,   -1.0,0.0,T3 /
C
C Initialize the plot (no labels).
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGENV(-4.,4.,-4.,4.,1,-2)
      CALL PGSCI(2)
      CALL PGSLS(1)
      CALL PGSLW(1)
C
C Write a heading.
C
      CALL PGLAB(' ',' ','PGPLOT Example 11:  Dodecahedron')
C
C Mark the vertices.
C
      DO 2 I=1,NVERT
          ZZ = VERT(3,I)
          CALL PGPT1(VERT(1,I)+0.2*ZZ,VERT(2,I)+0.3*ZZ,9)
    2 CONTINUE
C
C Draw the edges - test all vertex pairs to find the edges of the 
C correct length.
C
      CALL PGSLW(3)
      DO 20 I=2,NVERT
          DO 10 J=1,I-1
              R = 0.
              DO 5 K=1,3
                  R = R + (VERT(K,I)-VERT(K,J))**2
    5         CONTINUE
              R = SQRT(R)
              IF(ABS(R-2.0).GT.0.1) GOTO 10
              ZZ = VERT(3,I)
              X(1) = VERT(1,I)+0.2*ZZ
              Y(1) = VERT(2,I)+0.3*ZZ
              ZZ = VERT(3,J)
              X(2) = VERT(1,J)+0.2*ZZ
              Y(2) = VERT(2,J)+0.3*ZZ
              CALL PGLINE(2,X,Y)
   10     CONTINUE
   20 CONTINUE
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX12
C-----------------------------------------------------------------------
C Test routine for PGPLOT: draw arrows with PGARRO.
C-----------------------------------------------------------------------
      INTEGER NV, I, K
      REAL A, D, X, Y, XT, YT
C
C Number of arrows.
C
      NV =16
C
C Select a square viewport.
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGSCH(0.7)
      CALL PGSCI(2)
      CALL PGENV(-1.05,1.05,-1.05,1.05,1,-1)
      CALL PGLAB(' ', ' ', 'PGPLOT Example 12: PGARRO')
      CALL PGSCI(1)
C
C Draw the arrows
C
      K = 1
      D = 360.0/57.29577951/NV
      A = -D
      DO 20 I=1,NV
          A = A+D
          X = COS(A)
          Y = SIN(A)
          XT = 0.2*COS(A-D)
          YT = 0.2*SIN(A-D)
          CALL PGSAH(K, 80.0-3.0*I, 0.5*REAL(I)/REAL(NV))
          CALL PGSCH(0.25*I)
          CALL PGARRO(XT, YT, X, Y)
          K = K+1
          IF (K.GT.2) K=1
   20 CONTINUE
C
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX13
C----------------------------------------------------------------------
C This example illustrates the use of PGTBOX.
C----------------------------------------------------------------------
      INTEGER N
      PARAMETER (N=10)
      INTEGER I
      REAL X1(N), X2(N)
      CHARACTER*20 XOPT(N), BSL*1
      DATA X1 /   4*0.0, -8000.0, 100.3, 205.3, -45000.0, 2*0.0/
      DATA X2 /4*8000.0,  8000.0, 101.3, 201.1, 3*-100000.0/
      DATA XOPT / 'BSTN', 'BSTNZ', 'BSTNZH', 'BSTNZD', 'BSNTZHFO', 
     :      'BSTNZD', 'BSTNZHI', 'BSTNZHP', 'BSTNZDY', 'BSNTZHFOY'/
C
      BSL = CHAR(92)
      CALL PGPAGE
      CALL PGSAVE
      CALL PGBBUF
      CALL PGSCH(0.7)
      DO 100 I=1,N
        CALL PGSVP(0.15, 0.85, (0.7+REAL(N-I))/REAL(N), 
     :                         (0.7+REAL(N-I+1))/REAL(N)) 
        CALL PGSWIN(X1(I), X2(I), 0.0, 1.0)
        CALL PGTBOX(XOPT(I),0.0,0,' ',0.0,0)
        CALL PGLAB('Option = '//XOPT(I), ' ', ' ')
        IF (I.EQ.1) THEN
           CALL PGMTXT('B', -1.0, 0.5, 0.5, 
     :                 BSL//'fiAxes drawn with PGTBOX')
        END IF
  100 CONTINUE
      CALL PGEBUF
      CALL PGUNSA
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX14
C-----------------------------------------------------------------------
C Test routine for PGPLOT: polygon fill and color representation.
C-----------------------------------------------------------------------
      INTEGER I, J, N, M
      REAL PI, THINC, R, G, B, THETA
      REAL XI(100),YI(100),XO(100),YO(100),XT(3),YT(3)
      PARAMETER (PI=3.14159265359)
C
      N = 33
      M = 8
      THINC=2.0*PI/N
      DO 10 I=1,N
        XI(I) = 0.0
        YI(I) = 0.0
   10 CONTINUE
      CALL PGBBUF
      CALL PGSAVE
      CALL PGENV(-1.,1.,-1.,1.,1,-2)
      CALL PGLAB(' ', ' ', 'PGPLOT Example 14: PGPOLY and PGSCR')
      DO 50 J=1,M
        R = 1.0
        G = 1.0 - REAL(J)/REAL(M)
        B = G
        CALL PGSCR(J, R, G, B)
        THETA = -REAL(J)*PI/REAL(N)
        R = REAL(J)/REAL(M)
        DO 20 I=1,N
          THETA = THETA+THINC
          XO(I) = R*COS(THETA)
          YO(I) = R*SIN(THETA)
   20   CONTINUE
        DO 30 I=1,N
          XT(1) = XO(I)
          YT(1) = YO(I)
          XT(2) = XO(MOD(I,N)+1)
          YT(2) = YO(MOD(I,N)+1)
          XT(3) = XI(I)
          YT(3) = YI(I)
          CALL PGSCI(J)
          CALL PGSFS(1)
          CALL PGPOLY(3,XT,YT)
          CALL PGSFS(2)
          CALL PGSCI(1)
          CALL PGPOLY(3,XT,YT)
   30   CONTINUE
        DO 40 I=1,N
          XI(I) = XO(I)
          YI(I) = YO(I)
   40   CONTINUE
   50 CONTINUE
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      SUBROUTINE PGEX15
C----------------------------------------------------------------------
C This is a line-drawing test; it draws a regular n-gon joining
C each vertex to every other vertex. It is not optimized for pen
C plotters.
C----------------------------------------------------------------------
      INTEGER I, J, NV
      REAL A, D, X(100), Y(100)
C
C Set the number of vertices, and compute the 
C coordinates for unit circumradius.
C
      NV = 17
      D = 360.0/NV
      A = -D
      DO 20 I=1,NV
          A = A+D
          X(I) = COS(A/57.29577951)
          Y(I) = SIN(A/57.29577951)
   20 CONTINUE
C
C Select a square viewport.
C
      CALL PGBBUF
      CALL PGSAVE
      CALL PGSCH(0.5)
      CALL PGSCI(2)
      CALL PGENV(-1.05,1.05,-1.05,1.05,1,-1)
      CALL PGLAB(' ', ' ', 'PGPLOT Example 15: PGMOVE and PGDRAW')
      CALL PGSCR(0,0.2,0.3,0.3)
      CALL PGSCR(1,1.0,0.5,0.2)
      CALL PGSCR(2,0.2,0.5,1.0)
      CALL PGSCI(1)
C
C Draw the polygon.
C
      DO 40 I=1,NV-1
          DO 30 J=I+1,NV
            CALL PGMOVE(X(I),Y(I))
            CALL PGDRAW(X(J),Y(J))
   30     CONTINUE
   40 CONTINUE
C
C Flush the buffer.
C
      CALL PGUNSA
      CALL PGEBUF
C-----------------------------------------------------------------------
      END

      REAL FUNCTION PGBSJ0(XX)
      REAL XX
C-----------------------------------------------------------------------
C Bessel function of order 0 (approximate).
C Reference: Abramowitz and Stegun: Handbook of Mathematical Functions.
C-----------------------------------------------------------------------
      REAL X, XO3, T, F0, THETA0
C     
      X = ABS(XX)
      IF (X .LE. 3.0) THEN
         XO3 = X/3.0
         T   = XO3*XO3
         PGBSJ0 = 1.0 + T*(-2.2499997 +
     1                  T*( 1.2656208 +
     2                  T*(-0.3163866 +
     3                  T*( 0.0444479 +
     4                  T*(-0.0039444 +
     5                  T*( 0.0002100))))))
      ELSE
         T = 3.0/X
         F0 =     0.79788456 +
     1        T*(-0.00000077 + 
     2        T*(-0.00552740 +
     3        T*(-0.00009512 +
     4        T*( 0.00137237 +
     5        T*(-0.00072805 +
     6        T*( 0.00014476))))))
         THETA0 = X - 0.78539816 +
     1            T*(-0.04166397 +
     2            T*(-0.00003954 +
     3            T*( 0.00262573 +
     4            T*(-0.00054125 +
     5            T*(-0.00029333 +
     6            T*( 0.00013558))))))
         PGBSJ0 = F0*COS(THETA0)/SQRT(X)
      END IF
C-----------------------------------------------------------------------
      END

      REAL FUNCTION PGBSJ1(XX)
      REAL XX
C-----------------------------------------------------------------------
C Bessel function of order 1 (approximate).
C Reference: Abramowitz and Stegun: Handbook of Mathematical Functions.
C-----------------------------------------------------------------------
      REAL X, XO3, T, F1, THETA1
C
      X = ABS(XX)
      IF (X .LE. 3.0) THEN
         XO3 = X/3.0
         T = XO3*XO3
         PGBSJ1 = 0.5 + T*(-0.56249985 +
     1                  T*( 0.21093573 +
     2                  T*(-0.03954289 +
     3                  T*( 0.00443319 +
     4                  T*(-0.00031761 +
     5                  T*( 0.00001109))))))
         PGBSJ1 = PGBSJ1*XX
      ELSE
         T = 3.0/X
         F1 =    0.79788456 +
     1       T*( 0.00000156 +
     2       T*( 0.01659667 + 
     3       T*( 0.00017105 +
     4       T*(-0.00249511 +
     5       T*( 0.00113653 + 
     6       T*(-0.00020033))))))
         THETA1 = X   -2.35619449 + 
     1             T*( 0.12499612 +
     2             T*( 0.00005650 +
     3             T*(-0.00637879 +
     4             T*( 0.00074348 +
     5             T*( 0.00079824 +
     6             T*(-0.00029166))))))
         PGBSJ1 = F1*COS(THETA1)/SQRT(X)
      END IF
      IF (XX .LT. 0.0) PGBSJ1 = -PGBSJ1
C-----------------------------------------------------------------------
      END

      REAL FUNCTION PGRNRM (ISEED)
      INTEGER ISEED
C-----------------------------------------------------------------------
C Returns a normally distributed deviate with zero mean and unit 
C variance. The routine uses the Box-Muller transformation of uniform
C deviates. For a more efficient implementation of this algorithm,
C see Press et al., Numerical Recipes, Sec. 7.2.
C
C Arguments:
C  ISEED  (in/out) : seed used for PGRAND random-number generator.
C
C Subroutines required:
C  PGRAND -- return a uniform random deviate between 0 and 1.
C
C History:
C  1995 Dec 12 - TJP.
C-----------------------------------------------------------------------
      REAL R, X, Y, PGRAND
C
 10   X = 2.0*PGRAND(ISEED) - 1.0
      Y = 2.0*PGRAND(ISEED) - 1.0
      R = X**2 + Y**2
      IF (R.GE.1.0) GOTO 10
      PGRNRM = X*SQRT(-2.0*LOG(R)/R)
C-----------------------------------------------------------------------
      END

      REAL FUNCTION PGRAND(ISEED)
      INTEGER ISEED
C-----------------------------------------------------------------------
C Returns a uniform random deviate between 0.0 and 1.0.
C
C NOTE: this is not a good random-number generator; it is only
C intended for exercising the PGPLOT routines.
C
C Based on: Park and Miller's "Minimal Standard" random number
C   generator (Comm. ACM, 31, 1192, 1988)
C
C Arguments:
C  ISEED  (in/out) : seed.
C-----------------------------------------------------------------------
      INTEGER   IM, IA, IQ, IR
      PARAMETER (IM=2147483647)
      PARAMETER (IA=16807, IQ=127773, IR= 2836)
      REAL      AM
      PARAMETER (AM=128.0/IM)
      INTEGER   K
C-
      K = ISEED/IQ
      ISEED = IA*(ISEED-K*IQ) - IR*K
      IF (ISEED.LT.0) ISEED = ISEED+IM
      PGRAND = AM*(ISEED/128)
      RETURN
      END