klvnzo Subroutine

subroutine klvnzo ( nt, kd, zo )

  !*****************************************************************************80
  !
  !! KLVNZO computes zeros of the Kelvin functions.
  !
  !  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:
  !
  !    15 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, integer ( kind = 4 ) NT, the number of zeros.
  !
  !    Input, integer ( kind = 4 ) KD, the function code.
  !    1 for ber x, 
  !    2 for bei x,
  !    3 for ker x, 
  !    4 for kei x,
  !    5 for ber' x, 
  !    6 for bei' x,
  !    7 for ker' x, 
  !    8 for kei' x.
  !
  !    Output, real ( kind = 8 ) ZO(NT), the zeros of the given Kelvin function.
  !
  implicit none

  integer ( kind = 4 ) nt

  real ( kind = 8 ) bei
  real ( kind = 8 ) ber
  real ( kind = 8 ) ddi
  real ( kind = 8 ) ddr
  real ( kind = 8 ) dei
  real ( kind = 8 ) der
  real ( kind = 8 ) gdi
  real ( kind = 8 ) gdr
  real ( kind = 8 ) gei
  real ( kind = 8 ) ger
  real ( kind = 8 ) hei
  real ( kind = 8 ) her
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  real ( kind = 8 ) rt
  real ( kind = 8 ) rt0(8)
  real ( kind = 8 ) zo(nt)

  rt0(1) = 2.84891D+00
  rt0(2) = 5.02622D+00
  rt0(3) = 1.71854D+00
  rt0(4) = 3.91467D+00
  rt0(5) = 6.03871D+00
  rt0(6) = 3.77268D+00
  rt0(7) = 2.66584D+00
  rt0(8) = 4.93181D+00

  rt = rt0(kd)

  do m = 1, nt

     do

        call klvna ( rt, ber, bei, ger, gei, der, dei, her, hei )

        if ( kd == 1 ) then
           rt = rt - ber / der
        else if ( kd == 2 ) then
           rt = rt - bei / dei
        else if ( kd == 3 ) then
           rt = rt - ger / her
        else if ( kd == 4 ) then
           rt = rt - gei / hei
        else if ( kd == 5 ) then
           ddr = - bei - der / rt
           rt = rt - der / ddr
        else if ( kd == 6 ) then
           ddi = ber - dei / rt
           rt = rt - dei / ddi
        else if ( kd == 7 ) then
           gdr = - gei - her / rt
           rt = rt - her / gdr
        else
           gdi = ger - hei / rt
           rt = rt - hei / gdi
        end if

        if ( abs ( rt - rt0(kd) ) <= 5.0D-10 ) then
           exit
        end if

        rt0(kd) = rt

     end do

     zo(m) = rt
     rt = rt + 4.44D+00

  end do

  return
end subroutine klvnzo