Blame view
fvn_fnlib/aie.f
4.66 KB
38581db0c
git-svn-id: https...
|
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
|