Commit 83a7ea50 by Martin Schlipf

implement calculation of exponential integral

transfer implementation from FLEUR
use algorithm described in Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51

git-svn-id: http://qboxcode.org/svn/qb/branches/hse-dev@1373 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 1a50efab
////////////////////////////////////////////////////////////////////////////////
//
// 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;
}
}
////////////////////////////////////////////////////////////////////////////////
//
// 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.h
//
////////////////////////////////////////////////////////////////////////////////
//
// 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
//
#ifndef EXPONENTIAL_INTEGRAL_H
#define EXPONENTIAL_INTEGRAL_H
namespace util {
// constant at which we switch from series to gauss laguerre expansion
static const double series_cutoff = 4.0;
// 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);
// 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);
// 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);
}
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment