exponential_integral.f90 3.77 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 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 120 121
! 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(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

end module m_exponential_integral