ExponentialIntegral.C 4.04 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 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 <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// 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 <cmath>
#include <algorithm>
#include <vector>

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;

}

}