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 C (C) Copr. 1986-92 Numerical Recipes Software .