/************************************************************* * 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 #include #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