1 SUBROUTINE calcei(ARG,RESULT,INT) 2C---------------------------------------------------------------------- 3C 4C This Fortran 77 packet computes the exponential integrals Ei(x), 5C E1(x), and exp(-x)*Ei(x) for real arguments x where 6C 7C integral (from t=-infinity to t=x) (exp(t)/t), x > 0, 8C Ei(x) = 9C -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, 10C 11C and where the first integral is a principal value integral. 12C The packet contains three function type subprograms: EI, EONE, 13C and EXPEI; and one subroutine type subprogram: CALCEI. The 14C calling statements for the primary entries are 15C 16C Y = EI(X), where X .NE. 0, 17C 18C Y = EONE(X), where X .GT. 0, 19C and 20C Y = EXPEI(X), where X .NE. 0, 21C 22C and where the entry points correspond to the functions Ei(x), 23C E1(x), and exp(-x)*Ei(x), respectively. The routine CALCEI 24C is intended for internal packet use only, all computations within 25C the packet being concentrated in this routine. The function 26C subprograms invoke CALCEI with the Fortran statement 27C CALL CALCEI(ARG,RESULT,INT) 28C where the parameter usage is as follows 29C 30C Function Parameters for CALCEI 31C Call ARG RESULT INT 32C 33C EI(X) X .NE. 0 Ei(X) 1 34C EONE(X) X .GT. 0 -Ei(-X) 2 35C EXPEI(X) X .NE. 0 exp(-X)*Ei(X) 3 36C 37C The main computation involves evaluation of rational Chebyshev 38C approximations published in Math. Comp. 22, 641-649 (1968), and 39C Math. Comp. 23, 289-303 (1969) by Cody and Thacher. This 40C transportable program is patterned after the machine-dependent 41C FUNPACK packet NATSEI, but cannot match that version for 42C efficiency or accuracy. This version uses rational functions 43C that theoretically approximate the exponential integrals to 44C at least 18 significant decimal digits. The accuracy achieved 45C depends on the arithmetic system, the compiler, the intrinsic 46C functions, and proper selection of the machine-dependent 47C constants. 48C 49C 50C******************************************************************* 51C******************************************************************* 52C 53C Explanation of machine-dependent constants 54C 55C beta = radix for the floating-point system. 56C minexp = smallest representable power of beta. 57C maxexp = smallest power of beta that overflows. 58C XBIG = largest argument acceptable to EONE; solution to 59C equation: 60C exp(-x)/x * (1 + 1/x) = beta ** minexp. 61C XINF = largest positive machine number; approximately 62C beta ** maxexp 63C XMAX = largest argument acceptable to EI; solution to 64C equation: exp(x)/x * (1 + 1/x) = beta ** maxexp. 65C 66C Approximate values for some important machines are: 67C 68C beta minexp maxexp 69C 70C CRAY-1 (S.P.) 2 -8193 8191 71C Cyber 180/185 72C under NOS (S.P.) 2 -975 1070 73C IEEE (IBM/XT, 74C SUN, etc.) (S.P.) 2 -126 128 75C IEEE (IBM/XT, 76C SUN, etc.) (D.P.) 2 -1022 1024 77C IBM 3033 (D.P.) 16 -65 63 78C VAX D-Format (D.P.) 2 -128 127 79C VAX G-Format (D.P.) 2 -1024 1023 80C 81C XBIG XINF XMAX 82C 83C CRAY-1 (S.P.) 5670.31 5.45E+2465 5686.21 84C Cyber 180/185 85C under NOS (S.P.) 669.31 1.26E+322 748.28 86C IEEE (IBM/XT, 87C SUN, etc.) (S.P.) 82.93 3.40E+38 93.24 88C IEEE (IBM/XT, 89C SUN, etc.) (D.P.) 701.84 1.79D+308 716.35 90C IBM 3033 (D.P.) 175.05 7.23D+75 179.85 91C VAX D-Format (D.P.) 84.30 1.70D+38 92.54 92C VAX G-Format (D.P.) 703.22 8.98D+307 715.66 93C 94C******************************************************************* 95C******************************************************************* 96C 97C Error returns 98C 99C The following table shows the types of error that may be 100C encountered in this routine and the function value supplied 101C in each case. 102C 103C Error Argument Function values for 104C Range EI EXPEI EONE 105C 106C UNDERFLOW (-)X .GT. XBIG 0 - 0 107C OVERFLOW X .GE. XMAX XINF - - 108C ILLEGAL X X = 0 -XINF -XINF XINF 109C ILLEGAL X X .LT. 0 - - USE ABS(X) 110C 111C Intrinsic functions required are: 112C 113C ABS, SQRT, EXP 114C 115C 116C Author: W. J. Cody 117C Mathematics abd Computer Science Division 118C Argonne National Laboratory 119C Argonne, IL 60439 120C 121C Latest modification: September 9, 1988 122C 123C---------------------------------------------------------------------- 124 INTEGER I,INT 125CS REAL 126 DOUBLE PRECISION 127 1 A,ARG,B,C,D,EXP40,E,EI,F,FOUR,FOURTY,FRAC,HALF,ONE,P, 128 2 PLG,PX,P037,P1,P2,Q,QLG,QX,Q1,Q2,R,RESULT,S,SIX,SUMP, 129 3 SUMQ,T,THREE,TWELVE,TWO,TWO4,W,X,XBIG,XINF,XMAX,XMX0, 130 4 X0,X01,X02,X11,Y,YSQ,ZERO 131 DIMENSION A(7),B(6),C(9),D(9),E(10),F(10),P(10),Q(10),R(10), 132 1 S(9),P1(10),Q1(9),P2(10),Q2(9),PLG(4),QLG(4),PX(10),QX(10) 133C---------------------------------------------------------------------- 134C Mathematical constants 135C EXP40 = exp(40) 136C X0 = zero of Ei 137C X01/X11 + X02 = zero of Ei to extra precision 138C---------------------------------------------------------------------- 139CS DATA ZERO,P037,HALF,ONE,TWO/0.0E0,0.037E0,0.5E0,1.0E0,2.0E0/, 140CS 1 THREE,FOUR,SIX,TWELVE,TWO4/3.0E0,4.0E0,6.0E0,12.E0,24.0E0/, 141CS 2 FOURTY,EXP40/40.0E0,2.3538526683701998541E17/, 142CS 3 X01,X11,X02/381.5E0,1024.0E0,-5.1182968633365538008E-5/, 143CS 4 X0/3.7250741078136663466E-1/ 144 DATA ZERO,P037,HALF,ONE,TWO/0.0D0,0.037D0,0.5D0,1.0D0,2.0D0/, 145 1 THREE,FOUR,SIX,TWELVE,TWO4/3.0D0,4.0D0,6.0D0,12.D0,24.0D0/, 146 2 FOURTY,EXP40/40.0D0,2.3538526683701998541D17/, 147 3 X01,X11,X02/381.5D0,1024.0D0,-5.1182968633365538008D-5/, 148 4 X0/3.7250741078136663466D-1/ 149C---------------------------------------------------------------------- 150C Machine-dependent constants 151C---------------------------------------------------------------------- 152CS DATA XINF/3.40E+38/,XMAX/93.246E0/,XBIG/82.93E0/ 153 DATA XINF/1.79D+308/,XMAX/716.351D0/,XBIG/701.84D0/ 154C---------------------------------------------------------------------- 155C Coefficients for -1.0 <= X < 0.0 156C---------------------------------------------------------------------- 157CS DATA A/1.1669552669734461083368E2, 2.1500672908092918123209E3, 158CS 1 1.5924175980637303639884E4, 8.9904972007457256553251E4, 159CS 2 1.5026059476436982420737E5,-1.4815102102575750838086E5, 160CS 3 5.0196785185439843791020E0/ 161CS DATA B/4.0205465640027706061433E1, 7.5043163907103936624165E2, 162CS 1 8.1258035174768735759855E3, 5.2440529172056355429883E4, 163CS 2 1.8434070063353677359298E5, 2.5666493484897117319268E5/ 164 DATA A/1.1669552669734461083368D2, 2.1500672908092918123209D3, 165 1 1.5924175980637303639884D4, 8.9904972007457256553251D4, 166 2 1.5026059476436982420737D5,-1.4815102102575750838086D5, 167 3 5.0196785185439843791020D0/ 168 DATA B/4.0205465640027706061433D1, 7.5043163907103936624165D2, 169 1 8.1258035174768735759855D3, 5.2440529172056355429883D4, 170 2 1.8434070063353677359298D5, 2.5666493484897117319268D5/ 171C---------------------------------------------------------------------- 172C Coefficients for -4.0 <= X < -1.0 173C---------------------------------------------------------------------- 174CS DATA C/3.828573121022477169108E-1, 1.107326627786831743809E+1, 175CS 1 7.246689782858597021199E+1, 1.700632978311516129328E+2, 176CS 2 1.698106763764238382705E+2, 7.633628843705946890896E+1, 177CS 3 1.487967702840464066613E+1, 9.999989642347613068437E-1, 178CS 4 1.737331760720576030932E-8/ 179CS DATA D/8.258160008564488034698E-2, 4.344836335509282083360E+0, 180CS 1 4.662179610356861756812E+1, 1.775728186717289799677E+2, 181CS 2 2.953136335677908517423E+2, 2.342573504717625153053E+2, 182CS 3 9.021658450529372642314E+1, 1.587964570758947927903E+1, 183CS 4 1.000000000000000000000E+0/ 184 DATA C/3.828573121022477169108D-1, 1.107326627786831743809D+1, 185 1 7.246689782858597021199D+1, 1.700632978311516129328D+2, 186 2 1.698106763764238382705D+2, 7.633628843705946890896D+1, 187 3 1.487967702840464066613D+1, 9.999989642347613068437D-1, 188 4 1.737331760720576030932D-8/ 189 DATA D/8.258160008564488034698D-2, 4.344836335509282083360D+0, 190 1 4.662179610356861756812D+1, 1.775728186717289799677D+2, 191 2 2.953136335677908517423D+2, 2.342573504717625153053D+2, 192 3 9.021658450529372642314D+1, 1.587964570758947927903D+1, 193 4 1.000000000000000000000D+0/ 194C---------------------------------------------------------------------- 195C Coefficients for X < -4.0 196C---------------------------------------------------------------------- 197CS DATA E/1.3276881505637444622987E+2,3.5846198743996904308695E+4, 198CS 1 1.7283375773777593926828E+5,2.6181454937205639647381E+5, 199CS 2 1.7503273087497081314708E+5,5.9346841538837119172356E+4, 200CS 3 1.0816852399095915622498E+4,1.0611777263550331766871E03, 201CS 4 5.2199632588522572481039E+1,9.9999999999999999087819E-1/ 202CS DATA F/3.9147856245556345627078E+4,2.5989762083608489777411E+5, 203CS 1 5.5903756210022864003380E+5,5.4616842050691155735758E+5, 204CS 2 2.7858134710520842139357E+5,7.9231787945279043698718E+4, 205CS 3 1.2842808586627297365998E+4,1.1635769915320848035459E+3, 206CS 4 5.4199632588522559414924E+1,1.0E0/ 207 DATA E/1.3276881505637444622987D+2,3.5846198743996904308695D+4, 208 1 1.7283375773777593926828D+5,2.6181454937205639647381D+5, 209 2 1.7503273087497081314708D+5,5.9346841538837119172356D+4, 210 3 1.0816852399095915622498D+4,1.0611777263550331766871D03, 211 4 5.2199632588522572481039D+1,9.9999999999999999087819D-1/ 212 DATA F/3.9147856245556345627078D+4,2.5989762083608489777411D+5, 213 1 5.5903756210022864003380D+5,5.4616842050691155735758D+5, 214 2 2.7858134710520842139357D+5,7.9231787945279043698718D+4, 215 3 1.2842808586627297365998D+4,1.1635769915320848035459D+3, 216 4 5.4199632588522559414924D+1,1.0D0/ 217C---------------------------------------------------------------------- 218C Coefficients for rational approximation to ln(x/a), |1-x/a| < .1 219C---------------------------------------------------------------------- 220CS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02, 221CS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/ 222CS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02, 223CS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/ 224 DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, 225 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ 226 DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, 227 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ 228C---------------------------------------------------------------------- 229C Coefficients for 0.0 < X < 6.0, 230C ratio of Chebyshev polynomials 231C---------------------------------------------------------------------- 232CS DATA P/-1.2963702602474830028590E01,-1.2831220659262000678155E03, 233CS 1 -1.4287072500197005777376E04,-1.4299841572091610380064E06, 234CS 2 -3.1398660864247265862050E05,-3.5377809694431133484800E08, 235CS 3 3.1984354235237738511048E08,-2.5301823984599019348858E10, 236CS 4 1.2177698136199594677580E10,-2.0829040666802497120940E11/ 237CS DATA Q/ 7.6886718750000000000000E01,-5.5648470543369082846819E03, 238CS 1 1.9418469440759880361415E05,-4.2648434812177161405483E06, 239CS 2 6.4698830956576428587653E07,-7.0108568774215954065376E08, 240CS 3 5.4229617984472955011862E09,-2.8986272696554495342658E10, 241CS 4 9.8900934262481749439886E10,-8.9673749185755048616855E10/ 242 DATA P/-1.2963702602474830028590D01,-1.2831220659262000678155D03, 243 1 -1.4287072500197005777376D04,-1.4299841572091610380064D06, 244 2 -3.1398660864247265862050D05,-3.5377809694431133484800D08, 245 3 3.1984354235237738511048D08,-2.5301823984599019348858D10, 246 4 1.2177698136199594677580D10,-2.0829040666802497120940D11/ 247 DATA Q/ 7.6886718750000000000000D01,-5.5648470543369082846819D03, 248 1 1.9418469440759880361415D05,-4.2648434812177161405483D06, 249 2 6.4698830956576428587653D07,-7.0108568774215954065376D08, 250 3 5.4229617984472955011862D09,-2.8986272696554495342658D10, 251 4 9.8900934262481749439886D10,-8.9673749185755048616855D10/ 252C---------------------------------------------------------------------- 253C J-fraction coefficients for 6.0 <= X < 12.0 254C---------------------------------------------------------------------- 255CS DATA R/-2.645677793077147237806E00,-2.378372882815725244124E00, 256CS 1 -2.421106956980653511550E01, 1.052976392459015155422E01, 257CS 2 1.945603779539281810439E01,-3.015761863840593359165E01, 258CS 3 1.120011024227297451523E01,-3.988850730390541057912E00, 259CS 4 9.565134591978630774217E00, 9.981193787537396413219E-1/ 260CS DATA S/ 1.598517957704779356479E-4, 4.644185932583286942650E00, 261CS 1 3.697412299772985940785E02,-8.791401054875438925029E00, 262CS 2 7.608194509086645763123E02, 2.852397548119248700147E01, 263CS 3 4.731097187816050252967E02,-2.369210235636181001661E02, 264CS 4 1.249884822712447891440E00/ 265 DATA R/-2.645677793077147237806D00,-2.378372882815725244124D00, 266 1 -2.421106956980653511550D01, 1.052976392459015155422D01, 267 2 1.945603779539281810439D01,-3.015761863840593359165D01, 268 3 1.120011024227297451523D01,-3.988850730390541057912D00, 269 4 9.565134591978630774217D00, 9.981193787537396413219D-1/ 270 DATA S/ 1.598517957704779356479D-4, 4.644185932583286942650D00, 271 1 3.697412299772985940785D02,-8.791401054875438925029D00, 272 2 7.608194509086645763123D02, 2.852397548119248700147D01, 273 3 4.731097187816050252967D02,-2.369210235636181001661D02, 274 4 1.249884822712447891440D00/ 275C---------------------------------------------------------------------- 276C J-fraction coefficients for 12.0 <= X < 24.0 277C---------------------------------------------------------------------- 278CS DATA P1/-1.647721172463463140042E00,-1.860092121726437582253E01, 279CS 1 -1.000641913989284829961E01,-2.105740799548040450394E01, 280CS 2 -9.134835699998742552432E-1,-3.323612579343962284333E01, 281CS 3 2.495487730402059440626E01, 2.652575818452799819855E01, 282CS 4 -1.845086232391278674524E00, 9.999933106160568739091E-1/ 283CS DATA Q1/ 9.792403599217290296840E01, 6.403800405352415551324E01, 284CS 1 5.994932325667407355255E01, 2.538819315630708031713E02, 285CS 2 4.429413178337928401161E01, 1.192832423968601006985E03, 286CS 3 1.991004470817742470726E02,-1.093556195391091143924E01, 287CS 4 1.001533852045342697818E00/ 288 DATA P1/-1.647721172463463140042D00,-1.860092121726437582253D01, 289 1 -1.000641913989284829961D01,-2.105740799548040450394D01, 290 2 -9.134835699998742552432D-1,-3.323612579343962284333D01, 291 3 2.495487730402059440626D01, 2.652575818452799819855D01, 292 4 -1.845086232391278674524D00, 9.999933106160568739091D-1/ 293 DATA Q1/ 9.792403599217290296840D01, 6.403800405352415551324D01, 294 1 5.994932325667407355255D01, 2.538819315630708031713D02, 295 2 4.429413178337928401161D01, 1.192832423968601006985D03, 296 3 1.991004470817742470726D02,-1.093556195391091143924D01, 297 4 1.001533852045342697818D00/ 298C---------------------------------------------------------------------- 299C J-fraction coefficients for X .GE. 24.0 300C---------------------------------------------------------------------- 301CS DATA P2/ 1.75338801265465972390E02,-2.23127670777632409550E02, 302CS 1 -1.81949664929868906455E01,-2.79798528624305389340E01, 303CS 2 -7.63147701620253630855E00,-1.52856623636929636839E01, 304CS 3 -7.06810977895029358836E00,-5.00006640413131002475E00, 305CS 4 -3.00000000320981265753E00, 1.00000000000000485503E00/ 306CS DATA Q2/ 3.97845977167414720840E04, 3.97277109100414518365E00, 307CS 1 1.37790390235747998793E02, 1.17179220502086455287E02, 308CS 2 7.04831847180424675988E01,-1.20187763547154743238E01, 309CS 3 -7.99243595776339741065E00,-2.99999894040324959612E00, 310CS 4 1.99999999999048104167E00/ 311 DATA P2/ 1.75338801265465972390D02,-2.23127670777632409550D02, 312 1 -1.81949664929868906455D01,-2.79798528624305389340D01, 313 2 -7.63147701620253630855D00,-1.52856623636929636839D01, 314 3 -7.06810977895029358836D00,-5.00006640413131002475D00, 315 4 -3.00000000320981265753D00, 1.00000000000000485503D00/ 316 DATA Q2/ 3.97845977167414720840D04, 3.97277109100414518365D00, 317 1 1.37790390235747998793D02, 1.17179220502086455287D02, 318 2 7.04831847180424675988D01,-1.20187763547154743238D01, 319 3 -7.99243595776339741065D00,-2.99999894040324959612D00, 320 4 1.99999999999048104167D00/ 321C---------------------------------------------------------------------- 322 X = ARG 323 IF (X .EQ. ZERO) THEN 324 EI = -XINF 325 IF (INT .EQ. 2) EI = -EI 326 ELSE IF ((X .LT. ZERO) .OR. (INT .EQ. 2)) THEN 327C---------------------------------------------------------------------- 328C Calculate EI for negative argument or for E1. 329C---------------------------------------------------------------------- 330 Y = ABS(X) 331 IF (Y .LE. ONE) THEN 332 SUMP = A(7) * Y + A(1) 333 SUMQ = Y + B(1) 334 DO 110 I = 2, 6 335 SUMP = SUMP * Y + A(I) 336 SUMQ = SUMQ * Y + B(I) 337 110 CONTINUE 338 EI = LOG(Y) - SUMP / SUMQ 339 IF (INT .EQ. 3) EI = EI * EXP(Y) 340 ELSE IF (Y .LE. FOUR) THEN 341 W = ONE / Y 342 SUMP = C(1) 343 SUMQ = D(1) 344 DO 130 I = 2, 9 345 SUMP = SUMP * W + C(I) 346 SUMQ = SUMQ * W + D(I) 347 130 CONTINUE 348 EI = - SUMP / SUMQ 349 IF (INT .NE. 3) EI = EI * EXP(-Y) 350 ELSE 351 IF ((Y .GT. XBIG) .AND. (INT .LT. 3)) THEN 352 EI = ZERO 353 ELSE 354 W = ONE / Y 355 SUMP = E(1) 356 SUMQ = F(1) 357 DO 150 I = 2, 10 358 SUMP = SUMP * W + E(I) 359 SUMQ = SUMQ * W + F(I) 360 150 CONTINUE 361 EI = -W * (ONE - W * SUMP / SUMQ ) 362 IF (INT .NE. 3) EI = EI * EXP(-Y) 363 END IF 364 END IF 365 IF (INT .EQ. 2) EI = -EI 366 ELSE IF (X .LT. SIX) THEN 367C---------------------------------------------------------------------- 368C To improve conditioning, rational approximations are expressed 369C in terms of Chebyshev polynomials for 0 <= X < 6, and in 370C continued fraction form for larger X. 371C---------------------------------------------------------------------- 372 T = X + X 373 T = T / THREE - TWO 374 PX(1) = ZERO 375 QX(1) = ZERO 376 PX(2) = P(1) 377 QX(2) = Q(1) 378 DO 210 I = 2, 9 379 PX(I+1) = T * PX(I) - PX(I-1) + P(I) 380 QX(I+1) = T * QX(I) - QX(I-1) + Q(I) 381 210 CONTINUE 382 SUMP = HALF * T * PX(10) - PX(9) + P(10) 383 SUMQ = HALF * T * QX(10) - QX(9) + Q(10) 384 FRAC = SUMP / SUMQ 385 XMX0 = (X - X01/X11) - X02 386 IF (ABS(XMX0) .GE. P037) THEN 387 EI = LOG(X/X0) + XMX0 * FRAC 388 IF (INT .EQ. 3) EI = EXP(-X) * EI 389 ELSE 390C---------------------------------------------------------------------- 391C Special approximation to ln(X/X0) for X close to X0 392C---------------------------------------------------------------------- 393 Y = XMX0 / (X + X0) 394 YSQ = Y*Y 395 SUMP = PLG(1) 396 SUMQ = YSQ + QLG(1) 397 DO 220 I = 2, 4 398 SUMP = SUMP*YSQ + PLG(I) 399 SUMQ = SUMQ*YSQ + QLG(I) 400 220 CONTINUE 401 EI = (SUMP / (SUMQ*(X+X0)) + FRAC) * XMX0 402 IF (INT .EQ. 3) EI = EXP(-X) * EI 403 END IF 404 ELSE IF (X .LT. TWELVE) THEN 405 FRAC = ZERO 406 DO 230 I = 1, 9 407 FRAC = S(I) / (R(I) + X + FRAC) 408 230 CONTINUE 409 EI = (R(10) + FRAC) / X 410 IF (INT .NE. 3) EI = EI * EXP(X) 411 ELSE IF (X .LE. TWO4) THEN 412 FRAC = ZERO 413 DO 240 I = 1, 9 414 FRAC = Q1(I) / (P1(I) + X + FRAC) 415 240 CONTINUE 416 EI = (P1(10) + FRAC) / X 417 IF (INT .NE. 3) EI = EI * EXP(X) 418 ELSE 419 IF ((X .GE. XMAX) .AND. (INT .LT. 3)) THEN 420 EI = XINF 421 ELSE 422 Y = ONE / X 423 FRAC = ZERO 424 DO 250 I = 1, 9 425 FRAC = Q2(I) / (P2(I) + X + FRAC) 426 250 CONTINUE 427 FRAC = P2(10) + FRAC 428 EI = Y + Y * Y * FRAC 429 IF (INT .NE. 3) THEN 430 IF (X .LE. XMAX-TWO4) THEN 431 EI = EI * EXP(X) 432 ELSE 433C---------------------------------------------------------------------- 434C Calculation reformulated to avoid premature overflow 435C---------------------------------------------------------------------- 436 EI = (EI * EXP(X-FOURTY)) * EXP40 437 END IF 438 END IF 439 END IF 440 END IF 441 RESULT = EI 442 RETURN 443C---------- Last line of CALCEI ---------- 444 END 445 446 447 448 SUBROUTINE einlib(X, RESULT) 449C FUNCTION EINLIB(X) 450C-------------------------------------------------------------------- 451C 452C This function program computes approximate values for the 453C exponential integral Ei(x), where x is real. 454C 455C Author: W. J. Cody 456C 457C Latest modification: January 12, 1988 458C Latest modification: 20130629 by TWY 459C 460C-------------------------------------------------------------------- 461 INTEGER INT 462CS REAL EI 463CS REAL X 464CS REAL RESULT 465 DOUBLE PRECISION X 466CD DOUBLE PRECISION EI 467 DOUBLE PRECISION RESULT 468C-------------------------------------------------------------------- 469 INT = 1 470 CALL calcei(X,RESULT,INT) 471CD EI = RESULT 472 RETURN 473C---------- Last line of EI ---------- 474 END 475 476 477 478 SUBROUTINE expeinl(X, RESULT) 479C FUNCTION EXPEINL(X) 480C-------------------------------------------------------------------- 481C 482C This function program computes approximate values for the 483C function exp(-x) * Ei(x), where Ei(x) is the exponential 484C integral, and x is real. 485C 486C Author: W. J. Cody 487C 488C Latest modification: January 12, 1988 489C Latest modification: 20130629 by TWY 490C 491C-------------------------------------------------------------------- 492 INTEGER INT 493CS REAL EXPEI 494CS REAL X 495CS REAL RESULT 496CD DOUBLE PRECISION EXPEI 497 DOUBLE PRECISION X 498 DOUBLE PRECISION RESULT 499C-------------------------------------------------------------------- 500 INT = 3 501 CALL calcei(X,RESULT,INT) 502CD EXPEI = RESULT 503 RETURN 504C---------- Last line of EXPEI ---------- 505 END 506 507 508 509 SUBROUTINE eonenl(X, RESULT) 510C FUNCTION EONENL(X) 511C-------------------------------------------------------------------- 512C 513C This function program computes approximate values for the 514C exponential integral E1(x), where x is real. 515C 516C Author: W. J. Cody 517C 518C Latest modification: January 12, 1988 519C Latest modification: 20130629 by TWY 520C 521C-------------------------------------------------------------------- 522 INTEGER INT 523CS REAL EONE 524CS REAL X 525CS REAL RESULT 526CD DOUBLE PRECISION EONE 527 DOUBLE PRECISION X 528 DOUBLE PRECISION RESULT 529C-------------------------------------------------------------------- 530 INT = 2 531 CALL calcei(X,RESULT,INT) 532CD EONE = RESULT 533 RETURN 534C---------- Last line of EONE ---------- 535 END 536