exponential_integral.f90 4.02 KB
Newer Older
1 2 3 4
! Calculate the exponential integral using the algorithm of
! [1] Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
module m_exponential_integral

Matthias Redies's avatar
Matthias Redies committed
5 6 7
   implicit none

   real, parameter :: series_laguerre = 4.0
8 9 10

contains

Matthias Redies's avatar
Matthias Redies committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
   ! 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(out) :: 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
120 121

end module m_exponential_integral