airya Subroutine

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

  !*****************************************************************************80
  !
  !! AIRYA computes Airy functions and their derivatives.
  !
  !  Licensing:
  !
  !    The original FORTRAN77 version of this routine is copyrighted by 
  !    Shanjie Zhang and Jianming Jin.  However, they give permission to 
  !    incorporate this routine into a user program that the copyright 
  !    is acknowledged.
  !
  !  Modified:
  !
  !    30 June 2012
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 argument of the Airy function.
  !
  !    Output, real ( kind = 8 ) AI, BI, AD, BD, the values of Ai(x), Bi(x),
  !    Ai'(x), 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 ) pir
  real ( kind = 8 ) sr3
  real ( kind = 8 ) vi1
  real ( kind = 8 ) vi2
  real ( kind = 8 ) vj1
  real ( kind = 8 ) vj2
  real ( kind = 8 ) vk1
  real ( kind = 8 ) vk2
  real ( kind = 8 ) vy1
  real ( kind = 8 ) vy2
  real ( kind = 8 ) x
  real ( kind = 8 ) xa
  real ( kind = 8 ) xq
  real ( kind = 8 ) z

  xa = abs ( x )
  pir = 0.318309886183891D+00
  c1 = 0.355028053887817D+00
  c2 = 0.258819403792807D+00
  sr3 = 1.732050807568877D+00
  z = xa ** 1.5D+00 / 1.5D+00
  xq = sqrt ( xa )

  call ajyik ( z, vj1, vj2, vy1, vy2, vi1, vi2, vk1, vk2 )

  if ( x == 0.0D+00 ) then
     ai = c1
     bi = sr3 * c1
     ad = - c2
     bd = sr3 * c2
  else if ( 0.0D+00 < x ) then
     ai = pir * xq / sr3 * vk1
     bi = xq * ( pir * vk1 + 2.0D+00 / sr3 * vi1 )
     ad = - xa / sr3 * pir * vk2
     bd = xa * ( pir * vk2 + 2.0D+00 / sr3 * vi2 )
  else
     ai = 0.5D+00 * xq * ( vj1 - vy1 / sr3 )
     bi = - 0.5D+00 * xq * ( vj1 / sr3 + vy1 )
     ad = 0.5D+00 * xa * ( vj2 + vy2 / sr3 )
     bd = 0.5D+00 * xa * ( vj2 / sr3 - vy2 )
  end if

  return
end subroutine airya