ai.f
2.5 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
function ai (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
dimension aifcs(9), aigcs(8)
external aie, 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
data naif, naig, x3sml, xmax / 2*0, 2*0.0 /
c
if (naif.ne.0) go to 10
naif = inits (aifcs, 9, 0.1*r1mach(3))
naig = inits (aigcs, 8, 0.1*r1mach(3))
c
x3sml = r1mach(3)**0.3334
xmax = (-1.5*alog(r1mach(1)))**0.6667
xmax = xmax - xmax*alog(xmax)/(4.0*xmax*sqrt(xmax)+1.0) - 0.01
c
10 if (x.ge.(-1.0)) go to 20
call r9aimp (x, xm, theta)
ai = 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
ai = 0.375 + (csevl (z, aifcs, naif) - x*(0.25 +
1 csevl (z, aigcs, naig)) )
return
c
30 if (x.gt.xmax) go to 40
ai = aie(x) * exp(-2.0*x*sqrt(x)/3.0)
return
c
40 ai = 0.0
call seteru (30hai x so big ai underflows, 30, 1, 0)
return
c
end