itairy Subroutine

subroutine itairy ( x, apt, bpt, ant, bnt )

  !****************************************************************************80
  !
  !! ITAIRY computes the integrals of Airy functions.
  !
  !  Discussion:
  !
  !    Compute the integrals of Airy functions with respect to t,
  !    from 0 and x.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    19 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) X, the upper limit of the integral.
  !
  !    Output, real ( kind = 8 ) APT, BPT, ANT, BNT, the integrals, from 0 to x,
  !    of Ai(t), Bi(t), Ai(-t), and Bi(-t).
  !       
  implicit none

  real ( kind = 8 ), save, dimension ( 16 ) :: a = (/ &
       0.569444444444444D+00, 0.891300154320988D+00, &
       0.226624344493027D+01, 0.798950124766861D+01, &
       0.360688546785343D+02, 0.198670292131169D+03, &
       0.129223456582211D+04, 0.969483869669600D+04, &
       0.824184704952483D+05, 0.783031092490225D+06, &
       0.822210493622814D+07, 0.945557399360556D+08, &
       0.118195595640730D+10, 0.159564653040121D+11, &
       0.231369166433050D+12, 0.358622522796969D+13 /)
  real ( kind = 8 ) ant
  real ( kind = 8 ) apt
  real ( kind = 8 ) bnt
  real ( kind = 8 ) bpt
  real ( kind = 8 ) c1
  real ( kind = 8 ) c2
  real ( kind = 8 ) eps
  real ( kind = 8 ) fx
  real ( kind = 8 ) gx
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  real ( kind = 8 ) pi
  real ( kind = 8 ) q0
  real ( kind = 8 ) q1
  real ( kind = 8 ) q2
  real ( kind = 8 ) r
  real ( kind = 8 ) sr3
  real ( kind = 8 ) su1
  real ( kind = 8 ) su2
  real ( kind = 8 ) su3
  real ( kind = 8 ) su4
  real ( kind = 8 ) su5
  real ( kind = 8 ) su6
  real ( kind = 8 ) x
  real ( kind = 8 ) xe
  real ( kind = 8 ) xp6
  real ( kind = 8 ) xr1
  real ( kind = 8 ) xr2

  eps = 1.0D-15
  pi = 3.141592653589793D+00
  c1 = 0.355028053887817D+00
  c2 = 0.258819403792807D+00
  sr3 = 1.732050807568877D+00

  if ( x == 0.0D+00 ) then

     apt = 0.0D+00
     bpt = 0.0D+00
     ant = 0.0D+00
     bnt = 0.0D+00

  else

     if ( abs ( x ) <= 9.25D+00 ) then

        do l = 0, 1

           x = ( -1.0D+00 ) ** l * x
           fx = x
           r = x

           do k = 1, 40
              r = r * ( 3.0D+00 * k - 2.0D+00 ) &
                   / ( 3.0D+00 * k + 1.0D+00 ) * x / ( 3.0D+00 * k ) &
                   * x / ( 3.0D+00 * k - 1.0D+00 ) * x 
              fx = fx + r
              if ( abs ( r ) < abs ( fx ) * eps ) then
                 exit
              end if
           end do

           gx = 0.5D+00 * x * x
           r = gx

           do k = 1, 40
              r = r * ( 3.0D+00 * k - 1.0D+00 ) &
                   / ( 3.0D+00 * k + 2.0D+00 ) * x / ( 3.0D+00 * k ) * x &
                   / ( 3.0D+00 * k + 1.0D+00 ) * x
              gx = gx + r
              if ( abs ( r ) < abs ( gx ) * eps ) then
                 exit
              end if
           end do

           ant = c1 * fx - c2 * gx
           bnt = sr3 * ( c1 * fx + c2 * gx )

           if ( l == 0 ) then
              apt = ant
              bpt = bnt
           else
              ant = -ant
              bnt = -bnt
              x = -x
           end if

        end do

     else

        q2 = 1.414213562373095D+00
        q0 = 0.3333333333333333D+00
        q1 = 0.6666666666666667D+00
        xe = x * sqrt ( x ) / 1.5D+00
        xp6 = 1.0D+00 / sqrt ( 6.0D+00 * pi * xe )
        su1 = 1.0D+00
        r = 1.0D+00
        xr1 = 1.0D+00 / xe
        do k = 1, 16
           r = - r * xr1
           su1 = su1 + a(k) * r
        end do
        su2 = 1.0D+00
        r = 1.0D+00
        do k = 1, 16
           r = r * xr1
           su2 = su2 + a(k) * r
        end do

        apt = q0 - exp ( - xe ) * xp6 * su1
        bpt = 2.0D+00 * exp ( xe ) * xp6 * su2
        su3 = 1.0D+00
        r = 1.0D+00
        xr2 = 1.0D+00 / ( xe * xe )
        do k = 1, 8
           r = - r * xr2
           su3 = su3 + a(2*k) * r
        end do
        su4 = a(1) * xr1
        r = xr1
        do k = 1, 7
           r = -r * xr2
           su4 = su4 + a(2*k+1) * r
        end do
        su5 = su3 + su4
        su6 = su3 - su4
        ant = q1 - q2 * xp6 * ( su5 * cos ( xe ) - su6 * sin ( xe ) )
        bnt = q2 * xp6 * ( su5 * sin ( xe ) + su6 * cos ( xe ) )

     end if

  end if

  return
end subroutine itairy