jyzo Subroutine

subroutine jyzo ( n, nt, rj0, rj1, ry0, ry1 )

  !*****************************************************************************80
  !
  !! JYZO computes the zeros of Bessel functions Jn(x), Yn(x) and 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:
  !
  !    28 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 ) N, the order of the Bessel functions.
  !
  !    Input, integer ( kind = 4 ) NT, the number of zeros.
  !
  !    Output, real ( kind = 8 ) RJ0(NT), RJ1(NT), RY0(NT), RY1(NT), the zeros 
  !    of Jn(x), Jn'(x), Yn(x), Yn'(x).
  !
  implicit none

  integer ( kind = 4 ) nt

  real ( kind = 8 ) bjn
  real ( kind = 8 ) byn
  real ( kind = 8 ) djn
  real ( kind = 8 ) dyn
  real ( kind = 8 ) fjn
  real ( kind = 8 ) fyn
  integer ( kind = 4 ) l
  integer ( kind = 4 ) n
  real ( kind = 8 ) n_r8
  real ( kind = 8 ) rj0(nt)
  real ( kind = 8 ) rj1(nt)
  real ( kind = 8 ) ry0(nt)
  real ( kind = 8 ) ry1(nt)
  real ( kind = 8 ) x
  real ( kind = 8 ) x0

  n_r8 = real ( n, kind = 8 )

  if ( n <= 20 ) then
     x = 2.82141D+00 + 1.15859D+00 * n_r8 
  else
     x = n + 1.85576D+00 * n_r8 ** 0.33333D+00 &
          + 1.03315D+00 / n_r8 ** 0.33333D+00
  end if

  l = 0

  do

     x0 = x
     call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
     x = x - bjn / djn

     if ( 1.0D-09 < abs ( x - x0 ) ) then
        cycle
     end if

     l = l + 1
     rj0(l) = x
     x = x + 3.1416D+00 + ( 0.0972D+00 + 0.0679D+00 * n_r8 &
          - 0.000354D+00 * n_r8 ** 2 ) / l

     if ( nt <= l ) then
        exit
     end if

  end do

  if ( n <= 20 ) then
     x = 0.961587D+00 + 1.07703D+00 * n_r8 
  else
     x = n_r8 + 0.80861D+00 * n_r8 ** 0.33333D+00 &
          + 0.07249D+00 / n_r8 ** 0.33333D+00
  end if

  if ( n == 0 ) then
     x = 3.8317D+00
  end if

  l = 0

  do

     x0 = x
     call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
     x = x - djn / fjn
     if ( 1.0D-09 < abs ( x - x0 ) ) then
        cycle
     end if
     l = l + 1
     rj1(l) = x
     x = x + 3.1416D+00 + ( 0.4955D+00 + 0.0915D+00 * n_r8 &
          - 0.000435D+00 * n_r8 ** 2 ) / l

     if ( nt <= l ) then
        exit
     end if

  end do

  if ( n <= 20 ) then
     x = 1.19477D+00 + 1.08933D+00 * n_r8 
  else
     x = n_r8 + 0.93158D+00 * n_r8 ** 0.33333D+00 &
          + 0.26035D+00 / n_r8 ** 0.33333D+00
  end if

  l = 0

  do

     x0 = x
     call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
     x = x - byn / dyn

     if ( 1.0D-09 < abs ( x - x0 ) ) then
        cycle
     end if

     l = l + 1
     ry0(l) = x 
     x = x + 3.1416D+00 + ( 0.312D+00 + 0.0852D+00 * n_r8 &
          - 0.000403D+00 * n_r8 ** 2 ) / l

     if ( nt <= l ) then
        exit
     end if

  end do

  if ( n <= 20 ) then
     x = 2.67257D+00 + 1.16099D+00 * n_r8 
  else
     x = n_r8 + 1.8211D+00 * n_r8 ** 0.33333D+00 &
          + 0.94001D+00 / n_r8 ** 0.33333D+00
  end if

  l = 0

  do

     x0 = x
     call jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )
     x = x - dyn / fyn

     if ( 1.0D-09 < abs ( x - x0 ) ) then
        cycle
     end if

     l = l + 1
     ry1(l) = x
     x = x + 3.1416D+00 + ( 0.197D+00 + 0.0643D+00 * n_r8 &
          -0.000286D+00 * n_r8 ** 2 ) / l 

     if ( nt <= l ) then
        exit
     end if

  end do

  return
end subroutine jyzo