double precision function dai (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, aifcs(13), aigcs(13), theta, xm, xmax, x3sml,
     1  z, d1mach, dcsevl, daie, dcos, dexp, dlog, dsqrt
      external d1mach, daie, dcos, dcsevl, dexp, dlog, dsqrt, initds
c
c series for aif        on the interval -1.00000e+00 to  1.00000e+00
c                                        with weighted error   8.37e-33
c                                         log weighted error  32.08
c                               significant figures required  30.87
c                                    decimal places required  32.63
c
      data aif cs(  1) / -.3797135849 6669997496 1970894694 14 d-1     /
      data aif cs(  2) / +.5919188853 7263638574 3197280137 77 d-1     /
      data aif cs(  3) / +.9862928057 7279975365 6038910440 60 d-3     /
      data aif cs(  4) / +.6848843819 0765667554 8548301824 12 d-5     /
      data aif cs(  5) / +.2594202596 2194713019 4892790814 03 d-7     /
      data aif cs(  6) / +.6176612774 0813750329 4457496972 36 d-10    /
      data aif cs(  7) / +.1009245417 2466117901 4295562246 01 d-12    /
      data aif cs(  8) / +.1201479251 1179938141 2880332253 33 d-15    /
      data aif cs(  9) / +.1088294558 8716991878 5252954666 66 d-18    /
      data aif cs( 10) / +.7751377219 6684887039 2384000000 00 d-22    /
      data aif cs( 11) / +.4454811203 7175638391 4666666666 66 d-25    /
      data aif cs( 12) / +.2109284523 1692343466 6666666666 66 d-28    /
      data aif cs( 13) / +.8370173591 0741333333 3333333333 33 d-32    /
c
c series for aig        on the interval -1.00000e+00 to  1.00000e+00
c                                        with weighted error   7.47e-34
c                                         log weighted error  33.13
c                               significant figures required  31.50
c                                    decimal places required  33.68
c
      data aig cs(  1) / +.1815236558 1161273011 5562099578 64 d-1     /
      data aig cs(  2) / +.2157256316 6010755534 0306388199 68 d-1     /
      data aig cs(  3) / +.2567835698 7483249659 0524280901 33 d-3     /
      data aig cs(  4) / +.1426521411 9792403898 8294969217 21 d-5     /
      data aig cs(  5) / +.4572114920 0180426070 4340975581 91 d-8     /
      data aig cs(  6) / +.9525170843 5647098607 3922788405 92 d-11    /
      data aig cs(  7) / +.1392563460 5771399051 1504206861 90 d-13    /
      data aig cs(  8) / +.1507099914 2762379592 3069911386 66 d-16    /
      data aig cs(  9) / +.1255914831 2567778822 7032053333 33 d-19    /
      data aig cs( 10) / +.8306307377 0821340343 8293333333 33 d-23    /
      data aig cs( 11) / +.4465753849 3718567445 3333333333 33 d-26    /
      data aig cs( 12) / +.1990085503 4518869333 3333333333 33 d-29    /
      data aig cs( 13) / +.7470288525 6533333333 3333333333 33 d-33    /
c
      data naif, naig, x3sml, xmax / 2*0, 2*0.d0 /
c
      if (naif.ne.0) go to 10
      naif = initds (aifcs, 13, 0.1*sngl(d1mach(3)))
      naig = initds (aigcs, 13, 0.1*sngl(d1mach(3)))
c
      x3sml = d1mach(3)**0.3334d0
      xmax = (-1.5d0*dlog(d1mach(1)))**0.6667d0
      xmax = xmax - xmax*dlog(xmax)/(4.0d0*xmax*dsqrt(xmax)+1.0d0)
     1  - 0.01d0
c
 10   if (x.ge.(-1.d0)) go to 20
      call d9aimp (x, xm, theta)
      dai = xm * dcos(theta)
      return
c
 20   if (x.gt.1.0d0) go to 30
      z = 0.0d0
      if (dabs(x).gt.x3sml) z = x**3
      dai = 0.375d0 + (dcsevl (z, aifcs, naif) - x*(0.25d0 +
     1  dcsevl (z, aigcs, naig)) )
      return
c
 30   if (x.gt.xmax) go to 40
      dai = daie(x) * dexp(-2.0d0*x*dsqrt(x)/3.0d0)
      return
c
 40   dai = 0.0d0
      call seteru (30hdai     x so big ai underflows, 30, 1, 0)
      return
c
      end