r9knus.f
6.55 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
subroutine r9knus (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
dimension alpha(15), beta(15), a(15), c0kcs(16), znu1cs(12)
external csevl,gamma, inits, r1mach
c
c series for c0k on the interval 0. to 2.50000d-01
c with weighted error 1.60e-17
c log weighted error 16.79
c significant figures required 15.99
c decimal places required 17.40
c
data c0k cs( 1) / .0601830572 42626108e0 /
data c0k cs( 2) / -.1536487143 3017286e0 /
data c0k cs( 3) / -.0117511760 08210492e0 /
data c0k cs( 4) / -.0008524878 88919795e0 /
data c0k cs( 5) / -.0000613298 38767496e0 /
data c0k cs( 6) / -.0000044052 28124551e0 /
data c0k cs( 7) / -.0000003163 12467283e0 /
data c0k cs( 8) / -.0000000227 10719382e0 /
data c0k cs( 9) / -.0000000016 30564460e0 /
data c0k cs(10) / -.0000000001 17069392e0 /
data c0k cs(11) / -.0000000000 08405206e0 /
data c0k cs(12) / -.0000000000 00603466e0 /
data c0k cs(13) / -.0000000000 00043326e0 /
data c0k cs(14) / -.0000000000 00003110e0 /
data c0k cs(15) / -.0000000000 00000223e0 /
data c0k cs(16) / -.0000000000 00000016e0 /
c
c series for znu1 on the interval -7.00000d-01 to 0.
c with weighted error 1.43e-17
c log weighted error 16.85
c significant figures required 16.08
c decimal places required 17.38
c
data znu1cs( 1) / .2033067569 9419173e0 /
data znu1cs( 2) / .1400779334 1321977e0 /
data znu1cs( 3) / .0079167969 61001613e0 /
data znu1cs( 4) / .0003398011 82532104e0 /
data znu1cs( 5) / .0000117419 75688989e0 /
data znu1cs( 6) / .0000003393 57570612e0 /
data znu1cs( 7) / .0000000084 25941769e0 /
data znu1cs( 8) / .0000000001 83336677e0 /
data znu1cs( 9) / .0000000000 03549698e0 /
data znu1cs(10) / .0000000000 00061903e0 /
data znu1cs(11) / .0000000000 00000981e0 /
data znu1cs(12) / .0000000000 00000014e0 /
c
data euler / 0.5772156649 0153286e0 /
data sqpi2 / 1.253314137 3155003e0 /
data aln2 / 0.693147180 55994531e0 /
data ntc0k, ntznu1 / 2*0 /
data xnusml, xsml, alnsml, alnbig, alneps / 5*0.0 /
c
if (ntc0k.ne.0) go to 10
ntc0k = inits (c0kcs, 16, 0.1*r1mach(3))
ntznu1 = inits (znu1cs, 12, 0.1*r1mach(3))
c
xnusml = sqrt (r1mach(3)/8.0)
xsml = 0.1*r1mach(3)
alnsml = alog (r1mach(1))
alnbig = alog (r1mach(2))
alneps = alog (0.1*r1mach(3))
c
10 if (xnu.lt.0. .or. xnu.ge.1.0) call seteru (
1 33hr9knus xnu must be ge 0 and lt 1, 33, 1, 2)
if (x.le.0.) call seteru (22hr9knus x must be gt 0, 22, 2, 2)
c
iswtch = 0
if (x.gt.2.0) 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.5) v = 1.0 - xnu
c
c carefully find (x/2)**xnu and z**xnu where z = x*x/4.
alnz = 2.0 * (alog(x) - aln2)
c
if (x.gt.xnu) go to 20
if (-0.5*xnu*alnz-aln2-alog(xnu).gt.alnbig) call seteru (
1 45hr9knus x so small bessel k-sub-xnu overflows, 45, 3, 2)
c
20 vlnz = v*alnz
x2tov = exp (0.5*vlnz)
ztov = 0.0
if (vlnz.gt.alnsml) ztov = x2tov**2
c
a0 = 0.5*gamma(1.0+v)
b0 = 0.5*gamma(1.0-v)
c0 = -euler
if (ztov.gt.0.5 .and. v.gt.xnusml) c0 = -0.75 +
1 csevl ((8.0*v)*v-1., c0kcs, ntc0k)
c
if (ztov.le.0.5) alpha(1) = (a0-ztov*b0)/v
if (ztov.gt.0.5) alpha(1) = c0 - alnz*(0.75 +
1 csevl (vlnz/0.35+1.0, znu1cs, ntznu1))*b0
beta(1) = -0.5*(a0+ztov*b0)
c
z = 0.0
if (x.gt.xsml) z = 0.25*x*x
nterms = max1 (2.0, 11.0+(8.*alnz-25.19-alneps)/(4.28-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.0*xi*a0)/(xi*(xi+v))
beta(i) = (xi-0.5*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.5*(xnu+1.)*alnz-2.0*aln2.gt.alnbig) iswtch = 1
if (iswtch.eq.1) return
bknud = expx*bknud*2.0/(x2tov*x)
c
if (xnu.le.0.5) bknu1 = v*bknu/x - bknud
if (xnu.le.0.5) return
c
bknu0 = bknu
bknu = -v*bknu/x - bknud
bknu1 = 2.0*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.0/xsml) go to 90
an = -1.56 + 4.0/x
bn = -0.29 - 0.22/x
nterms = min0 (15, max1 (3., an+bn*alneps))
c
do 80 inu=1,2
xmu = 0.
if (inu.eq.1 .and. xnu.gt.xnusml) xmu = (4.0*xnu)*xnu
if (inu.eq.2) xmu = 4.0*(abs(xnu)+1.)**2
c
a(1) = 1.0 - xmu
a(2) = 9.0 - xmu
a(3) = 25.0 - xmu
if (a(2).eq.0.) result = sqpi2*(16.*x+xmu+7.)/(16.*x*sqrtx)
if (a(2).eq.0.) go to 70
c
alpha(1) = 1.0
alpha(2) = (16.*x+a(2))/a(2)
alpha(3) = ((768.*x+48.*a(3))*x + a(2)*a(3))/(a(2)*a(3))
c
beta(1) = 1.0
beta(2) = (16.*x+(xmu+7.))/a(2)
beta(3) = ((768.*x+48.*(xmu+23.))*x + ((xmu+62.)*xmu+129.))
1 / (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.)**2 - xmu
qq = 16.*x2n/a(i)
p1 = -x2n*(float(12*n*n-20*n)-a(1))/((x2n-2.)*a(i)) - qq*x
p2 = (float(12*n*n-28*n+8)-a(1))/a(i) - qq*x
p3 = -x2n*a(i-3)/((x2n-2.)*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