#------------------------------------------------------------------------------- # 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]