#-------------------------------------------------------------------------------
# Name: expifun-theis.py
# Purpose: Calculating Theis well function W(u) with two different formulas
# from Srivastava(1998), the standard Python function scipy.special.expn
# and compare with reference values from Huisman(1972)
# Source: Practical Approximations of the Well Function
# Rajesh Srivastava and Amado Guzman-Guzman
# Ground Water Vol.36 No.5 september-october 1998
# Reference values for W(u) from:
# Huisman, 1972. Groundwater Recovery. Macmillan. Table 4.1
# Created: 15-01-2012
# Copyright: (c) Grondwaterformules.nl
#-------------------------------------------------------------------------------
#!/usr/bin/env python
from scipy import *
from scipy.special import expn as ExpInt
gamma = 0.5772 # Euler-Macheroni constant
uref = array([1.0e-15,1.0e-10,1.0e-5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,9.9])
HuismanWu = array([33.9616,22.4486,10.9357,0.2194,0.04890,0.01305,0.003779,0.001148,0.003601,0.0001155,0.00003767,0.00001245,0.000004637])
def accurateWu(u):
# accurate approximation for Wu according to equation 4a and 4b from Srivastava(1998)
# Values of constants
a = array([-0.57722,0.99999,-0.24991,0.05519,-0.00976,0.00108])
b = array([0.2677,8.63476,18.05902,8.57333,1])
c = array([3.95850,21.09965,25.63296,9.57332,1])
# Wu for u<1
u0 = where(u<1.0,u,1) # u=1 is just a dummy to make log() work on all values
asum = 0
for i in range(len(a)):
asum += a[i]*u**i
Wu0 = -log(u) + asum
#Wu for u>=1
u1 = where(u>=1.0,u,1) # u=1 is just a dummy
bsum = 0
csum = 0
for i in range(len(b)):
bsum += b[i]*u**i
csum += c[i]*u**i
Wu1 = 1/(u*exp(u))*bsum/csum
# combine Wu0 and Wu1
Wu = where(u<1.0,Wu0,Wu1)
return Wu
def fastWu(u):
# Fast approximation for Wu according to equation 7a and 7b from Srivastava(1998)
# Wu for u<1
u0 = where(u<1.0,u,1) # u=1 is just a dummy to make ln() work on all values
Wu0 = log(exp(-gamma)/u0) + 0.9653*u - 0.1690*u0**2
#Wu for u>=1
u1 = where(u>=1.0,u,1) # u=1 is just a dummy to make ln() work on all values
Wu1 = 1/(u1*exp(u1))*(u1+0.3575)/(u1+1.280)
# combine Wu0 and Wu1
Wu = where(u<1.0,Wu0,Wu1)
return Wu
if __name__ == '__main__':
u = uref
wu1 = accurateWu(u)
wu2 = fastWu(u)
wu3 = ExpInt(1,u)
print " u Huisman Accurate W(u) Fast W(u) Scipy expn"
for i in range(len(u)):
print '%4.1E'%u[i],'%12.8f'%HuismanWu[i],'%12.8f'%wu1[i],'%12.8f'%wu2[i],'%12.8f'%wu3[i]