module rjk implicit none contains SUBROUTINE RKBESL(X,ALPHA,NB,IZE,BK,NCALC) C------------------------------------------------------------------- C C This FORTRAN 77 routine calculates modified Bessel functions C of the second kind, K SUB(N+ALPHA) (X), for non-negative C argument X, and non-negative order N+ALPHA, with or without C exponential scaling. C C Explanation of variables in the calling sequence C C Description of output values .. C C X - Working precision non-negative real argument for which C K's or exponentially scaled K's (K*EXP(X)) C are to be calculated. If K's are to be calculated, C X must not be greater than XMAX (see below). C ALPHA - Working precision fractional part of order for which C K's or exponentially scaled K's (K*EXP(X)) are C 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 IZE - Integer type. IZE = 1 if unscaled K's are to be calculated, C and 2 if exponentially scaled K's are to be calculated. C BK - Working precision output vector of length NB. If the C routine terminates normally (NCALC=NB), the vector BK C contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X), C or the corresponding exponentially scaled functions. C If (0 .LT. NCALC .LT. NB), BK(I) contains correct function C values for I .LE. NCALC, and contains the ratios C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. C NCALC - Integer output variable indicating possible errors. C Before using the vector BK, 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 minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C EPS = The smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C XMAX = Upper limit on the magnitude of X when IZE=1; Solution C to equation: C W(X) * (1-1/8X+9/128X**2) = beta**minexp C where W(X) = EXP(-X)*SQRT(PI/2X) C SQXMIN = Square root of beta**minexp C XINF = Largest positive machine number; approximately C beta**maxexp C XMIN = Smallest positive machine number; approximately C beta**minexp C C C Approximate values for some important machines are: C C beta minexp maxexp EPS C C CRAY-1 (S.P.) 2 -8193 8191 7.11E-15 C Cyber 180/185 C under NOS (S.P.) 2 -975 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 -126 128 1.19E-7 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 -1022 1024 2.22D-16 C IBM 3033 (D.P.) 16 -65 63 2.22D-16 C VAX (S.P.) 2 -128 127 5.96E-8 C VAX D-Format (D.P.) 2 -128 127 1.39D-17 C VAX G-Format (D.P.) 2 -1024 1023 1.11D-16 C C C SQXMIN XINF XMIN XMAX C C CRAY-1 (S.P.) 6.77E-1234 5.45E+2465 4.59E-2467 5674.858 C Cyber 180/855 C under NOS (S.P.) 1.77E-147 1.26E+322 3.14E-294 672.788 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.08E-19 3.40E+38 1.18E-38 85.337 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.49D-154 1.79D+308 2.23D-308 705.342 C IBM 3033 (D.P.) 7.35D-40 7.23D+75 5.40D-79 177.852 C VAX (S.P.) 5.42E-20 1.70E+38 2.94E-39 86.715 C VAX D-Format (D.P.) 5.42D-20 1.70D+38 2.94D-39 86.715 C VAX G-Format (D.P.) 7.46D-155 8.98D+307 5.57D-309 706.728 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all K'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, the B-vector is not calculated, C and NCALC is set to MIN0(NB,0)-2 so that NCALC .NE. NB. C NCALC = -1: Either K(ALPHA,X) .GE. XINF or C K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) .GE. XINF. In this case, C the B-vector is not calculated. Note that again C NCALC .NE. NB. C C 0 .LT. NCALC .LT. NB: Not all requested function values could C be calculated accurately. BK(I) contains correct function C values for I .LE. NCALC, and contains the ratios C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. C C C Intrinsic functions required are: C C ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT C C C Acknowledgement C C This program is based on a program written by J. B. Campbell C (2) that computes values of the Bessel functions K of real C argument and real order. Modifications include the addition C of non-scaled functions, parameterization of machine C dependencies, and the use of more accurate approximations C for SINH and SIN. C C References: "On Temme's Algorithm for the Modified Bessel C Functions of the Third Kind," Campbell, J. B., C TOMS 6(4), Dec. 1980, pp. 581-586. C C "A FORTRAN IV Subroutine for the Modified Bessel C Functions of the Third Kind of Real Order and Real C Argument," Campbell, J. B., Report NRC/ERB-925, C National Research Council, Canada. C C Latest modification: May 30, 1989 C C Modified by: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C------------------------------------------------------------------- INTEGER I,IEND,ITEMP,IZE,J,K,M,MPLUS1,NB,NCALC CS REAL DOUBLE PRECISION 1 A,ALPHA,BLPHA,BK,BK1,BK2,C,D,DM,D1,D2,D3,ENU,EPS,ESTF,ESTM, 2 EX,FOUR,F0,F1,F2,HALF,ONE,P,P0,Q,Q0,R,RATIO,S,SQXMIN,T,TINYX, 3 TWO,TWONU,TWOX,T1,T2,WMINF,X,XINF,XMAX,XMIN,X2BY4,ZERO DIMENSION BK(1),P(8),Q(7),R(5),S(4),T(6),ESTM(6),ESTF(7) C--------------------------------------------------------------------- C Mathematical constants C A = LOG(2.D0) - Euler's constant C D = SQRT(2.D0/PI) C--------------------------------------------------------------------- CS DATA HALF,ONE,TWO,ZERO/0.5E0,1.0E0,2.0E0,0.0E0/ CS DATA FOUR,TINYX/4.0E0,1.0E-10/ CS DATA A/ 0.11593151565841244881E0/,D/0.797884560802865364E0/ CD DATA HALF,ONE,TWO,ZERO/0.5D0,1.0D0,2.0D0,0.0D0/ CD DATA FOUR,TINYX/4.0D0,1.0D-10/ CD DATA A/ 0.11593151565841244881D0/,D/0.797884560802865364D0/ C--------------------------------------------------------------------- C Machine dependent parameters C--------------------------------------------------------------------- CS DATA EPS/1.19E-7/,SQXMIN/1.08E-19/,XINF/3.40E+38/ CS DATA XMIN/1.18E-38/,XMAX/85.337E0/ CD DATA EPS/2.22D-16/,SQXMIN/1.49D-154/,XINF/1.79D+308/ CD DATA XMIN/2.23D-308/,XMAX/705.342D0/ C--------------------------------------------------------------------- C P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA C + Euler's constant C Coefficients converted from hex to decimal and modified C by W. J. Cody, 2/26/82 C R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) C T - Approximation for SINH(Y)/Y C--------------------------------------------------------------------- CS DATA P/ 0.805629875690432845E00, 0.204045500205365151E02, CS 1 0.157705605106676174E03, 0.536671116469207504E03, CS 2 0.900382759291288778E03, 0.730923886650660393E03, CS 3 0.229299301509425145E03, 0.822467033424113231E00/ CS DATA Q/ 0.294601986247850434E02, 0.277577868510221208E03, CS 1 0.120670325591027438E04, 0.276291444159791519E04, CS 2 0.344374050506564618E04, 0.221063190113378647E04, CS 3 0.572267338359892221E03/ CS DATA R/-0.48672575865218401848E+0, 0.13079485869097804016E+2, CS 1 -0.10196490580880537526E+3, 0.34765409106507813131E+3, CS 2 0.34958981245219347820E-3/ CS DATA S/-0.25579105509976461286E+2, 0.21257260432226544008E+3, CS 1 -0.61069018684944109624E+3, 0.42269668805777760407E+3/ CS DATA T/ 0.16125990452916363814E-9, 0.25051878502858255354E-7, CS 1 0.27557319615147964774E-5, 0.19841269840928373686E-3, CS 2 0.83333333333334751799E-2, 0.16666666666666666446E+0/ CS DATA ESTM/5.20583E1, 5.7607E0, 2.7782E0, 1.44303E1, 1.853004E2, CS 1 9.3715E0/ CS DATA ESTF/4.18341E1, 7.1075E0, 6.4306E0, 4.25110E1, 1.35633E0, CS 1 8.45096E1, 2.0E1/ CD DATA P/ 0.805629875690432845D00, 0.204045500205365151D02, CD 1 0.157705605106676174D03, 0.536671116469207504D03, CD 2 0.900382759291288778D03, 0.730923886650660393D03, CD 3 0.229299301509425145D03, 0.822467033424113231D00/ CD DATA Q/ 0.294601986247850434D02, 0.277577868510221208D03, CD 1 0.120670325591027438D04, 0.276291444159791519D04, CD 2 0.344374050506564618D04, 0.221063190113378647D04, CD 3 0.572267338359892221D03/ CD DATA R/-0.48672575865218401848D+0, 0.13079485869097804016D+2, CD 1 -0.10196490580880537526D+3, 0.34765409106507813131D+3, CD 2 0.34958981245219347820D-3/ CD DATA S/-0.25579105509976461286D+2, 0.21257260432226544008D+3, CD 1 -0.61069018684944109624D+3, 0.42269668805777760407D+3/ CD DATA T/ 0.16125990452916363814D-9, 0.25051878502858255354D-7, CD 1 0.27557319615147964774D-5, 0.19841269840928373686D-3, CD 2 0.83333333333334751799D-2, 0.16666666666666666446D+0/ CD DATA ESTM/5.20583D1, 5.7607D0, 2.7782D0, 1.44303D1, 1.853004D2, CD 1 9.3715D0/ CD DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0, CD 1 8.45096D1, 2.0D1/ C--------------------------------------------------------------------- EX = X ENU = ALPHA NCALC = MIN(NB,0)-2 IF ((NB .GT. 0) .AND. ((ENU .GE. ZERO) .AND. (ENU .LT. ONE)) 1 .AND. ((IZE .GE. 1) .AND. (IZE .LE. 2)) .AND. 2 ((IZE .NE. 1) .OR. (EX .LE. XMAX)) .AND. 3 (EX .GT. ZERO)) THEN K = 0 IF (ENU .LT. SQXMIN) ENU = ZERO IF (ENU .GT. HALF) THEN K = 1 ENU = ENU - ONE END IF TWONU = ENU+ENU IEND = NB+K-1 C = ENU*ENU D3 = -C IF (EX .LE. ONE) THEN C--------------------------------------------------------------------- C Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA C Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA C--------------------------------------------------------------------- D1 = ZERO D2 = P(1) T1 = ONE T2 = Q(1) DO 10 I = 2,7,2 D1 = C*D1+P(I) D2 = C*D2+P(I+1) T1 = C*T1+Q(I) T2 = C*T2+Q(I+1) 10 CONTINUE D1 = ENU*D1 T1 = ENU*T1 F1 = LOG(EX) F0 = A+ENU*(P(8)-ENU*(D1+D2)/(T1+T2))-F1 Q0 = EXP(-ENU*(A-ENU*(P(8)+ENU*(D1-D2)/(T1-T2))-F1)) F1 = ENU*F0 P0 = EXP(F1) C--------------------------------------------------------------------- C Calculation of F0 = C--------------------------------------------------------------------- D1 = R(5) T1 = ONE DO 20 I = 1,4 D1 = C*D1+R(I) T1 = C*T1+S(I) 20 CONTINUE IF (ABS(F1) .LE. HALF) THEN F1 = F1*F1 D2 = ZERO DO 30 I = 1,6 D2 = F1*D2+T(I) 30 CONTINUE D2 = F0+F0*F1*D2 ELSE D2 = SINH(F1)/ENU END IF F0 = D2-ENU*D1/(T1*P0) IF (EX .LE. TINYX) THEN C-------------------------------------------------------------------- C X.LE.1.0E-10 C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) C-------------------------------------------------------------------- BK(1) = F0+EX*F0 IF (IZE .EQ. 1) BK(1) = BK(1)-EX*BK(1) RATIO = P0/F0 C = EX*XINF IF (K .NE. 0) THEN C-------------------------------------------------------------------- C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X), C ALPHA .GE. 1/2 C-------------------------------------------------------------------- NCALC = -1 IF (BK(1) .GE. C/RATIO) GO TO 500 BK(1) = RATIO*BK(1)/EX TWONU = TWONU+TWO RATIO = TWONU END IF NCALC = 1 IF (NB .EQ. 1) GO TO 500 C-------------------------------------------------------------------- C Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X), L = 1, 2, ... , NB-1 C-------------------------------------------------------------------- NCALC = -1 DO 80 I = 2,NB IF (RATIO .GE. C) GO TO 500 BK(I) = RATIO/EX TWONU = TWONU+TWO RATIO = TWONU 80 CONTINUE NCALC = 1 GO TO 420 ELSE C-------------------------------------------------------------------- C 1.0E-10 .LT. X .LE. 1.0 C-------------------------------------------------------------------- C = ONE X2BY4 = EX*EX/FOUR P0 = HALF*P0 Q0 = HALF*Q0 D1 = -ONE D2 = ZERO BK1 = ZERO BK2 = ZERO F1 = F0 F2 = P0 100 D1 = D1+TWO D2 = D2+ONE D3 = D1+D3 C = X2BY4*C/D2 F0 = (D2*F0+P0+Q0)/D3 P0 = P0/(D2-ENU) Q0 = Q0/(D2+ENU) T1 = C*F0 T2 = C*(P0-D2*F0) BK1 = BK1+T1 BK2 = BK2+T2 IF ((ABS(T1/(F1+BK1)) .GT. EPS) .OR. 1 (ABS(T2/(F2+BK2)) .GT. EPS)) GO TO 100 BK1 = F1+BK1 BK2 = TWO*(F2+BK2)/EX IF (IZE .EQ. 2) THEN D1 = EXP(EX) BK1 = BK1*D1 BK2 = BK2*D1 END IF WMINF = ESTF(1)*EX+ESTF(2) END IF ELSE IF (EPS*EX .GT. ONE) THEN C-------------------------------------------------------------------- C X .GT. ONE/EPS C-------------------------------------------------------------------- NCALC = NB BK1 = ONE / (D*SQRT(EX)) DO 110 I = 1, NB BK(I) = BK1 110 CONTINUE GO TO 500 ELSE C-------------------------------------------------------------------- C X .GT. 1.0 C-------------------------------------------------------------------- TWOX = EX+EX BLPHA = ZERO RATIO = ZERO IF (EX .LE. FOUR) THEN C-------------------------------------------------------------------- C Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 .LE. X .LE. 4.0 C-------------------------------------------------------------------- D2 = AINT(ESTM(1)/EX+ESTM(2)) M = INT(D2) D1 = D2+D2 D2 = D2-HALF D2 = D2*D2 DO 120 I = 2,M D1 = D1-TWO D2 = D2-D1 RATIO = (D3+D2)/(TWOX+D1-RATIO) 120 CONTINUE C-------------------------------------------------------------------- C Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward C recurrence and K(ALPHA,X) from the wronskian C-------------------------------------------------------------------- D2 = AINT(ESTM(3)*EX+ESTM(4)) M = INT(D2) C = ABS(ENU) D3 = C+C D1 = D3-ONE F1 = XMIN F0 = (TWO*(C+D2)/EX+HALF*EX/(C+D2+ONE))*XMIN DO 130 I = 3,M D2 = D2-ONE F2 = (D3+D2+D2)*F0 BLPHA = (ONE+D1/D2)*(F2+BLPHA) F2 = F2/EX+F1 F1 = F0 F0 = F2 130 CONTINUE F1 = (D3+TWO)*F0/EX+F1 D1 = ZERO T1 = ONE DO 140 I = 1,7 D1 = C*D1+P(I) T1 = C*T1+Q(I) 140 CONTINUE P0 = EXP(C*(A+C*(P(8)-C*D1/T1)-LOG(EX)))/EX F2 = (C+HALF-RATIO)*F1/EX BK1 = P0+(D3*F0-F2+F0+BLPHA)/(F2+F1+F0)*P0 IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) WMINF = ESTF(3)*EX+ESTF(4) ELSE C-------------------------------------------------------------------- C Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by backward C recurrence, for X .GT. 4.0 C-------------------------------------------------------------------- DM = AINT(ESTM(5)/EX+ESTM(6)) M = INT(DM) D2 = DM-HALF D2 = D2*D2 D1 = DM+DM DO 160 I = 2,M DM = DM-ONE D1 = D1-TWO D2 = D2-D1 RATIO = (D3+D2)/(TWOX+D1-RATIO) BLPHA = (RATIO+RATIO*BLPHA)/DM 160 CONTINUE BK1 = ONE/((D+D*BLPHA)*SQRT(EX)) IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) WMINF = ESTF(5)*(EX-ABS(EX-ESTF(7)))+ESTF(6) END IF C-------------------------------------------------------------------- C Calculation of K(ALPHA+1,X) from K(ALPHA,X) and C K(ALPHA+1,X)/K(ALPHA,X) C-------------------------------------------------------------------- BK2 = BK1+BK1*(ENU+HALF-RATIO)/EX END IF C-------------------------------------------------------------------- C Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, C K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 C-------------------------------------------------------------------- NCALC = NB BK(1) = BK1 IF (IEND .EQ. 0) GO TO 500 J = 2-K IF (J .GT. 0) BK(J) = BK2 IF (IEND .EQ. 1) GO TO 500 M = MIN(INT(WMINF-ENU),IEND) DO 190 I = 2,M T1 = BK1 BK1 = BK2 TWONU = TWONU+TWO IF (EX .LT. ONE) THEN IF (BK1 .GE. (XINF/TWONU)*EX) GO TO 195 GO TO 187 ELSE IF (BK1/EX .GE. XINF/TWONU) GO TO 195 END IF 187 CONTINUE BK2 = TWONU/EX*BK1+T1 ITEMP = I J = J+1 IF (J .GT. 0) BK(J) = BK2 190 CONTINUE 195 M = ITEMP IF (M .EQ. IEND) GO TO 500 RATIO = BK2/BK1 MPLUS1 = M+1 NCALC = -1 DO 410 I = MPLUS1,IEND TWONU = TWONU+TWO RATIO = TWONU/EX+ONE/RATIO J = J+1 IF (J .GT. 1) THEN BK(J) = RATIO ELSE IF (BK2 .GE. XINF/RATIO) GO TO 500 BK2 = RATIO*BK2 END IF 410 CONTINUE NCALC = MAX(MPLUS1-K,1) IF (NCALC .EQ. 1) BK(1) = BK2 IF (NB .EQ. 1) GO TO 500 420 J = NCALC+1 DO 430 I = J,NB IF (BK(NCALC) .GE. XINF/BK(I)) GO TO 500 BK(I) = BK(NCALC)*BK(I) NCALC = I 430 CONTINUE END IF 500 RETURN C---------- Last line of RKBESL ---------- END subroutine rkbesl end module rjk