d9knus.f
9.21 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
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, dgamma, 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