import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import pastas as ps

3. Gamma

A = 5
n = 1.5
a = 50
p = [A, n, a]

gamma = ps.Gamma()
tmax = gamma.get_tmax(p)
t = np.arange(0, tmax)

step = gamma.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(gamma.impulse, 0, t[i], args=(p))[0]

plt.plot(t[1:], step, label='analytic')
plt.plot(t, stepnum, '--', label='numerical')
plt.xlabel('time (d)')
plt.ylabel('step (m)')
plt.grid()
plt.legend();

4. Exponential

A = 5
a = 50
p = [A, a]

exponential = ps.Exponential()
tmax = exponential.get_tmax(p)
t = np.arange(0, tmax)

step = exponential.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(exponential.impulse, 0, t[i], args=(p))[0]

plt.plot(t[1:], step, label='analytic')
plt.plot(t, stepnum, '--', label='numerical')
plt.xlabel('time (d)')
plt.ylabel('step (m)')
plt.grid()
plt.legend();

5. Hantush

A = 5
a = 50
b = 2
p = [A, a, b]

hantush = ps.Hantush()
tmax = hantush.get_tmax(p)
t = np.arange(0, tmax)

step = hantush.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(hantush.impulse, 0, t[i], args=(p))[0]

plt.plot(t[1:], step, label='analytic')
plt.plot(t, stepnum, '--', label='numerical')
plt.xlabel('time (d)')
plt.ylabel('step (m)')
plt.grid()
plt.legend();

6. Polder - not correct yet

A = np.exp(-1)
a = 100
b = 0.25
p = [A, a, b]

polder = ps.Polder()
tmax = polder.get_tmax(p)
t = np.arange(0, tmax)

step = polder.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(polder.impulse, 0, t[i], args=(p))[0]

plt.plot(t[1:], step, label='analytic')
plt.plot(t, stepnum, '--', label='numerical')
plt.xlabel('time (d)')
plt.ylabel('step (m)')
plt.grid()
plt.legend();