airyzo Subroutine

subroutine airyzo ( nt, kf, xa, xb, xc, xd )

  !*****************************************************************************80
  !
  !! AIRYZO computes the first NT zeros of Ai(x) and Ai'(x).
  !
  !   Discussion:
  !
  !    Compute the first NT zeros of Airy functions Ai(x) and Ai'(x), 
  !    a and a', and the associated values of Ai(a') and Ai'(a); and 
  !    the first NT zeros of Airy functions Bi(x) and Bi'(x), b and
  !    b', and the associated values of Bi(b') and Bi'(b).
  !
  !  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:
  !
  !    14 March 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, integer ( kind = 4 ) NT, the number of zeros.
  !
  !    Input, integer ( kind = 4 ) KF, the function code.
  !    1 for Ai(x) and Ai'(x);
  !    2 for Bi(x) and Bi'(x).
  !
  !    Output, real ( kind = 8 ) XA(m), a, the m-th zero of Ai(x) or
  !    b, the m-th zero of Bi(x).
  !
  !    Output, real ( kind = 8 ) XB(m), a', the m-th zero of Ai'(x) or
  !    b', the m-th zero of Bi'(x).
  !
  !    Output, real ( kind = 8 ) XC(m), Ai(a') or Bi(b').
  !
  !    Output, real ( kind = 8 ) XD(m), Ai'(a) or Bi'(b)
  !
  implicit none

  integer ( kind = 4 ) nt

  real ( kind = 8 ) ad
  real ( kind = 8 ) ai
  real ( kind = 8 ) bd
  real ( kind = 8 ) bi
  integer ( kind = 4 ) i
  integer ( kind = 4 ) kf
  real ( kind = 8 ) pi
  real ( kind = 8 ) rt
  real ( kind = 8 ) rt0
  real ( kind = 8 ) u
  real ( kind = 8 ) u1
  real ( kind = 8 ) x
  real ( kind = 8 ) xa(nt)
  real ( kind = 8 ) xb(nt)
  real ( kind = 8 ) xc(nt)
  real ( kind = 8 ) xd(nt)

  pi = 3.141592653589793D+00

  do i = 1, nt

     if (kf == 1) then
        u = 3.0D+00 * pi * ( 4.0D+00 * i - 1 ) / 8.0D+00
        u1 = 1.0D+00 / ( u * u )
        rt0 = - ( u * u ) ** ( 1.0 / 3.0 ) &
             * (((( -15.5902D+00 * u1 + 0.929844D+00 ) * u1 &
             - 0.138889D+00 ) * u1 + 0.10416667D+00 ) * u1 + 1.0D+00 )
     else if ( kf == 2 ) then
        if ( i == 1 ) then
           rt0 = -1.17371D+00
        else
           u = 3.0D+00 * pi * ( 4.0D+00 * i - 3.0D+00 ) / 8.0D+00
           u1 = 1.0D+00 / ( u * u )
           rt0 = - ( u * u ) ** ( 1.0D+00 / 3.0D+00 ) &
                * (((( -15.5902D+00 * u1 + 0.929844D+00 ) * u1 &
                - 0.138889D+00 ) * u1 + 0.10416667D+00 ) * u1 + 1.0D+00 )
        end if
     end if

     do

        x = rt0
        call airyb ( x, ai, bi, ad, bd )

        if ( kf == 1 ) then
           rt = rt0 - ai / ad
        else
           rt = rt0 - bi / bd
        end if

        if ( abs ( ( rt - rt0 ) / rt ) <= 1.0D-09 ) then
           exit
        end if
        rt0 = rt

     end do

     xa(i) = rt
     if ( kf == 1 ) then
        xd(i) = ad
     else
        xd(i) = bd
     end if

  end do

  do i = 1, nt

     if ( kf == 1 ) then
        if ( i == 1 ) then
           rt0 = -1.01879D+00
        else
           u = 3.0D+00 * pi * ( 4.0D+00 * i - 3.0D+00 ) / 8.0D+00
           u1 = 1.0D+00 / ( u * u )
           rt0 = - ( u * u ) ** ( 1.0D+00 / 3.0D+00 ) &
                * (((( 15.0168D+00 * u1 - 0.873954D+00 ) &
                * u1 + 0.121528D+00 ) * u1 - 0.145833D+00 ) * u1 + 1.0D+00 )
        end if
     else if ( kf == 2 ) then
        if ( i == 1 ) then
           rt0 = -2.29444D+00
        else
           u = 3.0D+00 * pi * ( 4.0D+00 * i - 1.0D+00 ) / 8.0D+00
           u1 = 1.0D+00 / ( u * u )
           rt0 = - ( u * u ) ** ( 1.0D+00 / 3.0D+00 ) &
                * (((( 15.0168D+00 * u1 - 0.873954D+00 ) &
                * u1 + 0.121528D+00 ) * u1 - 0.145833D+00 ) * u1 + 1.0D+00 )
        end if
     end if

     do

        x = rt0
        call airyb ( x, ai, bi, ad, bd )

        if ( kf == 1 ) then
           rt = rt0 - ad / ( ai * x )
        else
           rt = rt0 - bd / ( bi * x )
        end if

        if ( abs ( ( rt - rt0 ) / rt ) <= 1.0D-09 ) then
           exit
        end if

        rt0 = rt

     end do

     xb(i) = rt
     if ( kf == 1 ) then
        xc(i) = ai
     else
        xc(i) = bi
     end if

  end do

  return
end subroutine airyzo