SUBROUTINE RYBESL(X,ALPHA,NB,BY,NCALC) C---------------------------------------------------------------------- C C This routine calculates Bessel functions Y SUB(N+ALPHA) (X) C for non-negative argument X, and non-negative order N+ALPHA. C C C Explanation of variables in the calling sequence C C X - Working precision non-negative real argument for which C Y's are to be calculated. C ALPHA - Working precision fractional part of order for which C Y's are to be calculated. 0 .LE. ALPHA .LT. 1.0. C NB - Integer number of functions to be calculated, NB .GT. 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C BY - Working precision output vector of length NB. If the C routine terminates normally (NCALC=NB), the vector BY C contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X), C If (0 .LT. NCALC .LT. NB), BY(I) contains correct function C values for I .LE. NCALC, and contains the ratios C Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array. C NCALC - Integer output variable indicating possible errors. C Before using the vector BY, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See error returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C beta = Radix for the floating-point system C p = Number of significant base-beta digits in the C significand of a floating-point number C minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C EPS = beta ** (-p) C DEL = Machine number below which sin(x)/x = 1; approximately C SQRT(EPS). C XMIN = Smallest acceptable argument for RBESY; approximately C max(2*beta**minexp,2/XINF), rounded up C XINF = Largest positive machine number; approximately C beta**maxexp C THRESH = Lower bound for use of the asymptotic form; approximately C AINT(-LOG10(EPS/2.0))+1.0 C XLARGE = Upper bound on X; approximately 1/DEL, because the sine C and cosine functions have lost about half of their C precision at that point. C C C Approximate values for some important machines are: C C beta p minexp maxexp EPS C C CRAY-1 (S.P.) 2 48 -8193 8191 3.55E-15 C Cyber 180/185 C under NOS (S.P.) 2 48 -975 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 24 -126 128 5.96E-8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 53 -1022 1024 1.11D-16 C IBM 3033 (D.P.) 16 14 -65 63 1.39D-17 C VAX (S.P.) 2 24 -128 127 5.96E-8 C VAX D-Format (D.P.) 2 56 -128 127 1.39D-17 C VAX G-Format (D.P.) 2 53 -1024 1023 1.11D-16 C C C DEL XMIN XINF THRESH XLARGE C C CRAY-1 (S.P.) 5.0E-8 3.67E-2466 5.45E+2465 15.0E0 2.0E7 C Cyber 180/855 C under NOS (S.P.) 5.0E-8 6.28E-294 1.26E+322 15.0E0 2.0E7 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.0E-4 2.36E-38 3.40E+38 8.0E0 1.0E4 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.0D-8 4.46D-308 1.79D+308 16.0D0 1.0D8 C IBM 3033 (D.P.) 1.0D-8 2.77D-76 7.23D+75 17.0D0 1.0D8 C VAX (S.P.) 1.0E-4 1.18E-38 1.70E+38 8.0E0 1.0E4 C VAX D-Format (D.P.) 1.0D-9 1.18D-38 1.70D+38 17.0D0 1.0D9 C VAX G-Format (D.P.) 1.0D-8 2.23D-308 8.98D+307 16.0D0 1.0D8 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all Y's are C calculated to the desired accuracy. C C NCALC .LT. -1: An argument is out of range. For example, C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. C XMAX. In this case, BY(1) = 0.0, the remainder of the C BY-vector is not calculated, and NCALC is set to C MIN0(NB,0)-2 so that NCALC .NE. NB. C NCALC = -1: Y(ALPHA,X) .GE. XINF. The requested function C values are set to 0.0. C 1 .LT. NCALC .LT. NB: Not all requested function values could C be calculated accurately. BY(I) contains correct function C values for I .LE. NCALC, and and the remaining NB-NCALC C array elements contain 0.0. C C C Intrinsic functions required are: C C DBLE, EXP, INT, MAX, MIN, REAL, SQRT C C C Acknowledgement C C This program draws heavily on Temme's Algol program for Y(a,x) C and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's C scheme is used for x < THRESH, and Campbell's scheme is used C in the asymptotic region. Segments of code from both sources C have been translated into Fortran 77, merged, and heavily modified. C Modifications include parameterization of machine dependencies, C use of a new approximation for ln(gamma(x)), and built-in C protection against over/underflow. C C References: "Bessel functions J_nu(x) and Y_nu(x) of real C order and real argument," Campbell, J. B., C Comp. Phy. Comm. 18, 1979, pp. 133-142. C C "On the numerical evaluation of the ordinary C Bessel function of the second kind," Temme, C N. M., J. Comput. Phys. 21, 1976, pp. 343-350. C C Latest modification: March 19, 1990 C C Modified by: W. J. Cody C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,K,NA,NB,NCALC CS REAL CD DOUBLE PRECISION 1 ALFA,ALPHA,AYE,B,BY,C,CH,COSMU,D,DEL,DEN,DDIV,DIV,DMU,D1,D2, 2 E,EIGHT,EN,ENU,EN1,EPS,EVEN,EX,F,FIVPI,G,GAMMA,H,HALF,ODD, 3 ONBPI,ONE,ONE5,P,PA,PA1,PI,PIBY2,PIM5,Q,QA,QA1,Q0,R,S,SINMU, 4 SQ2BPI,TEN9,TERM,THREE,THRESH,TWO,TWOBYX,X,XINF,XLARGE,XMIN, 5 XNA,X2,YA,YA1,ZERO DIMENSION BY(NB),CH(21) C---------------------------------------------------------------------- C Mathematical constants C FIVPI = 5*PI C PIM5 = 5*PI - 15 C ONBPI = 1/PI C PIBY2 = PI/2 C SQ2BPI = SQUARE ROOT OF 2/PI C---------------------------------------------------------------------- CS DATA ZERO,HALF,ONE,TWO,THREE/0.0E0,0.5E0,1.0E0,2.0E0,3.0E0/ CS DATA EIGHT,ONE5,TEN9/8.0E0,15.0E0,1.9E1/ CS DATA FIVPI,PIBY2/1.5707963267948966192E1,1.5707963267948966192E0/ CS DATA PI,SQ2BPI/3.1415926535897932385E0,7.9788456080286535588E-1/ CS DATA PIM5,ONBPI/7.0796326794896619231E-1,3.1830988618379067154E-1/ CD DATA ZERO,HALF,ONE,TWO,THREE/0.0D0,0.5D0,1.0D0,2.0D0,3.0D0/ CD DATA EIGHT,ONE5,TEN9/8.0D0,15.0D0,1.9D1/ CD DATA FIVPI,PIBY2/1.5707963267948966192D1,1.5707963267948966192D0/ CD DATA PI,SQ2BPI/3.1415926535897932385D0,7.9788456080286535588D-1/ CD DATA PIM5,ONBPI/7.0796326794896619231D-1,3.1830988618379067154D-1/ C---------------------------------------------------------------------- C Machine-dependent constants C---------------------------------------------------------------------- CS DATA DEL,XMIN,XINF,EPS/1.0E-4,2.36E-38,3.40E38,5.96E-8/ CS DATA THRESH,XLARGE/8.0E0,1.0E4/ CD DATA DEL,XMIN,XINF,EPS/1.0D-8,4.46D-308,1.79D308,1.11D-16/ CD DATA THRESH,XLARGE/16.0D0,1.0D8/ C---------------------------------------------------------------------- C Coefficients for Chebyshev polynomial expansion of C 1/gamma(1-x), abs(x) .le. .5 C---------------------------------------------------------------------- CS DATA CH/-0.67735241822398840964E-23,-0.61455180116049879894E-22, CS 1 0.29017595056104745456E-20, 0.13639417919073099464E-18, CS 2 0.23826220476859635824E-17,-0.90642907957550702534E-17, CS 3 -0.14943667065169001769E-14,-0.33919078305362211264E-13, CS 4 -0.17023776642512729175E-12, 0.91609750938768647911E-11, CS 5 0.24230957900482704055E-09, 0.17451364971382984243E-08, CS 6 -0.33126119768180852711E-07,-0.86592079961391259661E-06, CS 7 -0.49717367041957398581E-05, 0.76309597585908126618E-04, CS 8 0.12719271366545622927E-02, 0.17063050710955562222E-02, CS 9 -0.76852840844786673690E-01,-0.28387654227602353814E+00, CS A 0.92187029365045265648E+00/ CD DATA CH/-0.67735241822398840964D-23,-0.61455180116049879894D-22, CD 1 0.29017595056104745456D-20, 0.13639417919073099464D-18, CD 2 0.23826220476859635824D-17,-0.90642907957550702534D-17, CD 3 -0.14943667065169001769D-14,-0.33919078305362211264D-13, CD 4 -0.17023776642512729175D-12, 0.91609750938768647911D-11, CD 5 0.24230957900482704055D-09, 0.17451364971382984243D-08, CD 6 -0.33126119768180852711D-07,-0.86592079961391259661D-06, CD 7 -0.49717367041957398581D-05, 0.76309597585908126618D-04, CD 8 0.12719271366545622927D-02, 0.17063050710955562222D-02, CD 9 -0.76852840844786673690D-01,-0.28387654227602353814D+00, CD A 0.92187029365045265648D+00/ C---------------------------------------------------------------------- EX = X ENU = ALPHA IF ((NB .GT. 0) .AND. (X .GE. XMIN) .AND. (EX .LT. XLARGE) 1 .AND. (ENU .GE. ZERO) .AND. (ENU .LT. ONE)) THEN XNA = AINT(ENU+HALF) NA = INT(XNA) IF (NA .EQ. 1) ENU = ENU - XNA IF (ENU .EQ. -HALF) THEN P = SQ2BPI/SQRT(EX) YA = P * SIN(EX) YA1 = -P * COS(EX) ELSE IF (EX .LT. THREE) THEN C---------------------------------------------------------------------- C Use Temme's scheme for small X C---------------------------------------------------------------------- B = EX * HALF D = -LOG(B) F = ENU * D E = B**(-ENU) IF (ABS(ENU) .LT. DEL) THEN C = ONBPI ELSE C = ENU / SIN(ENU*PI) END IF C---------------------------------------------------------------------- C Computation of sinh(f)/f C---------------------------------------------------------------------- IF (ABS(F) .LT. ONE) THEN X2 = F*F EN = TEN9 S = ONE DO 80 I = 1, 9 S = S*X2/EN/(EN-ONE)+ONE EN = EN - TWO 80 CONTINUE ELSE S = (E - ONE/E) * HALF / F END IF C---------------------------------------------------------------------- C Computation of 1/gamma(1-a) using Chebyshev polynomials C---------------------------------------------------------------------- X2 = ENU*ENU*EIGHT AYE = CH(1) EVEN = ZERO ALFA = CH(2) ODD = ZERO DO 40 I = 3, 19, 2 EVEN = -(AYE+AYE+EVEN) AYE = -EVEN*X2 - AYE + CH(I) ODD = -(ALFA+ALFA+ODD) ALFA = -ODD*X2 - ALFA + CH(I+1) 40 CONTINUE EVEN = (EVEN*HALF+AYE)*X2 - AYE + CH(21) ODD = (ODD+ALFA)*TWO GAMMA = ODD*ENU + EVEN C---------------------------------------------------------------------- C End of computation of 1/gamma(1-a) C---------------------------------------------------------------------- G = E * GAMMA E = (E + ONE/E) * HALF F = TWO*C*(ODD*E+EVEN*S*D) E = ENU*ENU P = G*C Q = ONBPI / G C = ENU*PIBY2 IF (ABS(C) .LT. DEL) THEN R = ONE ELSE R = SIN(C)/C END IF R = PI*C*R*R C = ONE D = - B*B H = ZERO YA = F + R*Q YA1 = P EN = ZERO 100 EN = EN + ONE IF (ABS(G/(ONE+ABS(YA))) 1 + ABS(H/(ONE+ABS(YA1))) .GT. EPS) THEN F = (F*EN+P+Q)/(EN*EN-E) C = C * D/EN P = P/(EN-ENU) Q = Q/(EN+ENU) G = C*(F+R*Q) H = C*P - EN*G YA = YA + G YA1 = YA1+H GO TO 100 END IF YA = -YA YA1 = -YA1/B ELSE IF (EX .LT. THRESH) THEN C---------------------------------------------------------------------- C Use Temme's scheme for moderate X C---------------------------------------------------------------------- C = (HALF-ENU)*(HALF+ENU) B = EX + EX E = (EX*ONBPI*COS(ENU*PI)/EPS) E = E*E P = ONE Q = -EX R = ONE + EX*EX S = R EN = TWO 200 IF (R*EN*EN .LT. E) THEN EN1 = EN+ONE D = (EN-ONE+C/EN)/S P = (EN+EN-P*D)/EN1 Q = (-B+Q*D)/EN1 S = P*P + Q*Q R = R*S EN = EN1 GO TO 200 END IF F = P/S P = F G = -Q/S Q = G 220 EN = EN - ONE IF (EN .GT. ZERO) THEN R = EN1*(TWO-P)-TWO S = B + EN1*Q D = (EN-ONE+C/EN)/(R*R+S*S) P = D*R Q = D*S E = F + ONE F = P*E - G*Q G = Q*E + P*G EN1 = EN GO TO 220 END IF F = ONE + F D = F*F + G*G PA = F/D QA = -G/D D = ENU + HALF -P Q = Q + EX PA1 = (PA*Q-QA*D)/EX QA1 = (QA*Q+PA*D)/EX B = EX - PIBY2*(ENU+HALF) C = COS(B) S = SIN(B) D = SQ2BPI/SQRT(EX) YA = D*(PA*S+QA*C) YA1 = D*(QA1*S-PA1*C) ELSE C---------------------------------------------------------------------- C Use Campbell's asymptotic scheme. C---------------------------------------------------------------------- NA = 0 D1 = AINT(EX/FIVPI) I = INT(D1) DMU = ((EX-ONE5*D1)-D1*PIM5)-(ALPHA+HALF)*PIBY2 IF (I-2*(I/2) .EQ. 0) THEN COSMU = COS(DMU) SINMU = SIN(DMU) ELSE COSMU = -COS(DMU) SINMU = -SIN(DMU) END IF DDIV = EIGHT * EX DMU = ALPHA DEN = SQRT(EX) DO 350 K = 1, 2 P = COSMU COSMU = SINMU SINMU = -P D1 = (TWO*DMU-ONE)*(TWO*DMU+ONE) D2 = ZERO DIV = DDIV P = ZERO Q = ZERO Q0 = D1/DIV TERM = Q0 DO 310 I = 2, 20 D2 = D2 + EIGHT D1 = D1 - D2 DIV = DIV + DDIV TERM = -TERM*D1/DIV P = P + TERM D2 = D2 + EIGHT D1 = D1 - D2 DIV = DIV + DDIV TERM = TERM*D1/DIV Q = Q + TERM IF (ABS(TERM) .LE. EPS) GO TO 320 310 CONTINUE 320 P = P + ONE Q = Q + Q0 IF (K .EQ. 1) THEN YA = SQ2BPI * (P*COSMU-Q*SINMU) / DEN ELSE YA1 = SQ2BPI * (P*COSMU-Q*SINMU) / DEN END IF DMU = DMU + ONE 350 CONTINUE END IF IF (NA .EQ. 1) THEN H = TWO*(ENU+ONE)/EX IF (H .GT. ONE) THEN IF (ABS(YA1) .GT. XINF/H) THEN H = ZERO YA = ZERO END IF END IF H = H*YA1 - YA YA = YA1 YA1 = H END IF C---------------------------------------------------------------------- C Now have first one or two Y's C---------------------------------------------------------------------- BY(1) = YA BY(2) = YA1 IF (YA1 .EQ. ZERO) THEN NCALC = 1 ELSE AYE = ONE + ALPHA TWOBYX = TWO/EX NCALC = 2 DO 400 I = 3, NB IF (TWOBYX .LT. ONE) THEN IF (ABS(BY(I-1))*TWOBYX .GE. XINF/AYE) 1 GO TO 450 ELSE IF (ABS(BY(I-1)) .GE. XINF/AYE/TWOBYX ) 1 GO TO 450 END IF BY(I) = TWOBYX*AYE*BY(I-1) - BY(I-2) AYE = AYE + ONE NCALC = NCALC + 1 400 CONTINUE END IF 450 DO 460 I = NCALC+1, NB BY(I) = ZERO 460 CONTINUE ELSE BY(1) = ZERO NCALC = MIN(NB,0) - 1 END IF 900 RETURN C---------- Last line of RYBESL ---------- END