/****************************************************************************
 *              fresnel.c -
 *  Calculation of Fresnel integrals by expansion to Chebyshev series
 *  Expansions are taken from the book
 *  Y.L. Luke. Mathematical functions and their approximations.
 *  Moscow, "Mir", 1980. PP. 145-149 (Russian edition)
 ****************************************************************************
 */
/*
  Modified for Ruby bindings
          2006/Dec/24 Y. TSUNESADA
*/
#include <math.h>
#include "include/rb_gsl.h"
#include "include/rb_gsl_sf.h"

static const double sqrt_pi_2   = 1.2533141373155002512078826424; /* sqrt(pi/2) */
static const double sqrt_2_pi   = 0.7978845608028653558798921199; /* sqrt(2/pi) */
static const double _1_sqrt_2pi = 0.3989422804014326779399460599; /* 1/sqrt(2*pi) */
static const double pi_2        = 1.5707963267948966192313216916; /* pi/2 */

static double f_data_a[18] =
{
  0.76435138664186000189,
  -0.43135547547660179313,
  0.43288199979726653054,
  -0.26973310338387111029,
  0.08416045320876935378,
  -0.01546524484461381958,
  0.00187855423439822018,
  -0.00016264977618887547,
  0.00001057397656383260,
  -0.00000053609339889243,
  0.00000002181658454933,
  -0.00000000072901621186,
  0.00000000002037332546,
  -0.00000000000048344033,
  0.00000000000000986533,
  -0.00000000000000017502,
  0.00000000000000000272,
  -0.00000000000000000004
};

static double f_data_b[17] =
{
  0.63041404314570539241,
  -0.42344511405705333544,
  0.37617172643343656625,
  -0.16249489154509567415,
  0.03822255778633008694,
  -0.00564563477132190899,
  0.00057454951976897367,
  -0.00004287071532102004,
  0.00000245120749923299,
  -0.00000011098841840868,
  0.00000000408249731696,
  -0.00000000012449830219,
  0.00000000000320048425,
  -0.00000000000007032416,
  0.00000000000000133638,
  -0.00000000000000002219,
  0.00000000000000000032
};

static double fresnel_cos_0_8(double x)
{
  double x_8 = x/8.0;
  double xx = 2.0*x_8*x_8 - 1.0;

  double t0 = 1.0;
  double t1 = xx;
  double sumC = f_data_a[0] + f_data_a[1]*t1;
  double t2;
  int n;
  for (n = 2; n < 18; n++)
  {
    t2 = 2.0*xx*t1 - t0;
    sumC += f_data_a[n]*t2;
    t0 = t1; t1 = t2;
  }
  return _1_sqrt_2pi*sqrt(x)*sumC;
}

static double fresnel_sin_0_8(double x)
{
  double x_8 = x/8.0;
  double xx = 2.0*x_8*x_8 - 1.0;
  double t0 = 1.;
  double t1 = xx;
  double ot1 = x_8;
  double ot2 = 2.0*x_8*t1 - ot1;
  double sumS = f_data_b[0]*ot1 + f_data_b[1]*ot2;
  int n;
  double t2;
  for (n = 2; n < 17; n++)
  {
    t2 = 2.0*xx*t1 - t0;
    ot1 = ot2;
    ot2 = 2.0*x_8*t2 - ot1;
    sumS += f_data_b[n]*ot2;
    t0 = t1; t1 = t2;
  }
  return _1_sqrt_2pi*sqrt(x)*sumS;
}

static double f_data_e[41] =
{
  0.97462779093296822410,
  -0.02424701873969321371,
  0.00103400906842977317,
  -0.00008052450246908016,
  0.00000905962481966582,
  -0.00000131016996757743,
  0.00000022770820391497,
  -0.00000004558623552026,
  0.00000001021567537083,
  -0.00000000251114508133,
  0.00000000066704761275,
  -0.00000000018931512852,
  0.00000000005689898935,
  -0.00000000001798219359,
  0.00000000000594162963,
  -0.00000000000204285065,
  0.00000000000072797580,
  -0.00000000000026797428,
  0.00000000000010160694,
  -0.00000000000003958559,
  0.00000000000001581262,
  -0.00000000000000646411,
  0.00000000000000269981,
  -0.00000000000000115038,
  0.00000000000000049942,
  -0.00000000000000022064,
  0.00000000000000009910,
  -0.00000000000000004520,
  0.00000000000000002092,
  -0.00000000000000000982,
  0.00000000000000000467,
  -0.00000000000000000225,
  0.00000000000000000110,
  -0.00000000000000000054,
  0.00000000000000000027,
  -0.00000000000000000014,
  0.00000000000000000007,
  -0.00000000000000000004,
  0.00000000000000000002,
  -0.00000000000000000001,
  0.00000000000000000001
};

static double f_data_f[35] =
{
  0.99461545179407928910,
  -0.00524276766084297210,
  0.00013325864229883909,
  -0.00000770856452642713,
  0.00000070848077032045,
  -0.00000008812517411602,
  0.00000001359784717148,
  -0.00000000246858295747,
  0.00000000050925789921,
  -0.00000000011653400634,
  0.00000000002906578309,
  -0.00000000000779847361,
  0.00000000000222802542,
  -0.00000000000067239338,
  0.00000000000021296411,
  -0.00000000000007041482,
  0.00000000000002419805,
  -0.00000000000000861080,
  0.00000000000000316287,
  -0.00000000000000119596,
  0.00000000000000046444,
  -0.00000000000000018485,
  0.00000000000000007527,
  -0.00000000000000003131,
  0.00000000000000001328,
  -0.00000000000000000574,
  0.00000000000000000252,
  -0.00000000000000000113,
  0.00000000000000000051,
  -0.00000000000000000024,
  0.00000000000000000011,
  -0.00000000000000000005,
  0.00000000000000000002,
  -0.00000000000000000001,
  0.00000000000000000001
};

static double fresnel_cos_8_inf(double x)
{
  double xx = 128.0/(x*x) - 1.0;  /* 2.0*(8/x)^2 - 1 */
  double t0 = 1.0;
  double t1 = xx;
  double sumP = f_data_e[0] + f_data_e[1]*t1;
  double sumQ = f_data_f[0] + f_data_f[1]*t1;
  double t2;
  int n;
  for(n = 2; n < 35; n++)
  {
    t2 = 2.0*xx*t1 - t0;
    sumP += f_data_e[n]*t2; /*  sumP += f_data_e[n]*ChebyshevT(n,xx) */
    sumQ += f_data_f[n]*t2; /*  sumQ += f_data_f[n]*ChebyshevT(n,xx) */
    t0 = t1; t1 = t2;
  }
  for(n = 35; n < 41; n++)
  {
    t2 = 2.0*xx*t1 - t0;
    sumP += f_data_e[n]*t2; /*  sumP += f_data_e[n]*ChebyshevT(n,xx) */
    t0 = t1; t1 = t2;
  }
  return 0.5 - _1_sqrt_2pi*(0.5*sumP*cos(x)/x - sumQ*sin(x))/sqrt(x);
}

static double fresnel_sin_8_inf(double x)
{
  double xx = 128.0/(x*x) - 1.0;  /* 2.0*(8/x)^2 - 1 */
  double t0 = 1.0;
  double t1 = xx;
  double sumP = f_data_e[0] + f_data_e[1]*t1;
  double sumQ = f_data_f[0] + f_data_f[1]*t1;
  double t2;
  int n;
  for(n = 2; n < 35; n++)
  {
    t2 = 2.0*xx*t1 - t0;
    sumP += f_data_e[n]*t2; /*  sumP += f_data_e[n]*ChebyshevT(n,xx) */
    sumQ += f_data_f[n]*t2; /*  sumQ += f_data_f[n]*ChebyshevT(n,xx) */
    t0 = t1; t1 = t2;
  }
  for(n = 35; n < 41; n++)
  {
    t2 = 2.0*xx*t1 - t0;
    sumP += f_data_e[n]*t2; /*  sumQ += f_data_f[n]*ChebyshevT(n,xx) */
    t0 = t1; t1 = t2;
  }
  return 0.5 - _1_sqrt_2pi*(0.5*sumP*sin(x)/x + sumQ*cos(x))/sqrt(x);
}


double fresnel_c(double x)
{
  double xx = x*x*pi_2;
  double ret_val;
  if(xx<=8.0)
    ret_val = fresnel_cos_0_8(xx);
  else
    ret_val = fresnel_cos_8_inf(xx);
  return (x<0.0) ? -ret_val : ret_val;
}

double fresnel_s(double x)
{
  double xx = x*x*pi_2;
  double ret_val;
  if(xx<=8.0)
    ret_val = fresnel_sin_0_8(xx);
  else
    ret_val = fresnel_sin_8_inf(xx);
  return (x<0.0) ? -ret_val : ret_val;
}

double fresnel_c1(double x)
{
  return fresnel_c(x*sqrt_2_pi);
}

double fresnel_s1(double x)
{
  return fresnel_s(x*sqrt_2_pi);
}

static VALUE rb_fresnel_c(VALUE obj, VALUE x)
{
  return rb_gsl_sf_eval1(fresnel_c, x);
}
static VALUE rb_fresnel_s(VALUE obj, VALUE x)
{
  return rb_gsl_sf_eval1(fresnel_s, x);
}
static VALUE rb_fresnel_c1(VALUE obj, VALUE x)
{
  return rb_gsl_sf_eval1(fresnel_c1, x);
}
static VALUE rb_fresnel_s1(VALUE obj, VALUE x)
{
  return rb_gsl_sf_eval1(fresnel_s1, x);
}
void Init_fresnel(VALUE module)
{
  VALUE mfresnel;
  mfresnel = rb_define_module_under(module, "Fresnel");
  rb_define_module_function(module, "fresnel_c", rb_fresnel_c, 1);
  rb_define_module_function(module, "fresnel_s", rb_fresnel_s, 1);
  rb_define_module_function(module, "fresnel_c1", rb_fresnel_c1, 1);
  rb_define_module_function(module, "fresnel_s1", rb_fresnel_s1, 1);
  rb_define_module_function(mfresnel, "c", rb_fresnel_c, 1);
  rb_define_module_function(mfresnel, "s", rb_fresnel_s, 1);
  rb_define_module_function(mfresnel, "c1", rb_fresnel_c1, 1);
  rb_define_module_function(mfresnel, "s1", rb_fresnel_s1, 1);
}