cfs Subroutine

subroutine cfs ( z, zf, zd )

  !*****************************************************************************80
  !
  !! CFS computes the complex Fresnel integral S(z) and S'(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:
  !
  !    24 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 S(z) and S'(z).
  !
  implicit none

  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
  complex ( kind = 8 ) s
  real ( kind = 8 ) w0
  real ( kind = 8 ) wb
  real ( kind = 8 ) wb0
  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 == z0 ) then

     s = z0

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

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

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

     m = 85
     s = 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 /= int ( k / 2 ) * 2 ) then
           s = s + cf
        end if
        cf1 = cf0
        cf0 = cf
     end do
     s = sqrt ( 2.0D+00 / ( pi * zp ) ) * sin ( zp ) / cf * s

  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
     s = 0.5D+00 - ( cf * cos ( zp ) + cg * sin ( zp ) ) / ( pi * z )

  end if

  zf = s
  zd = sin ( 0.5D+00 * pi * z * z )

  return
end subroutine cfs