cfc Subroutine

subroutine cfc ( z, zf, zd )

  !*****************************************************************************80
  !
  !! CFC computes the complex Fresnel integral C(z) and C'(z).
  !
  !  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:
  !
  !    26 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, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, complex ( kind = 8 ) ZF, ZD, the values of C(z) and C'(z).
  !
  implicit none

  complex ( kind = 8 ) c
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf0
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cg
  complex ( kind = 8 ) cr
  real ( kind = 8 ) eps
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) pi
  real ( kind = 8 ) w0
  real ( kind = 8 ) wa
  real ( kind = 8 ) wa0
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z0
  complex ( kind = 8 ) zd
  complex ( kind = 8 ) zf
  complex ( kind = 8 ) zp
  complex ( kind = 8 ) zp2

  eps = 1.0D-14
  pi = 3.141592653589793D+00
  w0 = abs ( z )
  zp = 0.5D+00 * pi * z * z
  zp2 = zp * zp
  z0 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )

  if ( z .eq. z0 ) then

     c = z0

  else if ( w0 <= 2.5D+00 ) then

     cr = z
     c = cr
     do k = 1, 80
        cr = -0.5D+00 * cr * ( 4.0D+00 * k - 3.0D+00 ) &
             / k / ( 2.0D+00 * k - 1.0D+00 ) &
             / ( 4.0D+00 * k + 1.0D+00 ) * zp2
        c = c + cr
        wa = abs ( c )
        if ( abs ( ( wa - wa0 ) / wa ) < eps .and. 10 < k ) then
           exit
        end if
        wa0 = wa
     end do

  else if ( 2.5D+00 < w0 .and. w0 < 4.5D+00 ) then

     m = 85
     c = z0
     cf1 = z0
     cf0 = cmplx ( 1.0D-30, 0.0D+00, kind = 8 )
     do k = m, 0, -1
        cf = ( 2.0D+00 * k + 3.0D+00 ) * cf0 / zp - cf1
        if ( k .eq. int ( k / 2 ) * 2 ) then
           c = c + cf
        end if
        cf1 = cf0
        cf0 = cf
     end do
     c = sqrt ( 2.0D+00 / ( pi * zp ) ) * sin ( zp ) / cf * c

  else

     cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     cf = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, 20
        cr = - 0.25D+00 * cr * ( 4.0D+00 * k - 1.0D+00 ) &
             * ( 4.0D+00 * k - 3.0D+00 ) / zp2
        cf = cf + cr
     end do
     cr = 1.0D+00 / ( pi * z * z )
     cg = cr
     do k = 1, 12
        cr = - 0.25D+00 * cr * ( 4.0D+00 * k + 1.0D+00 ) &
             * ( 4.0D+00 * k - 1.0D+00 ) / zp2
        cg = cg + cr
     end do
     c = 0.5D+00 + ( cf * sin ( zp ) - cg * cos ( zp ) ) / ( pi * z )

  end if

  zf = c
  zd = cos ( 0.5D+00 * pi * z * z )

  return
end subroutine cfc