fcs Subroutine

subroutine fcs ( x, c, s )

  !*****************************************************************************80
  !
  !! FCS computes Fresnel integrals C(x) and 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, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) C, S, the function values.
  !
  implicit none

  real ( kind = 8 ) c
  real ( kind = 8 ) eps
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) g
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) pi
  real ( kind = 8 ) px
  real ( kind = 8 ) q
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) su
  real ( kind = 8 ) t
  real ( kind = 8 ) t0
  real ( kind = 8 ) t2
  real ( kind = 8 ) x
  real ( kind = 8 ) xa

  eps = 1.0D-15
  pi = 3.141592653589793D+00
  xa = abs ( x )
  px = pi * xa
  t = 0.5D+00 * px * xa
  t2 = t * t

  if ( xa == 0.0D+00 ) then

     c = 0.0D+00
     s = 0.0D+00

  else if ( xa < 2.5D+00 ) then

     r = xa
     c = r
     do k = 1, 50
        r = -0.5D+00 * r * ( 4.0D+00 * k - 3.0D+00 ) / k &
             / ( 2.0D+00 * k - 1.0D+00 ) / ( 4.0D+00 * k + 1.0D+00 ) * t2
        c = c + r
        if ( abs ( r ) < abs ( c ) * eps ) then
           exit
        end if
     end do

     s = xa * t / 3.0D+00
     r = s
     do k = 1, 50
        r = - 0.5D+00 * r * ( 4.0D+00 * k - 1.0D+00 ) / k &
             / ( 2.0D+00 * k + 1.0D+00 ) / ( 4.0D+00 * k + 3.0D+00 ) * t2
        s = s + r
        if ( abs ( r ) < abs ( s ) * eps ) then
           if ( x < 0.0D+00 ) then
              c = -c
              s = -s
           end if
           return
        end if
     end do

  else if ( xa < 4.5D+00 ) then

     m = int ( 42.0D+00 + 1.75D+00 * t )
     su = 0.0D+00
     c = 0.0D+00
     s = 0.0D+00
     f1 = 0.0D+00
     f0 = 1.0D-100

     do k = m, 0, -1
        f = ( 2.0D+00 * k + 3.0D+00 ) * f0 / t - f1
        if ( k == int ( k / 2 ) * 2 ) then
           c = c + f
        else
           s = s + f
        end if
        su = su + ( 2.0D+00 * k + 1.0D+00 ) * f * f
        f1 = f0
        f0 = f
     end do

     q = sqrt ( su )
     c = c * xa / q
     s = s * xa / q

  else

     r = 1.0D+00
     f = 1.0D+00
     do k = 1, 20
        r = -0.25D+00 * r * ( 4.0D+00 * k - 1.0D+00 ) &
             * ( 4.0D+00 * k - 3.0D+00 ) / t2
        f = f + r
     end do
     r = 1.0D+00 / ( px * xa )
     g = r
     do k = 1, 12
        r = -0.25D+00 * r * ( 4.0D+00 * k + 1.0D+00 ) &
             * ( 4.0D+00 * k - 1.0D+00 ) / t2
        g = g + r
     end do

     t0 = t - int ( t / ( 2.0D+00 * pi ) ) * 2.0D+00 * pi
     c = 0.5D+00 + ( f * sin ( t0 ) - g * cos ( t0 ) ) / px
     s = 0.5D+00 - ( f * cos ( t0 ) + g * sin ( t0 ) ) / px

  end if

  if ( x < 0.0D+00 ) then
     c = -c
     s = -s
  end if

  return
end subroutine fcs