Add Edelman response function

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

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

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)