enxb Subroutine

subroutine enxb ( n, x, en )

  !*****************************************************************************80
  !
  !! ENXB computes the exponential integral En(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:
  !
  !    10 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, integer ( kind = 4 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) EN(0:N), the function values.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) en(0:n)
  real ( kind = 8 ) ens
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  real ( kind = 8 ) ps
  real ( kind = 8 ) r
  real ( kind = 8 ) rp
  real ( kind = 8 ) s
  real ( kind = 8 ) s0
  real ( kind = 8 ) t
  real ( kind = 8 ) t0
  real ( kind = 8 ) x

  if ( x == 0.0D+00 ) then

     en(0) = 1.0D+300
     en(1) = 1.0D+300
     do k = 2, n
        en(k) = 1.0D+00 / ( k - 1.0D+00 )
     end do
     return

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

     en(0) = exp ( - x ) / x
     do l = 1, n
        rp = 1.0D+00
        do j = 1, l - 1
           rp = - rp * x / j
        end do
        ps = -0.5772156649015328D+00
        do m = 1, l - 1
           ps = ps + 1.0D+00 / m
        end do
        ens = rp * ( - log ( x ) + ps )
        s = 0.0D+00
        do m = 0, 20
           if ( m /= l - 1 ) then
              r = 1.0D+00
              do j = 1, m
                 r = - r * x / j
              end do
              s = s + r / ( m - l + 1.0D+00 )
              if ( abs ( s - s0 ) < abs ( s ) * 1.0D-15 ) then
                 exit
              end if
              s0 = s
           end if
        end do

        en(l) = ens - s

     end do

  else

     en(0) = exp ( - x ) / x
     m = 15 + int ( 100.0D+00 / x )
     do l = 1, n
        t0 = 0.0D+00
        do k = m, 1, -1
           t0 = ( l + k - 1.0D+00 ) / ( 1.0D+00 + k / ( x + t0 ) )
        end do
        t = 1.0D+00 / ( x + t0 )
        en(l) = exp ( - x ) * t
     end do

  end if

  return
end subroutine enxb