fcszo Subroutine

subroutine fcszo ( kf, nt, zo )

  !*****************************************************************************80
  !
  !! FCSZO computes complex zeros of Fresnel integrals C(x) or S(x).
  !
  !  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:
  !
  !    17 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 ) KF, the function code.
  !    1 for C(z);
  !    2 for S(z)
  !
  !    Input, integer ( kind = 4 ) NT, the total number of zeros desired.
  !
  !    Output, complex ( kind = 8 ) Z0(NT), the zeros.
  !
  implicit none

  integer ( kind = 4 ) nt

  integer ( kind = 4 ) i
  integer ( kind = 4 ) it
  integer ( kind = 4 ) j
  integer ( kind = 4 ) kf
  integer ( kind = 4 ) nr
  real ( kind = 8 ) pi
  real ( kind = 8 ) psq
  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

     if ( kf == 1 ) then
        psq = sqrt ( 4.0D+00 * nr - 1.0D+00 )
     else
        psq = 2.0D+00 * sqrt ( real ( nr, kind = 8 ) )
     end if

     px = psq - log ( pi * psq ) / ( pi * pi * psq ** 3.0D+00 )
     py = log ( pi * psq ) / ( pi * psq )
     z = cmplx ( px, py )

     if ( kf == 2 ) then
        if ( nr == 2 ) then
           z = cmplx ( 2.8334D+00, 0.2443D+00 )
        else if ( nr == 3 ) then
           z = cmplx ( 3.4674D+00, 0.2185D+00 )
        else if ( nr == 4 ) then
           z = cmplx ( 4.0025D+00, 0.2008D+00 )
        end if
     end if

     it = 0

     do

        it = it + 1

        if ( kf == 1 ) then
           call cfc ( z, zf, zd )
        else
           call cfs ( z, zf, zd )
        end if

        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 = cdabs ( z )

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

        if ( 50 < it ) then
           exit
        end if

     end do

     zo(nr) = z

  end do

  return
end subroutine fcszo