Blame view
fvn_fnlib/ai.f
2.5 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 |
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 |