dbesk0.f
2.45 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
double precision function dbesk0 (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, bk0cs(16), xmax, xsml, y,
1 d1mach, dcsevl, dbesi0, dbsk0e, dexp, dlog, dsqrt
external d1mach, dbesi0, dbsk0e, dcsevl, initds
c
c series for bk0 on the interval 0. to 4.00000e+00
c with weighted error 3.08e-33
c log weighted error 32.51
c significant figures required 32.05
c decimal places required 33.11
c
data bk0 cs( 1) / -.3532739323 3902768720 1140060063 153 d-1 /
data bk0 cs( 2) / +.3442898999 2462848688 6344927529 213 d+0 /
data bk0 cs( 3) / +.3597993651 5361501626 5721303687 231 d-1 /
data bk0 cs( 4) / +.1264615411 4469259233 8479508673 447 d-2 /
data bk0 cs( 5) / +.2286212103 1194517860 8269830297 585 d-4 /
data bk0 cs( 6) / +.2534791079 0261494573 0790013428 354 d-6 /
data bk0 cs( 7) / +.1904516377 2202088589 7214059381 366 d-8 /
data bk0 cs( 8) / +.1034969525 7633624585 1008317853 089 d-10 /
data bk0 cs( 9) / +.4259816142 7910825765 2445327170 133 d-13 /
data bk0 cs( 10) / +.1374465435 8807508969 4238325440 000 d-15 /
data bk0 cs( 11) / +.3570896528 5083735909 9688597333 333 d-18 /
data bk0 cs( 12) / +.7631643660 1164373766 7498666666 666 d-21 /
data bk0 cs( 13) / +.1365424988 4407818590 8053333333 333 d-23 /
data bk0 cs( 14) / +.2075275266 9066680831 9999999999 999 d-26 /
data bk0 cs( 15) / +.2712814218 0729856000 0000000000 000 d-29 /
data bk0 cs( 16) / +.3082593887 9146666666 6666666666 666 d-32 /
c
data ntk0, xsml, xmax / 0, 2*0.d0 /
c
if (ntk0.ne.0) go to 10
ntk0 = initds (bk0cs, 16, 0.1*sngl(d1mach(3)))
xsml = dsqrt (4.0d0*d1mach(3))
xmax = -dlog(d1mach(1))
xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0)
c
10 if (x.le.0.d0) call seteru (29hdbesk0 x is zero or negative,
1 29, 2, 2)
if (x.gt.2.0d0) go to 20
c
y = 0.d0
if (x.gt.xsml) y = x*x
dbesk0 = -dlog(0.5d0*x)*dbesi0(x) - 0.25d0 + dcsevl (.5d0*y-1.d0,
1 bk0cs, ntk0)
return
c
20 dbesk0 = 0.d0
if (x.gt.xmax) call seteru (30hdbesk0 x so big k0 underflows,
1 30, 1, 0)
if (x.gt.xmax) return
c
dbesk0 = dexp(-x) * dbsk0e(x)
c
return
end