cerzo Subroutine

subroutine cerzo ( nt, zo )

  !*****************************************************************************80
  !
  !! CERZO evaluates the complex zeros of the error function.
  !
  !  Discussion:
  !
  !    The modified Newton method is used.
  !
  !  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.
  !
  !    Output, complex ( kind = 8 ) ZO(NT), the zeros.
  !
  implicit none

  integer ( kind = 4 ) nt

  integer ( kind = 4 ) i
  integer ( kind = 4 ) it
  integer ( kind = 4 ) j
  integer ( kind = 4 ) nr
  real ( kind = 8 ) pi
  real ( kind = 8 ) pu
  real ( kind = 8 ) pv
  real ( kind = 8 ) px
  real ( kind = 8 ) py
  real ( kind = 8 ) w
  real ( kind = 8 ) w0
  complex ( kind = 8 ) z
  complex ( kind = 8 ) zd
  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 ) zw

  pi = 3.141592653589793D+00

  do nr = 1, nt

     pu = sqrt ( pi * ( 4.0D+00 * nr - 0.5D+00 ) )
     pv = pi * sqrt ( 2.0D+00 * nr - 0.25D+00 )
     px = 0.5D+00 * pu - 0.5D+00 * log ( pv ) / pu
     py = 0.5D+00 * pu + 0.5D+00 * log ( pv ) / pu
     z = cmplx ( px, py, kind = 8 )
     it = 0

     do

        it = it + 1
        call cerf ( 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-11 ) then
           exit
        end if

     end do

     zo(nr) = z

  end do

  return
end subroutine cerzo