dbesk1.f
2.63 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
double precision function dbesk1 (x)
c july 1980 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, bk1cs(16), xmax, xmin, xsml, y,
1 d1mach, dcsevl, dbesi1, dbsk1e, dexp, dlog, dsqrt
external d1mach, dbesi1, dbsk1e, dcsevl, initds
c
c series for bk1 on the interval 0. to 4.00000e+00
c with weighted error 9.16e-32
c log weighted error 31.04
c significant figures required 30.61
c decimal places required 31.64
c
data bk1 cs( 1) / +.2530022733 8947770532 5311208685 33 d-1 /
data bk1 cs( 2) / -.3531559607 7654487566 7238316918 01 d+0 /
data bk1 cs( 3) / -.1226111808 2265714823 4790679300 42 d+0 /
data bk1 cs( 4) / -.6975723859 6398643501 8129202960 83 d-2 /
data bk1 cs( 5) / -.1730288957 5130520630 1765073689 79 d-3 /
data bk1 cs( 6) / -.2433406141 5659682349 6007350301 64 d-5 /
data bk1 cs( 7) / -.2213387630 7347258558 3152525451 26 d-7 /
data bk1 cs( 8) / -.1411488392 6335277610 9583302126 08 d-9 /
data bk1 cs( 9) / -.6666901694 1993290060 8537512643 73 d-12 /
data bk1 cs( 10) / -.2427449850 5193659339 2631968648 53 d-14 /
data bk1 cs( 11) / -.7023863479 3862875971 7837971200 00 d-17 /
data bk1 cs( 12) / -.1654327515 5100994675 4910293333 33 d-19 /
data bk1 cs( 13) / -.3233834745 9944491991 8933333333 33 d-22 /
data bk1 cs( 14) / -.5331275052 9265274999 4666666666 66 d-25 /
data bk1 cs( 15) / -.7513040716 2157226666 6666666666 66 d-28 /
data bk1 cs( 16) / -.9155085717 6541866666 6666666666 66 d-31 /
c
data ntk1, xmin, xsml, xmax / 0, 3*0.d0 /
c
if (ntk1.ne.0) go to 10
ntk1 = initds (bk1cs, 16, 0.1*sngl(d1mach(3)))
xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0)
xsml = dsqrt (4.0d0*d1mach(3))
xmax = -dlog(d1mach(1))
xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0) - 0.01d0
c
10 if (x.le.0.d0) call seteru (29hdbesk1 x is zero or negative,
1 29, 2, 2)
if (x.gt.2.0d0) go to 20
c
if (x.lt.xmin) call seteru (31hdbesk1 x so small k1 overflows,
1 31, 3, 2)
y = 0.d0
if (x.gt.xsml) y = x*x
dbesk1 = dlog(0.5d0*x)*dbesi1(x) + (0.75d0 + dcsevl (.5d0*y-1.d0,
1 bk1cs, ntk1))/x
return
c
20 dbesk1 = 0.d0
if (x.gt.xmax) call seteru (30hdbesk1 x so big k1 underflows,
1 30, 1, 0)
if (x.gt.xmax) return
c
dbesk1 = dexp(-x) * dbsk1e(x)
c
return
end