dbsk1e.f
8.42 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
double precision function dbsk1e (x)
c july 1980 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, bk1cs(16), ak1cs(38), ak12cs(33), xmin,
1 xsml, y, d1mach, dcsevl, dbesi1, dexp, dlog, dsqrt
external d1mach, dbesi1, 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
c series for ak1 on the interval 1.25000e-01 to 5.00000e-01
c with weighted error 3.07e-32
c log weighted error 31.51
c significant figures required 30.71
c decimal places required 32.30
c
data ak1 cs( 1) / +.2744313406 9738829695 2576662272 66 d+0 /
data ak1 cs( 2) / +.7571989953 1993678170 8923781492 90 d-1 /
data ak1 cs( 3) / -.1441051556 4754061229 8531161756 25 d-2 /
data ak1 cs( 4) / +.6650116955 1257479394 2513854770 36 d-4 /
data ak1 cs( 5) / -.4369984709 5201407660 5808450891 67 d-5 /
data ak1 cs( 6) / +.3540277499 7630526799 4171390085 34 d-6 /
data ak1 cs( 7) / -.3311163779 2932920208 9826882457 04 d-7 /
data ak1 cs( 8) / +.3445977581 9010534532 3114997709 92 d-8 /
data ak1 cs( 9) / -.3898932347 4754271048 9819374927 58 d-9 /
data ak1 cs( 10) / +.4720819750 4658356400 9474493390 05 d-10 /
data ak1 cs( 11) / -.6047835662 8753562345 3735915628 90 d-11 /
data ak1 cs( 12) / +.8128494874 8658747888 1938379856 63 d-12 /
data ak1 cs( 13) / -.1138694574 7147891428 9239159510 42 d-12 /
data ak1 cs( 14) / +.1654035840 8462282325 9729482050 90 d-13 /
data ak1 cs( 15) / -.2480902567 7068848221 5160104405 33 d-14 /
data ak1 cs( 16) / +.3829237890 7024096948 4292272991 57 d-15 /
data ak1 cs( 17) / -.6064734104 0012418187 7682103773 86 d-16 /
data ak1 cs( 18) / +.9832425623 2648616038 1940046506 66 d-17 /
data ak1 cs( 19) / -.1628416873 8284380035 6666201156 26 d-17 /
data ak1 cs( 20) / +.2750153649 6752623718 2841203370 66 d-18 /
data ak1 cs( 21) / -.4728966646 3953250924 2810695680 00 d-19 /
data ak1 cs( 22) / +.8268150002 8109932722 3920503466 66 d-20 /
data ak1 cs( 23) / -.1468140513 6624956337 1939648853 33 d-20 /
data ak1 cs( 24) / +.2644763926 9208245978 0858948266 66 d-21 /
data ak1 cs( 25) / -.4829015756 4856387897 9698688000 00 d-22 /
data ak1 cs( 26) / +.8929302074 3610130180 6563327999 99 d-23 /
data ak1 cs( 27) / -.1670839716 8972517176 9977514666 66 d-23 /
data ak1 cs( 28) / +.3161645603 4040694931 3686186666 66 d-24 /
data ak1 cs( 29) / -.6046205531 2274989106 5064106666 66 d-25 /
data ak1 cs( 30) / +.1167879894 2042732700 7184213333 33 d-25 /
data ak1 cs( 31) / -.2277374158 2653996232 8678400000 00 d-26 /
data ak1 cs( 32) / +.4481109730 0773675795 3058133333 33 d-27 /
data ak1 cs( 33) / -.8893288476 9020194062 3360000000 00 d-28 /
data ak1 cs( 34) / +.1779468001 8850275131 3920000000 00 d-28 /
data ak1 cs( 35) / -.3588455596 7329095821 9946666666 66 d-29 /
data ak1 cs( 36) / +.7290629049 2694257991 6799999999 99 d-30 /
data ak1 cs( 37) / -.1491844984 5546227073 0240000000 00 d-30 /
data ak1 cs( 38) / +.3073657387 2934276300 7999999999 99 d-31 /
c
c series for ak12 on the interval 0. to 1.25000e-01
c with weighted error 2.41e-32
c log weighted error 31.62
c significant figures required 30.25
c decimal places required 32.38
c
data ak12cs( 1) / +.6379308343 7390010366 0048853410 2 d-1 /
data ak12cs( 2) / +.2832887813 0497209358 3503028470 8 d-1 /
data ak12cs( 3) / -.2475370673 9052503454 1454556673 2 d-3 /
data ak12cs( 4) / +.5771972451 6072488204 7097662576 3 d-5 /
data ak12cs( 5) / -.2068939219 5365483027 4553319655 2 d-6 /
data ak12cs( 6) / +.9739983441 3818041803 0921309788 7 d-8 /
data ak12cs( 7) / -.5585336140 3806249846 8889551112 9 d-9 /
data ak12cs( 8) / +.3732996634 0461852402 2121285473 1 d-10 /
data ak12cs( 9) / -.2825051961 0232254451 3506575492 8 d-11 /
data ak12cs( 10) / +.2372019002 4841441736 4349695548 6 d-12 /
data ak12cs( 11) / -.2176677387 9917539792 6830166793 8 d-13 /
data ak12cs( 12) / +.2157914161 6160324539 3956268970 6 d-14 /
data ak12cs( 13) / -.2290196930 7182692759 9155133815 4 d-15 /
data ak12cs( 14) / +.2582885729 8232749619 1993956522 6 d-16 /
data ak12cs( 15) / -.3076752641 2684631876 2109817344 0 d-17 /
data ak12cs( 16) / +.3851487721 2804915970 9489684479 9 d-18 /
data ak12cs( 17) / -.5044794897 6415289771 1728250880 0 d-19 /
data ak12cs( 18) / +.6888673850 4185442370 1829222399 9 d-20 /
data ak12cs( 19) / -.9775041541 9501183030 0213248000 0 d-21 /
data ak12cs( 20) / +.1437416218 5238364610 0165973333 3 d-21 /
data ak12cs( 21) / -.2185059497 3443473734 9973333333 3 d-22 /
data ak12cs( 22) / +.3426245621 8092206316 4538880000 0 d-23 /
data ak12cs( 23) / -.5531064394 2464082325 0124800000 0 d-24 /
data ak12cs( 24) / +.9176601505 6859954037 8282666666 6 d-25 /
data ak12cs( 25) / -.1562287203 6180249114 4874666666 6 d-25 /
data ak12cs( 26) / +.2725419375 4843331323 4943999999 9 d-26 /
data ak12cs( 27) / -.4865674910 0748279923 7802666666 6 d-27 /
data ak12cs( 28) / +.8879388552 7235025873 5786666666 6 d-28 /
data ak12cs( 29) / -.1654585918 0392575489 3653333333 3 d-28 /
data ak12cs( 30) / +.3145111321 3578486743 0399999999 9 d-29 /
data ak12cs( 31) / -.6092998312 1931276124 1600000000 0 d-30 /
data ak12cs( 32) / +.1202021939 3698158346 2399999999 9 d-30 /
data ak12cs( 33) / -.2412930801 4594088413 8666666666 6 d-31 /
c
data ntk1, ntak1, ntak12, xmin, xsml / 3*0, 2*0.d0 /
c
if (ntk1.ne.0) go to 10
eta = 0.1*sngl(d1mach(3))
ntk1 = initds (bk1cs, 16, eta)
ntak1 = initds (ak1cs, 38, eta)
ntak12 = initds (ak12cs, 33, eta)
c
xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0)
xsml = dsqrt (4.0d0*d1mach(3))
c
10 if (x.le.0.d0) call seteru (29hdbsk1e x is zero or negative, 29,
1 1, 2)
if (x.gt.2.0d0) go to 20
c
if (x.lt.xmin) call seteru (31hdbsk1e x so small k1 overflows,
1 31, 2, 2)
y = 0.d0
if (x.gt.xsml) y = x*x
dbsk1e = dexp(x)*(dlog(0.5d0*x)*dbesi1(x) + (0.75d0 +
1 dcsevl (0.5d0*y-1.d0, bk1cs, ntk1))/x )
return
c
20 if (x.le.8.d0) dbsk1e = (1.25d0 + dcsevl ((16.d0/x-5.d0)/3.d0,
1 ak1cs, ntak1))/dsqrt(x)
if (x.gt.8.d0) dbsk1e = (1.25d0 +
1 dcsevl (16.d0/x-1.d0, ak12cs, ntak12))/dsqrt(x)
c
return
end