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(inout) :: 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
|
