besi1e.f
4.78 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
function besi1e (x)
c oct 1983 version. w. fullerton, c3, los alamos scientific lab.
dimension bi1cs(11), ai1cs(21), ai12cs(22)
external csevl, inits, r1mach
c
c series for bi1 on the interval 0. to 9.00000d+00
c with weighted error 2.40e-17
c log weighted error 16.62
c significant figures required 16.23
c decimal places required 17.14
c
data bi1 cs( 1) / -.0019717132 61099859e0 /
data bi1 cs( 2) / .4073488766 7546481e0 /
data bi1 cs( 3) / .0348389942 99959456e0 /
data bi1 cs( 4) / .0015453945 56300123e0 /
data bi1 cs( 5) / .0000418885 21098377e0 /
data bi1 cs( 6) / .0000007649 02676483e0 /
data bi1 cs( 7) / .0000000100 42493924e0 /
data bi1 cs( 8) / .0000000000 99322077e0 /
data bi1 cs( 9) / .0000000000 00766380e0 /
data bi1 cs(10) / .0000000000 00004741e0 /
data bi1 cs(11) / .0000000000 00000024e0 /
c
c series for ai1 on the interval 1.25000d-01 to 3.33333d-01
c with weighted error 6.98e-17
c log weighted error 16.16
c significant figures required 14.53
c decimal places required 16.82
c
data ai1 cs( 1) / -.0284674418 1881479e0 /
data ai1 cs( 2) / -.0192295323 1443221e0 /
data ai1 cs( 3) / -.0006115185 8579437e0 /
data ai1 cs( 4) / -.0000206997 1253350e0 /
data ai1 cs( 5) / .0000085856 1914581e0 /
data ai1 cs( 6) / .0000010494 9824671e0 /
data ai1 cs( 7) / -.0000002918 3389184e0 /
data ai1 cs( 8) / -.0000000155 9378146e0 /
data ai1 cs( 9) / .0000000131 8012367e0 /
data ai1 cs(10) / -.0000000014 4842341e0 /
data ai1 cs(11) / -.0000000002 9085122e0 /
data ai1 cs(12) / .0000000001 2663889e0 /
data ai1 cs(13) / -.0000000000 1664947e0 /
data ai1 cs(14) / -.0000000000 0166665e0 /
data ai1 cs(15) / .0000000000 0124260e0 /
data ai1 cs(16) / -.0000000000 0027315e0 /
data ai1 cs(17) / .0000000000 0002023e0 /
data ai1 cs(18) / .0000000000 0000730e0 /
data ai1 cs(19) / -.0000000000 0000333e0 /
data ai1 cs(20) / .0000000000 0000071e0 /
data ai1 cs(21) / -.0000000000 0000006e0 /
c
c series for ai12 on the interval 0. to 1.25000d-01
c with weighted error 3.55e-17
c log weighted error 16.45
c significant figures required 14.69
c decimal places required 17.12
c
data ai12cs( 1) / .0285762350 1828014e0 /
data ai12cs( 2) / -.0097610974 9136147e0 /
data ai12cs( 3) / -.0001105889 3876263e0 /
data ai12cs( 4) / -.0000038825 6480887e0 /
data ai12cs( 5) / -.0000002512 2362377e0 /
data ai12cs( 6) / -.0000000263 1468847e0 /
data ai12cs( 7) / -.0000000038 3538039e0 /
data ai12cs( 8) / -.0000000005 5897433e0 /
data ai12cs( 9) / -.0000000000 1897495e0 /
data ai12cs(10) / .0000000000 3252602e0 /
data ai12cs(11) / .0000000000 1412580e0 /
data ai12cs(12) / .0000000000 0203564e0 /
data ai12cs(13) / -.0000000000 0071985e0 /
data ai12cs(14) / -.0000000000 0040836e0 /
data ai12cs(15) / -.0000000000 0002101e0 /
data ai12cs(16) / .0000000000 0004273e0 /
data ai12cs(17) / .0000000000 0001041e0 /
data ai12cs(18) / -.0000000000 0000382e0 /
data ai12cs(19) / -.0000000000 0000186e0 /
data ai12cs(20) / .0000000000 0000033e0 /
data ai12cs(21) / .0000000000 0000028e0 /
data ai12cs(22) / -.0000000000 0000003e0 /
c
data nti1, ntai1, ntai12, xmin, xsml / 3*0, 2*0. /
c
if (nti1.ne.0) go to 10
nti1 = inits (bi1cs, 11, 0.1*r1mach(3))
ntai1 = inits (ai1cs, 21, 0.1*r1mach(3))
ntai12 = inits (ai12cs, 22, 0.1*r1mach(3))
c
xmin = 2.0*r1mach(1)
xsml = sqrt (8.0*r1mach(3))
c
10 y = abs(x)
if (y.gt.3.0) go to 20
c
besi1e = 0.0
if (y.eq.0.0) return
c
if (y.le.xmin) call seteru (
1 37hbesi1e abs(x) so small i1 underflows, 37, 1, 0)
besi1e = 0.0
if (y.gt.xmin) besi1e = 0.5*x
if (y.gt.xsml) besi1e = x * (.875 + csevl(y*y/4.5-1., bi1cs,nti1))
besi1e = exp(-y) * besi1e
return
c
20 if (y.le.8.) besi1e = (.375 + csevl ((48./y-11.)/5., ai1cs, ntai1)
1 ) / sqrt(y)
if (y.gt.8.) besi1e = (.375 + csevl (16./y-1.0, ai12cs, ntai12))
1 / sqrt(y)
besi1e = sign (besi1e, x)
c
return
end