# include <cstdlib>
# include <cmath>

using namespace std;

# include "fresnel.hpp"

//****************************************************************************80

void fresnel ( double x, double &c, double &s )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    11 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 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:
//
//    double X, the argument.
//
//  Output:
//
//    double C, S, the function values.
//
{
  double eps;
  double f;
  double f0;
  double f1;
  double g;
  int k;
  int m;
  double pi;
  double px;
  double q;
  double r;
  double su;
  double t;
  double t0;
  double t2;
  double xa;

  eps = 1.0E-15;
  pi = 3.141592653589793;
  xa = fabs ( 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
//
  else if ( xa < 2.5 )
  {
    r = xa;
     c = r;
    for ( k = 1; k <= 50; k++ )
    {
      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 ( fabs ( r ) < fabs (  c ) * eps )
      {
        break;
      }
    }

     s = xa * t / 3.0;
    r =  s;
    for ( k = 1; k <= 50; k++ )
    {
      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 ( fabs ( r ) < fabs (  s ) * eps )
      {
        if ( x < 0.0 )
        {
           c = -  c;
           s = -  s;
        }
        return;
      }
    }
  }
//
//  2.5 <= x < 4.5
//
  else if ( xa < 4.5 )
  {
    m = ( int ) ( 42.0 + 1.75 * t );
    su = 0.0;
     c = 0.0;
     s = 0.0;
    f1 = 0.0;
    f0 = 1.0E-100;

    for ( k = m; 0 <= k; k-- )
    {
      f = ( 2.0 * k + 3.0 ) * f0 / t - f1;
      if ( k == ( int ) ( k / 2 ) * 2 )
      {
         c =  c + f;
      }
      else
      {
         s =  s + f;
      }
      su = su + ( 2.0 * k + 1.0 ) * f * f;
      f1 = f0;
      f0 = f;
    }

    q = sqrt ( su );
     c =  c * xa / q;
     s =  s * xa / q;
  }
//
//  4.5 <= x
//
  else
  {
    r = 1.0;
    f = 1.0;
    for ( k = 1; k <= 20; k++ )
    {
      r = -0.25 * r * ( 4.0 * k - 1.0 ) * ( 4.0 * k - 3.0 ) / t2;
      f = f + r;
    }
    r = 1.0 / ( px * xa );
    g = r;
    for ( k = 1; k <= 12; k++ )
    {
      r = -0.25 * r * ( 4.0 * k + 1.0 ) * ( 4.0 * k - 1.0 ) / t2;
      g = g + r;
    }

    t0 = t - ( int ) ( 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;
  }
//
//  Apply symmetry for x < 0.
//
  if ( x < 0.0 )
  {
     c = -  c;
     s = -  s;
  }

  return;
}
//****************************************************************************80

double fresnel_cos ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    fresnel_cos() evaluates the Fresnel cosine integral C(x).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 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.
//
{
  double c;
  double s;

  fresnel ( x, c, s );

  return c;
}
//****************************************************************************80

double fresnel_sin ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    fresnel_sin() evaluates the Fresnel sine integral S(x).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 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.
//
{
  double c;
  double s;

  fresnel ( x, c, s );

  return s;
}