beskes.f
1.3 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
subroutine beskes (xnu, x, nin, bke)
c july 1980 edition. w. fullerton, c3, los alamos scientific lab.
dimension bke(1)
external r1mach
data alnbig / 0. /
c
if (alnbig.eq.0.) alnbig = alog (r1mach(2))
c
v = abs(xnu)
n = iabs(nin)
c
if (v.ge.1.) call seteru (29hbeskes abs(xnu) must be lt 1, 29,
1 1, 2)
if (x.le.0.) call seteru (17hbeskes x is le 0, 17, 2, 2)
if (n.eq.0) call seteru (
1 41hbeskes n the number in the sequence is 0, 41, 3, 2)
c
call r9knus (v, x, bke(1), bknu1, iswtch)
if (n.eq.1) return
c
vincr = sign (1.0, float(nin))
direct = vincr
if (xnu.ne.0.) direct = vincr*sign(1.0,xnu)
if (iswtch.eq.1 .and. direct.gt.0.) call seteru (
1 47hbeskes x so small bessel k-sub-xnu+1 overflows, 47, 4, 2)
bke(2) = bknu1
c
if (direct.lt.0.) call r9knus (abs(xnu+vincr), x, bke(2), bknu1,
1 iswtch)
if (n.eq.2) return
c
vend = abs(xnu+float(nin)) - 1.0
if (x+(vend-0.5)*alog(vend)+.27 - vend*(alog(x)+.306).gt.alnbig)
1 call seteru (67hbeskes x so small or abs(nu) so big that bessel
2 k-sub-nu overflows, 67, 5, 2)
c
v = xnu
do 10 i=3,n
v = v + vincr
bke(i) = 2.0*v*bke(i-1)/x + bke(i-2)
10 continue
c
return
end