*DECK DEXINT SUBROUTINE DEXINT (X, N, KODE, M, TOL, EN, NZ, IERR) C***BEGIN PROLOGUE DEXINT C***PURPOSE Compute an M member sequence of exponential integrals C E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. C***LIBRARY SLATEC C***CATEGORY C5 C***TYPE DOUBLE PRECISION (EXINT-S, DEXINT-D) C***KEYWORDS EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS C***AUTHOR Amos, D. E., (SNLA) C***DESCRIPTION C C DEXINT computes M member sequences of exponential integrals C E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. The C exponential integral is defined by C C E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N C C where X=0.0 and N=1 cannot occur simultaneously. Formulas C and notation are found in the NBS Handbook of Mathematical C Functions (ref. 1). C C The power series is implemented for X .LE. XCUT and the C confluent hypergeometric representation C C E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X) C C is computed for X .GT. XCUT. Since sequences are computed in C a stable fashion by recurring away from X, A is selected as C the integer closest to X within the constraint N .LE. A .LE. C N+M-1. For the U computation, A is further modified to be the C nearest even integer. Indices are carried forward or C backward by the two term recursion relation C C K*E(K+1,X) + X*E(K,X) = EXP(-X) C C once E(A,X) is computed. The U function is computed by means C of the backward recursive Miller algorithm applied to the C three term contiguous relation for U(A+K,A,X), K=0,1,... C This produces accurate ratios and determines U(A+K,A,X), and C hence E(A,X), to within a multiplicative constant C. C Another contiguous relation applied to C*U(A,A,X) and C C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to C E(A+1,X). The normalizing constant C is obtained from the C two term recursion relation above with K=A. C C The maximum number of significant digits obtainable C is the smaller of 14 and the number of digits carried in C double precision arithmetic. C C Description of Arguments C C Input * X and TOL are double precision * C X X .GT. 0.0 for N=1 and X .GE. 0.0 for N .GE. 2 C N order of the first member of the sequence, N .GE. 1 C (X=0.0 and N=1 is an error) C KODE a selection parameter for scaled values C KODE=1 returns E(N+K,X), K=0,1,...,M-1. C =2 returns EXP(X)*E(N+K,X), K=0,1,...,M-1. C M number of exponential integrals in the sequence, C M .GE. 1 C TOL relative accuracy wanted, ETOL .LE. TOL .LE. 0.1 C ETOL is the larger of double precision unit C roundoff = D1MACH(4) and 1.0D-18 C C Output * EN is a double precision vector * C EN a vector of dimension at least M containing values C EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M C depending on KODE C NZ underflow indicator C NZ=0 a normal return C NZ=M X exceeds XLIM and an underflow occurs. C EN(K)=0.0D0 , K=1,M returned on KODE=1 C IERR error flag C IERR=0, normal return, computation completed C IERR=1, input error, no computation C IERR=2, error, no computation C algorithm termination condition not met C C***REFERENCES M. Abramowitz and I. A. Stegun, Handbook of C Mathematical Functions, NBS AMS Series 55, U.S. Dept. C of Commerce, 1955. C D. E. Amos, Computation of exponential integrals, ACM C Transactions on Mathematical Software 6, (1980), C pp. 365-377 and pp. 420-428. C***ROUTINES CALLED D1MACH, DPSIXN, I1MACH C***REVISION HISTORY (YYMMDD) C 800501 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900326 Removed duplicate information from DESCRIPTION section. C (WRB) C 910408 Updated the REFERENCES section. (WRB) C 920207 Updated with code with a revision date of 880811 from C D. Amos. Included correction of argument list. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DEXINT