subroutine fresnel ( x, c, s )

!*****************************************************************************80
!
!! fresnel() 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:
!
!    Original Fortran77 version by Shanjie Zhang, Jianming Jin.
!    This version by John Burkardt.
!
!  Reference:
!
!    John D Cook,
!    Cornu's spiral,
!    Posted 23 March 2016.
!    https://www.johndcook.com/blog/2016/03/23/cornus-spiral/
!
!    Shanjie Zhang, Jianming Jin,
!    Computation of Special Functions,
!    Wiley, 1996,
!    ISBN: 0-471-11963-6,
!    LC: QA351.C45.
!
!  Input:
!
!    real ( kind = rk8 ) X, the argument.
!
!  Output:
!
!    real ( kind = rk8 ) C, S, the function values.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

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

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

    c = 0.0D+00
    s = 0.0D+00
!
!  0 < x < 2.5
!
  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
!
!  2.5 <= x < 4.5
!
  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
!
!  4.5 <= x
!
  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
!
!  Apply symmetry for x < 0.
!
  if ( x < 0.0D+00 ) then
    c = -c
    s = -s
  end if

  return
end
function fresnel_cos ( x )

!*****************************************************************************80
!
!! fresnel_cos() evaluates the Fresnel cosine integral C(x).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 July 2025
!
!  Author:
!
!    John Burkardt.
!
!  Reference:
!
!    Shanjie Zhang, Jianming Jin,
!    Computation of Special Functions,
!    Wiley, 1996,
!    ISBN: 0-471-11963-6,
!    LC: QA351.C45.
!
!  Input:
!
!    real X: the argument.
!
!  Output:
!
!    real C: the Fresnel cosine integral value at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) c
  real ( kind = rk8 ) fresnel_cos
  real ( kind = rk8 ) s
  real ( kind = rk8 ) x

  call fresnel ( x, c, s )

  fresnel_cos = c

  return
end
function fresnel_sin ( x )

!*****************************************************************************80
!
!! fresnel_sin() evaluates the Fresnel sine integral S(x).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 July 2025
!
!  Author:
!
!    John Burkardt.
!
!  Reference:
!
!    Shanjie Zhang, Jianming Jin,
!    Computation of Special Functions,
!    Wiley, 1996,
!    ISBN: 0-471-11963-6,
!    LC: QA351.C45.
!
!  Input:
!
!    real X: the argument.
!
!  Output:
!
!    real S: the Fresnel sine integral value at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) c
  real ( kind = rk8 ) fresnel_sin
  real ( kind = rk8 ) s
  real ( kind = rk8 ) x

  call fresnel ( x, c, s )

  fresnel_sin = s

  return
end