diff --git a/pastas_plugins/responses/rfunc.py b/pastas_plugins/responses/rfunc.py index 874420d..f0c7d11 100644 --- a/pastas_plugins/responses/rfunc.py +++ b/pastas_plugins/responses/rfunc.py @@ -4,7 +4,7 @@ from pandas import DataFrame from pastas.rfunc import RfuncBase from pastas.typing import ArrayLike -from scipy.special import exp1 +from scipy.special import erfc, erfcinv, exp1 class Theis(RfuncBase): @@ -123,3 +123,74 @@ def to_dict(self): # "t": self.t, } return data + + +class Edelman(RfuncBase): + """The function of Edelman, describing the propagation of an instantaneous + water level change into an adjacent half-infinite aquifer. + + Parameters + ---------- + up: bool or None, optional + indicates whether a positive stress will cause the head to go up (True, + default) or down (False), if None the head can go both ways. + gain_scale_factor: float, optional + the scale factor is used to set the initial value and the bounds of the gain + parameter, computed as 1 / gain_scale_factor. + cutoff: float, optional + proportion after which the step function is cut off. + + Notes + ----- + The Edelman function is explained in :cite:t:`edelman_over_1947`. + + The impulse response function for this class can be viewed on the Documentation + website or using `latexify` by running the following code in a Jupyter notebook + environment:: + + ps.Edelman.impulse + + """ + + _name = "Edelman" + + def __init__( + self, + cutoff: float = 0.999, + **kwargs, + ) -> None: + RfuncBase.__init__(self, cutoff=cutoff, **kwargs) + self.nparam = 1 + + def get_init_parameters(self, name: str) -> DataFrame: + parameters = DataFrame( + columns=["initial", "pmin", "pmax", "vary", "name", "dist"] + ) + beta_init = 1.0 + parameters.loc[name + "_beta"] = (beta_init, 0, 1000, True, name, "uniform") + return parameters + + def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float: + if cutoff is None: + cutoff = self.cutoff + return 1.0 / (p[0] * erfcinv(cutoff)) ** 2 + + @staticmethod + def gain(p: ArrayLike) -> float: + return 1.0 + + @staticmethod + def impulse(t: ArrayLike, p: ArrayLike) -> ArrayLike: + (a,) = p + return 1 / (np.sqrt(np.pi) * a * t**1.5) * np.exp(-1 / (a**2 * t)) + + def step( + self, + p: ArrayLike, + dt: float = 1.0, + cutoff: Optional[float] = None, + maxtmax: Optional[int] = None, + ) -> ArrayLike: + t = self.get_t(p=p, dt=dt, cutoff=cutoff, maxtmax=maxtmax) + s = erfc(1 / (p[0] * np.sqrt(t))) + return s diff --git a/pastas_plugins/responses/rfunc_utils.py b/pastas_plugins/responses/rfunc_utils.py index e38535f..8450607 100644 --- a/pastas_plugins/responses/rfunc_utils.py +++ b/pastas_plugins/responses/rfunc_utils.py @@ -30,3 +30,25 @@ def kraijenhoff_parameters( A = -N * L**2 / (2 * Ks * D) * (b**2 - (1 / 4)) a = Sy * L**2 / (pi**2 * Ks * D) return A, a, b + + +def exponential_parameters(c: float, Sy: float) -> tuple[float]: + """Get Pastas parameters for an Exponential response for a linear + reservoir system. + + Parameters + ---------- + c : float + The drainage resistance in days. + Sy : float + The specific yield of the aquifer [-]. + + Returns + ------- + tuple[float] + A tuple containing the response parameters A and a. + """ + + A = c + a = c * Sy + return A, a diff --git a/tests/test_responses.py b/tests/test_responses.py index 20be5b2..e1ef4c2 100644 --- a/tests/test_responses.py +++ b/tests/test_responses.py @@ -1,6 +1,7 @@ import numpy as np +from pandas import DataFrame -from pastas_plugins.responses.rfunc import Theis +from pastas_plugins.responses.rfunc import Edelman, Theis def test_theis_init(): @@ -47,3 +48,51 @@ def test_theis_to_dict(): assert data["class"] == "Theis" assert data["cutoff"] == theis.cutoff # assert data["nterms"] == theis.nterms + + +def test_edelman_init(): + cutoff = 0.999 + rfunc = Edelman(cutoff=cutoff) + assert rfunc.cutoff == cutoff + + +def test_edelman_get_init_parameters(): + rfunc = Edelman() + params = rfunc.get_init_parameters("Edelman") + assert isinstance(params, DataFrame) + assert len(params) == 1 + assert params.loc["Edelman_beta"]["initial"] == 1.0 + + +def test_edelman_get_tmax(): + p = np.array([1.0]) + cutoff = 0.999 + rfunc = Edelman(cutoff=cutoff) + tmax = rfunc.get_tmax(p, cutoff=cutoff) + assert isinstance(tmax, float) + + +def test_edelman_gain(): + p = np.array([1.0]) + rfunc = Edelman() + gain = rfunc.gain(p) + assert isinstance(gain, float) + assert gain == 1.0 + + +def test_edelman_impulse(): + t = np.array([1.0, 2.0, 3.0]) + p = np.array([1.0]) + rfunc = Edelman() + impulse = rfunc.impulse(t, p) + assert isinstance(impulse, np.ndarray) + assert len(impulse) == len(t) + + +def test_edelman_step(): + p = np.array([1.0]) + dt = 1.0 + cutoff = 0.999 + rfunc = Edelman(cutoff=cutoff) + step = rfunc.step(p, dt=dt, cutoff=cutoff) + assert isinstance(step, np.ndarray)