function ai (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      dimension aifcs(9), aigcs(8)
      external aie, alog, cos, csevl, exp, inits, r1mach, sqrt
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