airyb Subroutine

subroutine airyb ( x, ai, bi, ad, bd )

  !*****************************************************************************80
  !
  !! AIRYB computes Airy functions and their derivatives.
  !
  !  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:
  !
  !    02 June 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, argument of Airy function.
  !
  !    Output, real ( kind = 8 ) AI, Ai(x).
  !
  !    Output, real ( kind = 8 ) BI, Bi(x).
  !
  !    Output, real ( kind = 8 ) AD, Ai'(x).
  !
  !    Output, real ( kind = 8 ) BD, Bi'(x).
  !
  implicit none

  real ( kind = 8 ) ad
  real ( kind = 8 ) ai
  real ( kind = 8 ) bd
  real ( kind = 8 ) bi
  real ( kind = 8 ) c1
  real ( kind = 8 ) c2
  real ( kind = 8 ) ck(41)
  real ( kind = 8 ) df
  real ( kind = 8 ) dg
  real ( kind = 8 ) dk(41)
  real ( kind = 8 ) eps
  real ( kind = 8 ) fx
  real ( kind = 8 ) gx
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) rp
  real ( kind = 8 ) sad
  real ( kind = 8 ) sai
  real ( kind = 8 ) sbd
  real ( kind = 8 ) sbi
  real ( kind = 8 ) sda
  real ( kind = 8 ) sdb
  real ( kind = 8 ) sr3
  real ( kind = 8 ) ssa
  real ( kind = 8 ) ssb
  real ( kind = 8 ) x
  real ( kind = 8 ) xa
  real ( kind = 8 ) xar
  real ( kind = 8 ) xcs
  real ( kind = 8 ) xe
  real ( kind = 8 ) xf
  real ( kind = 8 ) xm
  real ( kind = 8 ) xp1
  real ( kind = 8 ) xq
  real ( kind = 8 ) xr1
  real ( kind = 8 ) xr2
  real ( kind = 8 ) xss

  eps = 1.0D-15
  pi = 3.141592653589793D+00
  c1 = 0.355028053887817D+00
  c2 = 0.258819403792807D+00
  sr3 = 1.732050807568877D+00
  xa = abs ( x )
  xq = sqrt ( xa )

  if ( x <= 0.0D+00 ) then
     xm = 8.0D+00
  else
     xm = 5.0D+00
  end if

  if ( x == 0.0D+00 ) then
     ai = c1
     bi = sr3 * c1
     ad = -c2
     bd = sr3 * c2
     return
  end if

  if ( xa <= xm ) then

     fx = 1.0D+00
     r = 1.0D+00
     do k = 1, 40
        r = r * 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 = x
     r = x
     do k = 1, 40
        r = r * 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

     ai = c1 * fx - c2 * gx
     bi = sr3 * ( c1 * fx + c2 * gx )
     df = 0.5D+00 * x * x
     r = df
     do k = 1, 40
        r = r * x / ( 3.0D+00 * k ) * x / ( 3.0D+00 * k + 2.0D+00 ) * x
        df = df + r 
        if ( abs ( r ) < abs ( df ) * eps ) then
           exit
        end if
     end do

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

     ad = c1 * df - c2 * dg
     bd = sr3 * ( c1 * df + c2 * dg )

  else

     xe = xa * xq / 1.5D+00
     xr1 = 1.0D+00 / xe
     xar = 1.0D+00 / xq
     xf = sqrt ( xar )
     rp = 0.5641895835477563D+00
     r = 1.0D+00
     do k = 1, 40
        r = r * ( 6.0D+00 * k - 1.0D+00 ) &
             / 216.0D+00 * ( 6.0D+00 * k - 3.0D+00 ) &
             / k * ( 6.0D+00 * k - 5.0D+00 ) / ( 2.0D+00 * k - 1.0D+00 )
        ck(k) = r
        dk(k) = - ( 6.0D+00 * k + 1.0D+00 ) / ( 6.0D+00 * k - 1.0D+00 ) * ck(k)
     end do

     km = int ( 24.5D+00 - xa )

     if ( xa < 6.0D+00 ) then
        km = 14
     end if

     if ( 15.0D+00 < xa ) then
        km = 10
     end if

     if ( 0.0D+00 < x ) then
        sai = 1.0D+00
        sad = 1.0D+00
        r = 1.0D+00
        do k = 1, km
           r = - r * xr1
           sai = sai + ck(k) * r
           sad = sad + dk(k) * r
        end do
        sbi = 1.0D+00
        sbd = 1.0D+00
        r = 1.0D+00
        do k = 1, km
           r = r * xr1
           sbi = sbi + ck(k) * r
           sbd = sbd + dk(k) * r
        end do
        xp1 = exp ( - xe )
        ai = 0.5D+00 * rp * xf * xp1 * sai
        bi = rp * xf / xp1 * sbi
        ad = -0.5D+00 * rp / xf * xp1 * sad
        bd = rp / xf / xp1 * sbd
     else
        xcs = cos ( xe + pi / 4.0D+00 )
        xss = sin ( xe + pi / 4.0D+00 )
        ssa = 1.0D+00
        sda = 1.0D+00
        r = 1.0D+00
        xr2 = 1.0D+00 / ( xe * xe )
        do k = 1, km
           r = - r * xr2
           ssa = ssa + ck(2*k) * r
           sda = sda + dk(2*k) * r
        end do
        ssb = ck(1) * xr1
        sdb = dk(1) * xr1
        r = xr1
        do k = 1, km
           r = - r * xr2
           ssb = ssb + ck(2*k+1) * r
           sdb = sdb + dk(2*k+1) * r
        end do
        ai = rp * xf * ( xss * ssa - xcs * ssb )
        bi = rp * xf * ( xcs * ssa + xss * ssb )
        ad = -rp / xf * ( xcs * sda + xss * sdb )
        bd =  rp / xf * ( xss * sda - xcs * sdb )
     end if

  end if

  return
end subroutine airyb