/******************************************************************** * Calculate the complex exponential integral E1(z) with a complex * * argument * * ----------------------------------------------------------------- * * EXPLANATION: * * * * ========================================================= * * Purpose: This program computes the complex exponential * * integral E1(z) using subroutine E1Z * * Example: * * z Re[E1(z)] Im[E1(z)] * * ----------------------------------------------- * * 3.0 2.0 -.90959209D-02 -.69001793D-02 * * 3.0 -2.0 -.90959209D-02 .69001793D-02 * * -3.0 2.0 -.28074891D+01 .59603353D+01 * * -3.0 -2.0 -.28074891D+01 -.59603353D+01 * * 25.0 10.0 -.29302080D-12 .40391222D-12 * * 25.0 -10.0 -.29302080D-12 -.40391222D-12 * * -25.0 10.0 .27279957D+10 -.49430610D+09 * * -25.0 -10.0 .27279957D+10 .49430610D+09 * * ========================================================= * * * * ----------------------------------------------------------------- * * SAMPLE RUNS: * * * * Please enter x and y (z =x+iy): 3 2 * * * * z Re[E1(z)] Im[E1(z)] * * ----------------------------------------------- * * 3.0 2.0 -9.095921e-003 -6.900179e-003 * * * * Please enter x and y (z =x+iy): 25 10 * * * * z Re[E1(z)] Im[E1(z)] * * ----------------------------------------------- * * 25.0 10.0 -5.534833e-014 7.834278e-014 * * * * ----------------------------------------------------------------- * * 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 CE1[2], Z[2]; //complex numbers double X, Y; //ABSOLUTE VALUE OF A COMPLEX NUMBER Z=X+I*Y double CABS(double *Z) { double X,Y,W; X=fabs(Z[0]); Y=fabs(Z[1]); if (X == 0.0) W=Y; else if (Y == 0.0) W=X; else if (X > Y) W=X*sqrt(1.0+(Y/X)*(Y/X)); else W=Y*sqrt(1.0+(X/Y)*(Y/X)); return W; } // Z=Z1/Z2 void CDIV(double *Z1, double *Z2, double *Z) { double D; D=Z2[0]*Z2[0]+Z2[1]*Z2[1]; if (D<1e-12) return; Z[0]=(Z1[0]*Z2[0]+Z1[1]*Z2[1])/D; Z[1]=(Z1[1]*Z2[0]-Z1[0]*Z2[1])/D; } // Z1=exp(Z) void CEXP(double *Z, double *Z1) { double tmp; if (exp(Z[0]) > 1e16) tmp=1e16; else tmp=exp(Z[0]); Z1[0]=tmp*cos(Z[1]); Z1[1]=tmp*sin(Z[1]); } void CLOG(double *ZA, double *ZB, int *IERR) { /***BEGIN PROLOGUE ZLOG DOUBLE PRECISION COMPLEX LOGARITHM B:=CLOG(A) IERR:=0,NORMAL RETURN IERR:=1, Z:=CMPLX(0.0,0.0) ***ROUTINES CALLED ZABS ***END PROLOGUE ZLOG */ // Labels: e10,e20,e30,e40,e50,e60,e70 double AR,AI,BR,BI,ZM, DTHETA, DPI, DHPI; DPI = 3.141592653589793238462643383e+00; DHPI= 1.570796326794896619231321696e+00; *IERR=0; AR=ZA[0]; AI=ZA[1]; if (AR == 0.0) goto e10; if (AI == 0.0) goto e20; DTHETA = atan(AI/AR); if (DTHETA <= 0.0) goto e40; if (AR < 0.0) DTHETA -= DPI; goto e50; e10: if (AI == 0.0) goto e60; BI = DHPI; BR = log(fabs(AI)); if (AI < 0.0) BI = -BI; goto e70; e20: if (AR > 0.0) goto e30; BR = log(fabs(AR)); BI = DPI; goto e70; e30: BR = log(AR); BI = 0.0; goto e70; e40: if (AR < 0.0) DTHETA += DPI; e50: ZM = CABS(ZA); BR = log(ZM); BI = DTHETA; goto e70; e60: *IERR=1; e70: ZB[0]=BR; ZB[1]=BI; } //CLOG //Z=Z1*Z2 void CMUL(double *Z1, double *Z2, double *Z) { Z[0]=Z1[0]*Z2[0] - Z1[1]*Z2[1]; Z[1]=Z1[0]*Z2[1] + Z1[1]*Z2[0]; } double Sqr(double x) { return x*x; } void E1Z(double *Z, double *CE1) { /* ===================================================== Purpose: Compute complex exponential integral E1(z) Input : z --- Argument of E1(z) Output: CE1 --- E1(z) ====================================================== */ // Label e15 double A0,EL,X; double CK[2],CONE[2],CR[2],CT[2],CT0[2],TMP[2],TMP1[2]; int IERR, K; double PI = 4.0*atan(1.0); EL=0.5772156649015328; X=Z[0]; A0=CABS(Z); if (A0 == 0.0) { CE1[0]=1.0E+200; CE1[1]=0.0; } else if (A0 <= 10.0 || X < 0.0 && A0 <20.0) { CE1[0]=1.0; CE1[1]=0.0; CR[0]=1.0; CR[1]=0.0; for (K=1; K<=150; K++) { // CR=-CR*K*Z/Sqr(K+1.0) CMUL(CR,Z,TMP); CR[0]=TMP[0]; CR[1]=TMP[1]; CR[0]=-K*CR[0]/Sqr(K+1.0); CR[1]=-K*CR[1]/Sqr(K+1.0); CE1[0]=CE1[0]+CR[0]; CE1[1]=CE1[1]+CR[1]; if (CABS(CR) <= CABS(CE1)*1.0e-15) goto e15; } e15: // CE1=-EL-CLOG(Z)+Z*CE1 CLOG(Z,TMP,&IERR); CMUL(Z,CE1,TMP1); CE1[0]=-EL-TMP[0]+TMP1[0]; CE1[1]= -TMP[1]+TMP1[1]; } else { CT0[0]=0.0; CT0[1]=0.0; for (K=120; K>0; K--) { // CT0=K/(1.0+K/(Z+CT0)) CK[0]=1.0*K; CK[1]=0.0; TMP[0]=1.0+CK[0]; TMP[1]=0.0; TMP1[0]=Z[0]+CT0[0]; TMP1[1]=Z[1]+CT0[1]; CDIV(TMP,TMP1,CT); CDIV(CK,CT,CT0); } // CT=1.0/(Z+CT0) CONE[0]=1.0; CONE[1]=0.0; TMP[0]=Z[0]+CT0[0]; TMP[1]=Z[1]+CT0[1]; CDIV(CONE,TMP,CT); TMP[0]=-Z[0]; TMP[1]=-Z[1]; CEXP(TMP,TMP1); CMUL(TMP1,CT,CE1); if (X <= 0.0 && Z[1] == 0.0) { // CE1=CE1-PI*(0.0D0,1.0D0) TMP[0]=0.0; TMP[1]=1.0; CE1[0]=CE1[0]-PI*TMP[0]; CE1[1]=CE1[1]-PI*TMP[1]; } } } void main() { printf("\n Please enter x and y (z =x+iy): "); scanf("%lf %lf", &X, &Y); Z[0]=X; Z[1]=Y; E1Z(Z, CE1); printf("\n"); printf(" z Re[E1(z)] Im[E1(z)]\n"); printf(" -------------------------------------------------\n"); printf(" %5.1f %5.1f %15e %15e\n\n", X, Y, CE1[0], CE1[1]); } // end of file me1z.cpp