```
c-------------------------------------------------------------------------

function expint(ind,x) 2
c
c   ROUTINE:   expint
c   PURPOSE:   compute the exponential integral of x.
c              the exponential integral of index n is given by
c              integral( exp(-tx)/x^n dx), x ranges from 1 to infinity
c
c   USEAGE:    result=expint(ind,x,nz)
c
c   INPUT:
c     ind      type of exponential integral, for example use ind=1
c              to get first exponential integral, E1
c     x        argument to exponential integral ( 0 < x < infinity)
c
c   OUTPUT:
c     result   exponential integral
c
dimension a(4),b(4),c(6)

data a/0.2677737343,8.6347608925,18.0590169730,8.5733287401/
data b/3.9584969228,21.0996530827,25.6329561486,9.5733223454/
data c/-0.57721566,0.99999193,-0.24991055,0.05519968,
&       -0.00976004, 0.00107857/
c
if(ind .le. 0) then
write(*,*) 'illegal value of IND in EXPINT'
stop
endif
c
expx=exp(-x)
if (x.le.1.) then
expint=((((c(6)*x+c(5))*x+c(4))*x+c(3))*x+c(2))*x+c(1)-log(x)
else
a1=(((x+a(4))*x+a(3))*x+a(2))*x+a(1)
a2=(((x+b(4))*x+b(3))*x+b(2))*x+b(1)
expint=expx*a1/(a2*x)
endif
c
c use recursion to get higher order exponential integrals
c

if(ind.eq.1) return
do 10 n=1,ind-1
expint=(expx-x*expint)/n
10   continue
return
end
```