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