View source for Exponential Integral (a.k.a. Theis' Well Function)

; Language: C++, Visual Basic
; Purpose: Calculate the Exponential integral function, <math>\mathbf{E_1}(x)</math>, also known as the Theis Well function <math>W(u)</math>
; 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

===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 <math>10^{-9}</math>
* 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