DECLARE SUB E1Z (Z#(), CE1#()) DECLARE FUNCTION CABS# (Z#()) DECLARE SUB CDIV (Z1#(), Z2#(), Z#()) DECLARE SUB CEXP (Z#(), Z1#()) DECLARE SUB CLOG (ZA#(), ZB#(), IERR%) DECLARE SUB CMUL (Z1#(), Z2#(), Z#()) '********************************************************************* '* 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.095920874794272D-03 -6.900179262212203D-03 * '* * '* Please enter x and y (z =x+iy): 25,10 * '* * '* z Re[E1(z)] Im[E1(z)] * '* -------------------------------------------------------------- * '* 25.0 10.0 -5.534833138109094D-14 7.834277777448312D-14 * '* * '* ----------------------------------------------------------------- * '* REFERENCE:"Fortran Routines for Computation of Special Functions, * '* jin.ece.uiuc.edu/routines/routines.html". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************************* 'PROGRAM ME1Z DEFDBL A-H, O-Z DEFINT I-N DIM CE1(2), Z(2) 'Complex numbers CLS PRINT INPUT " Please enter x and y (z =x+iy): ", X, Y Z(1) = X: Z(2) = Y CALL E1Z(Z(), CE1()) F\$ = " ###.# ###.# " PRINT PRINT " z Re[E1(z)] Im[E1(z)]" PRINT " -----------------------------------------------------------------" PRINT USING F\$; X; Y; PRINT CE1(1); " "; CE1(2) PRINT END 'of main program 'ABSOLUTE VALUE OF A COMPLEX NUMBER Z=X+I*Y FUNCTION CABS (Z()) X = ABS(Z(1)) Y = ABS(Z(2)) IF X = 0# THEN W = Y ELSE IF Y = 0# THEN W = X ELSE IF X > Y THEN W = X * SQR(1# + (Y / X) * (Y / X)) ELSE W = Y * SQR(1# + (X / Y) * (Y / X)) END IF END IF END IF CABS = W END FUNCTION ' Z=Z1/Z2 SUB CDIV (Z1(), Z2(), Z()) D = Z2(1) * Z2(1) + Z2(2) * Z2(2) IF D < 1E-12 THEN EXIT SUB Z(1) = (Z1(1) * Z2(1) + Z1(2) * Z2(2)) / D Z(2) = (Z1(2) * Z2(1) - Z1(1) * Z2(2)) / D END SUB ' z1=exp(z) SUB CEXP (Z(), Z1()) IF EXP(Z(1)) > 1E+16 THEN TMP = 1E+16 ELSE TMP = EXP(Z(1)) END IF Z1(1) = TMP * COS(Z(2)) Z1(2) = TMP * SIN(Z(2)) END SUB SUB CLOG (ZA(), ZB(), 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 DPI = 3.141592653589793# DHPI = 1.570796326794897# IERR = 0 AR = ZA(1): AI = ZA(2) IF AR = 0# THEN GOTO 10 IF AI = 0# THEN GOTO 20 DTHETA = ATN(AI / AR) IF DTHETA <= 0# THEN GOTO 40 IF AR < 0! THEN DTHETA = DTHETA - DPI GOTO 50 10 IF AI = 0# THEN GOTO 60 BI = DHPI BR = LOG(ABS(AI)) IF AI < 0! THEN BI = -BI GOTO 70 20 IF AR > 0! THEN GOTO 30 BR = LOG(ABS(AR)) BI = DPI GOTO 70 30 BR = LOG(AR) BI = 0# GOTO 70 40 IF AR < 0# THEN DTHETA = DTHETA + DPI 50 ZM = CABS(ZA()) BR = LOG(ZM) BI = DTHETA GOTO 70 60 IERR = 1 70 ZB(1) = BR: ZB(2) = BI END SUB 'CLOG ' Z=Z1*Z2 SUB CMUL (Z1(), Z2(), Z()) Z(1) = Z1(1) * Z2(1) - Z1(2) * Z2(2) Z(2) = Z1(1) * Z2(2) + Z1(2) * Z2(1) END SUB SUB E1Z (Z(), CE1()) ' ===================================================== ' Purpose: Compute complex exponential integral E1(z) ' Input : z --- Argument of E1(z) ' Output: CE1 --- E1(z) ' ===================================================== DIM CK(2), CONE(2), CR(2), CT(2), CT0(2), TMP(2), TMP1(2)'Complex numbers EL = .5772156649015328# PI = 4# * ATN(1#) X = Z(1) A0 = CABS(Z()) IF A0 = 0# THEN CE1(1) = 1D+200: CE1(2) = 0# ELSEIF A0 <= 10# OR (X < 0# AND A0 < 20#) THEN CE1(1) = 1#: CE1(2) = 0# CR(1) = 1#: CR(2) = 0# FOR K = 1 TO 150 'CR=-CR*K*Z/(K+1.0)^2 CALL CMUL(CR(), Z(), TMP()) CR(1) = TMP(1): CR(2) = TMP(2) CR(1) = -K * CR(1) / (K + 1!) ^ 2 CR(2) = -K * CR(2) / (K + 1!) ^ 2 CE1(1) = CE1(1) + CR(1) CE1(2) = CE1(2) + CR(2) IF CABS(CR()) <= CABS(CE1()) * .000000000000001# THEN GOTO 15 NEXT K 15 'CE1=-EL-CLOG(Z)+Z*CE1 CALL CLOG(Z(), TMP(), IERR) CALL CMUL(Z(), CE1(), TMP1()) CE1(1) = -EL - TMP(1) + TMP1(1) CE1(2) = -TMP(2) + TMP1(2) ELSE CT0(1) = 0#: CT0(2) = 0# FOR K = 120 TO 1 STEP -1 'CT0=K/(1.0+K/(Z+CT0)) CK(1) = 1# * K: CK(2) = 0# TMP(1) = 1# + CK(1): TMP(2) = 0# TMP1(1) = Z(1) + CT0(1) TMP1(2) = Z(2) + CT0(2) CALL CDIV(TMP(), TMP1(), CT()) TMP(1) = CT(1): TMP(2) = CT(2) CALL CDIV(CK(), TMP(), CT0()) NEXT K 'CT=1.0/(Z+CT0) CONE(1) = 1#: CONE(2) = 0# TMP(1) = Z(1) + CT0(1): TMP(2) = Z(2) + CT0(2) CALL CDIV(CONE(), TMP(), CT()) TMP(1) = -Z(1): TMP(2) = -Z(2) CALL CEXP(TMP(), TMP1()) CALL CMUL(TMP1(), CT(), CE1()) IF X <= 0# AND Z(2) = 0# THEN 'CE1=CE1-PI*(0.0D0,1.0D0) TMP(1) = 0#: TMP(2) = 1# CE1(1) = CE1(1) - PI * TMP(1) CE1(2) = CE1(2) - PI * TMP(2) END IF END IF END SUB