1! Dalton, a molecular electronic structure program 2! Copyright (C) by the authors of Dalton. 3! 4! This program is free software; you can redistribute it and/or 5! modify it under the terms of the GNU Lesser General Public 6! License version 2.1 as published by the Free Software Foundation. 7! 8! This program is distributed in the hope that it will be useful, 9! but WITHOUT ANY WARRANTY; without even the implied warranty of 10! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 11! Lesser General Public License for more details. 12! 13! If a copy of the GNU LGPL v2.1 was not distributed with this 14! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 15C 16 17C***************************************************************************** 18C Utilities for J. Toulouse 19C***************************************************************************** 20 21c --------------------------------------------------------------------- 22c 23c Routines for Ei (from Netlib) 24c 25c --------------------------------------------------------------------- 26 27 SUBROUTINE CALCEI(ARG,RESULT,INT) 28C---------------------------------------------------------------------- 29C 30C This Fortran 77 packet computes the exponential integrals Ei(x), 31C E1(x), and exp(-x)*Ei(x) for real arguments x where 32C 33C integral (from t=-infinity to t=x) (exp(t)/t), x > 0, 34C Ei(x) = 35C -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, 36C 37C and where the first integral is a principal value integral. 38C The packet contains three function type subprograms: EI, EONE, 39C and EXPEI; and one subroutine type subprogram: CALCEI. The 40C calling statements for the primary entries are 41C 42C Y = EI(X), where X .NE. 0, 43C 44C Y = EONE(X), where X .GT. 0, 45C and 46C Y = EXPEI(X), where X .NE. 0, 47C 48C and where the entry points correspond to the functions Ei(x), 49C E1(x), and exp(-x)*Ei(x), respectively. The routine CALCEI 50C is intended for internal packet use only, all computations within 51C the packet being concentrated in this routine. The function 52C subprograms invoke CALCEI with the Fortran statement 53C CALL CALCEI(ARG,RESULT,INT) 54C where the parameter usage is as follows 55C 56C Function Parameters for CALCEI 57C Call ARG RESULT INT 58C 59C EI(X) X .NE. 0 Ei(X) 1 60C EONE(X) X .GT. 0 -Ei(-X) 2 61C EXPEI(X) X .NE. 0 exp(-X)*Ei(X) 3 62C 63C The main computation involves evaluation of rational Chebyshev 64C approximations published in Math. Comp. 22, 641-649 (1968), and 65C Math. Comp. 23, 289-303 (1969) by Cody and Thacher. This 66C transportable program is patterned after the machine-dependent 67C FUNPACK packet NATSEI, but cannot match that version for 68C efficiency or accuracy. This version uses rational functions 69C that theoretically approximate the exponential integrals to 70C at least 18 significant decimal digits. The accuracy achieved 71C depends on the arithmetic system, the compiler, the intrinsic 72C functions, and proper selection of the machine-dependent 73C constants. 74C 75C 76C******************************************************************* 77C******************************************************************* 78C 79C Explanation of machine-dependent constants 80C 81C beta = radix for the floating-point system. 82C minexp = smallest representable power of beta. 83C maxexp = smallest power of beta that overflows. 84C XBIG = largest argument acceptable to EONE; solution to 85C equation: 86C exp(-x)/x * (1 + 1/x) = beta ** minexp. 87C XINF = largest positive machine number; approximately 88C beta ** maxexp 89C XMAX = largest argument acceptable to EI; solution to 90C equation: exp(x)/x * (1 + 1/x) = beta ** maxexp. 91C 92C Approximate values for some important machines are: 93C 94C beta minexp maxexp 95C 96C CRAY-1 (S.P.) 2 -8193 8191 97C Cyber 180/185 98C under NOS (S.P.) 2 -975 1070 99C IEEE (IBM/XT, 100C SUN, etc.) (S.P.) 2 -126 128 101C IEEE (IBM/XT, 102C SUN, etc.) (D.P.) 2 -1022 1024 103C IBM 3033 (D.P.) 16 -65 63 104C VAX D-Format (D.P.) 2 -128 127 105C VAX G-Format (D.P.) 2 -1024 1023 106C 107C XBIG XINF XMAX 108C 109C CRAY-1 (S.P.) 5670.31 5.45E+2465 5686.21 110C Cyber 180/185 111C under NOS (S.P.) 669.31 1.26E+322 748.28 112C IEEE (IBM/XT, 113C SUN, etc.) (S.P.) 82.93 3.40E+38 93.24 114C IEEE (IBM/XT, 115C SUN, etc.) (D.P.) 701.84 1.79D+308 716.35 116C IBM 3033 (D.P.) 175.05 7.23D+75 179.85 117C VAX D-Format (D.P.) 84.30 1.70D+38 92.54 118C VAX G-Format (D.P.) 703.22 8.98D+307 715.66 119C 120C******************************************************************* 121C******************************************************************* 122C 123C Error returns 124C 125C The following table shows the types of error that may be 126C encountered in this routine and the function value supplied 127C in each case. 128C 129C Error Argument Function values for 130C Range EI EXPEI EONE 131C 132C UNDERFLOW (-)X .GT. XBIG 0 - 0 133C OVERFLOW X .GE. XMAX XINF - - 134C ILLEGAL X X = 0 -XINF -XINF XINF 135C ILLEGAL X X .LT. 0 - - USE ABS(X) 136C 137C Intrinsic functions required are: 138C 139C ABS, SQRT, EXP 140C 141C 142C Author: W. J. Cody 143C Mathematics abd Computer Science Division 144C Argonne National Laboratory 145C Argonne, IL 60439 146C 147C Latest modification: September 9, 1988 148C 149C---------------------------------------------------------------------- 150 INTEGER I,INT 151CS REAL 152 DOUBLE PRECISION 153 1 A,ARG,B,C,D,EXP40,E,EI,F,FOUR,FOURTY,FRAC,HALF,ONE,P, 154 2 PLG,PX,P037,P1,P2,Q,QLG,QX,Q1,Q2,R,RESULT,S,SIX,SUMP, 155 3 SUMQ,T,THREE,TWELVE,TWO,TWO4,W,X,XBIG,XINF,XMAX,XMX0, 156 4 X0,X01,X02,X11,Y,YSQ,ZERO 157 DIMENSION A(7),B(6),C(9),D(9),E(10),F(10),P(10),Q(10),R(10), 158 1 S(9),P1(10),Q1(9),P2(10),Q2(9),PLG(4),QLG(4),PX(10),QX(10) 159C---------------------------------------------------------------------- 160C Mathematical constants 161C EXP40 = exp(40) 162C X0 = zero of Ei 163C X01/X11 + X02 = zero of Ei to extra precision 164C---------------------------------------------------------------------- 165CS DATA ZERO,P037,HALF,ONE,TWO/0.0E0,0.037E0,0.5E0,1.0E0,2.0E0/, 166CS 1 THREE,FOUR,SIX,TWELVE,TWO4/3.0E0,4.0E0,6.0E0,12.E0,24.0E0/, 167CS 2 FOURTY,EXP40/40.0E0,2.3538526683701998541E17/, 168CS 3 X01,X11,X02/381.5E0,1024.0E0,-5.1182968633365538008E-5/, 169CS 4 X0/3.7250741078136663466E-1/ 170 DATA ZERO,P037,HALF,ONE,TWO/0.0D0,0.037D0,0.5D0,1.0D0,2.0D0/, 171 1 THREE,FOUR,SIX,TWELVE,TWO4/3.0D0,4.0D0,6.0D0,12.D0,24.0D0/, 172 2 FOURTY,EXP40/40.0D0,2.3538526683701998541D17/, 173 3 X01,X11,X02/381.5D0,1024.0D0,-5.1182968633365538008D-5/, 174 4 X0/3.7250741078136663466D-1/ 175C---------------------------------------------------------------------- 176C Machine-dependent constants 177C---------------------------------------------------------------------- 178CS DATA XINF/3.40E+38/,XMAX/93.246E0/,XBIG/82.93E0/ 179 DATA XINF/1.79D+308/,XMAX/716.351D0/,XBIG/701.84D0/ 180C---------------------------------------------------------------------- 181C Coefficients for -1.0 <= X < 0.0 182C---------------------------------------------------------------------- 183CS DATA A/1.1669552669734461083368E2, 2.1500672908092918123209E3, 184CS 1 1.5924175980637303639884E4, 8.9904972007457256553251E4, 185CS 2 1.5026059476436982420737E5,-1.4815102102575750838086E5, 186CS 3 5.0196785185439843791020E0/ 187CS DATA B/4.0205465640027706061433E1, 7.5043163907103936624165E2, 188CS 1 8.1258035174768735759855E3, 5.2440529172056355429883E4, 189CS 2 1.8434070063353677359298E5, 2.5666493484897117319268E5/ 190 DATA A/1.1669552669734461083368D2, 2.1500672908092918123209D3, 191 1 1.5924175980637303639884D4, 8.9904972007457256553251D4, 192 2 1.5026059476436982420737D5,-1.4815102102575750838086D5, 193 3 5.0196785185439843791020D0/ 194 DATA B/4.0205465640027706061433D1, 7.5043163907103936624165D2, 195 1 8.1258035174768735759855D3, 5.2440529172056355429883D4, 196 2 1.8434070063353677359298D5, 2.5666493484897117319268D5/ 197C---------------------------------------------------------------------- 198C Coefficients for -4.0 <= X < -1.0 199C---------------------------------------------------------------------- 200CS DATA C/3.828573121022477169108E-1, 1.107326627786831743809E+1, 201CS 1 7.246689782858597021199E+1, 1.700632978311516129328E+2, 202CS 2 1.698106763764238382705E+2, 7.633628843705946890896E+1, 203CS 3 1.487967702840464066613E+1, 9.999989642347613068437E-1, 204CS 4 1.737331760720576030932E-8/ 205CS DATA D/8.258160008564488034698E-2, 4.344836335509282083360E+0, 206CS 1 4.662179610356861756812E+1, 1.775728186717289799677E+2, 207CS 2 2.953136335677908517423E+2, 2.342573504717625153053E+2, 208CS 3 9.021658450529372642314E+1, 1.587964570758947927903E+1, 209CS 4 1.000000000000000000000E+0/ 210 DATA C/3.828573121022477169108D-1, 1.107326627786831743809D+1, 211 1 7.246689782858597021199D+1, 1.700632978311516129328D+2, 212 2 1.698106763764238382705D+2, 7.633628843705946890896D+1, 213 3 1.487967702840464066613D+1, 9.999989642347613068437D-1, 214 4 1.737331760720576030932D-8/ 215 DATA D/8.258160008564488034698D-2, 4.344836335509282083360D+0, 216 1 4.662179610356861756812D+1, 1.775728186717289799677D+2, 217 2 2.953136335677908517423D+2, 2.342573504717625153053D+2, 218 3 9.021658450529372642314D+1, 1.587964570758947927903D+1, 219 4 1.000000000000000000000D+0/ 220C---------------------------------------------------------------------- 221C Coefficients for X < -4.0 222C---------------------------------------------------------------------- 223CS DATA E/1.3276881505637444622987E+2,3.5846198743996904308695E+4, 224CS 1 1.7283375773777593926828E+5,2.6181454937205639647381E+5, 225CS 2 1.7503273087497081314708E+5,5.9346841538837119172356E+4, 226CS 3 1.0816852399095915622498E+4,1.0611777263550331766871E03, 227CS 4 5.2199632588522572481039E+1,9.9999999999999999087819E-1/ 228CS DATA F/3.9147856245556345627078E+4,2.5989762083608489777411E+5, 229CS 1 5.5903756210022864003380E+5,5.4616842050691155735758E+5, 230CS 2 2.7858134710520842139357E+5,7.9231787945279043698718E+4, 231CS 3 1.2842808586627297365998E+4,1.1635769915320848035459E+3, 232CS 4 5.4199632588522559414924E+1,1.0E0/ 233 DATA E/1.3276881505637444622987D+2,3.5846198743996904308695D+4, 234 1 1.7283375773777593926828D+5,2.6181454937205639647381D+5, 235 2 1.7503273087497081314708D+5,5.9346841538837119172356D+4, 236 3 1.0816852399095915622498D+4,1.0611777263550331766871D03, 237 4 5.2199632588522572481039D+1,9.9999999999999999087819D-1/ 238 DATA F/3.9147856245556345627078D+4,2.5989762083608489777411D+5, 239 1 5.5903756210022864003380D+5,5.4616842050691155735758D+5, 240 2 2.7858134710520842139357D+5,7.9231787945279043698718D+4, 241 3 1.2842808586627297365998D+4,1.1635769915320848035459D+3, 242 4 5.4199632588522559414924D+1,1.0D0/ 243C---------------------------------------------------------------------- 244C Coefficients for rational approximation to ln(x/a), |1-x/a| < .1 245C---------------------------------------------------------------------- 246CS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02, 247CS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/ 248CS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02, 249CS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/ 250 DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, 251 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ 252 DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, 253 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ 254C---------------------------------------------------------------------- 255C Coefficients for 0.0 < X < 6.0, 256C ratio of Chebyshev polynomials 257C---------------------------------------------------------------------- 258CS DATA P/-1.2963702602474830028590E01,-1.2831220659262000678155E03, 259CS 1 -1.4287072500197005777376E04,-1.4299841572091610380064E06, 260CS 2 -3.1398660864247265862050E05,-3.5377809694431133484800E08, 261CS 3 3.1984354235237738511048E08,-2.5301823984599019348858E10, 262CS 4 1.2177698136199594677580E10,-2.0829040666802497120940E11/ 263CS DATA Q/ 7.6886718750000000000000E01,-5.5648470543369082846819E03, 264CS 1 1.9418469440759880361415E05,-4.2648434812177161405483E06, 265CS 2 6.4698830956576428587653E07,-7.0108568774215954065376E08, 266CS 3 5.4229617984472955011862E09,-2.8986272696554495342658E10, 267CS 4 9.8900934262481749439886E10,-8.9673749185755048616855E10/ 268 DATA P/-1.2963702602474830028590D01,-1.2831220659262000678155D03, 269 1 -1.4287072500197005777376D04,-1.4299841572091610380064D06, 270 2 -3.1398660864247265862050D05,-3.5377809694431133484800D08, 271 3 3.1984354235237738511048D08,-2.5301823984599019348858D10, 272 4 1.2177698136199594677580D10,-2.0829040666802497120940D11/ 273 DATA Q/ 7.6886718750000000000000D01,-5.5648470543369082846819D03, 274 1 1.9418469440759880361415D05,-4.2648434812177161405483D06, 275 2 6.4698830956576428587653D07,-7.0108568774215954065376D08, 276 3 5.4229617984472955011862D09,-2.8986272696554495342658D10, 277 4 9.8900934262481749439886D10,-8.9673749185755048616855D10/ 278C---------------------------------------------------------------------- 279C J-fraction coefficients for 6.0 <= X < 12.0 280C---------------------------------------------------------------------- 281CS DATA R/-2.645677793077147237806E00,-2.378372882815725244124E00, 282CS 1 -2.421106956980653511550E01, 1.052976392459015155422E01, 283CS 2 1.945603779539281810439E01,-3.015761863840593359165E01, 284CS 3 1.120011024227297451523E01,-3.988850730390541057912E00, 285CS 4 9.565134591978630774217E00, 9.981193787537396413219E-1/ 286CS DATA S/ 1.598517957704779356479E-4, 4.644185932583286942650E00, 287CS 1 3.697412299772985940785E02,-8.791401054875438925029E00, 288CS 2 7.608194509086645763123E02, 2.852397548119248700147E01, 289CS 3 4.731097187816050252967E02,-2.369210235636181001661E02, 290CS 4 1.249884822712447891440E00/ 291 DATA R/-2.645677793077147237806D00,-2.378372882815725244124D00, 292 1 -2.421106956980653511550D01, 1.052976392459015155422D01, 293 2 1.945603779539281810439D01,-3.015761863840593359165D01, 294 3 1.120011024227297451523D01,-3.988850730390541057912D00, 295 4 9.565134591978630774217D00, 9.981193787537396413219D-1/ 296 DATA S/ 1.598517957704779356479D-4, 4.644185932583286942650D00, 297 1 3.697412299772985940785D02,-8.791401054875438925029D00, 298 2 7.608194509086645763123D02, 2.852397548119248700147D01, 299 3 4.731097187816050252967D02,-2.369210235636181001661D02, 300 4 1.249884822712447891440D00/ 301C---------------------------------------------------------------------- 302C J-fraction coefficients for 12.0 <= X < 24.0 303C---------------------------------------------------------------------- 304CS DATA P1/-1.647721172463463140042E00,-1.860092121726437582253E01, 305CS 1 -1.000641913989284829961E01,-2.105740799548040450394E01, 306CS 2 -9.134835699998742552432E-1,-3.323612579343962284333E01, 307CS 3 2.495487730402059440626E01, 2.652575818452799819855E01, 308CS 4 -1.845086232391278674524E00, 9.999933106160568739091E-1/ 309CS DATA Q1/ 9.792403599217290296840E01, 6.403800405352415551324E01, 310CS 1 5.994932325667407355255E01, 2.538819315630708031713E02, 311CS 2 4.429413178337928401161E01, 1.192832423968601006985E03, 312CS 3 1.991004470817742470726E02,-1.093556195391091143924E01, 313CS 4 1.001533852045342697818E00/ 314 DATA P1/-1.647721172463463140042D00,-1.860092121726437582253D01, 315 1 -1.000641913989284829961D01,-2.105740799548040450394D01, 316 2 -9.134835699998742552432D-1,-3.323612579343962284333D01, 317 3 2.495487730402059440626D01, 2.652575818452799819855D01, 318 4 -1.845086232391278674524D00, 9.999933106160568739091D-1/ 319 DATA Q1/ 9.792403599217290296840D01, 6.403800405352415551324D01, 320 1 5.994932325667407355255D01, 2.538819315630708031713D02, 321 2 4.429413178337928401161D01, 1.192832423968601006985D03, 322 3 1.991004470817742470726D02,-1.093556195391091143924D01, 323 4 1.001533852045342697818D00/ 324C---------------------------------------------------------------------- 325C J-fraction coefficients for X .GE. 24.0 326C---------------------------------------------------------------------- 327CS DATA P2/ 1.75338801265465972390E02,-2.23127670777632409550E02, 328CS 1 -1.81949664929868906455E01,-2.79798528624305389340E01, 329CS 2 -7.63147701620253630855E00,-1.52856623636929636839E01, 330CS 3 -7.06810977895029358836E00,-5.00006640413131002475E00, 331CS 4 -3.00000000320981265753E00, 1.00000000000000485503E00/ 332CS DATA Q2/ 3.97845977167414720840E04, 3.97277109100414518365E00, 333CS 1 1.37790390235747998793E02, 1.17179220502086455287E02, 334CS 2 7.04831847180424675988E01,-1.20187763547154743238E01, 335CS 3 -7.99243595776339741065E00,-2.99999894040324959612E00, 336CS 4 1.99999999999048104167E00/ 337 DATA P2/ 1.75338801265465972390D02,-2.23127670777632409550D02, 338 1 -1.81949664929868906455D01,-2.79798528624305389340D01, 339 2 -7.63147701620253630855D00,-1.52856623636929636839D01, 340 3 -7.06810977895029358836D00,-5.00006640413131002475D00, 341 4 -3.00000000320981265753D00, 1.00000000000000485503D00/ 342 DATA Q2/ 3.97845977167414720840D04, 3.97277109100414518365D00, 343 1 1.37790390235747998793D02, 1.17179220502086455287D02, 344 2 7.04831847180424675988D01,-1.20187763547154743238D01, 345 3 -7.99243595776339741065D00,-2.99999894040324959612D00, 346 4 1.99999999999048104167D00/ 347C---------------------------------------------------------------------- 348 X = ARG 349 IF (X .EQ. ZERO) THEN 350 EI = -XINF 351 IF (INT .EQ. 2) EI = -EI 352 ELSE IF ((X .LT. ZERO) .OR. (INT .EQ. 2)) THEN 353C---------------------------------------------------------------------- 354C Calculate EI for negative argument or for E1. 355C---------------------------------------------------------------------- 356 Y = ABS(X) 357 IF (Y .LE. ONE) THEN 358 SUMP = A(7) * Y + A(1) 359 SUMQ = Y + B(1) 360 DO 110 I = 2, 6 361 SUMP = SUMP * Y + A(I) 362 SUMQ = SUMQ * Y + B(I) 363 110 CONTINUE 364 EI = LOG(Y) - SUMP / SUMQ 365 IF (INT .EQ. 3) EI = EI * EXP(Y) 366 ELSE IF (Y .LE. FOUR) THEN 367 W = ONE / Y 368 SUMP = C(1) 369 SUMQ = D(1) 370 DO 130 I = 2, 9 371 SUMP = SUMP * W + C(I) 372 SUMQ = SUMQ * W + D(I) 373 130 CONTINUE 374 EI = - SUMP / SUMQ 375 IF (INT .NE. 3) EI = EI * EXP(-Y) 376 ELSE 377 IF ((Y .GT. XBIG) .AND. (INT .LT. 3)) THEN 378 EI = ZERO 379 ELSE 380 W = ONE / Y 381 SUMP = E(1) 382 SUMQ = F(1) 383 DO 150 I = 2, 10 384 SUMP = SUMP * W + E(I) 385 SUMQ = SUMQ * W + F(I) 386 150 CONTINUE 387 EI = -W * (ONE - W * SUMP / SUMQ ) 388 IF (INT .NE. 3) EI = EI * EXP(-Y) 389 END IF 390 END IF 391 IF (INT .EQ. 2) EI = -EI 392 ELSE IF (X .LT. SIX) THEN 393C---------------------------------------------------------------------- 394C To improve conditioning, rational approximations are expressed 395C in terms of Chebyshev polynomials for 0 <= X < 6, and in 396C continued fraction form for larger X. 397C---------------------------------------------------------------------- 398 T = X + X 399 T = T / THREE - TWO 400 PX(1) = ZERO 401 QX(1) = ZERO 402 PX(2) = P(1) 403 QX(2) = Q(1) 404 DO 210 I = 2, 9 405 PX(I+1) = T * PX(I) - PX(I-1) + P(I) 406 QX(I+1) = T * QX(I) - QX(I-1) + Q(I) 407 210 CONTINUE 408 SUMP = HALF * T * PX(10) - PX(9) + P(10) 409 SUMQ = HALF * T * QX(10) - QX(9) + Q(10) 410 FRAC = SUMP / SUMQ 411 XMX0 = (X - X01/X11) - X02 412 IF (ABS(XMX0) .GE. P037) THEN 413 EI = LOG(X/X0) + XMX0 * FRAC 414 IF (INT .EQ. 3) EI = EXP(-X) * EI 415 ELSE 416C---------------------------------------------------------------------- 417C Special approximation to ln(X/X0) for X close to X0 418C---------------------------------------------------------------------- 419 Y = XMX0 / (X + X0) 420 YSQ = Y*Y 421 SUMP = PLG(1) 422 SUMQ = YSQ + QLG(1) 423 DO 220 I = 2, 4 424 SUMP = SUMP*YSQ + PLG(I) 425 SUMQ = SUMQ*YSQ + QLG(I) 426 220 CONTINUE 427 EI = (SUMP / (SUMQ*(X+X0)) + FRAC) * XMX0 428 IF (INT .EQ. 3) EI = EXP(-X) * EI 429 END IF 430 ELSE IF (X .LT. TWELVE) THEN 431 FRAC = ZERO 432 DO 230 I = 1, 9 433 FRAC = S(I) / (R(I) + X + FRAC) 434 230 CONTINUE 435 EI = (R(10) + FRAC) / X 436 IF (INT .NE. 3) EI = EI * EXP(X) 437 ELSE IF (X .LE. TWO4) THEN 438 FRAC = ZERO 439 DO 240 I = 1, 9 440 FRAC = Q1(I) / (P1(I) + X + FRAC) 441 240 CONTINUE 442 EI = (P1(10) + FRAC) / X 443 IF (INT .NE. 3) EI = EI * EXP(X) 444 ELSE 445 IF ((X .GE. XMAX) .AND. (INT .LT. 3)) THEN 446 EI = XINF 447 ELSE 448 Y = ONE / X 449 FRAC = ZERO 450 DO 250 I = 1, 9 451 FRAC = Q2(I) / (P2(I) + X + FRAC) 452 250 CONTINUE 453 FRAC = P2(10) + FRAC 454 EI = Y + Y * Y * FRAC 455 IF (INT .NE. 3) THEN 456 IF (X .LE. XMAX-TWO4) THEN 457 EI = EI * EXP(X) 458 ELSE 459C---------------------------------------------------------------------- 460C Calculation reformulated to avoid premature overflow 461C---------------------------------------------------------------------- 462 EI = (EI * EXP(X-FOURTY)) * EXP40 463 END IF 464 END IF 465 END IF 466 END IF 467 RESULT = EI 468 RETURN 469C---------- Last line of CALCEI ---------- 470 END 471 FUNCTION EI(X) 472C-------------------------------------------------------------------- 473C 474C This function program computes approximate values for the 475C exponential integral Ei(x), where x is real. 476C 477C Author: W. J. Cody 478C 479C Latest modification: January 12, 1988 480C 481C-------------------------------------------------------------------- 482 INTEGER INT 483CS REAL EI, X, RESULT 484 DOUBLE PRECISION EI, X, RESULT 485C-------------------------------------------------------------------- 486 INT = 1 487 CALL CALCEI(X,RESULT,INT) 488 EI = RESULT 489 RETURN 490C---------- Last line of EI ---------- 491 END 492 FUNCTION fEXPEI(X) 493C-------------------------------------------------------------------- 494C 495C This function program computes approximate values for the 496C function exp(-x) * Ei(x), where Ei(x) is the exponential 497C integral, and x is real. 498C 499C Author: W. J. Cody 500C 501C Latest modification: January 12, 1988 502C 503C-------------------------------------------------------------------- 504 INTEGER INT 505CS REAL EXPEI, X, RESULT 506 Real*8 fEXPEI, X, RESULT 507C-------------------------------------------------------------------- 508 INT = 3 509 CALL CALCEI(X,RESULT,INT) 510 fEXPEI = RESULT 511 RETURN 512C---------- Last line of EXPEI ---------- 513 END 514 FUNCTION EONE(X) 515C-------------------------------------------------------------------- 516C 517C This function program computes approximate values for the 518C exponential integral E1(x), where x is real. 519C 520C Author: W. J. Cody 521C 522C Latest modification: January 12, 1988 523C 524C-------------------------------------------------------------------- 525 INTEGER INT 526CS REAL EONE, X, RESULT 527 DOUBLE PRECISION EONE, X, RESULT 528C-------------------------------------------------------------------- 529 INT = 2 530 CALL CALCEI(X,RESULT,INT) 531 EONE = RESULT 532 RETURN 533C---------- Last line of EONE ---------- 534 END 535 536*===================================================================== * 537 subroutine CUT_INTO_FIELDS ( 538 . string, !. input 539 . charac, !. input 540 . maxfld, !. input 541 . fields, !. output 542 . lenfld, !. output 543 . NF) !. output 544* -------------------------------------------------------------------- * 545* 546* Comments: divides a STRING into NF fields delimited by charac 547* 548* author:F.Colonna 549* date :21 May 92 550* -------------------------------------------------------------------- * 551 implicit none 552* includes: 553 554* input : 555 integer maxfld 556 character*(*) string 557 character*(1) charac 558* output arrays: 559 character*(*) fields(maxfld) 560 integer lenfld(maxfld) 561* output scalars: 562 integer NF 563* local scalars: 564 integer lenstr, indf, indl, cursor, ind2, i 565 integer lenwrd, fld 566 character*(250) temp, copy 567 character*(1) pad 568 character*(2) char2 569 logical found, max, maxl 570 571! fields = "" 572 call strim (string, lenstr) 573 574* write(6,'(3a)') 575* .' debug_cutflg> input string:"',string(1:lenstr),'"', 576* .' delimiter character:"',charac,'"' 577* write(6,'(a,i2)') 578* .' maximun field number:"',maxfld 579 580 pad = ' ' 581 if(charac.eq.' ') pad ='.' 582 583 if(lenstr.eq.0) then 584* write(6,'(a)')' warning_CUT_INTO_FIELDS> string zero length ' 585 NF = 0 586 return 587 end if 588 589 if(index(string, charac).eq.0) then 590 fields(1)= string 591 NF = 1 592 return 593 end if 594 595* delete double occurences of charac: 596 char2 = charac//charac 597 copy = string 598 ind2 = index(copy(1:lenstr),char2) 599 found = ind2.ne.0 600 do while(found) 601 602 temp = copy (1:ind2)//copy(ind2+2:) 603 copy = temp 604 lenstr = lenstr -1 605 606 ind2 = index(copy(1:lenstr),charac//charac) 607 found = ind2.ne.0 608 end do 609 610* delete occurence in first position: 611 if(copy(1:1).eq.charac) then 612 temp = copy(2:lenstr) 613 copy = temp(1:lenstr-1) 614 lenstr = lenstr-1 615 end if 616 617 NF = 0 618 indf = 0 619 found = .true. 620 max = .false. 621 maxl = .false. 622 temp = copy 623 624!hj do while (found.and..not.max.and..not.maxl) 625 do while (found.and..not.maxl) 626 NF = NF +1 627 indl = index(temp(1:lenstr),charac) 628 found = indl .gt.0 629 630 if(indf+1.le.indl-1) then 631 if (NF .le. maxfld) then 632 fields(NF) = copy(indf+1:indl-1) 633 lenfld(NF) = indl - indf - 1 634 end if 635 temp (indl:indl) = pad 636 indf = indl 637 else if (NF .le. maxfld) then !. end of string 638 fields(NF) = copy(indf+1:) 639 lenfld(NF) = lenstr - indf 640 end if 641 maxl = indl .eq. lenstr 642 max = NF .eq. maxfld 643 end do 644 645* write(*,*) ' found=',found 646* write(*,*) ' maxl=',maxl 647* write(*,*) ' lenstr=',lenstr 648!hj if(max) then 649 if(NF .gt. maxfld) then 650 write(6,'(a,i5,a)') 651 .' warning_cutfld> maximum number of fields =',maxfld,' reached ' 652 write(6,'(3a)') 653 .' string :"',copy(1:lenstr),'"', 654 .' delimiter :"',charac,'"' 655 end if 656* 657* write(6,'(3a)') 658* .' debug_cutflg> output fields ' 659* do fld =1, NF 660* write(6,'(a,i3,3a)')' field #',fld,' "',fields(fld),'"' 661* end do 662 663* end of routine CUT_INTO_FIELDS 664 return 665 end 666 667* ==================================================================== * 668 subroutine STRIM (String, strlen) 669* ---------------------------------------------------------------------- 670* 671* Comments: cuts out ending blanks, returns new length 672* 673* ---------------------------------------------------------------------- 674 implicit none 675 character*(*) string 676 integer strlen 677* begin: 678 strlen=len(string) 679 if(strlen.gt.0) then 680 10 continue 681 if(strlen.ge.1.and.string(strlen:strlen).eq.' ') then 682 strlen=strlen-1 683 go to 10 684 end if 685 else 686 write(6,'(9x,a)')' error_STRIM> zero length string' 687 end if 688 return 689 end 690