57 lines
1.4 KiB
Fortran
57 lines
1.4 KiB
Fortran
FUNCTION expint(n,x)
|
|
INTEGER n,MAXIT
|
|
REAL expint,x,EPS,FPMIN,EULER
|
|
PARAMETER (MAXIT=100,EPS=1.e-7,FPMIN=1.e-30,EULER=.5772156649)
|
|
INTEGER i,ii,nm1
|
|
REAL a,b,c,d,del,fact,h,psi
|
|
nm1=n-1
|
|
if(n.lt.0.or.x.lt.0..or.(x.eq.0..and.(n.eq.0.or.n.eq.1)))then
|
|
pause 'bad arguments in expint'
|
|
else if(n.eq.0)then
|
|
expint=exp(-x)/x
|
|
else if(x.eq.0.)then
|
|
expint=1./nm1
|
|
else if(x.gt.1.)then
|
|
b=x+n
|
|
c=1./FPMIN
|
|
d=1./b
|
|
h=d
|
|
do 11 i=1,MAXIT
|
|
a=-i*(nm1+i)
|
|
b=b+2.
|
|
d=1./(a*d+b)
|
|
c=b+a/c
|
|
del=c*d
|
|
h=h*del
|
|
if(abs(del-1.).lt.EPS)then
|
|
expint=h*exp(-x)
|
|
return
|
|
endif
|
|
11 continue
|
|
pause 'continued fraction failed in expint'
|
|
else
|
|
if(nm1.ne.0)then
|
|
expint=1./nm1
|
|
else
|
|
expint=-log(x)-EULER
|
|
endif
|
|
fact=1.
|
|
do 13 i=1,MAXIT
|
|
fact=-fact*x/i
|
|
if(i.ne.nm1)then
|
|
del=-fact/(i-nm1)
|
|
else
|
|
psi=-EULER
|
|
do 12 ii=1,nm1
|
|
psi=psi+1./ii
|
|
12 continue
|
|
del=fact*(-log(x)+psi)
|
|
endif
|
|
expint=expint+del
|
|
if(abs(del).lt.abs(expint)*EPS) return
|
|
13 continue
|
|
pause 'series failed in expint'
|
|
endif
|
|
return
|
|
END
|