      SUBROUTINE BESJ(X,ALPHA,N,Y,NZ) 
C***BEGIN PROLOGUE  BESJ
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C10A3
C***KEYWORDS  BESSEL FUNCTION,J BESSEL FUNCTION,SPECIAL FUNCTION
C***AUTHOR  AMOS, D. E., (SNLA)
C           DANIEL, S. L., (SNLA)
C           WESTON, M. K., (SNLA)
C***PURPOSE  BESJ computes an N member sequence of J Bessel functions 
C            J/SUB(ALPHA+K-1)/(X), K=1,..,N for non-negative ALPHA and X 
C***DESCRIPTION
C
C     Written by D. E. Amos, S. L. Daniel and M. K. Weston, January 1975
C
C     References
C         SAND-75-0147
C
C         CDC 6600 subroutines IBESS and JBESS for Bessel functions
C         I(NU,X) and J(NU,X), X .GE. 0, NU .GE. 0  by D. E. Amos, S. L.
C         Daniel, M. K. Weston,  ACM Transactions on Mathematical 
C         Software, Vol. 3, pp. 76-92 (1977).
C
C         Tables of Bessel Functions of Moderate or Large Orders,
C         NPL Mathematical Tables, Vol. 6, by F. W. J. Olver, Her
C         Majesty's Stationery Office, London, 1962.
C
C     Abstract
C         BESJ computes an N member sequence of J bessel functions
C         J/sub(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA and X.
C         A combination of the power series, the asymptotic expansion
C         for X to infinity and the uniform asymptotic expansion for
C         NU to infinity are applied over subdivisions of the (NU,X) 
C         plane.  For values of (NU,X) not covered by one of these
C         formulae, the order is incremented or decremented by integer
C         values into a region where one of the formulae apply. Backward
C         recursion is applied to reduce orders by integer values except
C         where the entire sequence lies in the oscillatory region.  In
C         this case forward recursion is stable and values from the
C         asymptotic expansion for X to infinity start the recursion
C         when it is efficient to do so.  Leading terms of the series
C         and uniform expansion are tested for underflow.  If a sequence
C         is requested and the last member would underflow, the result
C         is set to zero and the next lower order tried, etc., until a
C         member comes on scale or all members are set to zero.
C         Overflow cannot occur.
C
C         BESJ calls ASYJY, JAIRY, ALNGAM, R1MACH, I1MACH, XERROR
C
C     Description of Arguments
C
C         Input 
C           X      - X .GE. 0.0E0
C           ALPHA  - order of first member of the sequence,
C                    ALPHA .GE. 0.0E0
C           N      - number of members in the sequence, N .GE. 1
C
C         Output
C           Y      - a vector whose first  N components contain
C                    values for J/sub(ALPHA+K-1)/(X), K=1,...,N
C           NZ     - number of components of Y set to zero due to 
C                    underflow,
C                    NZ=0   , normal return, computation completed
C                    NZ .NE. 0, last NZ components of Y set to zero,
C                             Y(K)=0.0E0, K=N-NZ+1,...,N.
C
C     Error Conditions
C         Improper input arguments - a fatal error
C         Underflow  - a non-fatal error (NZ .NE. 0)
C         NOTE: BESJ assumes that an underflow does NOT generate
C              an abort. This is compiler dependent. Some compilers
C	       abort on underflow as a default, and require  special
C              user action at compile time or the insertion of a
C              special CALL to alter this.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  ALNGAM,ASYJY,I1MACH,JAIRY,R1MACH,XERROR
C***END PROLOGUE  BESJ
C
      EXTERNAL JAIRY
      INTEGER I,IALP,IDALP,IFLW,IN,INLIM,IS,I1,I2,K,KK,KM,KT,N,NN,
     1        NS,NZ 
      INTEGER I1MACH
      REAL       AK,AKM,ALPHA,ANS,AP,ARG,COEF,DALPHA,DFN,DTM,EARG,
     1           ELIM1,ETX,FIDAL,FLGJY,FN,FNF,FNI,FNP1,FNU,FNULIM,
     2           GLN,PDF,PIDT,PP,RDEN,RELB,RTTP,RTWO,RTX,RZDEN,
     3           S,SA,SB,SXO2,S1,S2,T,TA,TAU,TB,TEMP,TFN,TM,TOL,
     4           TOLLN,TRX,TX,T1,T2,WK,X,XO2,XO2L,Y
      REAL R1MACH, ALNGAM
      DIMENSION Y(1), TEMP(3), FNULIM(2), PP(4), WK(7)
      DATA RTWO,PDF,RTTP,PIDT                    / 1.34839972492648E+00,
     1 7.85398163397448E-01, 7.97884560802865E-01, 1.57079632679490E+00/
      DATA  PP(1),  PP(2),  PP(3),  PP(4)        / 8.72909153935547E+00,
     1 2.65693932265030E-01, 1.24578576865586E-01, 7.70133747430388E-04/
      DATA INLIM           /      150            /
      DATA FNULIM(1), FNULIM(2) /      100.0E0,     60.0E0     /
C***FIRST EXECUTABLE STATEMENT  BESJ
      NZ = 0
      KT = 1
C     I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE
C     I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE
      TA = R1MACH(3)
      TOL = AMAX1(TA,1.0E-15) 
      I1 = I1MACH(11) + 1
      I2 = I1MACH(12)
      TB = R1MACH(5)
      ELIM1 = 2.303E0*(FLOAT(-I2)*TB-3.0E0)
C     TOLLN = -LN(TOL)
      TOLLN = 2.303E0*TB*FLOAT(I1)
      TOLLN = AMIN1(TOLLN,34.5388E0)
      IF (N-1) 720, 10, 20
   10 KT = 2
   20 NN = N
      IF (X) 730, 30, 80
   30 IF (ALPHA) 710, 40, 50
   40 Y(1) = 1.0E0
      IF (N.EQ.1) RETURN
      I1 = 2
      GO TO 60
   50 I1 = 1
   60 DO 70 I=I1,N
        Y(I) = 0.0E0
   70 CONTINUE
      RETURN
   80 CONTINUE
      IF (ALPHA.LT.0.0E0) GO TO 710
C
      IALP = INT(ALPHA)
      FNI = FLOAT(IALP+N-1)
      FNF = ALPHA - FLOAT(IALP)
      DFN = FNI + FNF
      FNU = DFN
      XO2 = X*0.5E0 
      SXO2 = XO2*XO2
C
C     DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
C     TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
C     APPLIED.
C
      IF (SXO2.LE.(FNU+1.0E0)) GO TO 90 
      TA = AMAX1(20.0E0,FNU)
      IF (X.GT.TA) GO TO 120
      IF (X.GT.12.0E0) GO TO 110
      XO2L = ALOG(XO2)
      NS = INT(SXO2-FNU) + 1
      GO TO 100
   90 FN = FNU
      FNP1 = FN + 1.0E0
      XO2L = ALOG(XO2)
      IS = KT
      IF (X.LE.0.50E0) GO TO 330
      NS = 0
  100 FNI = FNI + FLOAT(NS)
      DFN = FNI + FNF
      FN = DFN
      FNP1 = FN + 1.0E0
      IS = KT
      IF (N-1+NS.GT.0) IS = 3 
      GO TO 330
  110 ANS = AMAX1(36.0E0-FNU,0.0E0)
      NS = INT(ANS) 
      FNI = FNI + FLOAT(NS)
      DFN = FNI + FNF
      FN = DFN
      IS = KT
      IF (N-1+NS.GT.0) IS = 3 
      GO TO 130
  120 CONTINUE
      RTX = SQRT(X) 
      TAU = RTWO*RTX
      TA = TAU + FNULIM(KT)
      IF (FNU.LE.TA) GO TO 480
      FN = FNU
      IS = KT
C
C     UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
C
  130 CONTINUE
      I1 = IABS(3-IS)
      I1 = MAX0(I1,1)
      FLGJY = 1.0E0 
      CALL ASYJY(JAIRY,X,FN,FLGJY,I1,TEMP(IS),WK,IFLW)
      IF(IFLW.NE.0) GO TO 380 
      GO TO (320, 450, 620), IS
  310 TEMP(1) = TEMP(3)
      KT = 1
  320 IS = 2
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IF(I1.EQ.2) GO TO 450
      GO TO 130
C
C     SERIES FOR (X/2)**2.LE.NU+1
C
  330 CONTINUE
      GLN = ALNGAM(FNP1)
      ARG = FN*XO2L - GLN
      IF (ARG.LT.(-ELIM1)) GO TO 400
      EARG = EXP(ARG)
  340 CONTINUE
      S = 1.0E0
      IF (X.LT.TOL) GO TO 360 
      AK = 3.0E0
      T2 = 1.0E0
      T = 1.0E0
      S1 = FN
      DO 350 K=1,17 
        S2 = T2 + S1
        T = -T*SXO2/S2
        S = S + T
        IF (ABS(T).LT.TOL) GO TO 360
        T2 = T2 + AK
        AK = AK + 2.0E0
        S1 = S1 + FN
  350 CONTINUE
  360 CONTINUE
      TEMP(IS) = S*EARG
      GO TO (370, 450, 610), IS
  370 EARG = EARG*FN/XO2
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IS = 2
      GO TO 340
C
C     SET UNDERFLOW VALUE AND UPDATE PARAMETERS
C
  380 Y(NN) = 0.0E0 
      NN = NN - 1
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IF (NN-1) 440, 390, 130 
  390 KT = 2
      IS = 2
      GO TO 130
  400 Y(NN) = 0.0E0 
      NN = NN - 1
      FNP1 = FN
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IF (NN-1) 440, 410, 420 
  410 KT = 2
      IS = 2
  420 IF (SXO2.LE.FNP1) GO TO 430
      GO TO 130
  430 ARG = ARG - XO2L + ALOG(FNP1)
      IF (ARG.LT.(-ELIM1)) GO TO 400
      GO TO 330
  440 NZ = N - NN
      RETURN
C
C     BACKWARD RECURSION SECTION
C
  450 CONTINUE
      NZ = N - NN
      IF (KT.EQ.2) GO TO 470
C     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
      Y(NN) = TEMP(1)
      Y(NN-1) = TEMP(2)
      IF (NN.EQ.2) RETURN
      TRX = 2.0E0/X 
      DTM = FNI
      TM = (DTM+FNF)*TRX
      K = NN + 1
      DO 460 I=3,NN 
        K = K - 1
        Y(K-2) = TM*Y(K-1) - Y(K)
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
  460 CONTINUE
      RETURN
  470 Y(1) = TEMP(2)
      RETURN
C
C     ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN
C     OSCILLATORY REGION X.GT.MAX(20, NU), PROVIDED THE LAST MEMBER
C     OF THE SEQUENCE IS ALSO IN THE REGION.
C
  480 CONTINUE
      IN = INT(ALPHA-TAU+2.0E0)
      IF (IN.LE.0) GO TO 490
      IDALP = IALP - IN - 1
      KT = 1
      GO TO 500
  490 CONTINUE
      IDALP = IALP
      IN = 0
  500 IS = KT
      FIDAL = FLOAT(IDALP)
      DALPHA = FIDAL + FNF
      ARG = X - PIDT*DALPHA - PDF
      SA = SIN(ARG) 
      SB = COS(ARG) 
      COEF = RTTP/RTX
      ETX = 8.0E0*X 
  510 CONTINUE
      DTM = FIDAL + FIDAL
      DTM = DTM*DTM 
      TM = 0.0E0
      IF (FIDAL.EQ.0.0E0 .AND. ABS(FNF).LT.TOL) GO TO 520
      TM = 4.0E0*FNF*(FIDAL+FIDAL+FNF)
  520 CONTINUE
      TRX = DTM - 1.0E0
      T2 = (TRX+TM)/ETX
      S2 = T2
      RELB = TOL*ABS(T2)
      T1 = ETX
      S1 = 1.0E0
      FN = 1.0E0
      AK = 8.0E0
      DO 530 K=1,13 
        T1 = T1 + ETX
        FN = FN + AK
        TRX = DTM - FN
        AP = TRX + TM
        T2 = -T2*AP/T1
        S1 = S1 + T2
        T1 = T1 + ETX
        AK = AK + 8.0E0
        FN = FN + AK
        TRX = DTM - FN
        AP = TRX + TM
        T2 = T2*AP/T1
        S2 = S2 + T2
        IF (ABS(T2).LE.RELB) GO TO 540
        AK = AK + 8.0E0
  530 CONTINUE
  540 TEMP(IS) = COEF*(S1*SB-S2*SA)
      IF(IS.EQ.2) GO TO 560
  550 FIDAL = FIDAL + 1.0E0
      DALPHA = FIDAL + FNF
      IS = 2
      TB = SA
      SA = -SB
      SB = TB
      GO TO 510
C
C     FORWARD RECURSION SECTION
C
  560 IF (KT.EQ.2) GO TO 470
      S1 = TEMP(1)
      S2 = TEMP(2)
      TX = 2.0E0/X
      TM = DALPHA*TX
      IF (IN.EQ.0) GO TO 580
C
C     FORWARD RECUR TO INDEX ALPHA
C
      DO 570 I=1,IN 
        S = S2
        S2 = TM*S2 - S1
        TM = TM + TX
        S1 = S
  570 CONTINUE
      IF (NN.EQ.1) GO TO 600
      S = S2
      S2 = TM*S2 - S1
      TM = TM + TX
      S1 = S
  580 CONTINUE
C
C     FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1 
C
      Y(1) = S1
      Y(2) = S2
      IF (NN.EQ.2) RETURN
      DO 590 I=3,NN 
        Y(I) = TM*Y(I-1) - Y(I-2)
        TM = TM + TX
  590 CONTINUE
      RETURN
  600 Y(1) = S2
      RETURN
C
C     BACKWARD RECURSION WITH NORMALIZATION BY
C     ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
C
  610 CONTINUE
C     COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
      AKM = AMAX1(3.0E0-FN,0.0E0)
      KM = INT(AKM) 
      TFN = FN + FLOAT(KM)
      TA = (GLN+TFN-0.9189385332E0-0.0833333333E0/TFN)/(TFN+0.5E0)
      TA = XO2L - TA
      TB = -(1.0E0-1.5E0/TFN)/TFN
      AKM = TOLLN/(-TA+SQRT(TA*TA-TOLLN*TB)) + 1.5E0
      IN = KM + INT(AKM)
      GO TO 660
  620 CONTINUE
C     COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
      GLN = WK(3) + WK(2)
      IF (WK(6).GT.30.0E0) GO TO 640
      RDEN = (PP(4)*WK(6)+PP(3))*WK(6) + 1.0E0
      RZDEN = PP(1) + PP(2)*WK(6)
      TA = RZDEN/RDEN
      IF (WK(1).LT.0.10E0) GO TO 630
      TB = GLN/WK(5)
      GO TO 650
  630 TB=(1.259921049E0+(0.1679894730E0+0.0887944358E0*WK(1))*WK(1))
     1 /WK(7)
      GO TO 650
  640 CONTINUE
      TA = 0.5E0*TOLLN/WK(4)
      TA=((0.0493827160E0*TA-0.1111111111E0)*TA+0.6666666667E0)*TA*WK(6)
      IF (WK(1).LT.0.10E0) GO TO 630
      TB = GLN/WK(5)
  650 IN = INT(TA/TB+1.5E0)
      IF (IN.GT.INLIM) GO TO 310
  660 CONTINUE
      DTM = FNI + FLOAT(IN)
      TRX = 2.0E0/X 
      TM = (DTM+FNF)*TRX
      TA = 0.0E0
      TB = TOL
      KK = 1
  670 CONTINUE
C
C     BACKWARD RECUR UNINDEXED
C
      DO 680 I=1,IN 
        S = TB
        TB = TM*TB - TA
        TA = S
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
  680 CONTINUE
C     NORMALIZATION 
      IF (KK.NE.1) GO TO 690
      TA = (TA/TB)*TEMP(3)
      TB = TEMP(3)
      KK = 2
      IN = NS
      IF (NS.NE.0) GO TO 670
  690 Y(NN) = TB
      NZ = N - NN
      IF (NN.EQ.1) RETURN
      K = NN - 1
      Y(K) = TM*TB - TA
      IF (NN.EQ.2) RETURN
      DTM = DTM - 1.0E0
      TM = (DTM+FNF)*TRX
      KM = K - 1
C
C     BACKWARD RECUR INDEXED
C
      DO 700 I=1,KM 
        Y(K-1) = TM*Y(K) - Y(K+1)
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
        K = K - 1
  700 CONTINUE
      RETURN
C
C
C
  710 CONTINUE
      CALL XERROR( 'BESJ - ORDER, ALPHA, LESS THAN ZERO.', 36, 2, 1)
      RETURN
  720 CONTINUE
      CALL XERROR( 'BESJ - N LESS THAN ONE.', 23, 2, 1)
      RETURN
  730 CONTINUE
      CALL XERROR( 'BESJ - X LESS THAN ZERO.', 24, 2, 1)
      RETURN
      END 
      SUBROUTINE ASYJY(FUNJY,X,FNU,FLGJY,IN,Y,WK,IFLW)
C***BEGIN PROLOGUE  ASYJY
C***REFER TO  BESJ,BESY
C
C                 ASYJY computes Bessel functions J and Y
C               for arguments X.GT.0.0 and orders FNU.GE.35.0
C               on FLGJY = 1 and FLGJY = -1 respectively
C
C                                  INPUT
C
C      FUNJY - external function JAIRY or YAIRY
C          X - argument, X.GT.0.0E0
C        FNU - order of the first Bessel function
C      FLGJY - selection flag
C              FLGJY =  1.0E0 gives the J function
C              FLGJY = -1.0E0 gives the Y function
C         IN - number of functions desired, IN = 1 or 2
C
C                                  OUTPUT
C
C         Y  - a vector whose first in components contain the sequence
C       IFLW - a flag indicating underflow or overflow
C                    return variables for BESJ only
C      WK(1) = 1 - (X/FNU)**2 = W**2
C      WK(2) = SQRT(ABS(WK(1)))
C      WK(3) = ABS(WK(2) - ATAN(WK(2)))  or
C              ABS(LN((1 + WK(2))/(X/FNU)) - WK(2))
C            = ABS((2/3)*ZETA**(3/2))
C      WK(4) = FNU*WK(3)
C      WK(5) = (1.5*WK(3)*FNU)**(1/3) = SQRT(ZETA)*FNU**(1/3)
C      WK(6) = SIGN(1.,W**2)*WK(5)**2 = SIGN(1.,W**2)*ZETA*FNU**(2/3) 
C      WK(7) = FNU**(1/3)
C
C                                  Written by
C                                  D. E. Amos
C
C     Abstract
C         ASYJY implements the uniform asymptotic expansion of
C         the J and Y Bessel functions for FNU.GE.35 and real
C         X.GT.0.0E0. The forms are identical except for a change
C         in sign of some of the terms. This change in sign is
C         accomplished by means of the flag FLGJY = 1 or -1. On
C         FLGJY = 1 the AIRY functions AI(X) and DAI(X) are
C         supplied by the external function JAIRY, and on
C         FLGJY = -1 the AIRY functions BI(X) and DBI(X) are
C         supplied by the external funtion YAIRY.
C***ROUTINES CALLED  I1MACH,R1MACH
C***END PROLOGUE  ASYJY
C
      INTEGER I, IFLW, IN, J, JN,JR,JU,K, KB,KLAST,KMAX,KP1, KS, KSP1,
     1 KSTEMP, L, LR, LRP1
      INTEGER I1MACH
      REAL ABW2, AKM, ALFA, ALFA1, ALFA2, AP, AR, ASUM, AZ, 
     1 BETA, BETA1, BETA2, BETA3, BR, BSUM, C, CON1, CON2,
     2 CON548,CR,CRZ32, DFI,ELIM, DR,FI, FLGJY, FN, FNU,
     3 FN2, GAMA, PHI,  RCZ, RDEN, RELB, RFN2,  RTZ, RZDEN, 
     4 SA, SB, SUMA, SUMB, S1, TA, TAU, TB, TFN, TOL, TOLS, T2, UPOL, 
     5  WK, X, XX, Y, Z, Z32
      REAL R1MACH
      DIMENSION Y(1), WK(1), C(65)
      DIMENSION ALFA(26,4), BETA(26,5)
      DIMENSION ALFA1(26,2), ALFA2(26,2)
      DIMENSION BETA1(26,2), BETA2(26,2), BETA3(26,1)
      DIMENSION GAMA(26), KMAX(5), AR(8), BR(10), UPOL(10)
      DIMENSION CR(10), DR(10)
      EQUIVALENCE (ALFA(1,1),ALFA1(1,1))
      EQUIVALENCE (ALFA(1,3),ALFA2(1,1))
      EQUIVALENCE (BETA(1,1),BETA1(1,1))
      EQUIVALENCE (BETA(1,3),BETA2(1,1))
      EQUIVALENCE (BETA(1,5),BETA3(1,1))
      DATA TOLS            /-6.90775527898214E+00/
      DATA CON1,CON2,CON548/
     1 6.66666666666667E-01,3.33333333333333E-01,1.04166666666667E-01/
      DATA  AR(1),  AR(2),  AR(3),  AR(4),  AR(5),  AR(6),  AR(7),
     1      AR(8)          / 8.35503472222222E-02, 1.28226574556327E-01,
     2 2.91849026464140E-01, 8.81627267443758E-01, 3.32140828186277E+00,
     3 1.49957629868626E+01, 7.89230130115865E+01, 4.74451538868264E+02/
      DATA  BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8),
     1      BR(9), BR(10)  /-1.45833333333333E-01,-9.87413194444444E-02,
     2-1.43312053915895E-01,-3.17227202678414E-01,-9.42429147957120E-01,
     3-3.51120304082635E+00,-1.57272636203680E+01,-8.22814390971859E+01,
     4-4.92355370523671E+02,-3.31621856854797E+03/
      DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
     1     C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
     2     C(19), C(20), C(21), C(22), C(23), C(24)/
     3       -2.08333333333333E-01,        1.25000000000000E-01,
     4        3.34201388888889E-01,       -4.01041666666667E-01,
     5        7.03125000000000E-02,       -1.02581259645062E+00,
     6        1.84646267361111E+00,       -8.91210937500000E-01,
     7        7.32421875000000E-02,        4.66958442342625E+00,
     8       -1.12070026162230E+01,        8.78912353515625E+00,
     9       -2.36408691406250E+00,        1.12152099609375E-01,
     1       -2.82120725582002E+01,        8.46362176746007E+01,
     2       -9.18182415432400E+01,        4.25349987453885E+01,
     3       -7.36879435947963E+00,        2.27108001708984E-01,
     4        2.12570130039217E+02,       -7.65252468141182E+02,
     5        1.05999045252800E+03,       -6.99579627376133E+02/
      DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
     1     C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
     2     C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
     3        2.18190511744212E+02,       -2.64914304869516E+01,
     4        5.72501420974731E-01,       -1.91945766231841E+03,
     5        8.06172218173731E+03,       -1.35865500064341E+04,
     6        1.16553933368645E+04,       -5.30564697861340E+03,
     7        1.20090291321635E+03,       -1.08090919788395E+02,
     8        1.72772750258446E+00,        2.02042913309661E+04,
     9       -9.69805983886375E+04,        1.92547001232532E+05,
     1       -2.03400177280416E+05,        1.22200464983017E+05,
     2       -4.11926549688976E+04,        7.10951430248936E+03,
     3       -4.93915304773088E+02,        6.07404200127348E+00,
     4       -2.42919187900551E+05,        1.31176361466298E+06,
     5       -2.99801591853811E+06,        3.76327129765640E+06/
      DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
     1     C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
     2     C(65)/
     3       -2.81356322658653E+06,        1.26836527332162E+06,
     4       -3.31645172484564E+05,        4.52187689813627E+04,
     5       -2.49983048181121E+03,        2.43805296995561E+01,
     6        3.28446985307204E+06,       -1.97068191184322E+07,
     7        5.09526024926646E+07,       -7.41051482115327E+07,
     8        6.63445122747290E+07,       -3.75671766607634E+07,
     9        1.32887671664218E+07,       -2.78561812808645E+06,
     1        3.08186404612662E+05,       -1.38860897537170E+04,
     2        1.10017140269247E+02/
      DATA ALFA1(1,1), ALFA1(2,1), ALFA1(3,1), ALFA1(4,1), ALFA1(5,1),
     1     ALFA1(6,1), ALFA1(7,1), ALFA1(8,1), ALFA1(9,1), ALFA1(10,1),
     2     ALFA1(11,1),ALFA1(12,1),ALFA1(13,1),ALFA1(14,1),ALFA1(15,1),
     3     ALFA1(16,1),ALFA1(17,1),ALFA1(18,1),ALFA1(19,1),ALFA1(20,1),
     4     ALFA1(21,1),ALFA1(22,1),ALFA1(23,1),ALFA1(24,1),ALFA1(25,1),
     5     ALFA1(26,1)     /-4.44444444444444E-03,-9.22077922077922E-04,
     6-8.84892884892885E-05, 1.65927687832450E-04, 2.46691372741793E-04,
     7 2.65995589346255E-04, 2.61824297061501E-04, 2.48730437344656E-04,
     8 2.32721040083232E-04, 2.16362485712365E-04, 2.00738858762752E-04,
     9 1.86267636637545E-04, 1.73060775917876E-04, 1.61091705929016E-04,
     1 1.50274774160908E-04, 1.40503497391270E-04, 1.31668816545923E-04,
     2 1.23667445598253E-04, 1.16405271474738E-04, 1.09798298372713E-04,
     3 1.03772410422993E-04, 9.82626078369363E-05, 9.32120517249503E-05,
     4 8.85710852478712E-05, 8.42963105715700E-05, 8.03497548407791E-05/
      DATA ALFA1(1,2), ALFA1(2,2), ALFA1(3,2), ALFA1(4,2), ALFA1(5,2),
     1     ALFA1(6,2), ALFA1(7,2), ALFA1(8,2), ALFA1(9,2), ALFA1(10,2),
     2     ALFA1(11,2),ALFA1(12,2),ALFA1(13,2),ALFA1(14,2),ALFA1(15,2),
     3     ALFA1(16,2),ALFA1(17,2),ALFA1(18,2),ALFA1(19,2),ALFA1(20,2),
     4     ALFA1(21,2),ALFA1(22,2),ALFA1(23,2),ALFA1(24,2),ALFA1(25,2),
     5     ALFA1(26,2)     / 6.93735541354589E-04, 2.32241745182922E-04,
     6-1.41986273556691E-05,-1.16444931672049E-04,-1.50803558053049E-04,
     7-1.55121924918096E-04,-1.46809756646466E-04,-1.33815503867491E-04,
     8-1.19744975684254E-04,-1.06184319207974E-04,-9.37699549891194E-05,
     9-8.26923045588193E-05,-7.29374348155221E-05,-6.44042357721016E-05,
     1-5.69611566009369E-05,-5.04731044303562E-05,-4.48134868008883E-05,
     2-3.98688727717599E-05,-3.55400532972042E-05,-3.17414256609022E-05,
     3-2.83996793904175E-05,-2.54522720634871E-05,-2.28459297164725E-05,
     4-2.05352753106481E-05,-1.84816217627666E-05,-1.66519330021394E-05/
      DATA ALFA2(1,1), ALFA2(2,1), ALFA2(3,1), ALFA2(4,1), ALFA2(5,1),
     1     ALFA2(6,1), ALFA2(7,1), ALFA2(8,1), ALFA2(9,1), ALFA2(10,1),
     2     ALFA2(11,1),ALFA2(12,1),ALFA2(13,1),ALFA2(14,1),ALFA2(15,1),
     3     ALFA2(16,1),ALFA2(17,1),ALFA2(18,1),ALFA2(19,1),ALFA2(20,1),
     4     ALFA2(21,1),ALFA2(22,1),ALFA2(23,1),ALFA2(24,1),ALFA2(25,1),
     5     ALFA2(26,1)     /-3.54211971457744E-04,-1.56161263945159E-04,
     6 3.04465503594936E-05, 1.30198655773243E-04, 1.67471106699712E-04,
     7 1.70222587683593E-04, 1.56501427608595E-04, 1.36339170977445E-04,
     8 1.14886692029825E-04, 9.45869093034688E-05, 7.64498419250898E-05,
     9 6.07570334965197E-05, 4.74394299290509E-05, 3.62757512005344E-05,
     1 2.69939714979225E-05, 1.93210938247939E-05, 1.30056674793963E-05,
     2 7.82620866744497E-06, 3.59257485819352E-06, 1.44040049814252E-07,
     3-2.65396769697939E-06,-4.91346867098486E-06,-6.72739296091248E-06,
     4-8.17269379678658E-06,-9.31304715093561E-06,-1.02011418798016E-05/
      DATA ALFA2(1,2), ALFA2(2,2), ALFA2(3,2), ALFA2(4,2), ALFA2(5,2),
     1     ALFA2(6,2), ALFA2(7,2), ALFA2(8,2), ALFA2(9,2), ALFA2(10,2),
     2     ALFA2(11,2),ALFA2(12,2),ALFA2(13,2),ALFA2(14,2),ALFA2(15,2),
     3     ALFA2(16,2),ALFA2(17,2),ALFA2(18,2),ALFA2(19,2),ALFA2(20,2),
     4     ALFA2(21,2),ALFA2(22,2),ALFA2(23,2),ALFA2(24,2),ALFA2(25,2),
     5     ALFA2(26,2)     / 3.78194199201773E-04, 2.02471952761816E-04,
     6-6.37938506318862E-05,-2.38598230603006E-04,-3.10916256027362E-04,
     7-3.13680115247576E-04,-2.78950273791323E-04,-2.28564082619141E-04,
     8-1.75245280340847E-04,-1.25544063060690E-04,-8.22982872820208E-05,
     9-4.62860730588116E-05,-1.72334302366962E-05, 5.60690482304602E-06,
     1 2.31395443148287E-05, 3.62642745856794E-05, 4.58006124490189E-05,
     2 5.24595294959114E-05, 5.68396208545815E-05, 5.94349820393104E-05,
     3 6.06478527578422E-05, 6.08023907788436E-05, 6.01577894539460E-05,
     4 5.89199657344698E-05, 5.72515823777593E-05, 5.52804375585853E-05/
      DATA BETA1(1,1), BETA1(2,1), BETA1(3,1), BETA1(4,1), BETA1(5,1),
     1     BETA1(6,1), BETA1(7,1), BETA1(8,1), BETA1(9,1), BETA1(10,1),
     2     BETA1(11,1),BETA1(12,1),BETA1(13,1),BETA1(14,1),BETA1(15,1),
     3     BETA1(16,1),BETA1(17,1),BETA1(18,1),BETA1(19,1),BETA1(20,1),
     4     BETA1(21,1),BETA1(22,1),BETA1(23,1),BETA1(24,1),BETA1(25,1),
     5     BETA1(26,1)     / 1.79988721413553E-02, 5.59964911064388E-03,
     6 2.88501402231133E-03, 1.80096606761054E-03, 1.24753110589199E-03,
     7 9.22878876572938E-04, 7.14430421727287E-04, 5.71787281789705E-04,
     8 4.69431007606482E-04, 3.93232835462917E-04, 3.34818889318298E-04,
     9 2.88952148495752E-04, 2.52211615549573E-04, 2.22280580798883E-04,
     1 1.97541838033063E-04, 1.76836855019718E-04, 1.59316899661821E-04,
     2 1.44347930197334E-04, 1.31448068119965E-04, 1.20245444949303E-04,
     3 1.10449144504599E-04, 1.01828770740567E-04, 9.41998224204238E-05,
     4 8.74130545753834E-05, 8.13466262162801E-05, 7.59002269646219E-05/
      DATA BETA1(1,2), BETA1(2,2), BETA1(3,2), BETA1(4,2), BETA1(5,2),
     1     BETA1(6,2), BETA1(7,2), BETA1(8,2), BETA1(9,2), BETA1(10,2),
     2     BETA1(11,2),BETA1(12,2),BETA1(13,2),BETA1(14,2),BETA1(15,2),
     3     BETA1(16,2),BETA1(17,2),BETA1(18,2),BETA1(19,2),BETA1(20,2),
     4     BETA1(21,2),BETA1(22,2),BETA1(23,2),BETA1(24,2),BETA1(25,2),
     5     BETA1(26,2)     /-1.49282953213429E-03,-8.78204709546389E-04,
     6-5.02916549572035E-04,-2.94822138512746E-04,-1.75463996970783E-04,
     7-1.04008550460816E-04,-5.96141953046458E-05,-3.12038929076098E-05,
     8-1.26089735980230E-05,-2.42892608575730E-07, 8.05996165414274E-06,
     9 1.36507009262147E-05, 1.73964125472926E-05, 1.98672978842134E-05,
     1 2.14463263790823E-05, 2.23954659232457E-05, 2.28967783814713E-05,
     2 2.30785389811178E-05, 2.30321976080909E-05, 2.28236073720349E-05,
     3 2.25005881105292E-05, 2.20981015361991E-05, 2.16418427448104E-05,
     4 2.11507649256221E-05, 2.06388749782171E-05, 2.01165241997082E-05/
      DATA BETA2(1,1), BETA2(2,1), BETA2(3,1), BETA2(4,1), BETA2(5,1),
     1     BETA2(6,1), BETA2(7,1), BETA2(8,1), BETA2(9,1), BETA2(10,1),
     2     BETA2(11,1),BETA2(12,1),BETA2(13,1),BETA2(14,1),BETA2(15,1),
     3     BETA2(16,1),BETA2(17,1),BETA2(18,1),BETA2(19,1),BETA2(20,1),
     4     BETA2(21,1),BETA2(22,1),BETA2(23,1),BETA2(24,1),BETA2(25,1),
     5     BETA2(26,1)     / 5.52213076721293E-04, 4.47932581552385E-04,
     6 2.79520653992021E-04, 1.52468156198447E-04, 6.93271105657044E-05,
     7 1.76258683069991E-05,-1.35744996343269E-05,-3.17972413350427E-05,
     8-4.18861861696693E-05,-4.69004889379141E-05,-4.87665447413787E-05,
     9-4.87010031186735E-05,-4.74755620890087E-05,-4.55813058138628E-05,
     1-4.33309644511266E-05,-4.09230193157750E-05,-3.84822638603221E-05,
     2-3.60857167535411E-05,-3.37793306123367E-05,-3.15888560772110E-05,
     3-2.95269561750807E-05,-2.75978914828336E-05,-2.58006174666884E-05,
     4-2.41308356761280E-05,-2.25823509518346E-05,-2.11479656768913E-05/
      DATA BETA2(1,2), BETA2(2,2), BETA2(3,2), BETA2(4,2), BETA2(5,2),
     1     BETA2(6,2), BETA2(7,2), BETA2(8,2), BETA2(9,2), BETA2(10,2),
     2     BETA2(11,2),BETA2(12,2),BETA2(13,2),BETA2(14,2),BETA2(15,2),
     3     BETA2(16,2),BETA2(17,2),BETA2(18,2),BETA2(19,2),BETA2(20,2),
     4     BETA2(21,2),BETA2(22,2),BETA2(23,2),BETA2(24,2),BETA2(25,2),
     5     BETA2(26,2)     /-4.74617796559960E-04,-4.77864567147321E-04,
     6-3.20390228067038E-04,-1.61105016119962E-04,-4.25778101285435E-05,
     7 3.44571294294968E-05, 7.97092684075675E-05, 1.03138236708272E-04,
     8 1.12466775262204E-04, 1.13103642108481E-04, 1.08651634848774E-04,
     9 1.01437951597662E-04, 9.29298396593364E-05, 8.40293133016090E-05,
     1 7.52727991349134E-05, 6.69632521975731E-05, 5.92564547323195E-05,
     2 5.22169308826976E-05, 4.58539485165361E-05, 4.01445513891487E-05,
     3 3.50481730031328E-05, 3.05157995034347E-05, 2.64956119950516E-05,
     4 2.29363633690998E-05, 1.97893056664022E-05, 1.70091984636413E-05/
      DATA BETA3(1,1), BETA3(2,1), BETA3(3,1), BETA3(4,1), BETA3(5,1),
     1     BETA3(6,1), BETA3(7,1), BETA3(8,1), BETA3(9,1), BETA3(10,1),
     2     BETA3(11,1),BETA3(12,1),BETA3(13,1),BETA3(14,1),BETA3(15,1),
     3     BETA3(16,1),BETA3(17,1),BETA3(18,1),BETA3(19,1),BETA3(20,1),
     4     BETA3(21,1),BETA3(22,1),BETA3(23,1),BETA3(24,1),BETA3(25,1),
     5     BETA3(26,1)     / 7.36465810572578E-04, 8.72790805146194E-04,
     6 6.22614862573135E-04, 2.85998154194304E-04, 3.84737672879366E-06,
     7-1.87906003636972E-04,-2.97603646594555E-04,-3.45998126832656E-04,
     8-3.53382470916038E-04,-3.35715635775049E-04,-3.04321124789040E-04,
     9-2.66722723047613E-04,-2.27654214122820E-04,-1.89922611854562E-04,
     1-1.55058918599094E-04,-1.23778240761874E-04,-9.62926147717644E-05,
     2-7.25178327714425E-05,-5.22070028895634E-05,-3.50347750511901E-05,
     3-2.06489761035552E-05,-8.70106096849767E-06, 1.13698686675100E-06,
     4 9.16426474122779E-06, 1.56477785428873E-05, 2.08223629482467E-05/
      DATA GAMA(1),   GAMA(2),   GAMA(3),   GAMA(4),   GAMA(5),
     1     GAMA(6),   GAMA(7),   GAMA(8),   GAMA(9),   GAMA(10),
     2     GAMA(11),  GAMA(12),  GAMA(13),  GAMA(14),  GAMA(15),
     3     GAMA(16),  GAMA(17),  GAMA(18),  GAMA(19),  GAMA(20),
     4     GAMA(21),  GAMA(22),  GAMA(23),  GAMA(24),  GAMA(25),
     5     GAMA(26)        / 6.29960524947437E-01, 2.51984209978975E-01,
     6 1.54790300415656E-01, 1.10713062416159E-01, 8.57309395527395E-02,
     7 6.97161316958684E-02, 5.86085671893714E-02, 5.04698873536311E-02,
     8 4.42600580689155E-02, 3.93720661543510E-02, 3.54283195924455E-02,
     9 3.21818857502098E-02, 2.94646240791158E-02, 2.71581677112934E-02,
     1 2.51768272973862E-02, 2.34570755306079E-02, 2.19508390134907E-02,
     2 2.06210828235646E-02, 1.94388240897881E-02, 1.83810633800683E-02,
     3 1.74293213231963E-02, 1.65685837786612E-02, 1.57865285987918E-02,
     4 1.50729501494096E-02, 1.44193250839955E-02, 1.38184805735342E-02/
C***FIRST EXECUTABLE STATEMENT  ASYJY
      TA = R1MACH(3)
      TOL = AMAX1(TA,1.0E-15) 
      TB = R1MACH(5)
      JU = I1MACH(12)
      IF(FLGJY.EQ.1.0E0) GO TO 6
      JR = I1MACH(11)
      ELIM = 2.303E0*TB*(FLOAT(-JU)-FLOAT(JR))
      GO TO 7
    6 CONTINUE
      ELIM = 2.303E0*(TB*FLOAT(-JU)-3.0E0)
    7 CONTINUE
      FN = FNU
      IFLW = 0
      DO 170 JN=1,IN
        XX = X/FN
        WK(1) = 1.0E0 - XX*XX 
        ABW2 = ABS(WK(1))
        WK(2) = SQRT(ABW2)
        WK(7) = FN**CON2
        IF (ABW2.GT.0.27750E0) GO TO 80 
C
C     ASYMPTOTIC EXPANSION
C     CASES NEAR X=FN, ABS(1.-(X/FN)**2).LE.0.2775
C     COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES
C
C     ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES
C
C     KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=MAX(2,SA)
C
        SA = 0.0E0
        IF (ABW2.EQ.0.0E0) GO TO 10
        SA = TOLS/ALOG(ABW2)
   10   SB = SA
        DO 20 I=1,5 
          AKM = AMAX1(SA,2.0E0)
          KMAX(I) = INT(AKM)
          SA = SA + SB
   20   CONTINUE
        KB = KMAX(5)
        KLAST = KB - 1
        SA = GAMA(KB)
        DO 30 K=1,KLAST
          KB = KB - 1
          SA = SA*WK(1) + GAMA(KB)
   30   CONTINUE
        Z = WK(1)*SA
        AZ = ABS(Z) 
        RTZ = SQRT(AZ)
        WK(3) = CON1*AZ*RTZ
        WK(4) = WK(3)*FN
        WK(5) = RTZ*WK(7)
        WK(6) = -WK(5)*WK(5)
        IF(Z.LE.0.0E0) GO TO 35
        IF(WK(4).GT.ELIM) GO TO 75
        WK(6) = -WK(6)
   35   CONTINUE
        PHI = SQRT(SQRT(SA+SA+SA+SA))
C
C     B(ZETA) FOR S=0
C
        KB = KMAX(5)
        KLAST = KB - 1
        SB = BETA(KB,1)
        DO 40 K=1,KLAST
          KB = KB - 1
          SB = SB*WK(1) + BETA(KB,1)
   40   CONTINUE
        KSP1 = 1
        FN2 = FN*FN 
        RFN2 = 1.0E0/FN2
        RDEN = 1.0E0
        ASUM = 1.0E0
        RELB = TOL*ABS(SB)
        BSUM = SB
        DO 60 KS=1,4
          KSP1 = KSP1 + 1
          RDEN = RDEN*RFN2
C
C     A(ZETA) AND B(ZETA) FOR S=1,2,3,4 
C
          KSTEMP = 5 - KS
          KB = KMAX(KSTEMP)
          KLAST = KB - 1
          SA = ALFA(KB,KS)
          SB = BETA(KB,KSP1)
          DO 50 K=1,KLAST
            KB = KB - 1
            SA = SA*WK(1) + ALFA(KB,KS) 
            SB = SB*WK(1) + BETA(KB,KSP1)
   50     CONTINUE
          TA = SA*RDEN
          TB = SB*RDEN
          ASUM = ASUM + TA
          BSUM = BSUM + TB
          IF (ABS(TA).LE.TOL .AND. ABS(TB).LE.RELB) GO TO 70
   60   CONTINUE
   70   CONTINUE
        BSUM = BSUM/(FN*WK(7))
        GO TO 160
C
   75   CONTINUE
        IFLW = 1
        RETURN
C
   80   CONTINUE
        UPOL(1) = 1.0E0
        TAU = 1.0E0/WK(2)
        T2 = 1.0E0/WK(1)
        IF (WK(1).GE.0.0E0) GO TO 90
C
C     CASES FOR (X/FN).GT.SQRT(1.2775)
C
        WK(3) = ABS(WK(2)-ATAN(WK(2)))
        WK(4) = WK(3)*FN
        RCZ = -CON1/WK(4)
        Z32 = 1.5E0*WK(3)
        RTZ = Z32**CON2
        WK(5) = RTZ*WK(7)
        WK(6) = -WK(5)*WK(5)
        GO TO 100
   90   CONTINUE
C
C     CASES FOR (X/FN).LT.SQRT(0.7225)
C
        WK(3) = ABS(ALOG((1.0E0+WK(2))/XX)-WK(2)) 
        WK(4) = WK(3)*FN
        RCZ = CON1/WK(4)
        IF(WK(4).GT.ELIM) GO TO 75
        Z32 = 1.5E0*WK(3)
        RTZ = Z32**CON2
        WK(7) = FN**CON2
        WK(5) = RTZ*WK(7)
        WK(6) = WK(5)*WK(5)
  100   CONTINUE
        PHI = SQRT((RTZ+RTZ)*TAU)
        TB = 1.0E0
        ASUM = 1.0E0
        TFN = TAU/FN
        UPOL(2) = (C(1)*T2+C(2))*TFN
        CRZ32 = CON548*RCZ
        BSUM = UPOL(2) + CRZ32
        RELB = TOL*ABS(BSUM)
        AP = TFN
        KS = 0
        KP1 = 2
        RZDEN = RCZ 
        L = 2
        DO 140 LR=2,8,2
C
C     COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA)
C
          LRP1 = LR + 1
          DO 120 K=LR,LRP1
            KS = KS + 1
            KP1 = KP1 + 1
            L = L + 1
            S1 = C(L)
            DO 110 J=2,KP1
              L = L + 1
              S1 = S1*T2 + C(L)
  110       CONTINUE
            AP = AP*TFN
            UPOL(KP1) = AP*S1 
            CR(KS) = BR(KS)*RZDEN
            RZDEN = RZDEN*RCZ 
            DR(KS) = AR(KS)*RZDEN
  120     CONTINUE
          SUMA = UPOL(LRP1)
          SUMB = UPOL(LR+2) + UPOL(LRP1)*CRZ32
          JU = LRP1 
          DO 130 JR=1,LR
            JU = JU - 1
            SUMA = SUMA + CR(JR)*UPOL(JU)
            SUMB = SUMB + DR(JR)*UPOL(JU)
  130     CONTINUE
          TB = -TB
          IF (WK(1).GT.0.0E0) TB = ABS(TB)
          ASUM = ASUM + SUMA*TB
          BSUM = BSUM + SUMB*TB
          IF (ABS(SUMA).LE.TOL .AND. ABS(SUMB).LE.RELB) GO TO 150
  140   CONTINUE
  150   TB = WK(5)
        IF (WK(1).GT.0.0E0) TB = -TB
        BSUM = BSUM/TB
C
  160   CONTINUE
        CALL FUNJY(WK(6), WK(5), WK(4), FI, DFI)
        Y(JN) = FLGJY*PHI*(FI*ASUM+DFI*BSUM)/WK(7)
        FN = FN - FLGJY
  170 CONTINUE
      RETURN
      END 
      SUBROUTINE JAIRY(X,RX,C,AI,DAI)
C***BEGIN PROLOGUE  JAIRY
C***REFER TO  BESJ,BESY
C
C     1-2-74
C                  JAIRY computes the Airy function AI(X) 
C                   and its derivative DAI(X) for ASYJY
C
C                                   INPUT
C
C         X - Argument, computed by ASYJY, X unrestricted
C        RX - RX=SQRT(ABS(X)), computed by ASYJY
C         C - C=2.*(ABS(X)**1.5)/3., computed by ASYJY
C
C                                  OUTPUT
C
C        AI - Value of function AI(X)
C       DAI - Value of the derivative DAI(X)
C
C                                Written by
C
C                                D. E. Amos
C                               S. L. Daniel 
C                               M. K. Weston 
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  JAIRY
C
      INTEGER I, J, M1, M1D, M2, M2D, M3, M3D, M4, M4D, N1, N1D, N2,
     1 N2D, N3, N3D, N4, N4D
      REAL A, AI, AJN, AJP, AK1, AK2, AK3, B, C, CCV, CON2, CON3,
     1 CON4, CON5, CV, DA, DAI, DAJN, DAJP, DAK1, DAK2, DAK3, DB, EC, 
     2 E1, E2, FPI12, F1, F2, RTRX, RX, SCV, T, TEMP1, TEMP2, TT, X
      DIMENSION AJP(19), AJN(19), A(15), B(15)
      DIMENSION AK1(14), AK2(23), AK3(14)
      DIMENSION DAJP(19), DAJN(19), DA(15), DB(15)
      DIMENSION DAK1(14), DAK2(24), DAK3(14)
      DATA N1,N2,N3,N4/14,23,19,15/
      DATA M1,M2,M3,M4/12,21,17,13/
      DATA FPI12,CON2,CON3,CON4,CON5/
     1 1.30899693899575E+00, 5.03154716196777E+00,
     2 3.80004589867293E-01, 8.33333333333333E-01, 8.66025403784439E-01/
      DATA AK1(1), AK1(2), AK1(3), AK1(4), AK1(5), AK1(6), AK1(7),
     1     AK1(8), AK1(9), AK1(10),AK1(11),AK1(12),AK1(13), 
     2     AK1(14)         / 2.20423090987793E-01,-1.25290242787700E-01,
     3 1.03881163359194E-02, 8.22844152006343E-04,-2.34614345891226E-04,
     4 1.63824280172116E-05, 3.06902589573189E-07,-1.29621999359332E-07,
     5 8.22908158823668E-09, 1.53963968623298E-11,-3.39165465615682E-11,
     6 2.03253257423626E-12,-1.10679546097884E-14,-5.16169497785080E-15/
      DATA AK2(1), AK2(2), AK2(3), AK2(4), AK2(5), AK2(6), AK2(7),
     1     AK2(8), AK2(9), AK2(10),AK2(11),AK2(12),AK2(13),AK2(14),
     2     AK2(15),AK2(16),AK2(17),AK2(18),AK2(19),AK2(20),AK2(21),
     3     AK2(22),AK2(23) / 2.74366150869598E-01, 5.39790969736903E-03,
     4-1.57339220621190E-03, 4.27427528248750E-04,-1.12124917399925E-04,
     5 2.88763171318904E-05,-7.36804225370554E-06, 1.87290209741024E-06,
     6-4.75892793962291E-07, 1.21130416955909E-07,-3.09245374270614E-08,
     7 7.92454705282654E-09,-2.03902447167914E-09, 5.26863056595742E-10,
     8-1.36704767639569E-10, 3.56141039013708E-11,-9.31388296548430E-12,
     9 2.44464450473635E-12,-6.43840261990955E-13, 1.70106030559349E-13,
     1-4.50760104503281E-14, 1.19774799164811E-14,-3.19077040865066E-15/
      DATA AK3(1), AK3(2), AK3(3), AK3(4), AK3(5), AK3(6), AK3(7),
     1     AK3(8), AK3(9), AK3(10),AK3(11),AK3(12),AK3(13), 
     2     AK3(14)         / 2.80271447340791E-01,-1.78127042844379E-03,
     3 4.03422579628999E-05,-1.63249965269003E-06, 9.21181482476768E-08,
     4-6.52294330229155E-09, 5.47138404576546E-10,-5.24408251800260E-11,
     5 5.60477904117209E-12,-6.56375244639313E-13, 8.31285761966247E-14,
     6-1.12705134691063E-14, 1.62267976598129E-15,-2.46480324312426E-16/
      DATA AJP(1), AJP(2), AJP(3), AJP(4), AJP(5), AJP(6), AJP(7),
     1     AJP(8), AJP(9), AJP(10),AJP(11),AJP(12),AJP(13),AJP(14),
     2     AJP(15),AJP(16),AJP(17),AJP(18),
     3     AJP(19)         / 7.78952966437581E-02,-1.84356363456801E-01,
     4 3.01412605216174E-02, 3.05342724277608E-02,-4.95424702513079E-03,
     5-1.72749552563952E-03, 2.43137637839190E-04, 5.04564777517082E-05,
     6-6.16316582695208E-06,-9.03986745510768E-07, 9.70243778355884E-08,
     7 1.09639453305205E-08,-1.04716330588766E-09,-9.60359441344646E-11,
     8 8.25358789454134E-12, 6.36123439018768E-13,-4.96629614116015E-14,
     9-3.29810288929615E-15, 2.35798252031104E-16/
      DATA AJN(1), AJN(2), AJN(3), AJN(4), AJN(5), AJN(6), AJN(7),
     1     AJN(8), AJN(9), AJN(10),AJN(11),AJN(12),AJN(13),AJN(14),
     2     AJN(15),AJN(16),AJN(17),AJN(18),
     3     AJN(19)         / 3.80497887617242E-02,-2.45319541845546E-01,
     4 1.65820623702696E-01, 7.49330045818789E-02,-2.63476288106641E-02,
     5-5.92535597304981E-03, 1.44744409589804E-03, 2.18311831322215E-04,
     6-4.10662077680304E-05,-4.66874994171766E-06, 7.15218807277160E-07,
     7 6.52964770854633E-08,-8.44284027565946E-09,-6.44186158976978E-10,
     8 7.20802286505285E-11, 4.72465431717846E-12,-4.66022632547045E-13,
     9-2.67762710389189E-14, 2.36161316570019E-15/
      DATA A(1),   A(2),   A(3),   A(4),   A(5),   A(6),   A(7),
     1     A(8),   A(9),   A(10),  A(11),  A(12),  A(13),  A(14),
     2     A(15)           / 4.90275424742791E-01, 1.57647277946204E-03,
     3-9.66195963140306E-05, 1.35916080268815E-07, 2.98157342654859E-07,
     4-1.86824767559979E-08,-1.03685737667141E-09, 3.28660818434328E-10,
     5-2.57091410632780E-11,-2.32357655300677E-12, 9.57523279048255E-13,
     6-1.20340828049719E-13,-2.90907716770715E-15, 4.55656454580149E-15,
     7-9.99003874810259E-16/
      DATA B(1),   B(2),   B(3),   B(4),   B(5),   B(6),   B(7),
     1     B(8),   B(9),   B(10),  B(11),  B(12),  B(13),  B(14),
     2     B(15)           / 2.78593552803079E-01,-3.52915691882584E-03,
     3-2.31149677384994E-05, 4.71317842263560E-06,-1.12415907931333E-07,
     4-2.00100301184339E-08, 2.60948075302193E-09,-3.55098136101216E-11,
     5-3.50849978423875E-11, 5.83007187954202E-12,-2.04644828753326E-13,
     6-1.10529179476742E-13, 2.87724778038775E-14,-2.88205111009939E-15,
     7-3.32656311696166E-16/
      DATA N1D,N2D,N3D,N4D/14,24,19,15/ 
      DATA M1D,M2D,M3D,M4D/12,22,17,13/ 
      DATA DAK1(1), DAK1(2), DAK1(3), DAK1(4), DAK1(5), DAK1(6),
     1     DAK1(7), DAK1(8), DAK1(9), DAK1(10),DAK1(11),DAK1(12),
     2    DAK1(13),DAK1(14)/ 2.04567842307887E-01,-6.61322739905664E-02,
     3-8.49845800989287E-03, 3.12183491556289E-03,-2.70016489829432E-04,
     4-6.35636298679387E-06, 3.02397712409509E-06,-2.18311195330088E-07,
     5-5.36194289332826E-10, 1.13098035622310E-09,-7.43023834629073E-11,
     6 4.28804170826891E-13, 2.23810925754539E-13,-1.39140135641182E-14/
      DATA DAK2(1), DAK2(2), DAK2(3), DAK2(4), DAK2(5), DAK2(6),
     1     DAK2(7), DAK2(8), DAK2(9), DAK2(10),DAK2(11),DAK2(12),
     2     DAK2(13),DAK2(14),DAK2(15),DAK2(16),DAK2(17),DAK2(18),
     3     DAK2(19),DAK2(20),DAK2(21),DAK2(22),DAK2(23),
     4     DAK2(24)        / 2.93332343883230E-01,-8.06196784743112E-03,
     5 2.42540172333140E-03,-6.82297548850235E-04, 1.85786427751181E-04,
     6-4.97457447684059E-05, 1.32090681239497E-05,-3.49528240444943E-06,
     7 9.24362451078835E-07,-2.44732671521867E-07, 6.49307837648910E-08,
     8-1.72717621501538E-08, 4.60725763604656E-09,-1.23249055291550E-09,
     9 3.30620409488102E-10,-8.89252099772401E-11, 2.39773319878298E-11,
     1-6.48013921153450E-12, 1.75510132023731E-12,-4.76303829833637E-13,
     2 1.29498241100810E-13,-3.52679622210430E-14, 9.62005151585923E-15,
     3-2.62786914342292E-15/
      DATA DAK3(1), DAK3(2), DAK3(3), DAK3(4), DAK3(5), DAK3(6),
     1     DAK3(7), DAK3(8), DAK3(9), DAK3(10),DAK3(11),DAK3(12),
     2    DAK3(13),DAK3(14)/ 2.84675828811349E-01, 2.53073072619080E-03,
     3-4.83481130337976E-05, 1.84907283946343E-06,-1.01418491178576E-07,
     4 7.05925634457153E-09,-5.85325291400382E-10, 5.56357688831339E-11,
     5-5.90889094779500E-12, 6.88574353784436E-13,-8.68588256452194E-14,
     6 1.17374762617213E-14,-1.68523146510923E-15, 2.55374773097056E-16/
      DATA DAJP(1), DAJP(2), DAJP(3), DAJP(4), DAJP(5), DAJP(6),
     1     DAJP(7), DAJP(8), DAJP(9), DAJP(10),DAJP(11),DAJP(12),
     2     DAJP(13),DAJP(14),DAJP(15),DAJP(16),DAJP(17),DAJP(18),
     3     DAJP(19)        / 6.53219131311457E-02,-1.20262933688823E-01,
     4 9.78010236263823E-03, 1.67948429230505E-02,-1.97146140182132E-03,
     5-8.45560295098867E-04, 9.42889620701976E-05, 2.25827860945475E-05,
     6-2.29067870915987E-06,-3.76343991136919E-07, 3.45663933559565E-08,
     7 4.29611332003007E-09,-3.58673691214989E-10,-3.57245881361895E-11,
     8 2.72696091066336E-12, 2.26120653095771E-13,-1.58763205238303E-14,
     9-1.12604374485125E-15, 7.31327529515367E-17/
      DATA DAJN(1), DAJN(2), DAJN(3), DAJN(4), DAJN(5), DAJN(6),
     1     DAJN(7), DAJN(8), DAJN(9), DAJN(10),DAJN(11),DAJN(12),
     2     DAJN(13),DAJN(14),DAJN(15),DAJN(16),DAJN(17),DAJN(18),
     3     DAJN(19)        / 1.08594539632967E-02, 8.53313194857091E-02,
     4-3.15277068113058E-01,-8.78420725294257E-02, 5.53251906976048E-02,
     5 9.41674060503241E-03,-3.32187026018996E-03,-4.11157343156826E-04,
     6 1.01297326891346E-04, 9.87633682208396E-06,-1.87312969812393E-06,
     7-1.50798500131468E-07, 2.32687669525394E-08, 1.59599917419225E-09,
     8-2.07665922668385E-10,-1.24103350500302E-11, 1.39631765331043E-12,
     9 7.39400971155740E-14,-7.32887475627500E-15/
      DATA DA(1),  DA(2),  DA(3),  DA(4),  DA(5),  DA(6),  DA(7),
     1     DA(8),  DA(9),  DA(10), DA(11), DA(12), DA(13), DA(14),
     2     DA(15)          / 4.91627321104601E-01, 3.11164930427489E-03,
     3 8.23140762854081E-05,-4.61769776172142E-06,-6.13158880534626E-08,
     4 2.87295804656520E-08,-1.81959715372117E-09,-1.44752826642035E-10,
     5 4.53724043420422E-11,-3.99655065847223E-12,-3.24089119830323E-13,
     6 1.62098952568741E-13,-2.40765247974057E-14, 1.69384811284491E-16,
     7 8.17900786477396E-16/
      DATA DB(1),  DB(2),  DB(3),  DB(4),  DB(5),  DB(6),  DB(7),
     1     DB(8),  DB(9),  DB(10), DB(11), DB(12), DB(13), DB(14),
     2     DB(15)          /-2.77571356944231E-01, 4.44212833419920E-03,
     3-8.42328522190089E-05,-2.58040318418710E-06, 3.42389720217621E-07,
     4-6.24286894709776E-09,-2.36377836844577E-09, 3.16991042656673E-10,
     5-4.40995691658191E-12,-5.18674221093575E-12, 9.64874015137022E-13,
     6-4.90190576608710E-14,-1.77253430678112E-14, 5.55950610442662E-15,
     7-7.11793337579530E-16/
C***FIRST EXECUTABLE STATEMENT  JAIRY
      IF (X.LT.0.0E0) GO TO 90
      IF (C.GT.5.0E0) GO TO 60
      IF (X.GT.1.20E0) GO TO 30
      T = (X+X-1.2E0)*CON4
      TT = T + T
      J = N1
      F1 = AK1(J)
      F2 = 0.0E0
      DO 10 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + AK1(J)
        F2 = TEMP1
   10 CONTINUE
      AI = T*F1 - F2 + AK1(1) 
C
      J = N1D
      F1 = DAK1(J)
      F2 = 0.0E0
      DO 20 I=1,M1D 
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DAK1(J)
        F2 = TEMP1
   20 CONTINUE
      DAI = -(T*F1-F2+DAK1(1))
      RETURN
C
   30 CONTINUE
      T = (X+X-CON2)*CON3
      TT = T + T
      J = N2
      F1 = AK2(J)
      F2 = 0.0E0
      DO 40 I=1,M2
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + AK2(J)
        F2 = TEMP1
   40 CONTINUE
      RTRX = SQRT(RX)
      EC = EXP(-C)
      AI = EC*(T*F1-F2+AK2(1))/RTRX
      J = N2D
      F1 = DAK2(J)
      F2 = 0.0E0
      DO 50 I=1,M2D 
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DAK2(J)
        F2 = TEMP1
   50 CONTINUE
      DAI = -EC*(T*F1-F2+DAK2(1))*RTRX
      RETURN
C
   60 CONTINUE
      T = 10.0E0/C - 1.0E0
      TT = T + T
      J = N1
      F1 = AK3(J)
      F2 = 0.0E0
      DO 70 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + AK3(J)
        F2 = TEMP1
   70 CONTINUE
      RTRX = SQRT(RX)
      EC = EXP(-C)
      AI = EC*(T*F1-F2+AK3(1))/RTRX
      J = N1D
      F1 = DAK3(J)
      F2 = 0.0E0
      DO 80 I=1,M1D 
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DAK3(J)
        F2 = TEMP1
   80 CONTINUE
      DAI = -RTRX*EC*(T*F1-F2+DAK3(1))
      RETURN
C
   90 CONTINUE
      IF (C.GT.5.0E0) GO TO 120
      T = 0.4E0*C - 1.0E0
      TT = T + T
      J = N3
      F1 = AJP(J)
      E1 = AJN(J)
      F2 = 0.0E0
      E2 = 0.0E0
      DO 100 I=1,M3 
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + AJP(J)
        E1 = TT*E1 - E2 + AJN(J)
        F2 = TEMP1
        E2 = TEMP2
  100 CONTINUE
      AI = (T*E1-E2+AJN(1)) - X*(T*F1-F2+AJP(1))
      J = N3D
      F1 = DAJP(J)
      E1 = DAJN(J)
      F2 = 0.0E0
      E2 = 0.0E0
      DO 110 I=1,M3D
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + DAJP(J)
        E1 = TT*E1 - E2 + DAJN(J)
        F2 = TEMP1
        E2 = TEMP2
  110 CONTINUE
      DAI = X*X*(T*F1-F2+DAJP(1)) + (T*E1-E2+DAJN(1))
      RETURN
C
  120 CONTINUE
      T = 10.0E0/C - 1.0E0
      TT = T + T
      J = N4
      F1 = A(J)
      E1 = B(J)
      F2 = 0.0E0
      E2 = 0.0E0
      DO 130 I=1,M4 
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + A(J)
        E1 = TT*E1 - E2 + B(J)
        F2 = TEMP1
        E2 = TEMP2
  130 CONTINUE
      TEMP1 = T*F1 - F2 + A(1)
      TEMP2 = T*E1 - E2 + B(1)
      RTRX = SQRT(RX)
      CV = C - FPI12
      CCV = COS(CV) 
      SCV = SIN(CV) 
      AI = (TEMP1*CCV-TEMP2*SCV)/RTRX
      J = N4D
      F1 = DA(J)
      E1 = DB(J)
      F2 = 0.0E0
      E2 = 0.0E0
      DO 140 I=1,M4D
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + DA(J)
        E1 = TT*E1 - E2 + DB(J)
        F2 = TEMP1
        E2 = TEMP2
  140 CONTINUE
      TEMP1 = T*F1 - F2 + DA(1)
      TEMP2 = T*E1 - E2 + DB(1)
      E1 = CCV*CON5 + 0.5E0*SCV
      E2 = SCV*CON5 - 0.5E0*CCV
      DAI = (TEMP1*E1-TEMP2*E2)*RTRX
      RETURN
      END 
      REAL FUNCTION ALNGAM(X)
C***BEGIN PROLOGUE  ALNGAM
C***DATE WRITTEN   770601   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7A
C***KEYWORDS  GAMMA FUNCTION,LOGARITHM,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the log of the absolute value of the Gamma
C            function 
C***DESCRIPTION
C
C ALNGAM(X) computes the logarithm of the absolute value of the
C gamma function at X.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  GAMMA,R1MACH,R9LGMC,XERROR
C***END PROLOGUE  ALNGAM
      REAL SQ2PIL,SQPI2L,PI,XMAX,DXREL,R1MACH,X,Y,GAMMA,SINPIY
      REAL R9LGMC
      DATA SQ2PIL / 0.9189385332 0467274E0/
      DATA SQPI2L / 0.2257913526 4472743E0/
      DATA PI     / 3.1415926535 8979324E0/
      DATA XMAX, DXREL / 0., 0. /
C***FIRST EXECUTABLE STATEMENT  ALNGAM
      IF (XMAX.NE.0.) GO TO 10
      XMAX = R1MACH(2)/ALOG(R1MACH(2))
      DXREL = SQRT (R1MACH(4))
C
 10   Y = ABS(X)
      IF (Y.GT.10.0) GO TO 20 
C
C ALOG (ABS (GAMMA(X))) FOR  ABS(X) .LE. 10.0
C
      ALNGAM = ALOG (ABS (GAMMA(X)))
      RETURN
C
C ALOG (ABS (GAMMA(X))) FOR ABS(X) .GT. 10.0
C
 20   IF (Y.GT.XMAX) CALL XERROR ( 'ALNGAM  ABS(X) SO BIG ALNGAM OVERFLO
     1WS', 38, 2, 2)
C
      IF (X.GT.0.) ALNGAM = SQ2PIL + (X-0.5)*ALOG(X) - X + R9LGMC(Y)
      IF (X.GT.0.) RETURN
C
      SINPIY = ABS (SIN(PI*Y))
      IF (SINPIY.EQ.0.) CALL XERROR ( 'ALNGAM  X IS A NEGATIVE INTEGER',
     1  31, 3, 2)
C
      IF (ABS((X-AINT(X-0.5))/X).LT.DXREL) CALL XERROR ( 'ALNGAM  ANSWER
     1 LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER', 68, 1, 1)
C
      ALNGAM = SQPI2L + (X-0.5)*ALOG(Y) - X - ALOG(SINPIY) - R9LGMC(Y)
      RETURN
C
      END 
      REAL FUNCTION CSEVL(X,CS,N)
C***BEGIN PROLOGUE  CSEVL
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C3A2
C***KEYWORDS  CHEBYSHEV,FNLIB,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Evaluate the N-term Chebyshev series CS at X.
C***DESCRIPTION
C
C Evaluate the N-term Chebyshev series CS at X.  Adapted from
C R. Broucke, Algorithm 446, C.A.C.M., 16, 254 (1973). Also see Fox 
C and Parker, Chebyshev Polynomials in Numerical Analysis, Oxford Press,
C page 56.
C
C       Input Arguments --
C X    value at which the series is to be evaluated.
C CS   array of N terms of a Chebyshev series.  In eval-
C      uating CS, only half the first coefficient is summed.
C N    number of terms in array CS.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  CSEVL
C
       INTEGER N,I,NI
       REAL X,B0,CS,TWOX,B2,B1
       DIMENSION CS(*)
C***FIRST EXECUTABLE STATEMENT  CSEVL
       IF(N.LT.1) CALL XERROR( 'CSEVL   NUMBER OF TERMS LE 0', 28, 2,2)
       IF(N.GT.1000) CALL XERROR ( 'CSEVL   NUMBER OF TERMS GT 1000', 
     1   31,3,2)
       IF (X.LT. -1.0 .OR. X.GT. 1.0) CALL XERROR( 'CSEVL   X OUTSIDE (-
     11,+1)', 25, 1, 1)
C
       B1=0.
       B0=0.
       TWOX=2.*X
       DO 10 I=1,N
       B2=B1
       B1=B0
       NI=N+1-I
       B0=TWOX*B1-B2+CS(NI)
 10    CONTINUE
C
       CSEVL = 0.5 * (B0-B2)
C
       RETURN
      END 
      SUBROUTINE GAMLIM(XMIN,XMAX)
C***BEGIN PROLOGUE  GAMLIM
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7A,R2
C***KEYWORDS  GAMMA FUNCTION,LIMITS,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the minimum and maximum bounds for X in GAMMA(X). 
C***DESCRIPTION
C
C Calculate the minimum and maximum legal bounds for X in GAMMA(X).
C XMIN and XMAX are not the only bounds, but they are the only non- 
C trivial ones to calculate.
C
C             Output Arguments --
C XMIN   minimum legal value of X in GAMMA(X).  Any smaller value of
C        X might result in underflow. 
C XMAX   maximum legal value of X in GAMMA(X).  Any larger value will
C        cause overflow.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH,XERROR
C***END PROLOGUE  GAMLIM
      INTEGER I
      REAL ALNSML,R1MACH,XMIN,XOLD,XLN,ALNBIG,XMAX
C***FIRST EXECUTABLE STATEMENT  GAMLIM
      ALNSML = ALOG(R1MACH(1))
      XMIN = -ALNSML
      DO 10 I=1,10
        XOLD = XMIN 
        XLN = ALOG(XMIN)
        XMIN = XMIN - XMIN*((XMIN+0.5)*XLN - XMIN - 0.2258 + ALNSML)
     1    / (XMIN*XLN + 0.5)
        IF (ABS(XMIN-XOLD).LT.0.005) GO TO 20
 10   CONTINUE
      CALL XERROR ( 'GAMLIM  UNABLE TO FIND XMIN', 27, 1, 2)
C
 20   XMIN = -XMIN + 0.01
C
      ALNBIG = ALOG(R1MACH(2))
      XMAX = ALNBIG 
      DO 30 I=1,10
        XOLD = XMAX 
        XLN = ALOG(XMAX)
        XMAX = XMAX - XMAX*((XMAX-0.5)*XLN - XMAX + 0.9189 - ALNBIG)
     1    / (XMAX*XLN - 0.5)
        IF (ABS(XMAX-XOLD).LT.0.005) GO TO 40
 30   CONTINUE
      CALL XERROR ( 'GAMLIM  UNABLE TO FIND XMAX', 27, 2, 2)
C
 40   XMAX = XMAX - 0.01
      XMIN = AMAX1 (XMIN, -XMAX+1.)
C
      RETURN
      END 
      REAL FUNCTION GAMMA(X)
C***BEGIN PROLOGUE  GAMMA
C***DATE WRITTEN   770601   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7A
C***KEYWORDS  GAMMA FUNCTION,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the Gamma function.
C***DESCRIPTION
C
C GAMMA computes the gamma function at X, where X is not 0, -1, -2, ....
C GAMMA and X are single precision.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL,GAMLIM,INITS,R1MACH,R9LGMC,XERROR
C***END PROLOGUE  GAMMA
      INTEGER NGCS,INITS,N,I
      REAL GCS,PI,SQ2PIL,XMIN,XMAX,DXREL,R1MACH,X,Y,CSEVL
      REAL R9LGMC,SINPIY
      DIMENSION GCS(23)
      DATA GCS   ( 1) / .0085711955 90989331E0/
      DATA GCS   ( 2) / .0044153813 24841007E0/
      DATA GCS   ( 3) / .0568504368 1599363E0/
      DATA GCS   ( 4) /-.0042198353 96418561E0/
      DATA GCS   ( 5) / .0013268081 81212460E0/
      DATA GCS   ( 6) /-.0001893024 529798880E0/
      DATA GCS   ( 7) / .0000360692 532744124E0/
      DATA GCS   ( 8) /-.0000060567 619044608E0/
      DATA GCS   ( 9) / .0000010558 295463022E0/
      DATA GCS   (10) /-.0000001811 967365542E0/
      DATA GCS   (11) / .0000000311 772496471E0/
      DATA GCS   (12) /-.0000000053 542196390E0/
      DATA GCS   (13) / .0000000009 193275519E0/
      DATA GCS   (14) /-.0000000001 577941280E0/
      DATA GCS   (15) / .0000000000 270798062E0/
      DATA GCS   (16) /-.0000000000 046468186E0/
      DATA GCS   (17) / .0000000000 007973350E0/
      DATA GCS   (18) /-.0000000000 001368078E0/
      DATA GCS   (19) / .0000000000 000234731E0/
      DATA GCS   (20) /-.0000000000 000040274E0/
      DATA GCS   (21) / .0000000000 000006910E0/
      DATA GCS   (22) /-.0000000000 000001185E0/
      DATA GCS   (23) / .0000000000 000000203E0/
      DATA PI /3.14159 26535 89793 24E0/
C SQ2PIL IS ALOG (SQRT (2.*PI) )
      DATA SQ2PIL /0.91893 85332 04672 74E0/
      DATA NGCS, XMIN, XMAX, DXREL /0, 3*0.0 /
C
C
C***FIRST EXECUTABLE STATEMENT  GAMMA
      IF (NGCS.NE.0) GO TO 10 
C
C ---------------------------------------------------------------------
C INITIALIZE.  FIND LEGAL BOUNDS FOR X, AND DETERMINE THE NUMBER OF
C TERMS IN THE SERIES REQUIRED TO ATTAIN AN ACCURACY TEN TIMES BETTER 
C THAN MACHINE PRECISION.
C
      NGCS = INITS (GCS, 23, 0.1*R1MACH(3))
C
      CALL GAMLIM (XMIN, XMAX)
      DXREL = SQRT (R1MACH(4))
C
C ---------------------------------------------------------------------
C FINISH INITIALIZATION.  START EVALUATING GAMMA(X).
C
 10   Y = ABS(X)
      IF (Y.GT.10.0) GO TO 50 
C
C COMPUTE GAMMA(X) FOR ABS(X) .LE. 10.0.  REDUCE INTERVAL AND
C FIND GAMMA(1+Y) FOR 0. .LE. Y .LT. 1. FIRST OF ALL.
C
      N = X
      IF (X.LT.0.) N = N - 1
      Y = X - FLOAT(N)
      N = N - 1
      GAMMA = 0.9375 + CSEVL(2.*Y-1., GCS, NGCS)
      IF (N.EQ.0) RETURN
C
      IF (N.GT.0) GO TO 30
C
C COMPUTE GAMMA(X) FOR X .LT. 1.
C
      N = -N
      IF (X.EQ.0.) CALL XERROR ( 'GAMMA   X IS 0', 14, 4, 2)
      IF (X.LT.0. .AND. X+FLOAT(N-2).EQ.0.) CALL XERROR (  'GAMMA   X IS
     1 A NEGATIVE INTEGER', 31, 4, 2)
      IF (X.LT.(-0.5) .AND. ABS((X-AINT(X-0.5))/X).LT.DXREL) CALL
     1  XERROR ( 'GAMMA   ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NE
     2GATIVE INTEGER', 68, 1, 1)
C
      DO 20 I=1,N
        GAMMA = GAMMA / (X+FLOAT(I-1))
 20   CONTINUE
      RETURN
C
C GAMMA(X) FOR X .GE. 2.
C
 30   DO 40 I=1,N
        GAMMA = (Y+FLOAT(I))*GAMMA
 40   CONTINUE
      RETURN
C
C COMPUTE GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
C
 50   IF (X.GT.XMAX) CALL XERROR ( 'GAMMA   X SO BIG GAMMA OVERFLOWS',
     1  32, 3, 2)
C
      GAMMA = 0.
      IF (X.LT.XMIN) CALL XERROR ( 'GAMMA   X SO SMALL GAMMA UNDERFLOWS'
     1  , 35, 2, 1) 
      IF (X.LT.XMIN) RETURN
C
      GAMMA = EXP((Y-0.5)*ALOG(Y) - Y + SQ2PIL + R9LGMC(Y) )
      IF (X.GT.0.) RETURN
C
      IF (ABS((X-AINT(X-0.5))/X).LT.DXREL) CALL XERROR ( 'GAMMA   ANSWER
     1 LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER'  , 61, 1, 1)
C
      SINPIY = SIN (PI*Y)
      IF (SINPIY.EQ.0.) CALL XERROR ( 'GAMMA   X IS A NEGATIVE INTEGER',
     1 31, 4, 2)
C
      GAMMA = -PI / (Y*SINPIY*GAMMA)
C
      RETURN
      END 
      INTEGER FUNCTION INITS(OS,NOS,ETA)
C***BEGIN PROLOGUE  INITS
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C3A2
C***KEYWORDS  INITIALIZE,ORTHOGONAL SERIES,SPECIAL FUNCTION 
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Initializes an orthogonal series so that it defines the
C            number of terms to carry in the series to meet a specified
C            error.
C***DESCRIPTION
C
C Initialize the orthogonal series so that INITS is the number of terms
C needed to insure the error is no larger than ETA.  Ordinarily, ETA
C will be chosen to be one-tenth machine precision.
C
C             Input Arguments --
C OS     array of NOS coefficients in an orthogonal series.
C NOS    number of coefficients in OS.
C ETA    requested accuracy of series.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  INITS
      INTEGER NOS,II,I
      REAL OS,ERR,ETA
      DIMENSION OS(NOS)
C***FIRST EXECUTABLE STATEMENT  INITS
      IF (NOS.LT.1) CALL XERROR ( 'INITS   NUMBER OF COEFFICIENTS LT 1',
     1 35, 2, 2)
C
      ERR = 0.
      DO 10 II=1,NOS
        I = NOS + 1 - II
        ERR = ERR + ABS(OS(I))
        IF (ERR.GT.ETA) GO TO 20
 10   CONTINUE
C
 20   IF (I.EQ.NOS) CALL XERROR ( 'INITS   ETA MAY BE TOO SMALL', 28, 
     1  1, 2)
      INITS = I
C
      RETURN
      END 
      REAL FUNCTION R9LGMC(X)
C***BEGIN PROLOGUE  R9LGMC
C***DATE WRITTEN   770801   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7E
C***KEYWORDS  CORRECTION FACTOR,LOG GAMMA,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the log Gamma correction factor so that
C            ALOG(GAMMA(X)) = ALOG(SQRT(2*PI)) + (X-.5)*ALOG(X) - X
C            + R9LGMC(X)
C***DESCRIPTION
C
C Compute the log gamma correction factor for X .GE. 10.0 so that
C  ALOG (GAMMA(X)) = ALOG(SQRT(2*PI)) + (X-.5)*ALOG(X) - X + R9LGMC(X)
C
C Series for ALGM       on the interval  0.          to  1.00000D-02
C                                        with weighted error   3.40E-16
C                                         log weighted error  15.47
C                               significant figures required  14.39
C                                    decimal places required  15.86
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL,INITS,R1MACH,XERROR
C***END PROLOGUE  R9LGMC
      INTEGER NALGM,INITS
      REAL ALGMCS,XBIG,XMAX,R1MACH,X,CSEVL
      DIMENSION ALGMCS(6)
      DATA ALGMCS( 1) /    .1666389480 45186E0 /
      DATA ALGMCS( 2) /   -.0000138494 817606E0 / 
      DATA ALGMCS( 3) /    .0000000098 108256E0 / 
      DATA ALGMCS( 4) /   -.0000000000 180912E0 / 
      DATA ALGMCS( 5) /    .0000000000 000622E0 / 
      DATA ALGMCS( 6) /   -.0000000000 000003E0 / 
      DATA NALGM, XBIG, XMAX / 0, 2*0.0 /
C***FIRST EXECUTABLE STATEMENT  R9LGMC
      IF (NALGM.NE.0) GO TO 10
      NALGM = INITS (ALGMCS, 6, R1MACH(3))
      XBIG = 1.0/SQRT(R1MACH(3))
      XMAX = EXP (AMIN1(ALOG(R1MACH(2)/12.0), -ALOG(12.0*R1MACH(1))) )
C
 10   IF (X.LT.10.0) CALL XERROR ( 'R9LGMC  X MUST BE GE 10', 23, 1, 2)
      IF (X.GE.XMAX) GO TO 20 
C
      R9LGMC = 1.0/(12.0*X)
      IF (X.LT.XBIG) R9LGMC = CSEVL (2.0*(10./X)**2-1., ALGMCS, NALGM)/X
      RETURN
C
 20   R9LGMC = 0.0
      CALL XERROR ( 'R9LGMC  X SO BIG R9LGMC UNDERFLOWS', 34, 2, 1)
      RETURN
C
      END
