subroutine d9knus (xnu, x, bknu, bknu1, iswtch) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. 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 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, an, bn double precision d1mach, dcsevl, dexp, dgamma, dlog, dsqrt external d1mach, dcsevl, dexp, dgamma, dlog, dsqrt, initds 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 data c0k cs( 1) / +.6018305724 2626108387 5774451803 29 d-1 / data c0k cs( 2) / -.1536487143 3017286092 9597559431 24 d+0 / data c0k cs( 3) / -.1175117600 8210492040 0682292262 13 d-1 / data c0k cs( 4) / -.8524878889 1979509827 0484015509 87 d-3 / data c0k cs( 5) / -.6132983876 7496791874 0981769221 11 d-4 / data c0k cs( 6) / -.4405228124 5510444562 6798895485 05 d-5 / data c0k cs( 7) / -.3163124672 8384488192 9154458921 99 d-6 / data c0k cs( 8) / -.2271071938 2899588330 6737717933 96 d-7 / data c0k cs( 9) / -.1630564460 8077609552 2746205153 60 d-8 / data c0k cs( 10) / -.1170693929 9414776568 7560440431 30 d-9 / data c0k cs( 11) / -.8405206378 6464437174 5465934137 92 d-11 / data c0k cs( 12) / -.6034667011 8979991487 0960507371 98 d-12 / data c0k cs( 13) / -.4332696033 5681371952 0459973669 03 d-13 / data c0k cs( 14) / -.3110735803 0203546214 6346977722 37 d-14 / data c0k cs( 15) / -.2233407822 6736982254 4861334098 40 d-15 / data c0k cs( 16) / -.1603514671 6864226300 6357915286 10 d-16 / data c0k cs( 17) / -.1151271736 3666556196 0356977053 05 d-17 / data c0k cs( 18) / -.8265759174 6836959105 1694790892 58 d-19 / data c0k cs( 19) / -.5934548080 6383948172 3334366959 84 d-20 / data c0k cs( 20) / -.4260813819 6467143926 4996130239 76 d-21 / data c0k cs( 21) / -.3059126686 4812876299 2636983705 42 d-22 / data c0k cs( 22) / -.2196354142 6734575224 9755018155 16 d-23 / data c0k cs( 23) / -.1576911326 1495836071 1057506847 60 d-24 / data c0k cs( 24) / -.1132171393 5950320948 7577310480 56 d-25 / data c0k cs( 25) / -.8128624883 4598404082 7923497144 33 d-27 / data c0k cs( 26) / -.5836090089 3453226552 8293493159 49 d-28 / data c0k cs( 27) / -.4190124162 3610922519 4523377809 05 d-29 / data c0k cs( 28) / -.3008373796 0206435069 5305042128 62 d-30 / data c0k cs( 29) / -.2159915206 7808647728 3421680898 32 d-31 / 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 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 / c data euler / 0.5772156649 0153286060 6512090082 40d0 / data sqpi2 / +1.253314137 3155002512 0788264240 55 d0 / data aln2 / 0.6931471805 5994530941 7232121458 18d0 / data ntc0k, ntznu1 / 2*0 / data xnusml, xsml, alnsml, alnbig, alneps / 4*0.d0, 0.0 / c if (ntc0k.ne.0) go to 10 eta = 0.1d0*d1mach(3) ntc0k = initds (c0kcs, 29, eta) ntznu1 = initds (znu1cs, 20, eta) c xnusml = dsqrt (d1mach(3)/8.d0) xsml = 0.1d0*d1mach(3) alnsml = dlog (d1mach(1)) alnbig = dlog (d1mach(2)) alneps = dlog (0.1d0*d1mach(3)) c 10 if (xnu.lt.0.d0 .or. xnu.ge.1.d0) call seteru ( 1 33hd9knus xnu must be ge 0 and lt 1, 33, 1, 2) if (x.le.0.) call seteru (22hd9knus x must be gt 0, 22, 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 * (dlog(x) - aln2) c if (x.gt.xnu) go to 20 if (-0.5d0*xnu*alnz-aln2-dlog(xnu).gt.alnbig) call seteru ( 1 45hd9knus x so small bessel k-sub-xnu overflows, 45, 3, 2) c 20 vlnz = v*alnz x2tov = dexp (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 = max1 (2.0, 11.0+(8.*sngl(alnz)-25.19-alneps) 1 /(4.28-sngl(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 = dexp(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 = dsqrt(x) if (x.gt.1.0d0/xsml) go to 90 an = -0.60 - 1.02/sngl(x) bn = -0.27 - 0.53/sngl(x) nterms = min0 (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*(dabs(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*(dble(float(12*n*n-20*n))-a(1))/((x2n-2.d0)*a(i)) 1 - qq*x p2 = (dble(float(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