/******************************************************************
!*      Purpose: This program computes the cosine and sine        *
!*               integrals using subroutine CISIA                 *
!*      Input :  x  --- Argument of Ci(x) and Si(x)               *
!*      Output:  CI --- Ci(x)                                     *
!*               SI --- Si(x)                                     *
!*      Example:                                                  * 
!*                   x         Ci(x)          Si(x)               *
!*                ------------------------------------            *
!*                  0.0     - ì             .00000000             *
!*                  5.0     -.19002975     1.54993124             *
!*                 10.0     -.04545643     1.65834759             *
!*                 20.0      .04441982     1.54824170             * 
!*                 30.0     -.03303242     1.56675654             *
!*                 40.0      .01902001     1.58698512             *
!* -------------------------------------------------------------- *
!* REFERENCE: "Fortran Routines for Computation of Special        *
!*             Functions jin.ece.uiuc.edu/routines/routines.html" *
!*                                                                *
!*                              C++ Release By J-P Moreau, Paris. *
!*                                      (www.jpmoreau.fr)         *
!*****************************************************************/
#include <stdio.h>
#include <math.h>

      double CI,SI,X;


      void CISIA(double X, double *CI, double *SI) {

/*      =============================================
!       Purpose: Compute cosine and sine integrals
!                Si(x) and Ci(x)  ( x ò 0 )
!       Input :  x  --- Argument of Ci(x) and Si(x)
!       Output:  CI --- Ci(x)
!                SI --- Si(x)
!       ============================================= */
        double EL, EPS, P2, X2, XA, XA0, XA1, XCS, XF, XG, XG1, XG2, XR, XS, XSS;
        double BJ[101];
	int K, M;
        P2=1.570796326794897;
        EL=.5772156649015329;
        EPS=1.0e-15;
        X2=X*X;
        if (X == 0.0) { 
           *CI=-1.0e+100;
           *SI=0.0;
        }
        else if (X <= 16.0) {
           XR=-0.25*X2;
           *CI=EL+log(X)+XR;
           for (K=2; K<41; K++) {
              XR=-0.5*XR*(K-1)/(K*K*(2*K-1))*X2;
              *CI=*CI+XR;
              if (fabs(XR) < fabs(*CI)*EPS) goto e15;
	   }
e15:       XR=X;
           *SI=X;
           for (K=1; K<41; K++) { 
              XR=-0.5*XR*(2*K-1)/K/(4*K*K+4*K+1)*X2;
              *SI=*SI+XR;
              if (fabs(XR) < fabs(*SI)*EPS) return;
	   }
        }
        else if (X <= 32.0) {
           M=int(47.2+.82*X);
           XA1=0.0;
           XA0=1.0e-100;
	   for (K=M; K>0; K--) {
              XA=4.0*K*XA0/X-XA1;
              BJ[K]=XA;
              XA1=XA0;
              XA0=XA;
           }
           XS=BJ[1];
           for (K=3; K<=M; K+=2)  XS=XS+2.0*BJ[K];
           BJ[1]=BJ[1]/XS;
           for (K=2; K<=M; K++)   BJ[K]=BJ[K]/XS;
           XR=1.0;
           XG1=BJ[1];
	   for (K=2; K<=M; K++) {
              XR=0.25*XR*pow(2.0*K-3.0,2.0)/((K-1.0)*pow(2.0*K-1.0,2.0))*X;
              XG1=XG1+BJ[K]*XR;
           }
           XR=1.0;
           XG2=BJ[1];
           for (K=2; K<=M; K++) { 
              XR=0.25*XR*pow(2.0*K-5.0,2.0)/((K-1.0)*pow(2.0*K-3.0,2.0))*X;
              XG2=XG2+BJ[K]*XR;
	   }
           XCS=cos(X/2.0);
           XSS=sin(X/2.0);
           *CI=EL+log(X)-X*XSS*XG1+2*XCS*XG2-2*XCS*XCS;
           *SI=X*XCS*XG1+2*XSS*XG2-sin(X);
        }
        else {
           XR=1.0;
           XF=1.0;
	   for (K=1; K<10; K++) {
              XR=-2.0*XR*K*(2*K-1)/X2;
              XF=XF+XR;
           }
           XR=1.0/X;
           XG=XR;
           for (K=1; K<9; K++) {
              XR=-2.0*XR*(2*K+1)*K/X2;
              XG=XG+XR;
           }
           *CI=XF*sin(X)/X-XG*cos(X)/X;
           *SI=P2-XF*cos(X)/X-XG*sin(X)/X;
        }
      }


void main() {

        printf("\n   x         Ci(x)          Si(x)\n");
        printf(" -----------------------------------\n");
       
	X=0.0;
	CISIA(X, &CI, &SI);
        printf("   0.0      -i             0\n");

	X=5.0;
        CISIA(X, &CI, &SI);
        printf(" %5.1f      %f      %f\n", X, CI, SI); 

	for (X=10; X<=40; X+=10) {
          CISIA(X, &CI, &SI);
	  if (CI < 0.0)
            printf(" %5.1f      %f      %f\n", X, CI, SI);
          else
            printf(" %5.1f       %f      %f\n", X, CI, SI);
        }
	printf("\n");
}

// end of file mcisia.cpp