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