PROGRAM 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)
      INTEGER   PGOPEN
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 PGOPEN to initiate PGPLOT and open the output device; PGOPEN
C will prompt the user to supply the device name and type. Always
C check the return code from PGOPEN.
C
      IF (PGOPEN('?') .LE. 0) STOP
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
C Finally, call PGCLOS to terminate things properly.
C
      CALL PGCLOS
C
      END


C-----------------------------------------------------------------------


      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


C-----------------------------------------------------------------------


      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