1*DECK EXINT 2 SUBROUTINE EXINT (X, N, KODE, M, TOL, EN, NZ, IERR) 3C***BEGIN PROLOGUE EXINT 4C***PURPOSE Compute an M member sequence of exponential integrals 5C E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. 6C***LIBRARY SLATEC 7C***CATEGORY C5 8C***TYPE SINGLE PRECISION (EXINT-S, DEXINT-D) 9C***KEYWORDS EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS 10C***AUTHOR Amos, D. E., (SNLA) 11C***DESCRIPTION 12C 13C EXINT computes M member sequences of exponential integrals 14C E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. The 15C exponential integral is defined by 16C 17C E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N 18C 19C where X=0.0 and N=1 cannot occur simultaneously. Formulas 20C and notation are found in the NBS Handbook of Mathematical 21C Functions (ref. 1). 22C 23C The power series is implemented for X .LE. XCUT and the 24C confluent hypergeometric representation 25C 26C E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X) 27C 28C is computed for X .GT. XCUT. Since sequences are computed in 29C a stable fashion by recurring away from X, A is selected as 30C the integer closest to X within the constraint N .LE. A .LE. 31C N+M-1. For the U computation, A is further modified to be the 32C nearest even integer. Indices are carried forward or 33C backward by the two term recursion relation 34C 35C K*E(K+1,X) + X*E(K,X) = EXP(-X) 36C 37C once E(A,X) is computed. The U function is computed by means 38C of the backward recursive Miller algorithm applied to the 39C three term contiguous relation for U(A+K,A,X), K=0,1,... 40C This produces accurate ratios and determines U(A+K,A,X), and 41C hence E(A,X), to within a multiplicative constant C. 42C Another contiguous relation applied to C*U(A,A,X) and 43C C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to 44C E(A+1,X). The normalizing constant C is obtained from the 45C two term recursion relation above with K=A. 46C 47C Description of Arguments 48C 49C Input 50C X X .GT. 0.0 for N=1 and X .GE. 0.0 for N .GE. 2 51C N order of the first member of the sequence, N .GE. 1 52C (X=0.0 and N=1 is an error) 53C KODE a selection parameter for scaled values 54C KODE=1 returns E(N+K,X), K=0,1,...,M-1. 55C =2 returns EXP(X)*E(N+K,X), K=0,1,...,M-1. 56C M number of exponential integrals in the sequence, 57C M .GE. 1 58C TOL relative accuracy wanted, ETOL .LE. TOL .LE. 0.1 59C ETOL = single precision unit roundoff = R1MACH(4) 60C 61C Output 62C EN a vector of dimension at least M containing values 63C EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M 64C depending on KODE 65C NZ underflow indicator 66C NZ=0 a normal return 67C NZ=M X exceeds XLIM and an underflow occurs. 68C EN(K)=0.0E0 , K=1,M returned on KODE=1 69C IERR error flag 70C IERR=0, normal return, computation completed 71C IERR=1, input error, no computation 72C IERR=2, error, no computation 73C algorithm termination condition not met 74C 75C***REFERENCES M. Abramowitz and I. A. Stegun, Handbook of 76C Mathematical Functions, NBS AMS Series 55, U.S. Dept. 77C of Commerce, 1955. 78C D. E. Amos, Computation of exponential integrals, ACM 79C Transactions on Mathematical Software 6, (1980), 80C pp. 365-377 and pp. 420-428. 81C***ROUTINES CALLED I1MACH, PSIXN, R1MACH 82C***REVISION HISTORY (YYMMDD) 83C 800501 DATE WRITTEN 84C 890531 Changed all specific intrinsics to generic. (WRB) 85C 890531 REVISION DATE from Version 3.2 86C 891214 Prologue converted to Version 4.0 format. (BAB) 87C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 88C 900326 Removed duplicate information from DESCRIPTION section. 89C (WRB) 90C 910408 Updated the REFERENCES section. (WRB) 91C 920207 Updated with code with a revision date of 880811 from 92C D. Amos. Included correction of argument list. (WRB) 93C 920501 Reformatted the REFERENCES section. (WRB) 94C***END PROLOGUE EXINT 95 REAL A,AA,AAMS,AH,AK,AT,B,BK,BT,CC,CNORM,CT,EM,EMX,EN, 96 1 ETOL,FNM,FX,PT,P1,P2,S,TOL,TX,X,XCUT,XLIM,XTOL,Y, 97 2 YT,Y1,Y2 98 REAL R1MACH,PSIXN 99 INTEGER I,IC,ICASE,ICT,IERR,IK,IND,IX,I1M,JSET,K,KK,KN,KODE,KS,M, 100 1 ML,MU,N,ND,NM,NZ 101 INTEGER I1MACH 102 DIMENSION EN(*), A(99), B(99), Y(2) 103C***FIRST EXECUTABLE STATEMENT EXINT 104 IERR = 0 105 NZ = 0 106 ETOL = MAX(R1MACH(4),0.5E-18) 107 IF (X.LT.0.0E0) IERR = 1 108 IF (N.LT.1) IERR = 1 109 IF (KODE.LT.1 .OR. KODE.GT.2) IERR = 1 110 IF (M.LT.1) IERR = 1 111 IF (TOL.LT.ETOL .OR. TOL.GT.0.1E0) IERR = 1 112 IF (X.EQ.0.0E0 .AND. N.EQ.1) IERR = 1 113 IF (IERR.NE.0) RETURN 114 I1M = -I1MACH(12) 115 PT = 2.3026E0*R1MACH(5)*I1M 116 XLIM = PT - 6.907755E0 117 BT = PT + (N+M-1) 118 IF (BT.GT.1000.0E0) XLIM = PT - LOG(BT) 119C 120 XCUT = 2.0E0 121 IF (ETOL.GT.2.0E-7) XCUT = 1.0E0 122 IF (X.GT.XCUT) GO TO 100 123 IF (X.EQ.0.0E0 .AND. N.GT.1) GO TO 80 124C----------------------------------------------------------------------- 125C SERIES FOR E(N,X) FOR X.LE.XCUT 126C----------------------------------------------------------------------- 127 TX = X + 0.5E0 128 IX = TX 129C----------------------------------------------------------------------- 130C ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1 131C ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2 132C----------------------------------------------------------------------- 133 ICASE = 2 134 IF (IX.GT.N) ICASE = 1 135 NM = N - ICASE + 1 136 ND = NM + 1 137 IND = 3 - ICASE 138 MU = M - IND 139 ML = 1 140 KS = ND 141 FNM = NM 142 S = 0.0E0 143 XTOL = 3.0E0*TOL 144 IF (ND.EQ.1) GO TO 10 145 XTOL = 0.3333E0*TOL 146 S = 1.0E0/FNM 147 10 CONTINUE 148 AA = 1.0E0 149 AK = 1.0E0 150 IC = 35 151 IF (X.LT.ETOL) IC = 1 152 DO 50 I=1,IC 153 AA = -AA*X/AK 154 IF (I.EQ.NM) GO TO 30 155 S = S - AA/(AK-FNM) 156 IF (ABS(AA).LE.XTOL*ABS(S)) GO TO 20 157 AK = AK + 1.0E0 158 GO TO 50 159 20 CONTINUE 160 IF (I.LT.2) GO TO 40 161 IF (ND-2.GT.I .OR. I.GT.ND-1) GO TO 60 162 AK = AK + 1.0E0 163 GO TO 50 164 30 S = S + AA*(-LOG(X)+PSIXN(ND)) 165 XTOL = 3.0E0*TOL 166 40 AK = AK + 1.0E0 167 50 CONTINUE 168 IF (IC.NE.1) GO TO 340 169 60 IF (ND.EQ.1) S = S + (-LOG(X)+PSIXN(1)) 170 IF (KODE.EQ.2) S = S*EXP(X) 171 EN(1) = S 172 EMX = 1.0E0 173 IF (M.EQ.1) GO TO 70 174 EN(IND) = S 175 AA = KS 176 IF (KODE.EQ.1) EMX = EXP(-X) 177 GO TO (220, 240), ICASE 178 70 IF (ICASE.EQ.2) RETURN 179 IF (KODE.EQ.1) EMX = EXP(-X) 180 EN(1) = (EMX-S)/X 181 RETURN 182 80 CONTINUE 183 DO 90 I=1,M 184 EN(I) = 1.0E0/(N+I-2) 185 90 CONTINUE 186 RETURN 187C----------------------------------------------------------------------- 188C BACKWARD RECURSIVE MILLER ALGORITHM FOR 189C E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X) 190C WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X. 191C U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION 192C----------------------------------------------------------------------- 193 100 CONTINUE 194 EMX = 1.0E0 195 IF (KODE.EQ.2) GO TO 130 196 IF (X.LE.XLIM) GO TO 120 197 NZ = M 198 DO 110 I=1,M 199 EN(I) = 0.0E0 200 110 CONTINUE 201 RETURN 202 120 EMX = EXP(-X) 203 130 CONTINUE 204 IX = X+0.5E0 205 KN = N + M - 1 206 IF (KN.LE.IX) GO TO 140 207 IF (N.LT.IX .AND. IX.LT.KN) GO TO 170 208 IF (N.GE.IX) GO TO 160 209 GO TO 340 210 140 ICASE = 1 211 KS = KN 212 ML = M - 1 213 MU = -1 214 IND = M 215 IF (KN.GT.1) GO TO 180 216 150 KS = 2 217 ICASE = 3 218 GO TO 180 219 160 ICASE = 2 220 IND = 1 221 KS = N 222 MU = M - 1 223 IF (N.GT.1) GO TO 180 224 IF (KN.EQ.1) GO TO 150 225 IX = 2 226 170 ICASE = 1 227 KS = IX 228 ML = IX - N 229 IND = ML + 1 230 MU = KN - IX 231 180 CONTINUE 232 IK = KS/2 233 AH = IK 234 JSET = 1 + KS - (IK+IK) 235C----------------------------------------------------------------------- 236C START COMPUTATION FOR 237C EN(IND) = C*U( A , A ,X) JSET=1 238C EN(IND) = C*U(A+1,A+1,X) JSET=2 239C FOR AN EVEN INTEGER A. 240C----------------------------------------------------------------------- 241 IC = 0 242 AA = AH + AH 243 AAMS = AA - 1.0E0 244 AAMS = AAMS*AAMS 245 TX = X + X 246 FX = TX + TX 247 AK = AH 248 XTOL = TOL 249 IF (TOL.LE.1.0E-3) XTOL = 20.0E0*TOL 250 CT = AAMS + FX*AH 251 EM = (AH+1.0E0)/((X+AA)*XTOL*SQRT(CT)) 252 BK = AA 253 CC = AH*AH 254C----------------------------------------------------------------------- 255C FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD 256C RECURSION 257C----------------------------------------------------------------------- 258 P1 = 0.0E0 259 P2 = 1.0E0 260 190 CONTINUE 261 IF (IC.EQ.99) GO TO 340 262 IC = IC + 1 263 AK = AK + 1.0E0 264 AT = BK/(BK+AK+CC+IC) 265 BK = BK + AK + AK 266 A(IC) = AT 267 BT = (AK+AK+X)/(AK+1.0E0) 268 B(IC) = BT 269 PT = P2 270 P2 = BT*P2 - AT*P1 271 P1 = PT 272 CT = CT + FX 273 EM = EM*AT*(1.0E0-TX/CT) 274 IF (EM*(AK+1.0E0).GT.P1*P1) GO TO 190 275 ICT = IC 276 KK = IC + 1 277 BT = TX/(CT+FX) 278 Y2 = (BK/(BK+CC+KK))*(P1/P2)*(1.0E0-BT+0.375E0*BT*BT) 279 Y1 = 1.0E0 280C----------------------------------------------------------------------- 281C BACKWARD RECURRENCE FOR 282C Y1= C*U( A ,A,X) 283C Y2= C*(A/(1+A/2))*U(A+1,A,X) 284C----------------------------------------------------------------------- 285 DO 200 K=1,ICT 286 KK = KK - 1 287 YT = Y1 288 Y1 = (B(KK)*Y1-Y2)/A(KK) 289 Y2 = YT 290 200 CONTINUE 291C----------------------------------------------------------------------- 292C THE CONTIGUOUS RELATION 293C X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X) 294C WITH B=A+1 , C=A IS USED FOR 295C Y(2) = C * U(A+1,A+1,X) 296C X IS INCORPORATED INTO THE NORMALIZING RELATION 297C----------------------------------------------------------------------- 298 PT = Y2/Y1 299 CNORM = 1.0E0 - PT*(AH+1.0E0)/AA 300 Y(1) = 1.0E0/(CNORM*AA+X) 301 Y(2) = CNORM*Y(1) 302 IF (ICASE.EQ.3) GO TO 210 303 EN(IND) = EMX*Y(JSET) 304 IF (M.EQ.1) RETURN 305 AA = KS 306 GO TO (220, 240), ICASE 307C----------------------------------------------------------------------- 308C RECURSION SECTION N*E(N+1,X) + X*E(N,X)=EMX 309C----------------------------------------------------------------------- 310 210 EN(1) = EMX*(1.0E0-Y(1))/X 311 RETURN 312 220 K = IND - 1 313 DO 230 I=1,ML 314 AA = AA - 1.0E0 315 EN(K) = (EMX-AA*EN(K+1))/X 316 K = K - 1 317 230 CONTINUE 318 IF (MU.LE.0) RETURN 319 AA = KS 320 240 K = IND 321 DO 250 I=1,MU 322 EN(K+1) = (EMX-X*EN(K))/AA 323 AA = AA + 1.0E0 324 K = K + 1 325 250 CONTINUE 326 RETURN 327 340 CONTINUE 328 IERR = 2 329 RETURN 330 END 331