LCOV - code coverage report
Current view: top level - hybrid - exponential_integral.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 15 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 3 0.0 %

          Line data    Source code
       1             : ! Calculate the exponential integral using the algorithm of
       2             : ! [1] Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
       3             : module m_exponential_integral
       4             : 
       5             :    implicit none
       6             : 
       7             :    real, parameter :: series_laguerre = 4.0
       8             : 
       9             : contains
      10             : 
      11             :    ! Calculate the exponential integral E_1(x):
      12             :    !
      13             :    !        inf
      14             :    !          /     -t
      15             :    !         |     e
      16             :    ! E (x) = | dt -----
      17             :    !  1      |      t
      18             :    !        /
      19             :    !         x
      20             :    !
      21             :    ! Input:  arg - position at which exponential integral is evaluated (arg > 0)
      22             :    ! Output: res - E_1(arg)
      23           0 :    pure subroutine calculateExponentialIntegral(arg, res)
      24             : 
      25             :       implicit none
      26             : 
      27             :       real, intent(in)  :: arg
      28             :       real, intent(out) :: res
      29             : 
      30             :       ! For arguments smaller than 4 the series expansion is used
      31           0 :       if (arg < series_laguerre) then
      32           0 :          res = seriesExpansion(arg)
      33             : 
      34             :          ! otherwise a Gauss-Laguerre expansion is better
      35             :       else
      36           0 :          res = exp(-arg)*gauss_laguerre(arg)
      37             :       endif
      38             : 
      39           0 :    end subroutine calculateExponentialIntegral
      40             : 
      41             :    ! Series expansion of the exponential integral
      42             :    !
      43             :    !                          n_cut
      44             :    !                          -----     n  n
      45             :    !                           \    (-1)  x
      46             :    ! E (x) = -gamma - ln(x) -   )   --------
      47             :    !  1                        /     n * n!
      48             :    !                          -----
      49             :    !                          n = 1
      50             :    !
      51             :    ! where gamma is the Euler constant.
      52             :    ! n_cut is set to 25
      53             :    ! Input: arg - argument for which the exponential integral is approximated
      54             :    ! Return: approximation by series expansion for E_1(arg)
      55           0 :    pure real function seriesExpansion(arg)
      56             : 
      57             :       implicit none
      58             : 
      59             :       real, intent(in) :: arg
      60             : 
      61             :       real    :: res, fact  ! result of the summation, 1 / n
      62             :       integer :: i          ! counter variable
      63             : 
      64             :       real, parameter :: EULER_GAMMA = 0.57721566490153286060651209008241 ! Euler constant
      65             :       integer, parameter :: ITERATION = 25 ! Cutoff for series expansion
      66             : 
      67             :       ! initialize summation result
      68           0 :       res = 0.0
      69             : 
      70             :       ! perform the summation
      71           0 :       do i = ITERATION, 2, -1
      72             :          ! calculate 1/n
      73           0 :          fact = 1.0/i
      74             :          ! add next term of summation
      75           0 :          res = arg*fact*(fact - res)
      76             :       end do
      77             : 
      78             :       ! calculate the final result
      79           0 :       seriesExpansion = -EULER_GAMMA - log(arg) + arg*(1.0 - res)
      80             : 
      81           0 :    end function seriesExpansion
      82             : 
      83             :    ! The Gauss Laguerre expansion of the exponential integral can be written as
      84             :    !
      85             :    !              N
      86             :    ! E (arg)    -----    a
      87             :    !  1          \        n
      88             :    ! ------- =    )   --------
      89             :    !   -arg      /    x  + arg
      90             :    !  e         -----  n
      91             :    !             n=1
      92             :    !
      93             :    ! where the a_n and x_n are determined by least quadrature and are given in [1]
      94             :    ! Input: arg - point at which Gaussian Laguerre quadrature is calculated
      95             :    ! Return: E_1(arg) in this approximation
      96           0 :    pure real function gauss_laguerre(arg)
      97             : 
      98             :       implicit none
      99             : 
     100             :       real, intent(in) :: arg
     101             : 
     102             :       ! the quadrature constants a_n and x_n from [1]
     103             :       real, parameter :: a(1:15) = (/ &
     104             :                          0.2182348859400869e+00, 0.3422101779228833e+00, 0.2630275779416801e+00, &
     105             :                          0.1264258181059305e+00, 0.4020686492100091e-01, 0.8563877803611838e-02, &
     106             :                          0.1212436147214252e-02, 0.1116743923442519e-03, 0.6459926762022901e-05, &
     107             :                          0.2226316907096273e-06, 0.4227430384979365e-08, 0.3921897267041089e-10, &
     108             :                          0.1456515264073126e-12, 0.1483027051113301e-15, 0.1600594906211133e-19/)
     109             :       real, parameter :: x(1:15) = (/ &
     110             :                          0.9330781201728180e-01, 0.4926917403018839e+00, 0.1215595412070949e+01, &
     111             :                          0.2269949526203743e+01, 0.3667622721751437e+01, 0.5425336627413553e+01, &
     112             :                          0.7565916226613068e+01, 0.1012022856801911e+02, 0.1313028248217572e+02, &
     113             :                          0.1665440770832996e+02, 0.2077647889944877e+02, 0.2562389422672878e+02, &
     114             :                          0.3140751916975394e+02, 0.3853068330648601e+02, 0.4802608557268579e+02/)
     115             : 
     116             :       ! Calculate the summation
     117           0 :       gauss_laguerre = sum(a/(x + arg))
     118             : 
     119           0 :    end function gauss_laguerre
     120             : 
     121             : end module m_exponential_integral

Generated by: LCOV version 1.13