/*************************************************************
*   Function sinintegral(x) to integrate the real function   *
*   FUNC(t) = sin(t)/t from t=0 to t=x by Simpson's method   *
* ---------------------------------------------------------- *
*   REFERENCE:  (for Simpson's method)                       * 
*              "Mathematiques en Turbo-Pascal (Part 1) By    *
*               Marc Ducamp and Alain Reverchon, Eyrolles,   *
*               Paris, 1987" [BIBLI 03].                     *
* ---------------------------------------------------------- *
*                                                            *
*                                C++ Version By J-P Moreau.  *
*                                    (www.jpmoreau.fr)       *
*************************************************************/
#include <stdio.h>
#include <math.h>

#define  gamma   0.57721566490153286   // Euler's constant

// Given function to integrate (2 kinds)
double FUNC(int kind, double t) {
  switch(kind) {
    case 1:
      if (fabs(t)<1e-10)
        return 1.0;
      else
        return sin(t)/t;        // for sinintegral
	  break;
    case 2:
      if (fabs(t)<1e-10)
        return 0.0;
       else
         return (cos(t)-1.0)/t; // for cosintegral
  }
  return 0;
}

/******************************************************
* Integral of a function FUNC(X) by Simpson's method  *
* --------------------------------------------------- *
* INPUTS:                                             *
*          kind   =1 for sinintegral                  *
*                 =2 for cosintegral                  * 
*          a      begin value of x variable           *
*          b      end value of x variable             *
*          n      number of integration steps         *
*                                                     *
* OUTPUT:                                             *
*          res    the integral of FUNC(X) from a to b *
*                                                     *
******************************************************/
void Integral_Simpson(int kind, double a, double b, int n, double *res) {
  int i; double step,r;
  step=(b-a)/2/n;
  r=FUNC(kind,a);
  *res=(r+FUNC(kind,b))/2;
  for (i=1; i<2*n; i++) {
    r=FUNC(kind,a+i*step);
    if ((i%2) != 0) *res += r+r;
    else *res += r;
  }
  *res *= step*2/3;
}

/******************************************************
*     Si(x) = Integral of sin(u)/u from u=0 to u=x    *
******************************************************/           
double sinintegral(double x)  {
  double  x0=0.0; //begin x value
//        x;        end   x value
  int     nstep;  //number of integration steps
  double  result; //result of integral
  
  nstep = (int) (200*fabs(x)); //this should ensure about 14 exact digits
  
  Integral_Simpson(1,x0,x,nstep,&result);   // kind=1

  return result;

}

/******************************************************
* Calcualte Ci(x)= gamma + Ln(x) + I(x) with I(x) =   *
* Integral of (cos(u)-1)/u from u=0 to u=x (x > 0).   *
* Gamma = 0.57721566... (Euler's constant).           *
******************************************************/          
double cosintegral(double x)  {
  double x0=0.0;  //begin x value
//       x          end x value
  int    nstep;   //number of integration steps
  double result;  //result of integral

  nstep = (int) (200*fabs(x)); //this should ensure about 14 exact digits
  
  Integral_Simpson(2,x0,x,nstep,&result);   // kind=2

  return (gamma + log(x) + result);

}

// End of file sinint.cpp