Exponential Integral (a.k.a. Theis' Well Function)

From Aemwiki

Language

C++, Visual Basic

Purpose

Calculate the Exponential integral function, , also known as the Theis Well function W(u)

Authors

Theo Olsthoorn (VB), James R. Craig (c++)

C++ Version:

#include <math.h>
const double EI_TOLERANCE=1e-9;
const double EULER=0.57721566490153;
/***********************************************************************
	 Exponential Integral / Well Function
************************************************************************
Exponential Integral, E_1(x)=Real(-Ei(-x))
Author: James R. Craig
translated from the VB version published in T. Olsthoorn, "Do a bit more 
with convolution", Groundwater 46(1), 2008
-----------------------------------------------------------------------*/
double Ei(const double &x)
{ 
 int    n;
 double Ei,term;
 Ei=0.0;
 if (x<15.0) 
 {
   n=1;
   term=-x;
   Ei =-EULER-log(x)-term;
   while(fabs(term)>EI_TOLERANCE)
   {
     term*=-x*n/(n+1)/(n+1);
     Ei-=term;
     n++;
   }
 }
 return Ei;
} 

Visual Basic Version:

'***********************************************************************
'	 Theis Well Function
'************************************************************************
'Exponential Integral
'from T. Olsthoorn, "Do a bit more with convolution", Groundwater 46(1), 2008
'-----------------------------------------------------------------------*/
Function W (u As Double) As Double
Dim n AS Integer, term As Double
W = 0
If u < 15 Then
 n = 1
 term = -u
 W =-0.5772156649 - Log(u) - term
 Do While Abs(term) > 0.000000001
  term = -u * n/(n + 1) ^ 2 * term
  W = W- term
  n = n + 1
  Loop
 End If
End Function

[edit] Notes

• This algorithm implements eq. 5.1.11 from Abramowitz and Stegun's handbook of mathematical functions
• The approximation used yields an error bound of 10 ^{− 9}
• These implementations do not address some important issues (i.e., x<0 shouldn't be allowed). A more robust implementation is available from Press et. al, Numerical Recipes