! Calculate the exponential integral using the algorithm of
! [1] Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
module m_exponential_integral

   implicit none

   real, parameter :: series_laguerre = 4.0

contains

   ! Calculate the exponential integral E_1(x):
   !
   !        inf
   !          /     -t
   !         |     e
   ! E (x) = | dt -----
   !  1      |      t
   !        /
   !         x
   !
   ! Input:  arg - position at which exponential integral is evaluated (arg > 0)
   ! Output: res - E_1(arg)
   pure subroutine calculateExponentialIntegral(arg, res)

      implicit none

      real, intent(in)  :: arg
      real, intent(inout) :: res

      ! For arguments smaller than 4 the series expansion is used
      if (arg < series_laguerre) then
         res = seriesExpansion(arg)

         ! otherwise a Gauss-Laguerre expansion is better
      else
         res = exp(-arg)*gauss_laguerre(arg)
      endif

   end subroutine calculateExponentialIntegral

   ! Series expansion of the exponential integral
   !
   !                          n_cut
   !                          -----     n  n
   !                           \    (-1)  x
   ! E (x) = -gamma - ln(x) -   )   --------
   !  1                        /     n * n!
   !                          -----
   !                          n = 1
   !
   ! where gamma is the Euler constant.
   ! n_cut is set to 25
   ! Input: arg - argument for which the exponential integral is approximated
   ! Return: approximation by series expansion for E_1(arg)
   pure real function seriesExpansion(arg)

      implicit none

      real, intent(in) :: arg

      real    :: res, fact  ! result of the summation, 1 / n
      integer :: i          ! counter variable

      real, parameter :: EULER_GAMMA = 0.57721566490153286060651209008241 ! Euler constant
      integer, parameter :: ITERATION = 25 ! Cutoff for series expansion

      ! initialize summation result
      res = 0.0

      ! perform the summation
      do i = ITERATION, 2, -1
         ! calculate 1/n
         fact = 1.0/i
         ! add next term of summation
         res = arg*fact*(fact - res)
      end do

      ! calculate the final result
      seriesExpansion = -EULER_GAMMA - log(arg) + arg*(1.0 - res)

   end function seriesExpansion

   ! The Gauss Laguerre expansion of the exponential integral can be written as
   !
   !              N
   ! E (arg)    -----    a
   !  1          \        n
   ! ------- =    )   --------
   !   -arg      /    x  + arg
   !  e         -----  n
   !             n=1
   !
   ! where the a_n and x_n are determined by least quadrature and are given in [1]
   ! Input: arg - point at which Gaussian Laguerre quadrature is calculated
   ! Return: E_1(arg) in this approximation
   pure real function gauss_laguerre(arg)

      implicit none

      real, intent(in) :: arg

      ! the quadrature constants a_n and x_n from [1]
      real, parameter :: a(1:15) = (/ &
                         0.2182348859400869e+00, 0.3422101779228833e+00, 0.2630275779416801e+00, &
                         0.1264258181059305e+00, 0.4020686492100091e-01, 0.8563877803611838e-02, &
                         0.1212436147214252e-02, 0.1116743923442519e-03, 0.6459926762022901e-05, &
                         0.2226316907096273e-06, 0.4227430384979365e-08, 0.3921897267041089e-10, &
                         0.1456515264073126e-12, 0.1483027051113301e-15, 0.1600594906211133e-19/)
      real, parameter :: x(1:15) = (/ &
                         0.9330781201728180e-01, 0.4926917403018839e+00, 0.1215595412070949e+01, &
                         0.2269949526203743e+01, 0.3667622721751437e+01, 0.5425336627413553e+01, &
                         0.7565916226613068e+01, 0.1012022856801911e+02, 0.1313028248217572e+02, &
                         0.1665440770832996e+02, 0.2077647889944877e+02, 0.2562389422672878e+02, &
                         0.3140751916975394e+02, 0.3853068330648601e+02, 0.4802608557268579e+02/)

      ! Calculate the summation
      gauss_laguerre = sum(a/(x + arg))

   end function gauss_laguerre

end module m_exponential_integral
