import scipy.special as sp

import numpy as np
from scipy.special import expi

def W(u):  return  -expi(-u)

def W1(u):
    """Returns Theis' well function axpproximation by numerical intergration
    
    Works only for scalar u
    """
    if not np.isscalar(u):
        raise ValueError("","u must be a scalar")
        
    LOG10INF = 2 # sufficient as exp(-100) is in the order of 1e-50
    y = np.logspace(np.log10(u), LOG10INF, 1000)  # we use thousand intermediate values
    ym = 0.5 * (y[:-1]  + y[1:])
    Dy = np.diff(y)
    w = np.sum( np.exp(-ym) / ym * Dy )
    return w

"""
Explanation: Investigating the character of the Theis well function
Introduction
In the previous section the Theis well function was introduced. The function, which is in fact the function known as exponential integral by mathematicians, proved available in the standard library of Ptyhon module scipy.special. We modified it a little to make it match the Theis well function exactly and gave it the name "W" like it has in groundwater hydrology books. Then we used it in some examples.
In this chapter we will investigate the Theis well funchtion character a more accurately.
Instead of looking for the function in the available library we could have computed the function ourselfs, for instance by numerical integration.
$$ W(u) = \intop_u^{-\infty} \frac {e^{-y}} y dy \approx \sum_0^N \frac {e^{-y_i}} {y_i} \Delta y_i $$
where $y_0 = u_0$ and $N$ has a a sufficiently large value.
End of explanation
"""


U = 4 * 10** -np.arange(11.)   # generates values 4, 4e-1, 4e-2 .. 4e-10
print("{:>10s}   {:>10s}  {:>10s}".format('u   ', 'W(u)','W1(u)  '))
for u in U:
    print("{0:10.1e}   {1:10.4e}   {2:10.4e}".format(u, W(u), W1(u)))

"""
Explanation: Try it out
End of explanation
"""


def W2(u):
    """Returns Theis well function computed as a power series"""
    tol = 1e-5
    w = -0.5772 -np.log(u) + u
    a = u
    for n in range(1, 100):
        a = -a * u * n / (n+1)**2 # new term (next term)
        w += a
        if np.all(a) < tol:
            return w

"""
Explanation: Is seems that our numerical integration is a fair approximation to four significant digits, but not better, even when computed with 1000 steps as we did. So it is relatively easy to create one's own numerically computed value of an analytical expression like the exponential integral
Theis well function as a power series
The theis well function can be expressed also as a power series. This expression has certain advanages as it gives insight into the behavior of its character and allows important simplifications and deductions.
$$ W(u) = -0.5773 - \ln(u) + u - \frac {u^2} {2 . 2!} + \frac {u^3} {3 . 3!} - \frac {u^4} {4 . 4!} + ... $$
This series too can be readily numerially comptuted by first defining a function for it. The sum will be computed in a loop. To prevent having to compute faculties, it is easiest to compute each successive term from the previous one.
So to get from term m to term n+1:
$$ \frac {u^{n+1}} {(n+1) . (n+1)!} = \frac {u^n} { n . n!} \times \frac {u \, n} {(n+1)^2} $$
This series is implemented below.
End of explanation
"""


U = 4.0 * 10** -np.arange(11.)   # generates values 4, 4e-1, 4e-2 .. 4e-10
print("{:>10s}  {:>10s}  {:>10s}  {:>10s}".format('u   ', 'W(u)  ','W1(u)  ', 'W2(u)  '))
for u in U:
    print("{0:10.1e}  {1:10.4e}   {2:10.4e}  {2:10.4e}".format(u, W(u), W1(u), W2(u)))

"""
Explanation: Compare the three methods of computing the well function.
End of explanation
"""


u = np.logspace(-7, 1, 71)

import matplotlib.pylab as plt
fig1= plt.figure()
ax1 = fig1.add_subplot(111)
ax1.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve versus u', yscale='log', xscale='log')
ax1.grid(True)
ax1.plot(u, W(u), 'b', label='-expi(-u)')
#ax1.plot(u, W1(u), 'rx', label='integal') # works only for scalars
ax1.plot(u, W2(u), 'g+', label='power series')
ax1.legend(loc='best')
plt.show()

"""
Explanation: We see that all three methods yiedld the same results.
Next we show the well function as it shown in groundwater hydrology books.
End of explanation
"""


fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
ax2.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve versus 1/u', yscale='log', xscale='log')
ax2.grid(True)
ax2.plot(1/u, W(u))
plt.show()

"""
Explanation: The curve W(u) versus u runs counter intuitively which and is, therefore, confusing. Therefore, it generally presented as W(u) versus 1/u instead as shown below
End of explanation
"""


fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
ax2.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve versus 1/u', yscale='linear', xscale='log')
ax2.grid(True)
ax2.plot(1/u, W(u))
ax2.invert_yaxis()
plt.show()

"""
Explanation: Now W(u) resembles the actual drawdown, which increases with time.
The reason that this is so, becomes clear from the fact that
$$ u = \frac {r^2 S} {4 kD t} $$
and that
$$ \frac 1 u = \frac {4 kDt} {r^2 S} = \frac {4 kD} S \frac t {r^2} $$
which shows that $\frac 1 u$ increases with time, so that the values of $\frac 1 u$ on the $\frac 1 u$ axis are propotional with time and so the drawdown, i.e., the well function $W(u)$ increases with $\frac 1 u$, which is less confusing.
The graph of $W(u)$ versus $\frac 1 u$ is called the Theis type curve. It's vertical axis is proportional to the drawdown and its horizontal axis proportional to time.
The same curve is shown below but now on linear vertical scale and a logarithmic horizontal scale. The vertical scale was reversed (see values on y-axis) to obtain a curve that illustrates the decline of groundwater head with time caused by the extraction. This way of presending is probably least confusing when reading the curve.
End of explanation
"""


U = np.logspace(-2, 0, 21)

Wa = lambda u : -0.5772 - np.log(u)

print("{:>12s}  {:>12s}   {:>12s}   {:>12s}".format('u','W(u)','Wa(u)','1-Wa(u)/W(u)'))
print("{:>12s}  {:>12s}   {:>12s}   {:>12s}".format(' ','    ','     ','the error'))
for u in U:
    print("{:12.3g}  {:12.3g}  {:12.3g}  {:12.1%}".format(u, W(u), Wa(u), 1-Wa(u)/W(u)))

U = np.logspace(-7, 1, 81)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve and its logarithmic approximation', yscale='linear', xscale='log')
ax.grid(True)
ax.plot(1/U, W(U), 'b',  linewidth = 2., label='Theis type curve')
ax.plot(1/U, Wa(U), 'r', linewidth = 0.25, label='log approximation')
ax.invert_yaxis()
plt.legend(loc='best')
plt.show()

"""
Explanation: Logarithmic approximaion of the Theis type curve
We see that after some time, the drawdown is linear when only the time-axis is logarithmic. This suggests that a logarithmic approximation of time-drawdown curve is accurate after some time.
That this is indeed the case can be deduced from the power series description of the type curve:
$$ W(u) = -0.5773 - \ln(u) + u - \frac {u^2} {2 . 2!} + \frac {u^3} {3 . 3!} - \frac {u^4} {4 . 4!} + ... $$
It is clear that all terms to the right of u will be smaller than u when $u<1$. Hence when u is so small that it can be neglected relative to $\ln(u)$, then all the terms to the right of $\ln(u)$ can be neglected. Therefore we have the following spproximation
$$ W(u) \approx -0.5772 -\ln(u) + O(u) $$
for
$$ -\ln(u)>>u \,\,\,\rightarrow \,\,\, \ln(u)<<-u \,\,\, \rightarrow \,\,\, u<<e^{-u} \, \approx  \,1 $$
which is practically the case for $u<0.01$, as can also be seen in the graph for $1/u = 10^2 $. From the graph one may conclude that even for 1/u>10 or u<0.1, the logarithmic type curve is straight and therefore can be accurately computed using a logarithmic approximation of the type curve.
Below the error between the full Theis curve $W(u)$ and the approximation $Wu(u) = -0.5772 - \ln(u)$ are computed and shown. This reveals that at $u=0.01$ the error is 5.4% and at $u=0.001$ it has come down to only 0.2%.
End of explanation
"""


U = np.logspace(-7, 1, 81)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set(xlabel='u', ylabel='W(u)', title='Theis type curve and its logarithmic approximation', yscale='linear', xscale='log')
ax.grid(True)
ax.plot(U, W(U), 'b',  linewidth = 2., label='Theis type curve')
ax.plot(U, Wa(U), 'r', linewidth = 0.25, label='log approximation')
ax.invert_yaxis()
plt.legend(loc='best')
plt.show()

"""
Explanation: Hence, in any practical situation, the logarithmic approximation is accurate enough when $u<0.01$.
The approximatin of the Theis type curve can no be elaborated:
$$ Wa (u) \approx -0.5772 - \ln(u) = \ln(e^{-0.5772}) - \ln(u) = \ln(0.5615) - \ln(u) = \ln \frac {0.5615} {u} $$
Because $u = \frac {r^2 S} {4 kD t}$ we have, with 4\times 0.5615 \approx 2.25
$$ W(u) \approx \ln \frac {2.25 kD t} {r^2 S} $$
and so the drawdown approximation becomes
$$ s \approx \frac Q {4 \pi kD} \ln \frac {2.25 kD t} {r^2 S} $$
The condition u<0.1 can be translated to $\frac {r^2 S} {4 kD t} < 0.1$ or
$$\frac t {r^2} > 2.5 \frac {S} {kD}$$
Radius of influence
The previous logarithmic drawdown type curve versus $1/u$ can be seen an image of the drawdown for a fixed distance and varying time. This is because $1/u$ is proportional to the real time. On the other hand, the drawdown type curve versus u may be regarded as the drawdown at a fixed time for varying distance. This follows from
s versus u is
$$ W(u)\approx \ln \frac {2.25 kD t} { r^2 S} \,\,\,\, versus\,\,\,\, u = \ln \frac {r^2 S} {4 kD t} = 2 \ln \left( \frac {S} {4 kD t} r\right) $$
That is, proportional r on log scale. The plot reveals this:
End of explanation
"""


#       t[min],  s[m]
H30 = [ [0.0, 0.0],
        [0.1, 0.04],
        [0.25, 0.08],
        [0.50, 0.13],
        [0.70, 0.18],
        [1.00, 0.23],
        [1.40, 0.28],
        [1.90, 0.33],
        [2.33, 0.36],
        [2.80, 0.39],
        [3.36, 0.42],
        [4.00, 0.45],
        [5.35, 0.50],
        [6.80, 0.54],
        [8.30, 0.57],
        [8.70, 0.58],
        [10.0, 0.60],
        [13.1, 0.64]]

#    t[min],  s[m]
H90= [[0.0, 0.0],
      [1.5, 0.015],
      [2.0, 0.021],
      [2.16, 0.23],
      [2.66, 0.044],
      [3.00, 0.054],
      [3.50, 0.075],
      [4.00, 0.090],
      [4.33, 0.104],
      [5.50, 0.133],
      [6.0, 0.154],
      [7.5, 0.178],
      [9.0, 0.206],
      [13.0, 0.250],
      [15.0, 0.275],
      [18.0, 0.305],
      [25.0, 0.348],
      [30.0, 0.364]]

#       t[min],  s[m]
H215=[[0.0, 0.0],
      [66.0, 0.089],
      [127., 0.138],
      [185., 0.165],
      [251., 0.186]]

"""
Explanation: This shows that the radius of influence is limited. We can now approximate this radius of influence by saying that the radius is where the appoximated Theis curve, that is the straight red line in the graph intersects the zero drawdown, i.e. $W(u) = 0$.
Hence, for the radius of influence, R, we have
$$ \ln \frac {2.25 kD t} {R^2 S} = 0 $$
impying that
$$ \frac {2.25 kD t } { R^2 S } = 1 $$
$$ R =\sqrt { \frac {2.25 kD t} D} $$
with R the radius of influence. Computing the radius of influence is an easy way to determine how far out the drawdown affects the groundwater heads.
Pumping test
Introduction
Below are the data given that were obtained from a pumping test carried out on the site "Oude Korendijk" south of Rotterdam in the Netherlands (See Kruseman and De Ridder, p56, 59). The piezometers are all open at 20 m below ground surface. The groundwater head is shallow, within a m from ground surface. The first18 m below ground surface consist of clay,peat and clayey fine sand. These layers form a practially impermeable confining unit. Below this, between 18 and25 m below ground surface are 7 m of sand an some gravel, that form the aquifer. Fine sandy and clayey sediments thereunder from the base of the aquifer, which is considered impermeable.
Piezometers wer installed at 30, 90 and 215 m from the well, open at 20 m below ground surface. The well has its screen installed over the whole thickness of the aquifer. We consider the aquifer as confined with no leakage. But we should look with a critical eye that the drawdown curves to verify to what extent this assumption holds true.
The drawdown data for the three piezometers is given below. The first column is time after the start of the pump in minutes; the second column is the drawdown in m.
The well extracts 788 m3/d
The objective of the pumping test is to determine the properties kD and S of the aquifer.
The data:
End of explanation
"""


