function dexpi(x) c*********************************************************************** c this subroutine calculates the exponential integral quantities. c the exponential integrals are calculates in double precision, c and then truncated to single precision values. c computations of exponential integral with fifteen significant c figure accuracy. c c for negative values of the argument c range 1 greater than -1.0d-20. gamma + log( qabs of x ) c range 2 -1.d-20 to -1.5, 3 point gaussian quadrature. c range 3 -1.5 to -4.65, ratio of 2 polynomials each with 7 terms. c range 4 -4.65 to -12.0, ratio of 2 polynomials each with 6 terms. c range 5 -12.0 to -170.0, 12 point gauss-laguerre quadrature. c c for positive values of the argument c range 1 less than 1.0d-20, gamma + log(x) c range 2 1.0d-20 to 40.0, 12 point gaussian quadrature. c range 3 40.0 to 173.0, 12 point gauss-laguerre quadrature. c c c*********************************************************************** implicit integer*4(i-n),real*8 (a-h,o-z) real * 8 a1(3),b1(3),a2(7),b2(7),a3(6),b3(6),a4(12),b4(12), 1 a5(12),b5(12) data a1 / 0.1193095930415965d+0, 0.3306046932331323d+0, 1 0.4662347571015760d+0/ data b1 / 0.4679139345726910d+0, 0.3607615730461336d+0, 1 0.1713244923791703d+0/ data a2 / 1 0.2823912701457739d-1, 0.3052042817823067d+1, 1 0.2158885931211323d+2, 0.4104611319636983d+2, 2 0.2785527592726121d+2, 0.7133086969436196d+1, 3 0.5758931590224375d+0/ data b2 / 1 0.1000000000000000d+1, 0.1383869728490638d+2, 1 0.4880858183073600d+2, 0.6348804630786363d+2, 2 0.3441289899236299d+2, 0.7708964199043784d+1, 3 0.5758934565014882d+0/ data a3 / 1 0.7630772325814641d-1, 0.2123699219410890d+1, 1 0.4745350785776186d+1, 0.2966421696379266d+1, 2 0.6444800036068992d+0, 0.4295808082119383d-1/ data b3 / 1 0.1000000000000000d+1, 0.5278950049492932d+1, 1 0.7196111390658517d+1, 0.3567945294128107d+1, 2 0.6874380519301884d+0, 0.4295808112146861d-1/ data a4 / 1 0.1157221173580207d+0, 0.6117574845151307d+0, 1 0.1512610269776419d+1, 0.2833751337743507d+1, 2 0.4599227639418348d+1, 0.6844525453115177d+1, 3 0.9621316842456867d+1, 0.1300605499330635d+2, 4 0.1711685518746226d+2, 0.2215109037939701d+2, 5 0.2848796725098400d+2, 0.3709912104446692d+2/ data b4 / 1 0.2647313710554432d+0, 0.3777592758731380d+0, 1 0.2440820113198776d+0, 0.9044922221168093d-1, 2 0.2010238115463410d-1, 0.2663973541865316d-2, 3 0.2032315926629994d-3, 0.8365055856819799d-5, 4 0.1668493876540910d-6, 0.1342391030515004d-8, 5 0.3061601635035021d-11, 0.8148077467426242d-15/ data a5 / 1 0.3202844643130281d-1, 0.9555943373680816d-1, 1 0.1575213398480817d+0, 0.2168967538130226d+0, 2 0.2727107356944198d+0, 0.3240468259684878d+0, 3 0.3700620957892772d+0, 0.4100009929869515d+0, 4 0.4432077635022005d+0, 0.4691372760013664d+0, 5 0.4873642779856547d+0, 0.4975936099985107d+0/ data b5 / 1 0.1279381953467522d+0, 0.1258374563468283d+0, 1 0.1216704729278034d+0, 0.1155056680537256d+0, 2 0.1074442701159656d+0, 0.9761865210411389d-1, 3 0.8619016153195328d-1, 0.7334648141108031d-1, 4 0.5929858491543678d-1, 0.4427743881741981d-1, 5 0.2853138862893366d-1, 0.1234122979998720d-1/ 10 format(10x,'the argument of expi is very close to zero. it is', 1e25.16//) 20 format(10x,'the argument of expi is very large. it is',e25.16//) c c dexpi = 0.0d+00 if (x) 200,100,300 c 100 write (6,10) x 100 continue return 200 ax = dabs(x) if ( x .gt. -1.0d-20 ) go to 201 if ( x .gt. -1.5 ) go to 205 if ( x .gt. -4.65 ) go to 215 if ( x .gt. -12.0 ) go to 225 if ( x .gt. -170.0 ) go to 235 return 201 dexpi = dlog(ax) +0.57721566490153d+00 return 205 yy = dexp(-0.5 * ax) yz = dexp(a1(1)*ax) sumn=(1.0 -yy/yz)/(0.5 +a1(1))+(1.0 -yy*yz)/(0.5 -a1(1)) yz = dexp(a1(2)*ax) sumd=(1.0 -yy/yz)/(0.5 +a1(2))+(1.0 -yy*yz)/(0.5 -a1(2)) yz = dexp(a1(3)*ax) sumt=(1.0 -yy/yz)/(0.5 +a1(3))+(1.0 -yy*yz)/(0.5 -a1(3)) dexpi= -0.5 *(b1(1)*sumn+b1(2)*sumd+b1(3)*sumt) dexpi = dexpi + dlog(ax) + 0.57721566490153d+00 return 215 sumn =(((((a2(7)*ax+a2(6))*ax+a2(5))*ax+a2(4))*ax+a2(3))*ax+a2(2)) 1*ax+a2(1) sumd =(((((b2(7)*ax+b2(6))*ax+b2(5))*ax+b2(4))*ax+b2(3))*ax+b2(2)) 1*ax+b2(1) 218 dexpi= sumn/(sumd*x) dexpi = dexpi * dexp(x) return 225 sumn=((((a3(6)*ax+a3(5))*ax+a3(4))*ax+a3(3))*ax+a3(2))*ax+a3(1) sumd=((((b3(6)*ax+b3(5))*ax+b3(4))*ax+b3(3))*ax+b3(2))*ax+b3(1) go to 218 235 do 238 j = 1,12 238 dexpi=dexpi + b4(j)/(1.0 + a4(j)/ax) dexpi = (dexp(x)/ax)*(-dexpi) return 300 if ( x .le. 1.0d-20 ) go to 301 if ( x .le. 40.0 ) go to 305 if ( x .le. 173.0 ) go to 335 c write (6,20 ) x return 301 dexpi= dlog(x) + 0.57721566490153d+00 return 305 yy = dexp( 0.5 * x) do 310 j = 1,12 yz = dexp(-a5(j)* x ) dexpi =(( 1.0 - yy/yz)/(0.5 + a5(j)) + ( 1.0 - yy*yz)/ 1 ( 0.5 - a5(j)) )*b5(j) + dexpi 310 continue dexpi = -0.5 * dexpi + dlog(x) + 0.57721566490153d+00 return 335 do 338 j = 1,12 338 dexpi = dexpi + b4(j)/(1.0-a4(j)/x) dexpi = (dexp(x)/x) * dexpi return end