;+
; NAME:
;  expintei
; PURPOSE: (one line)
;  Return the exponential integral, Ei
; DESCRIPTION:
;  Return the exponential integral, Ei
; CATEGORY:
;  Math
; CALLING SEQUENCE:
;  y = expintei(x)
; INPUTS:
;  x - argument, real number > 0
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  Ei(x)
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  x is a scalar
; PROCEDURE:
;  Based on ei.c, Numerical Recipes 2nd Ed, pp 225.
; MODIFICATION HISTORY:
;  Written 2006 Sept, by Leslie Young, SwRI
;-

function expintei, xin

    eulerg = 0.57721566         ; Euler's constant gamma
    maxit = 100                 ; Maximum number of iterations allowed
    fpmin = (machar()).xmin
    eps = 6.0e-8

    x = double(xin)

    y = x

    indx = where(x le 0.d, n)
    if n ne 0 then y[indx] = -!values.f_infinity

    indx = where(x lt fpmin, n)    
    if n ne 0 then y[indx] = alog(x[indx]) + eulerg

    indx = where(x le -alog(eps), n)    
    if n ne 0 then begin
        for i = 0, n-1 do begin
            xx = x[indx[i]]
            k = 1.
            sum = 0.
            fact = 1.
            term = 100.  ; larger than eps*sum to enter loop
            while k lt maxit and term ge eps * sum do begin
                fact = fact * xx/k
                term = fact / k
                sum = sum + term
                k = k + 1.
            end
            y[indx[i]] = sum + alog(xx) + eulerg
        end
    endif
    
    indx = where(x gt -alog(eps), n)    
    if n ne 0 then begin
        for i = 0, n-1 do begin
            xx = x[indx[i]]
            sum = 0.
            term = 1.
            k = 1.
            repeat begin
                prev = term
                term = term * k/xx
                sum = sum + term
                k = k + 1.
            endrep until term ge prev or term lt eps or k gt maxit
            if term lt prev then begin
                sum = sum - term - prev
            end
            y[indx[i]] = exp(xx) * (1.0 + sum) / xx
        endfor
    endif
    
    return, y

end