!********************************************************************* !* 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 -0.90959209D-02 -0.69001793D-02 * !* * !* Please enter x and y (z =x+iy): 25 10 * !* * !* z Re[E1(z)] Im[E1(z)] * !* ----------------------------------------------- * !* 25.0 10.0 -0.29302080D-12 0.40391222D-12 * !* * !* ----------------------------------------------------------------- * !* REFERENCE:"Fortran Routines for Computation of Special Functions, * !* jin.ece.uiuc.edu/routines/routines.html". * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !********************************************************************* PROGRAM ME1Z IMPLICIT COMPLEX*16 (C,Z) IMPLICIT DOUBLE PRECISION (D-H,O-Y) WRITE(*,10,advance='no'); READ *, X, Y Z=CMPLX(X,Y) CALL E1Z(Z,CE1) WRITE(*,*) WRITE(*,*)' z Re[E1(z)] Im[E1(z)]' WRITE(*,*)' -----------------------------------------------' WRITE(*,20) X, Y, CE1 WRITE(*,*) 10 FORMAT(/' Please enter x and y (z =x+iy): ') 20 FORMAT(1X,F5.1,2X,F5.1,1X,2D17.8) END SUBROUTINE E1Z(Z,CE1) ! ==================================================== ! Purpose: Compute complex exponential integral E1(z) ! Input : z --- Argument of E1(z) ! Output: CE1 --- E1(z) ! ==================================================== IMPLICIT COMPLEX*16 (C,Z) IMPLICIT DOUBLE PRECISION (D-H,O-Y) PI=3.141592653589793D0 EL=0.5772156649015328D0 X=REAL(Z) A0=CDABS(Z) IF (A0.EQ.0.0D0) THEN CE1=(1.0D+300,0.0D0) ELSE IF (A0.LE.10.0.OR.X.LT.0.0.AND.A0.LT.20.0) THEN CE1=(1.0D0,0.0D0) CR=(1.0D0,0.0D0) DO 10 K=1,150 CR=-CR*K*Z/(K+1.0D0)**2 CE1=CE1+CR IF (CDABS(CR).LE.CDABS(CE1)*1.0D-15) GO TO 15 10 CONTINUE 15 CE1=-EL-CDLOG(Z)+Z*CE1 ELSE CT0=(0.0D0,0.0D0) DO 20 K=120,1,-1 CT0=K/(1.0D0+K/(Z+CT0)) 20 CONTINUE CT=1.0D0/(Z+CT0) CE1=CDEXP(-Z)*CT IF (X.LE.0.0.AND.DIMAG(Z).EQ.0.0) CE1=CE1-PI*(0.0D0,1.0D0) ENDIF RETURN END ! end of file me1z.f90