aie.f
4.66 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
function aie (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate ai(x) for x .le. 0.0 and ai(x)*exp(zeta) where
c zeta = 2/3 * x**(3/2) for x .ge. 0.0
c
dimension aifcs(9), aigcs(8), aipcs(34)
external csevl, inits, r1mach
c
c series for aif on the interval -1.00000d+00 to 1.00000d+00
c with weighted error 1.09e-19
c log weighted error 18.96
c significant figures required 17.76
c decimal places required 19.44
c
data aif cs( 1) / -.0379713584 9666999750e0 /
data aif cs( 2) / .0591918885 3726363857e0 /
data aif cs( 3) / .0009862928 0577279975e0 /
data aif cs( 4) / .0000068488 4381907656e0 /
data aif cs( 5) / .0000000259 4202596219e0 /
data aif cs( 6) / .0000000000 6176612774e0 /
data aif cs( 7) / .0000000000 0010092454e0 /
data aif cs( 8) / .0000000000 0000012014e0 /
data aif cs( 9) / .0000000000 0000000010e0 /
c
c series for aig on the interval -1.00000d+00 to 1.00000d+00
c with weighted error 1.51e-17
c log weighted error 16.82
c significant figures required 15.19
c decimal places required 17.27
c
data aig cs( 1) / .0181523655 8116127e0 /
data aig cs( 2) / .0215725631 6601076e0 /
data aig cs( 3) / .0002567835 6987483e0 /
data aig cs( 4) / .0000014265 2141197e0 /
data aig cs( 5) / .0000000045 7211492e0 /
data aig cs( 6) / .0000000000 0952517e0 /
data aig cs( 7) / .0000000000 0001392e0 /
data aig cs( 8) / .0000000000 0000001e0 /
c
c series for aip on the interval 0. to 1.00000d+00
c with weighted error 5.10e-17
c log weighted error 16.29
c significant figures required 14.41
c decimal places required 17.06
c
data aip cs( 1) / -.0187519297 793868e0 /
data aip cs( 2) / -.0091443848 250055e0 /
data aip cs( 3) / .0009010457 337825e0 /
data aip cs( 4) / -.0001394184 127221e0 /
data aip cs( 5) / .0000273815 815785e0 /
data aip cs( 6) / -.0000062750 421119e0 /
data aip cs( 7) / .0000016064 844184e0 /
data aip cs( 8) / -.0000004476 392158e0 /
data aip cs( 9) / .0000001334 635874e0 /
data aip cs(10) / -.0000000420 735334e0 /
data aip cs(11) / .0000000139 021990e0 /
data aip cs(12) / -.0000000047 831848e0 /
data aip cs(13) / .0000000017 047897e0 /
data aip cs(14) / -.0000000006 268389e0 /
data aip cs(15) / .0000000002 369824e0 /
data aip cs(16) / -.0000000000 918641e0 /
data aip cs(17) / .0000000000 364278e0 /
data aip cs(18) / -.0000000000 147475e0 /
data aip cs(19) / .0000000000 060851e0 /
data aip cs(20) / -.0000000000 025552e0 /
data aip cs(21) / .0000000000 010906e0 /
data aip cs(22) / -.0000000000 004725e0 /
data aip cs(23) / .0000000000 002076e0 /
data aip cs(24) / -.0000000000 000924e0 /
data aip cs(25) / .0000000000 000417e0 /
data aip cs(26) / -.0000000000 000190e0 /
data aip cs(27) / .0000000000 000087e0 /
data aip cs(28) / -.0000000000 000040e0 /
data aip cs(29) / .0000000000 000019e0 /
data aip cs(30) / -.0000000000 000009e0 /
data aip cs(31) / .0000000000 000004e0 /
data aip cs(32) / -.0000000000 000002e0 /
data aip cs(33) / .0000000000 000001e0 /
data aip cs(34) / -.0000000000 000000e0 /
c
data naif, naig, naip / 3*0 /
data x3sml, x32sml, xbig / 3*0.0 /
c
if (naif.ne.0) go to 10
eta = 0.1*r1mach(3)
naif = inits (aifcs , 9, eta)
naig = inits (aigcs , 8, eta)
naip = inits (aipcs , 34, eta)
c
x3sml = eta**0.3333
x32sml = 1.3104*x3sml**2
xbig = r1mach(2)**0.6666
c
10 if (x.ge.(-1.0)) go to 20
call r9aimp (x, xm, theta)
aie = xm * cos(theta)
return
c
20 if (x.gt.1.0) go to 30
z = 0.0
if (abs(x).gt.x3sml) z = x**3
aie = 0.375 + (csevl (z, aifcs, naif) - x*(0.25 +
1 csevl (z, aigcs, naig)) )
if (x.gt.x32sml) aie = aie * exp(2.0*x*sqrt(x)/3.0)
return
c
30 sqrtx = sqrt(x)
z = -1.0
if (x.lt.xbig) z = 2.0/(x*sqrtx) - 1.0
aie = (.28125 + csevl (z, aipcs, naip))/sqrt(sqrtx)
return
c
end