*DECK D9KNUS SUBROUTINE D9KNUS (XNU, X, BKNU, BKNU1, ISWTCH) C***BEGIN PROLOGUE D9KNUS C***SUBSIDIARY C***PURPOSE Compute Bessel functions EXP(X)*K-SUB-XNU(X) and EXP(X)* C K-SUB-XNU+1(X) for 0.0 .LE. XNU .LT. 1.0. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C10B3 C***TYPE DOUBLE PRECISION (R9KNUS-S, D9KNUS-D) C***KEYWORDS BESSEL FUNCTION, FNLIB, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C Compute Bessel functions EXP(X) * K-sub-XNU (X) and C EXP(X) * K-sub-XNU+1 (X) for 0.0 .LE. XNU .LT. 1.0 . C C Series for C0K on the interval 0. to 2.50000E-01 C with weighted error 2.16E-32 C log weighted error 31.67 C significant figures required 30.86 C decimal places required 32.40 C C Series for ZNU1 on the interval -7.00000E-01 to 0. C with weighted error 2.45E-33 C log weighted error 32.61 C significant figures required 31.85 C decimal places required 33.26 C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, DCSEVL, DGAMMA, INITDS, XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 890911 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900720 Routine changed from user-callable to subsidiary. (WRB) C 900727 Added EXTERNAL statement. (WRB) C 920618 Removed space from variable names. (RWC, WRB) C***END PROLOGUE D9KNUS DOUBLE PRECISION XNU, X, BKNU, BKNU1, ALPHA(32), BETA(32), A(32), 1 C0KCS(29), ZNU1CS(20), ALNZ, ALN2, A0, BKNUD, BKNU0, 2 B0, C0, EULER, EXPX, P1, P2, P3, QQ, RESULT, SQPI2, SQRTX, V, 3 VLNZ, XI, XMU, XNUSML, XSML, X2N, X2TOV, Z, ZTOV, ALNSML, 4 ALNBIG REAL ALNEPS DOUBLE PRECISION D1MACH, DCSEVL, DGAMMA LOGICAL FIRST EXTERNAL DGAMMA SAVE C0KCS, ZNU1CS, EULER, SQPI2, ALN2, NTC0K, 1 NTZNU1, XNUSML, XSML, ALNSML, ALNBIG, ALNEPS, FIRST DATA C0KCS( 1) / +.6018305724 2626108387 5774451803 29 D-1 / DATA C0KCS( 2) / -.1536487143 3017286092 9597559431 24 D+0 / DATA C0KCS( 3) / -.1175117600 8210492040 0682292262 13 D-1 / DATA C0KCS( 4) / -.8524878889 1979509827 0484015509 87 D-3 / DATA C0KCS( 5) / -.6132983876 7496791874 0981769221 11 D-4 / DATA C0KCS( 6) / -.4405228124 5510444562 6798895485 05 D-5 / DATA C0KCS( 7) / -.3163124672 8384488192 9154458921 99 D-6 / DATA C0KCS( 8) / -.2271071938 2899588330 6737717933 96 D-7 / DATA C0KCS( 9) / -.1630564460 8077609552 2746205153 60 D-8 / DATA C0KCS( 10) / -.1170693929 9414776568 7560440431 30 D-9 / DATA C0KCS( 11) / -.8405206378 6464437174 5465934137 92 D-11 / DATA C0KCS( 12) / -.6034667011 8979991487 0960507371 98 D-12 / DATA C0KCS( 13) / -.4332696033 5681371952 0459973669 03 D-13 / DATA C0KCS( 14) / -.3110735803 0203546214 6346977722 37 D-14 / DATA C0KCS( 15) / -.2233407822 6736982254 4861334098 40 D-15 / DATA C0KCS( 16) / -.1603514671 6864226300 6357915286 10 D-16 / DATA C0KCS( 17) / -.1151271736 3666556196 0356977053 05 D-17 / DATA C0KCS( 18) / -.8265759174 6836959105 1694790892 58 D-19 / DATA C0KCS( 19) / -.5934548080 6383948172 3334366959 84 D-20 / DATA C0KCS( 20) / -.4260813819 6467143926 4996130239 76 D-21 / DATA C0KCS( 21) / -.3059126686 4812876299 2636983705 42 D-22 / DATA C0KCS( 22) / -.2196354142 6734575224 9755018155 16 D-23 / DATA C0KCS( 23) / -.1576911326 1495836071 1057506847 60 D-24 / DATA C0KCS( 24) / -.1132171393 5950320948 7577310480 56 D-25 / DATA C0KCS( 25) / -.8128624883 4598404082 7923497144 33 D-27 / DATA C0KCS( 26) / -.5836090089 3453226552 8293493159 49 D-28 / DATA C0KCS( 27) / -.4190124162 3610922519 4523377809 05 D-29 / DATA C0KCS( 28) / -.3008373796 0206435069 5305042128 62 D-30 / DATA C0KCS( 29) / -.2159915206 7808647728 3421680898 32 D-31 / DATA ZNU1CS( 1) / +.2033067569 9419172967 4444001216 911 D+0 / DATA ZNU1CS( 2) / +.1400779334 1321977106 2943670790 563 D+0 / DATA ZNU1CS( 3) / +.7916796961 0016135284 0972241972 320 D-2 / DATA ZNU1CS( 4) / +.3398011825 3210404535 2930092205 750 D-3 / DATA ZNU1CS( 5) / +.1174197568 8989336666 4507228352 690 D-4 / DATA ZNU1CS( 6) / +.3393575706 1226168033 3825865475 121 D-6 / DATA ZNU1CS( 7) / +.8425941769 7621991019 4629891264 803 D-8 / DATA ZNU1CS( 8) / +.1833366770 2485008918 4748150900 090 D-9 / DATA ZNU1CS( 9) / +.3549698447 0441631086 3007064469 557 D-11 / DATA ZNU1CS( 10) / +.6190324964 6988733220 5244342078 407 D-13 / DATA ZNU1CS( 11) / +.9819645356 8043942496 0346115456 527 D-15 / DATA ZNU1CS( 12) / +.1428513143 9649047421 1473563005 985 D-16 / DATA ZNU1CS( 13) / +.1918949218 8782529896 6162467488 436 D-18 / DATA ZNU1CS( 14) / +.2394309797 3949891416 2313140597 128 D-20 / DATA ZNU1CS( 15) / +.2788902468 1534735483 5870465474 995 D-22 / DATA ZNU1CS( 16) / +.3046066506 3303344258 2845214092 865 D-24 / DATA ZNU1CS( 17) / +.3131732370 4219181577 1564260932 089 D-26 / DATA ZNU1CS( 18) / +.3041330989 8785495164 5174908005 034 D-28 / DATA ZNU1CS( 19) / +.2798403846 3683308434 3185097659 733 D-30 / DATA ZNU1CS( 20) / +.2446371862 7449759648 5238794922 666 D-32 / DATA EULER / 0.5772156649 0153286060 6512090082 40D0 / DATA SQPI2 / +1.253314137 3155002512 0788264240 55 D0 / DATA ALN2 / 0.6931471805 5994530941 7232121458 18D0 / DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT D9KNUS IF (FIRST) THEN ETA = 0.1D0*D1MACH(3) NTC0K = INITDS (C0KCS, 29, ETA) NTZNU1 = INITDS (ZNU1CS, 20, ETA) C XNUSML = SQRT(D1MACH(3)/8.D0) XSML = 0.1D0*D1MACH(3) ALNSML = LOG (D1MACH(1)) ALNBIG = LOG (D1MACH(2)) ALNEPS = LOG (0.1D0*D1MACH(3)) ENDIF FIRST = .FALSE. C IF (XNU .LT. 0.D0 .OR. XNU .GE. 1.D0) CALL XERMSG ('SLATEC', + 'D9KNUS', 'XNU MUST BE GE 0 AND LT 1', 1, 2) IF (X .LE. 0.) CALL XERMSG ('SLATEC', 'D9KNUS', 'X MUST BE GT 0', + 2, 2) C ISWTCH = 0 IF (X.GT.2.0D0) GO TO 50 C C X IS SMALL. COMPUTE K-SUB-XNU (X) AND THE DERIVATIVE OF K-SUB-XNU (X) C THEN FIND K-SUB-XNU+1 (X). XNU IS REDUCED TO THE INTERVAL (-.5,+.5) C THEN TO (0., .5), BECAUSE K OF NEGATIVE ORDER (-NU) = K OF POSITIVE C ORDER (+NU). C V = XNU IF (XNU.GT.0.5D0) V = 1.0D0 - XNU C C CAREFULLY FIND (X/2)**XNU AND Z**XNU WHERE Z = X*X/4. ALNZ = 2.D0 * (LOG(X) - ALN2) C IF (X.GT.XNU) GO TO 20 IF (-0.5D0*XNU*ALNZ-ALN2-LOG(XNU) .GT. ALNBIG) CALL XERMSG + ('SLATEC', 'D9KNUS', 'X SO SMALL BESSEL K-SUB-XNU OVERFLOWS', + 3, 2) C 20 VLNZ = V*ALNZ X2TOV = EXP (0.5D0*VLNZ) ZTOV = 0.0D0 IF (VLNZ.GT.ALNSML) ZTOV = X2TOV**2 C A0 = 0.5D0*DGAMMA(1.0D0+V) B0 = 0.5D0*DGAMMA(1.0D0-V) C0 = -EULER IF (ZTOV.GT.0.5D0 .AND. V.GT.XNUSML) C0 = -0.75D0 + 1 DCSEVL ((8.0D0*V)*V-1.0D0, C0KCS, NTC0K) C IF (ZTOV.LE.0.5D0) ALPHA(1) = (A0-ZTOV*B0)/V IF (ZTOV.GT.0.5D0) ALPHA(1) = C0 - ALNZ*(0.75D0 + 1 DCSEVL (VLNZ/0.35D0+1.0D0, ZNU1CS, NTZNU1))*B0 BETA(1) = -0.5D0*(A0+ZTOV*B0) C Z = 0.0D0 IF (X.GT.XSML) Z = 0.25D0*X*X NTERMS = MAX (2.0, 11.0+(8.*REAL(ALNZ)-25.19-ALNEPS) 1 /(4.28-REAL(ALNZ))) DO 30 I=2,NTERMS XI = I - 1 A0 = A0/(XI*(XI-V)) B0 = B0/(XI*(XI+V)) ALPHA(I) = (ALPHA(I-1)+2.0D0*XI*A0)/(XI*(XI+V)) BETA(I) = (XI-0.5D0*V)*ALPHA(I) - ZTOV*B0 30 CONTINUE C BKNU = ALPHA(NTERMS) BKNUD = BETA(NTERMS) DO 40 II=2,NTERMS I = NTERMS + 1 - II BKNU = ALPHA(I) + BKNU*Z BKNUD = BETA(I) + BKNUD*Z 40 CONTINUE C EXPX = EXP(X) BKNU = EXPX*BKNU/X2TOV C IF (-0.5D0*(XNU+1.D0)*ALNZ-2.0D0*ALN2.GT.ALNBIG) ISWTCH = 1 IF (ISWTCH.EQ.1) RETURN BKNUD = EXPX*BKNUD*2.0D0/(X2TOV*X) C IF (XNU.LE.0.5D0) BKNU1 = V*BKNU/X - BKNUD IF (XNU.LE.0.5D0) RETURN C BKNU0 = BKNU BKNU = -V*BKNU/X - BKNUD BKNU1 = 2.0D0*XNU*BKNU/X + BKNU0 RETURN C C X IS LARGE. FIND K-SUB-XNU (X) AND K-SUB-XNU+1 (X) WITH Y. L. LUKE-S C RATIONAL EXPANSION. C 50 SQRTX = SQRT(X) IF (X.GT.1.0D0/XSML) GO TO 90 AN = -0.60 - 1.02/REAL(X) BN = -0.27 - 0.53/REAL(X) NTERMS = MIN (32, MAX1 (3.0, AN+BN*ALNEPS)) C DO 80 INU=1,2 XMU = 0.D0 IF (INU.EQ.1 .AND. XNU.GT.XNUSML) XMU = (4.0D0*XNU)*XNU IF (INU.EQ.2) XMU = 4.0D0*(ABS(XNU)+1.D0)**2 C A(1) = 1.0D0 - XMU A(2) = 9.0D0 - XMU A(3) = 25.0D0 - XMU IF (A(2).EQ.0.D0) RESULT = SQPI2*(16.D0*X+XMU+7.D0) / 1 (16.D0*X*SQRTX) IF (A(2).EQ.0.D0) GO TO 70 C ALPHA(1) = 1.0D0 ALPHA(2) = (16.D0*X+A(2))/A(2) ALPHA(3) = ((768.D0*X+48.D0*A(3))*X + A(2)*A(3))/(A(2)*A(3)) C BETA(1) = 1.0D0 BETA(2) = (16.D0*X+(XMU+7.D0))/A(2) BETA(3) = ((768.D0*X+48.D0*(XMU+23.D0))*X + 1 ((XMU+62.D0)*XMU+129.D0))/(A(2)*A(3)) C IF (NTERMS.LT.4) GO TO 65 DO 60 I=4,NTERMS N = I - 1 X2N = 2*N - 1 C A(I) = (X2N+2.D0)**2 - XMU QQ = 16.D0*X2N/A(I) P1 = -X2N*((12*N*N-20*N)-A(1))/((X2N-2.D0)*A(I)) 1 - QQ*X P2 = ((12*N*N-28*N+8)-A(1))/A(I) - QQ*X P3 = -X2N*A(I-3)/((X2N-2.D0)*A(I)) C ALPHA(I) = -P1*ALPHA(I-1) - P2*ALPHA(I-2) - P3*ALPHA(I-3) BETA(I) = -P1*BETA(I-1) - P2*BETA(I-2) - P3*BETA(I-3) 60 CONTINUE C 65 RESULT = SQPI2*BETA(NTERMS)/(SQRTX*ALPHA(NTERMS)) C 70 IF (INU.EQ.1) BKNU = RESULT IF (INU.EQ.2) BKNU1 = RESULT 80 CONTINUE RETURN C 90 BKNU = SQPI2/SQRTX BKNU1 = BKNU RETURN C END