How to use the pastas.rfunc.RfuncBase function in pastas

else:
            self.tmax = max(self.get_tmax(p, cutoff), 10 * dt)
            t = np.arange(dt, self.tmax, dt)
        tau = t / cS
        tau1 = tau[tau < rho / 2]
        tau2 = tau[tau >= rho / 2]
        w = (exp1(rho) - k0rho) / (exp1(rho) - exp1(rho / 2))
        F = np.zeros_like(tau)
        F[tau < rho / 2] = w * exp1(rho ** 2 / (4 * tau1)) - (w - 1) * exp1(
            tau1 + rho ** 2 / (4 * tau1))
        F[tau >= rho / 2] = 2 * k0rho - w * exp1(tau2) + (w - 1) * exp1(
            tau2 + rho ** 2 / (4 * tau2))
        return p[0] * F / (2 * k0rho)


class HantushWellModel(RfuncBase):
    """ A special implementation of the Hantush well function for
    multiple wells.

    Note: The parameter r (distance from the well to the observation point)
    is passed as a known value, and is used to scale the response function.
    The optimized parameters are slightly different from the original
    Hantush implementation:
    - A: To get the same A as the original Hantush:
        A_orig = A * 2 * k0(r / lambda) or use the gain() method
    - lab: lambda, the r parameter is passed separately to calculate
        rho = r / lambda internally
    - cS: stays the same
    - r: distance, used to calculate rho, see lab.

    Parameters
    ----------

self.up = -1
        self.meanstress = meanstress
        self.cutoff = cutoff
        self.tmax = 0

    def set_parameters(self, name):
        pass

    def step(self, p, dt=1):
        pass

    def block(self, p, dt=1):
        pass


class Gamma(RfuncBase):
    __doc__ = """
    Gamma response function with 3 parameters A, a, and n.

    step(t) = A * Gammainc(n, t / a)

    %(doc)s
    """ % {'doc': _class_doc}

    def __init__(self, up=True, meanstress=1, cutoff=0.99):
        RfuncBase.__init__(self, up, meanstress, cutoff)
        self.nparam = 3

    def set_parameters(self, name):
        parameters = pd.DataFrame(
            columns=['initial', 'pmin', 'pmax', 'vary', 'name'])
        parameters.loc[name + '_A'] = (

def gain(self, p):
        return p[0]

    def step(self, p, dt=1, cutoff=None):
        if isinstance(dt, np.ndarray):
            t = dt
        else:
            self.tmax = max(self.get_tmax(p, cutoff), 3 * dt)
            t = np.arange(dt, self.tmax, dt)

        s = p[0] * gammainc(p[1], t / p[2])
        return s


class Exponential(RfuncBase):
    """Exponential response function with 2 parameters: A and a.

    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.
    meanstress: float
        mean value of the stress, used to set the initial value such that
        the final step times the mean stress equals 1
    cutoff: float
        percentage after which the step function is cut off. default=0.99.

    Notes
    -----
    .. math::

return parameters

    def gain(self, p):
        return p[0]

    def step(self, p, dt=1, cutoff=None):
        if isinstance(dt, np.ndarray):
            return p[0] * np.ones(len(dt))
        else:
            return p[0] * np.ones(1)

    def block(self, p, dt=1, cutoff=None):
        return p[0] * np.ones(1)


class FourParam(RfuncBase):
    """Four Parameter response function with 4 parameters A, a, b, and n.

    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.
    meanstress: float
        mean value of the stress, used to set the initial value such that
        the final step times the mean stress equals 1
    cutoff: float
        percentage after which the step function is cut off. default=0.99.

    Notes
    -----

def step(self, p, dt=1):
        if isinstance(dt, np.ndarray):
            t = dt
        else:
            self.tmax = -np.log(1.0 / p[1]) * p[1]
            t = np.arange(dt, self.tmax, dt)
        s = self.up * p[0] * (1.0 - np.exp(-t / p[1]))
        return s

    def block(self, p, dt=1):
        s = self.step(p, dt)
        return s[1:] - s[:-1]


class Hantush(RfuncBase):
    """ The Hantush well function
    Parameters are rho = r / lambda and cS

    References
    ----------
    [1] Hantush, M. S., & Jacob, C. E. (1955). Non‐steady radial flow in an
    infinite leaky aquifer. Eos, Transactions American Geophysical Union,
    36(1), 95-100.

    [2] Veling, E. J. M., & Maas, C. (2010). Hantush well function revisited.
    Journal of hydrology, 393(3), 381-388.

    [3] Von Asmuth, J. R., Maas, K., Bakker, M., & Petersen, J. (2008). Modeling
    time series of ground water head fluctuations subjected to multiple
    stresses. Ground Water, 46(1), 30-40.

def gain(self, p):
        return p[0]

    def step(self, p, dt=1, cutoff=0.99):
        if isinstance(dt, np.ndarray):
            t = dt
        else:
            self.tmax = max(self.calc_tmax(p, cutoff), 3 * dt)
            t = np.arange(dt, self.tmax, dt)

        s = p[0] * (1 - ((1 - p[1]) * np.exp(-t / p[2]) +
                         p[1] * np.exp(-t / p[3])))
        return s


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.
    meanstress: float
        mean value of the stress, used to set the initial value such that
        the final step times the mean stress equals 1
    cutoff: float
        percentage after which the step function is cut off. default=0.99.

    Notes
    -----

else:
            self.tmax = max(self.get_tmax(p, cutoff), 3 * dt)
            t = np.arange(dt, self.tmax, dt)
        s = self.polder_function(p[0], p[1] * np.sqrt(t))
        if not self.up:
            s = -s
        return s

    @staticmethod
    def polder_function(x, y):
        s = .5 * np.exp(2 * x) * erfc(x / y + y) + \
            .5 * np.exp(-2 * x) * erfc(x / y - y)
        return s


class One(RfuncBase):
    """Dummy class for Constant. Returns 1

    """
    _name = "One"

    def __init__(self, up=None, meanstress=1, cutoff=0.999):
        RfuncBase.__init__(self, up, meanstress, cutoff)
        self.nparam = 1

    def get_init_parameters(self, name):
        parameters = DataFrame(
            columns=['initial', 'pmin', 'pmax', 'vary', 'name'])
        if self.up:
            parameters.loc[name + '_d'] = (
                self.meanstress, 0, np.nan, True, name)
        elif self.up is False:

return -p[1] * np.log(1 - cutoff)

    def gain(self, p):
        return p[0]

    def step(self, p, dt=1, cutoff=None):
        if isinstance(dt, np.ndarray):
            t = dt
        else:
            self.tmax = max(self.get_tmax(p, cutoff), 10 * dt)
            t = np.arange(dt, self.tmax, dt)
        s = p[0] * (1.0 - np.exp(-t / p[1]))
        return s


class Hantush(RfuncBase):
    """ The Hantush well function.

    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.
    meanstress: float
        mean value of the stress, used to set the initial value such that
        the final step times the mean stress equals 1
    cutoff: float
        percentage after which the step function is cut off. default=0.99.

    Notes
    -----
    The Hantush well function is explained in [1]_, [2]_ and [3]_.

function can be used for testing purposes.

    .. math::
        step(t) = \\frac{A}{quad(t^n*e^{-\\frac{t}{a} - \\frac{b}{t}},0,inf)} *
                            quad(t^n*e^{-\\frac{t}{a} - \\frac{b}{t}},0,t)

    """
    _name = "FourParamQuad"

    def __init__(self, up=True, meanstress=1, cutoff=0.999):
        FourParam.__init__(self, up, meanstress, cutoff)
        self.nparam = 4
        self.quad = True


class DoubleExponential(RfuncBase):
    """Gamma response function with 3 parameters A, a, and n.

    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.
    meanstress: float
        mean value of the stress, used to set the initial value such that
        the final step times the mean stress equals 1
    cutoff: float
        percentage after which the step function is cut off. default=0.99.

    Notes
    -----