# include <stdlib.h>
# include <math.h>

# include "fresnel.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

double fresnel_cos ( double x )

/******************************************************************************/
/*
  Purpose:

    fresnel_cos() evaluates the Fresnel cosine integral C(x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 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;
}
/******************************************************************************/

double fresnel_sin ( double x )

/******************************************************************************/
/*
  Purpose:

    fresnel_sin() evaluates the Fresnel sine integral S(x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 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;
}