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