ExponentialIntegral.C 4.04 KB
 Martin Schlipf committed Sep 27, 2013 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 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 //////////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2008 The Regents of the University of California // // This file is part of Qbox // // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // //////////////////////////////////////////////////////////////////////////////// // // ExponentialIntegral.C // //////////////////////////////////////////////////////////////////////////////// // // Calculate exponential integral using the algorithm of // Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51 // //////////////////////////////////////////////////////////////////////////////// // // Author: Martin Schlipf (2013) // Contact: martin.schlipf@gmail.com // #include "ExponentialIntegral.h" #include #include #include using namespace std; namespace util { // Calculate the exponential integral E_1(x): // // inf // / -t // | e // E (x) = | dt ----- // 1 | t // / // x // // Input: x - position at which exponential integral is evaluated (x > 0) // Return: E1(x) double E1(const double x) { if (x < series_cutoff) { // use series expansion return E1_series(x); } else { // use gauss_laguerre expression return exp(-x) * gauss_laguerre(x); } } // 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: x - position at which exponential integral is evaluated (x > 0) // Return: approximation by series expansion for E_1(x) double E1_series(const double x) { // Euler constant const double EULER_GAMMA = 0.57721566490153286060651209008241; // Cutoff for series expansion const int itmax = 25; // initialize summation result double res = 0.0; // perform the summation for (int it = itmax; it > 1; it--) { // calculate 1/n const double fact = 1.0 / it; // add next term of summation res = x * fact * (fact - res); } // add everything up return -EULER_GAMMA - log(x) + x * (1.0 - res); } // The Gauss Laguerre expansion of the exponential integral can be written as // // N // E (x0) ----- a // 1 \ n // ------ = ) --------- // -x0 / x + x0 // e ----- n // n=1 // // where the a_n and x_n are determined by least quadrature (see reference) // Input: x0 - point at which Gaussian Laguerre quadrature is calculated // Return: E_1(x0) / exp(-x0) in this approximation double gauss_laguerre(const double x0) { // initialize constants a_n and x_n const double size = 15; const double a[] = { 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 }; const double x[] = { 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 }; // initialize double res = 0.0; // evaluate a_n / ( x_n + x0 ) for ( int i = 0; i < size; i++ ) { res += a[i] / (x[i] + x0); } return res; } }