cyzo Subroutine

subroutine cyzo ( nt, kf, kc, zo, zv )

  !*****************************************************************************80
  !
  !! CYZO computes zeros of complex Bessel functions Y0(z) and Y1(z) and Y1'(z).
  !
  !  Parameters:
  !
  !    Ths procedure computes the complex zeros of Y0(z), Y1(z) and Y1'(z), 
  !    and their associated values at the zeros using the modified Newton's 
  !    iteration method.
  !
  !  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:
  !
  !    22 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 ) KF, the function choice.
  !    0 for Y0(z) and Y1(z0);
  !    1 for Y1(z) and Y0(z1);
  !    2 for Y1'(z) and Y1(z1').
  !
  !    Input, integer ( kind = 4 ) KC, complex/real choice.
  !    0, for complex roots;
  !    1, for real roots.
  !
  !    Output, real ( kind = 8 ) ZO(NT), ZV(NT), the zeros of Y0(z) or Y1(z) 
  !    or Y1'(z), and the value of Y0'(z) or Y1'(z) or Y1(z) at the L-th zero.
  !
  implicit none

  integer ( kind = 4 ) nt

  real ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ) it
  integer ( kind = 4 ) j
  integer ( kind = 4 ) kc
  integer ( kind = 4 ) kf
  integer ( kind = 4 ) nr
  real ( kind = 8 ) w
  real ( kind = 8 ) w0
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  complex ( kind = 8 ) z
  complex ( kind = 8 ) zd
  complex ( kind = 8 ) zero
  complex ( kind = 8 ) zf
  complex ( kind = 8 ) zfd
  complex ( kind = 8 ) zgd
  complex ( kind = 8 ) zo(nt)
  complex ( kind = 8 ) zp
  complex ( kind = 8 ) zq
  complex ( kind = 8 ) zv(nt)
  complex ( kind = 8 ) zw

  if ( kc == 0 ) then
     x = -2.4D+00
     y = 0.54D+00
     h = 3.14D+00
  else if ( kc == 1 ) then
     x = 0.89D+00
     y = 0.0D+00
     h = -3.14D+00
  end if

  if ( kf == 1 ) then
     x = -0.503D+00
  else if ( kf == 2 ) then
     x = 0.577D+00
  end if

  zero = cmplx ( x, y, kind = 8 )

  do nr = 1, nt

     if ( nr == 1 ) then
        z = zero
     else
        z = zo(nr-1) - h
     end if

     it = 0

     do

        it = it + 1
        call cy01 ( kf, z, zf, zd )

        zp = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do i = 1, nr - 1
           zp = zp * ( z - zo(i) )
        end do

        zfd = zf / zp

        zq = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        do i = 1, nr - 1
           zw = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
           do j = 1, nr - 1
              if ( j /= i ) then
                 zw = zw * ( z - zo(j) )
              end if
           end do
           zq = zq + zw
        end do

        zgd = ( zd - zq * zfd ) / zp
        z = z - zfd / zgd
        w0 = w
        w = abs ( z )

        if ( 50 < it .or. abs ( ( w - w0 ) / w ) <= 1.0D-12 ) then
           exit
        end if

     end do

     zo(nr) = z

  end do

  do i = 1, nt
     z = zo(i)
     if ( kf == 0 .or. kf == 2 ) then
        call cy01 ( 1, z, zf, zd )
        zv(i) = zf
     else if ( kf == 1 ) then
        call cy01 ( 0, z, zf, zd )
        zv(i) = zf
     end if
  end do

  return
end subroutine cyzo