e1xb Subroutine

subroutine e1xb ( x, e1 )

  !*****************************************************************************80
  !
  !! E1XB computes the exponential integral E1(x).
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    06 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) E1, the function value.
  !
  implicit none

  real ( kind = 8 ) e1
  real ( kind = 8 ) ga
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) r
  real ( kind = 8 ) t
  real ( kind = 8 ) t0
  real ( kind = 8 ) x

  if ( x == 0.0D+00 ) then

     e1 = 1.0D+300

  else if ( x <= 1.0D+00 ) then

     e1 = 1.0D+00
     r = 1.0D+00

     do k = 1, 25
        r = -r * k * x / ( k + 1.0D+00 )**2
        e1 = e1 + r
        if ( abs ( r ) <= abs ( e1 ) * 1.0D-15 ) then
           exit
        end if
     end do

     ga = 0.5772156649015328D+00
     e1 = - ga - log ( x ) + x * e1

  else

     m = 20 + int ( 80.0D+00 / x )
     t0 = 0.0D+00
     do k = m, 1, -1
        t0 = k / ( 1.0D+00 + k / ( x + t0 ) )
     end do
     t = 1.0D+00 / ( x + t0 )
     e1 = exp ( -x ) * t

  end if

  return
end subroutine e1xb