function [ c, s ] = fresnel ( x )

%*****************************************************************************80
%
%% fresnel() computes Fresnel integrals C(x) and S(x).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 July 2025
%
%  Author:
%
%    Original Fortran77 version by Shanjie Zhang, Jianming Jin.
%    This version by John Burkardt.
%
%  Reference:
%
%    John D Cook,
%    Cornu's spiral,
%    Posted 23 March 2026.
%    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 X: the argument.
%
%  Output:
%
%    real C, S: the values of the Fresnel cosine and sine integrals.
%
% eps = 1.0D-15
% pi = 3.141592653589793
  xa = abs ( x );
  px = pi * xa;
  t = 0.5 * px * xa;
  t2 = t * t;
%
%  x == 0
%
  if ( xa == 0.0 )

    c = 0.0;
    s = 0.0;
%
%  0 < x < 2.5
%
  elseif ( xa < 2.5 )

    r = xa;
    c = r;
    for k = 1 : 50
      r = -0.5 * r * ( 4.0 * k - 3.0 ) / k ...
        / ( 2.0 * k - 1.0 ) / ( 4.0 * k + 1.0 ) * t2;
      c = c + r;
      if ( abs ( r ) < abs ( c ) * eps )
        break;
      end
    end

    s = xa * t / 3.0;
    r = s;
    for k = 1 : 50
      r = - 0.5 * r * ( 4.0 * k - 1.0 ) / k ...
        / ( 2.0 * k + 1.0 ) / ( 4.0 * k + 3.0 ) * t2;
      s = s + r;
      if ( abs ( r ) < abs ( s ) * eps )
        if ( x < 0.0 )
          c = -c;
          s = -s;
        end
        return
      end
    end
%
%  2.5 <= x < 4.5
%
  elseif ( xa < 4.5 )

    m = floor ( 42.0 + 1.75 * t );
    su = 0.0;
    c = 0.0;
    s = 0.0;
    f1 = 0.0;
    f0 = 1.0E-100;

    for k = m : -1 : 0
      f = ( 2.0 * k + 3.0 ) * f0 / t - f1;
      if ( k == floor ( k / 2 ) * 2 )
        c = c + f;
      else
        s = s + f;
      end
      su = su + ( 2.0 * k + 1.0 ) * f * f;
      f1 = f0;
      f0 = f;
    end

    q = sqrt ( su );
    c = c * xa / q;
    s = s * xa / q;
%
%  4.5 <= x
%
  else

    r = 1.0;
    f = 1.0;
    for k = 1 : 20
      r = -0.25 * r * ( 4.0 * k - 1.0 ) * ( 4.0 * k - 3.0 ) / t2;
      f = f + r;
    end
    r = 1.0 / ( px * xa );
    g = r;
    for k = 1 : 12
      r = -0.25 * r * ( 4.0 * k + 1.0 ) * ( 4.0 * k - 1.0 ) / t2;
      g = g + r;
    end

    t0 = t - floor ( t / ( 2.0 * pi ) ) * 2.0 * pi;
    c = 0.5 + ( f * sin ( t0 ) - g * cos ( t0 ) ) / px;
    s = 0.5 - ( f * cos ( t0 ) + g * sin ( t0 ) ) / px;

  end
%
%  Apply symmetry for x < 0.
%
  if ( x < 0.0 )
    c = -c;
    s = -s;
  end

  return
end
