subroutine d9knus(xnu,x,bknu,bknu1,iswtch)
c***begin prologue d9knus
c***date written 770601 (yymmdd)
c***revision date 820801 (yymmdd)
c***category no. c10b3
c***keywords bessel function,double precision,fnlib,special function
c***author fullerton, w., (lanl)
c***purpose computes bessel functions exp(x)*k-sub-xnu(x) and exp(x)*
c k-sub-xnu+1(x) for 0. .le. xnu .lt. 1.
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***references (none)
c***routines called d1mach,dcsevl,dgamma,initds,xerror
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
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 /
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 ntc0k, ntznu1 / 2*0 /
data xnusml, xsml, alnsml, alnbig / 4*0.d0 /
data alneps / 0.0 /
c***first executable statement d9knus
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 xerror ( 'd9knus xnu must
1 be ge 0 and lt 1', 33, 1, 2)
if (x.le.0.) call xerror ( 'd9knus 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 xerror ( 'd9knu
1s 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
