39 lines
1.0 KiB
Fortran
39 lines
1.0 KiB
Fortran
subroutine expinx(x,em1)
|
|
c ========================
|
|
c
|
|
c expinx is a subroutine to calculate the value of e1, the exponential
|
|
c integral or em1=x*expo(x)*e1 at the point x. the polynomial
|
|
c expressions that are used come from abromowitz and stegen
|
|
c
|
|
c a modification of Tim Kallman's XSTAR routine
|
|
c
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
if(x.gt.1.) then
|
|
b1=9.5733223454
|
|
b2=25.6329561486
|
|
b3=21.0996530827
|
|
b4=3.9584969228
|
|
c1=8.5733287401
|
|
c2=18.0590169730
|
|
c3=8.6347608925
|
|
c4=0.2677737343
|
|
em1=x**4+c1*x**3+c2*x*x+c3*x+c4
|
|
em1=em1/(x**4+b1*x*x*x+b2*x*x+b3*x+b4)
|
|
else
|
|
a0=-0.57721566
|
|
a1=0.99999193
|
|
a2=-0.24991055
|
|
a3=0.05519968
|
|
a4=-0.00976004
|
|
a5=0.00107857
|
|
if(x.gt.0) then
|
|
e1= a0+a1*x+a2*x*x+a3*x**3+a4*x**4+a5*x**5-log(x)
|
|
else
|
|
e1=-a0+a1*x+a2*x*x+a3*x**3+a4*x**4+a5*x**5-log(-x)
|
|
end if
|
|
em1=e1*x*expo(x)
|
|
end if
|
|
return
|
|
end
|