CALGO 683 Exponential integral .
The description of all parameters and outputs can be found in the source code docstring.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| complex(kind=wp), | intent(in) | :: | z | |||
| integer, | intent(in) | :: | n | |||
| integer, | intent(in) | :: | kode | |||
| real(kind=wp), | intent(in) | :: | tol | |||
| integer, | intent(in) | :: | m | |||
| complex(kind=wp), | intent(out) | :: | cy(m) | |||
| integer, | intent(out) | :: | ierr |
SUBROUTINE cexint(z, n, kode, tol, m, cy, ierr) !! CALGO 683 Exponential integral \(\mathrm{E}_n(z)\). ! !! \(n \geq 1,\thinspace !! \lbrace z \in \mathbb{C} \mid -\pi \lt \arg(z) \leq \pi \rbrace \) ! !! The description of all parameters and outputs can be found in the source code !! docstring. ! ! ON KODE=1, CEXINT COMPUTES AN M MEMBER SEQUENCE OF COMPLEX(wp) ! EXPONENTIAL INTEGRALS CY(J)=E(N+J-1,Z), J=1,...,M, FOR ! POSITIVE ORDERS N,...,N+M-1 AND COMPLEX(wp) Z IN THE CUT PLANE ! -PI < ARG(Z) <= PI (N=1 AND Z=CMPLX(0.0,0.0) CANNOT HOLD AT ! THE SAME TIME). ON KODE=2, CEXINT COMPUTES SCALED FUNCTIONS ! ! CY(J)=E(N+J-1,Z)*CEXP(Z), J=1,...,M, ! ! WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND ! RIGHT HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE ! NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). ! ! INPUT ! Z - Z=CMPLX(X,Y), -PI < ARG(Z) <= PI ! N - INTEGER ORDER OF INITIAL E FUNCTION, N=1,2,... ! (N=1 AND Z=CMPLX(0.0,0.0) IS AN ERROR) ! KODE - A PARAMETER TO INDICATE THE SCALING OPTION ! KODE= 1 RETURNS ! CY(J)=E(N+J-1,Z), J=1,...,M ! = 2 RETURNS ! CY(J)=E(N+J-1,Z)*CEXP(Z), J=1,...,M ! TOL - PRECISION (ACCURACY) DESIRED FOR THE SEQUENCE, ! URND <= TOL <= 1.0E-3, WHERE URND IS LIMITED BY ! URND = MAX(UNIT ROUNDOFF,1.0E-18) AND UNIT ! ROUNDOFF = EPSILON(0.0_wp) ! M - NUMBER OF E FUNCTIONS IN THE SEQUENCE, M >= 1 ! ! OUTPUT ! CY - A COMPLEX(wp) VECTOR WHOSE FIRST M COMPONENTS CONTAIN ! VALUES FOR THE SEQUENCE ! CY(J)=E(N+J-1,Z) OR ! CY(J)=E(N+J-1,Z)*CEXP(Z), J=1,...,M ! DEPENDING ON KODE. ! IERR - ERROR FLAG ! IERR=0, NORMAL RETURN - COMPUTATION COMPLETED ! IERR=1, INPUT ERROR - NO COMPUTATION ! IERR=2, UNDERFLOW - FIRST M COMPONENTS OF CY ! SET TO ZERO, CY(J)=CMPLX(0.0,0.0), J=1,M, ! REAL(Z) > 0.0 TOO LARGE ON KODE=1 ! IERR=3, OVERFLOW - NO COMPUTATION, ! REAL(Z) < 0.0 TOO SMALL ON KODE=1 ! IERR=4, CABS(Z) OR N+M-1 LARGE - COMPUTATION DONE ! BUT LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION ! PRODUCE LESS THAN HALF OF MACHINE ACCURACY ! IERR=5, CABS(Z) OR N+M-1 TOO LARGE - NO COMPUTATION BECAUSE ! OF COMPLETE LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION ! IERR=6, ERROR - NO COMPUTATION, ! ALGORITHM TERMINATION CONDITION NOT MET. ! SEE LONG DESCRIPTION ABOUT PARAMETER ICMAX. ! IERR=7, ERROR - NO COMPUTATION, ! DISCRIMINATION ERROR. THIS CONDITION SHOULD NEVER OCCUR. ! ! LONG DESCRIPTION ! ! CEXINT USES A COMBINATION OF POWER SERIES AND BACKWARD RECURRENCE ! DESCRIBED IN REF. 2 FOR THE COMPLEX(wp) Z PLANE EXCEPT FOR A STRIP ! 2*YB WIDE ABOUT THE NEGATIVE REAL AXIS, WHERE ANALYTIC CONTINUATION IS ! CARRIED OUT BY LIMITED USE OF TAYLOR SERIES. ! THE SWITCH FROM BACKWARD RECURRENCE TO TAYLOR SERIES IS NECESSARY BECAUSE ! BACKWARD RECURRENCE IS SLOWLY CONVERGENT NEAR THE NEGATIVE REAL AXIS. ! THE BOUNDARIES Y=-YB AND Y=YB WERE DETERMINED SO THAT BACKWARD RECURRENCE ! WOULD CONVERGE EASILY WITH N AS LARGE AS 100 AND TOL AS SMALL AS 1.0e-18. ! SUBROUTINE CEXENZ DOES THE BACKWARD RECURRENCE AND SUBROUTINE CACEXI DOES ! THE ANALYTIC CONTINUATION. TO START THE CONTINUATION, CACEXI CALLS CEXENZ ! WITH ZB=CMPLX(X,YB). ! IF CEXENZ RETURNS IERR=6, THEN YB IS INCREASED BY 0.5 UNTIL CEXENZ RETURNS ! IERR=0 OR 10 TRIES, WHICHEVER COMES FIRST. ! WHEN IERR=0, THEN THE ANALYTIC CONTINUATION PROCEEDS VERTICALLY DOWN FROM ! ZB=CMPLX(X,YB) TO Z=CMPLX(X,Y), 0 <= Y < YB. ! CONJUGATION IS USED FOR Y < 0. YB INCREASES AS TOL DECREASES TO KEEP ! CONVERGENCE RATES UP AND RECURRENCE DOWN. ! ! PARAMETER ICDIM=250 ALLOCATES STORAGE FOR THE COEFFICIENTS OF THE BACKWARD ! RECURRENCE ALGORITHM. IF THE ALGORITHM TERMINATION CONDITION IS NOT MET ! IN ICDIM STEPS, THEN RECURRENCE PROCEEDS WITH NO ADDITIONAL STORAGE UNTIL ! THE TERMINATION CONDITION IS MET OR THE LIMIT ICMAX=2000 IS EXCEEDED. ! THE PURPOSE OF STORAGE IS TO MAKE THE ALGORITHM MORE EFFICIENT. ! THE TERMINATION CONDITION IS MET IN LESS THAN 250 STEPS OVER MOST OF ! THE COMPLEX(wp) PLANE EXCLUDING THE STRIP ABS(Y) < YB, X < 0. ! EXCEPTIONS TO THIS RULE ARE GENERATED NEAR STRIP BOUNDARIES WHEN N+M-1 ! AND ABS(Z) ARE LARGE AND NEARLY EQUAL. IN THESE CASES, THE CONVERGENCE ! IS VERY SLOW AND ADDITIONAL RECURRENCE (UP TO ICMAX) MUST BE USED. ! ON THE OTHERHAND, THESE REGIONS OF SLOW CONVERGENCE ARE KEPT SMALL BY ! ADJUSTING YB AS A FUNCTION OF TOL. THESE REGIONS COULD BE ELIMINATED ! ENTIRELY BY MAKING YB SUFFICIENTLY LARGE, BUT THE EXPENSE AND INSTABILITY ! OF CONTINUATION BY TAYLOR SERIES NOT ONLY GOES UP, BUT THE COMPUTATIONAL ! EXPENSE BECOMES EXCESSIVELY LARGE IN OTHER PARTS OF THE LEFT HALF PLANE ! (Y < YB) WHERE THE BACKWARD RECURRENCE ALGORITHM WOULD CONVERGE RAPIDLY. ! ! DERIVATIVES FOR SUCCESSIVE POWER SERIES ARE NOT COMPUTED BY EVALUATING ! DERIVATIVES OF A PREVIOUS POWER SERIES. BECAUSE OF THE RELATION ! ! (1) DE(N,Z)/DZ = - E(N-1,Z), ! ! SUCCESSIVE DERIVATIVES AT Z ARE GIVEN BY LOWER ORDER FUNCTIONS AND CAN BE ! COMPUTED IN A STABLE FASHION BY BACKWARD RECURRENCE USING (2) PROVIDED ! THAT THE BEGINNING ORDER NUB IS SMALLER THAN THE ARGUMENT. ! TO ACHIEVE THIS FOR ALL INTERMEDIATE VALUES ZZ BETWEEN ZB AND Z, WE TAKE ! NUB=MINO(N+M-1,INT(CABS(Z)+0.5)). ! TO START, E(NUB,ZB) IS EVALUATED BY THE BACKWARD RECURRENCE ALGORITHM OF ! REF. 3. TO CONTINUE THE FUNCTION FROM ZB TO Z VIA INTERMEDIATE VALUES ZZ, ! DERIVATIVES OF E(NUB,ZB) ARE COMPUTED BY BACKWARD RECURRENCE ON (2). ! THIS ALLOWS A STEP (NO LARGER THAN 0.5) TO ZZ FOR E(NUB,ZZ) USING THE ! TAYLOR SERIES ABOUT ZB. NOW, WE APPLY (2) AGAIN STARTING AT E(NUB,ZZ) FOR ! THE DERIVATIVES AT ZZ AND TAKE ANOTHER STEP, ETC. NOTICE THAT THE ! STABILITY CONDITION FOR BACKWARD RECURRENCE, NUB <= ABS(Z) <= ABS(ZZ) ! <= ABS(ZB), IS SATISFIED FOR ALL INTERMEDIATE VALUES ZZ. THE FINAL ! SEQUENCE FOR ORDERS N,...,N+M-1 IS GENERATED FROM (2) BY BACKWARD ! RECURRENCE, FORWARD RECURRENCE OR BOTH ONCE E(NUB,Z) HAS BEEN COMPUTED. ! ! RECURRENCE WITH THE RELATION ! ! (2) N*E(N+1,Z) + Z*E(N,Z) = CEXP(-Z) ! ! IN A DIRECTION AWAY FROM THE INTEGER CLOSEST TO ABS(Z) IS STABLE. ! FOR NEGATIVE ORDERS, THE RECURRENCE ! ! E( 0,Z) = CEXP(-Z)/Z ! E(-N,Z) = ( CEXP(-Z)+N*E(-N+1,Z) )/Z ,N=1,2,... ! ! IS NUMERICALLY STABLE FOR ALL Z. ! ! THE (CAPITAL) SINE AND COSINE INTEGRALS CAN BE COMPUTED FROM ! ! SI(Z) = (E(1,I*Z)-E(1,-I*Z))/(2*I) + PI/2 ! CI(Z) = -(E(1,I*Z)+E(1,-I*Z))/2 ! ! IN -PI/2 < ARG(Z) <= PI/2, (I**2=-1), WHILE THE PRINCIPAL ! VALUED EXPONENTIAL INTEGRAL EI(X) CAN BE COMPUTED FROM ! ! EI( X) = -(E(1,-X+I*0)+E(1,-X-I*0))/2 = -REAL(E(1,-X)) ! EI(-X) = -REAL(E(1,X)) ! ! FOR X > 0.0 TO AN ACCURACY TOL. IF Z = X > 0 THEN THE REAL SINE AND ! COSINE INTEGRALS ARE GIVEN BY ! ! SI(X) = AIMAG(E(1,I*X)) + PI/2 ! CI(X) = -REAL(E(1,I*X)) . ! ! THE ANALYTIC CONTINUATION TO OTHER SHEETS CAN BE DONE BY THE RELATIONS ! ! E(N,Z*CEXP(2*PI*M*I)) = E(N,Z) - 2*PI*M*I*(-Z)**(N-1)/(N-1)! ! ! WHERE M=+1 OR M=-1 AND I**2=-1. ! ! REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND ! I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF COMMERCE, 1955. ! ! COMPUTATION OF EXPONENTIAL INTEGRALS OF COMPLEX(wp) ARGUMENT ! BY D. E. AMOS, ACM TRANS. MATH. SOFTWARE ! ! COMPUTATION OF EXPONENTIAL INTEGRALS BY D. E. AMOS, ACM TRANS. ! MATH. SOFTWARE, VOL 6, NO. 3 SEPTEMBER 1980, PP. 365-377; ! ALGORITHM 556, EXPONENTIAL INTEGRALS, PP. 420-428. ! ! REMARK ON ALGORITHM 556 ! BY D. E. AMOS, ACM TRANS. MATH. SOFTWARE, VOL 9, NO. 4 ! DECEMBER 1983, P. 525. ! ! UNIFORM ASYMPTOTIC EXPANSIONS FOR EXPONENTIAL INTEGRALS E(N,X) ! AND BICKLEY FUNCTIONS KI(N,X) BY D. E. AMOS, ACM TRANS. MATH. ! SOFTWARE, VOL 9, NO. 4, DECEMBER. 1983, PP. 467-479; ! ALGORITHM 609, A PORTABLE FORTRAN SUBROUTINE FOR BICKLEY ! FUNCTIONS KI(N,X), PP. 480-493. ! ! ROUTINES CALLED CACEXI, CEXENZ COMPLEX(wp), INTENT(IN) :: z INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: kode REAL(wp), INTENT(IN) :: tol INTEGER, INTENT(IN) :: m COMPLEX(wp), INTENT(OUT) :: cy(m) INTEGER, INTENT(OUT) :: ierr INTEGER :: i, k, k1, k2 REAL(wp) :: aa, alim, az, bb, ca(250), d, elim, & fn, rbry, r1m5, urnd, x, y, yb, htol COMPLEX(wp) :: cb(250) !----------------------------------------------------------------------- ! DIMENSION CA(ICDIM),CB(ICDIM) !----------------------------------------------------------------------- INTEGER, PARAMETER :: icdim = 250 ! FIRST EXECUTABLE STATEMENT CEXINT ierr = 0 x = REAL(z, KIND=wp) y = AIMAG(z) IF (x == 0.0_wp .AND. y == 0.0_wp .AND. n == 1) ierr = 1 IF (n < 1) ierr = 1 IF (kode < 1 .OR. kode > 2) ierr = 1 IF (m < 1) ierr = 1 urnd = MAX(EPSILON(0.0_wp), 1.0E-18_wp) IF (tol < urnd .OR. tol > 1.0E-3_wp) ierr = 1 IF (ierr /= 0) RETURN IF (x /= 0.0_wp .OR. y /= 0.0_wp .OR. n <= 1) THEN !----------------------------------------------------------------------- ! SET PARAMETERS RELATED TO MACHINE CONSTANTS. ! URND IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. ! ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. ! EXP(-ELIM) < EXP(-ALIM)=EXP(-ELIM)/URND AND ! EXP(ELIM) > EXP(ALIM)=EXP(ELIM)*URND ARE INTERVALS NEAR ! UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. !----------------------------------------------------------------------- k1 = MINEXPONENT(0.0_wp) k2 = MAXEXPONENT(0.0_wp) r1m5 = LOG10(DBLE(RADIX(0.0_wp))) k = MIN(ABS(k1),ABS(k2)) elim = 2.303_wp * (k*r1m5 - 3.0_wp) k1 = DIGITS(0.0_wp) - 1 aa = r1m5 * k1 aa = aa * 2.303 alim = elim + MAX(-aa, -41.45_wp) rbry = 2.0 IF (urnd > 1.0E-8_wp) rbry = 1.0_wp !----------------------------------------------------------------------- ! TEST VARIABLES FOR RANGE. ABS(Z) CANNOT BE LARGER THAN THE ARGUMENT ! OF THE INT( ) FUNCTION. !----------------------------------------------------------------------- az = ABS(z) fn = n+m-1 aa = 0.5_wp / urnd bb = HUGE(0.0_wp) * 0.5_wp aa = MIN(aa,bb) IF (az > aa) GO TO 20 IF (fn > aa) GO TO 20 aa = SQRT(aa) IF (az > aa) ierr = 4 IF (fn > aa) ierr = 4 IF (x >= 0.0_wp) THEN !----------------------------------------------------------------------- ! BACKWARD RECURRENCE FOR THE RIGHT HALF PLANE, X >= 0.0E0 !----------------------------------------------------------------------- CALL cexenz(z, n, kode, m, cy, ierr, rbry, tol, elim, alim, icdim, ca, cb) RETURN END IF IF (az <= rbry) THEN !----------------------------------------------------------------------- ! POWER SERIES FOR ABS(Z) <= RBRY AND X < 0.0E0 !----------------------------------------------------------------------- CALL cexenz(z, n, kode, m, cy, ierr, rbry, tol, elim, alim, icdim, ca, cb) RETURN END IF d = -0.4342945_wp * LOG(tol) yb = 10.5_wp - 0.538460_wp * (18.0_wp-d) IF (ABS(y) >= yb) THEN !----------------------------------------------------------------------- ! BACKWARD RECURRENCE EXTERIOR TO THE STRIP ABS(Y) < YB, X < 0.0 !----------------------------------------------------------------------- htol = 0.125_wp * tol CALL cexenz(z, n, kode, m, cy, ierr, rbry, htol, elim, alim, icdim, ca, cb) RETURN END IF !----------------------------------------------------------------------- ! TAYLOR SERIES IN CACEXI FOR ANALYTIC CONTINUATION !----------------------------------------------------------------------- CALL cacexi(z, n, kode, m, cy, ierr, yb, rbry, tol, elim, alim, icdim, ca, cb ) RETURN END IF DO i = 1, m cy(i) = 1.0_wp / (n+i-2) END DO RETURN 20 ierr = 5 RETURN END SUBROUTINE cexint