dbi.f
6.65 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
double precision function dbi (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, bifcs(13), bigcs(13), bif2cs(15), big2cs(15),
1 theta, xm, xmax, x3sml, z, d1mach, dcsevl, dbie, dexp,
2 dlog, dsin, dsqrt
external d1mach, dbie, dcsevl, initds
c
c series for bif on the interval -1.00000e+00 to 1.00000e+00
c with weighted error 1.45e-32
c log weighted error 31.84
c significant figures required 30.85
c decimal places required 32.40
c
data bif cs( 1) / -.1673021647 1986649483 5374239281 76 d-1 /
data bif cs( 2) / +.1025233583 4249445611 4263627777 57 d+0 /
data bif cs( 3) / +.1708309250 7381516539 4296502420 13 d-2 /
data bif cs( 4) / +.1186254546 7744681179 2164592100 40 d-4 /
data bif cs( 5) / +.4493290701 7792133694 5318879272 42 d-7 /
data bif cs( 6) / +.1069820714 3387889067 5677676636 28 d-9 /
data bif cs( 7) / +.1748064339 9771824706 0105176285 73 d-12 /
data bif cs( 8) / +.2081023107 1761711025 8818918343 99 d-15 /
data bif cs( 9) / +.1884981469 5665416509 9279717333 33 d-18 /
data bif cs( 10) / +.1342577917 3097804625 8826666666 66 d-21 /
data bif cs( 11) / +.7715959342 9658887893 3333333333 33 d-25 /
data bif cs( 12) / +.3653387961 7478566399 9999999999 99 d-28 /
data bif cs( 13) / +.1449756592 7953066666 6666666666 66 d-31 /
c
c series for big on the interval -1.00000e+00 to 1.00000e+00
c with weighted error 1.29e-33
c log weighted error 32.89
c significant figures required 31.48
c decimal places required 33.45
c
data big cs( 1) / +.2246622324 8574522283 4682201390 24 d-1 /
data big cs( 2) / +.3736477545 3019545441 7275616667 52 d-1 /
data big cs( 3) / +.4447621895 7212285696 2152943266 39 d-3 /
data big cs( 4) / +.2470807563 6329384245 4945919488 82 d-5 /
data big cs( 5) / +.7919135339 5149635134 8624262855 96 d-8 /
data big cs( 6) / +.1649807985 1827779880 8878724027 06 d-10 /
data big cs( 7) / +.2411990666 4835455909 2475011228 41 d-13 /
data big cs( 8) / +.2610373623 6091436985 1847812693 33 d-16 /
data big cs( 9) / +.2175308297 7160323853 1237920000 00 d-19 /
data big cs( 10) / +.1438694640 0390433219 4837333333 33 d-22 /
data big cs( 11) / +.7734912561 2083468629 3333333333 33 d-26 /
data big cs( 12) / +.3446929203 3849002666 6666666666 66 d-29 /
data big cs( 13) / +.1293891927 3216000000 0000000000 00 d-32 /
c
c series for bif2 on the interval 1.00000e+00 to 8.00000e+00
c with weighted error 6.08e-32
c log weighted error 31.22
c approx significant figures required 30.8
c decimal places required 31.80
c
data bif2cs( 1) / +.0998457269 3816041044 6828425799 3 d+0 /
data bif2cs( 2) / +.4786249778 6300553772 2114673182 31 d+0 /
data bif2cs( 3) / +.2515521196 0433011771 3244154366 75 d-1 /
data bif2cs( 4) / +.5820693885 2326456396 5156978722 16 d-3 /
data bif2cs( 5) / +.7499765964 4377865943 8614573782 17 d-5 /
data bif2cs( 6) / +.6134602870 3493836681 4030103564 74 d-7 /
data bif2cs( 7) / +.3462753885 1480632900 4342687333 59 d-9 /
data bif2cs( 8) / +.1428891008 0270254287 7708467489 31 d-11 /
data bif2cs( 9) / +.4496270429 8334641895 0564721792 00 d-14 /
data bif2cs( 10) / +.1114232306 5833011708 4283001066 66 d-16 /
data bif2cs( 11) / +.2230479106 6175002081 5178666666 66 d-19 /
data bif2cs( 12) / +.3681577873 6393142842 9226666666 66 d-22 /
data bif2cs( 13) / +.5096086844 9338261333 3333333333 33 d-25 /
data bif2cs( 14) / +.6000338692 6288554666 6666666666 66 d-28 /
data bif2cs( 15) / +.6082749744 6570666666 6666666666 66 d-31 /
c
c series for big2 on the interval 1.00000e+00 to 8.00000e+00
c with weighted error 4.91e-33
c log weighted error 32.31
c approx significant figures required 31.6
c decimal places required 32.90
c
data big2cs( 1) / +.0333056621 4551434046 5176188111 647 d+0 /
data big2cs( 2) / +.1613092151 2319706761 3287532084 943 d+0 /
data big2cs( 3) / +.6319007309 6134286912 1615634921 173 d-2 /
data big2cs( 4) / +.1187904568 1625173638 9780192304 567 d-3 /
data big2cs( 5) / +.1304534588 6200265614 7116485012 843 d-5 /
data big2cs( 6) / +.9374125995 5352172954 6809615508 936 d-8 /
data big2cs( 7) / +.4745801886 7472515378 8510169834 595 d-10 /
data big2cs( 8) / +.1783107265 0948139980 0065667560 946 d-12 /
data big2cs( 9) / +.5167591927 8495818037 4276356640 000 d-15 /
data big2cs( 10) / +.1190045083 8682712512 9496251733 333 d-17 /
data big2cs( 11) / +.2229828806 6640351727 7063466666 666 d-20 /
data big2cs( 12) / +.3465519230 2768941972 2666666666 666 d-23 /
data big2cs( 13) / +.4539263363 2050451413 3333333333 333 d-26 /
data big2cs( 14) / +.5078849965 1352234666 6666666666 666 d-29 /
data big2cs( 15) / +.4910206746 9653333333 3333333333 333 d-32 /
c
data nbif, nbig, nbif2, nbig2, x3sml, xmax / 4*0, 2*0.d0 /
c
if (nbif.ne.0) go to 10
eta = 0.1*sngl(d1mach(3))
nbif = initds (bifcs, 13, eta)
nbig = initds (bigcs, 13, eta)
nbif2 = initds (bif2cs, 15, eta)
nbig2 = initds (big2cs, 15, eta)
c
x3sml = eta**0.3333
xmax = (1.5*dlog(d1mach(2)))**0.6666d0
c
10 if (x.ge.(-1.0d0)) go to 20
call d9aimp (x, xm, theta)
dbi = xm * dsin(theta)
return
c
20 if (x.gt.1.0d0) go to 30
z = 0.d0
if (dabs(x).gt.x3sml) z = x**3
dbi = 0.625 + dcsevl (z, bifcs, nbif) + x*(0.4375d0 +
1 dcsevl (z, bigcs, nbig))
return
c
30 if (x.gt.2.0d0) go to 40
z = (2.0d0*x**3 - 9.0d0)/7.d0
dbi = 1.125d0 + dcsevl (z, bif2cs, nbif2) + x*(0.625d0 +
1 dcsevl (z, big2cs, nbig2))
return
c
40 if (x.gt.xmax) call seteru (
1 34hdbi x so big that bi overflows, 34, 1, 2)
c
dbi = dbie(x) * dexp(2.0d0*x*dsqrt(x)/3.0d0)
return
c
end