<commit_before>"""Infiltration module.

"""
import numpy as np
from scipy.optimize import fsolve

def _green_ampt_cum_eq(F, psi, dtheta, K, t):
    """The Green-Ampt cumulative infiltration equation

    """
    tmp = psi*dtheta

    # np.log(x) computes ln(x)
    return F - tmp*np.log(1 + F/tmp) - K*t

def green_ampt_cum_infiltration(psi, eff_theta, eff_sat, K, t):
    """Compute the Green-Ampt cumulative infiltration

    Compute the potential cumulative infiltration up to time `t`,
    using Green-Ampt.

    Parameters
    ----------
    psi : scalar
        Soil suction head at wetting front.
    eff_theta : scalar
        Effective porosity.
    eff_sat : scalar
        Effective saturation.
    K : scalar
        Saturated hydraulic conductivity.
    t : scalar
        Time since beginning of event

    Returns
    -------
    soln : scalar
        Cumulative infiltration up to time `t`.

    Raises
    ------
    ValueError - If no solution can be found.

    """

    dtheta = (1 - eff_sat)*eff_theta

    F = K*t # initial guess

    soln, infodict, ierr, mesg = fsolve(_green_ampt_cum_eq, F,
                                        args=(psi, dtheta, K, t),
                                        full_output=True)

    if ierr == 1:
        return soln
    else:
        raise ValueError(mesg)

def green_ampt_infiltration_rate(psi, eff_theta, eff_sat, K, t):
    """Compute the Green-Ampt infiltration rate

    Compute the potential infiltration rate at time `t`, using
    Green-Ampt.

    Parameters
    ----------
    psi : scalar
        Soil suction head at wetting front.
    eff_theta : scalar
        Effective porosity.
    eff_sat : scalar
        Effective saturation.
    K : scalar
        Saturated hydraulic conductivity.
    t : scalar
        Time since beginning of event

    Returns
    -------
    ft : scalar
        Infiltration rate at time `t`.

    Raises
    ------
    ValueError - If no solution can be found.

    """

    F = green_ampt_cum_infiltration(psi, eff_theta, eff_sat, K, t)

    dtheta = (1 - eff_sat)*eff_theta

    return K*((psi*dtheta)/F + 1)

def test_basic_green_ampt():
    """Test the Green-Ampt infiltration solution

    This test implements Example 4.3.1 pg 116 from Chow et al. 1988,
    'Applied Hydrology'

    """

    psi = 16.7
    eff_theta = 0.486
    eff_sat = 0.3
    K = 0.65
    t = 1

    cum_infil = green_ampt_cum_infiltration(psi, eff_theta, eff_sat, K, t)

    assert np.allclose(cum_infil, [3.16721373])

    infil_rate = green_ampt_infiltration_rate(psi, eff_theta, eff_sat, K, t)

    assert np.allclose(infil_rate, [1.81596836])<commit_msg>ENH: Implement Green-Ampt with ponding taken into account<commit_after>"""Infiltration module.

"""
import numpy as np
from scipy.optimize import fsolve

def _green_ampt_cum_eq(F_t1, F_t0, psi, dtheta, K, dt):
    """The Green-Ampt cumulative infiltration equation

    """
    tmp = psi*dtheta

    # np.log(x) computes ln(x)
    return F_t1 - F_t0 - tmp*np.log((F_t1 + tmp)/(F_t0 + tmp)) - K*dt

def _green_ampt_infiltration_rate(F, psi, eff_theta, eff_sat, K):
    """Compute the Green-Ampt infiltration rate

    Compute the infiltration rate using the Green-Ampt cumulative
    infiltration.

    Parameters
    ----------
    F : scalar
        The cumulative infiltration for the time-period.
    psi : scalar
        Soil suction head at wetting front.
    eff_theta : scalar
        Effective porosity.
    eff_sat : scalar
        Effective saturation.
    K : scalar
        Saturated hydraulic conductivity.

    Returns
    -------
    ft : scalar
        Infiltration rate for the given `F`.

    Raises
    ------
    ValueError - If `F` is zero or negative.

    """
    if F <= 0:
        raise ValueError('F must be greater than zero.')

    dtheta = (1 - eff_sat)*eff_theta

    return K*((psi*dtheta)/F + 1)

def green_ampt_cum_infiltration(rain, psi, eff_theta,
                                eff_sat, K, dt, F_t0=0.0):
    """Compute the Green-Ampt cumulative infiltration

    Compute the Green-Ampt cumulative infiltration depth for a given
    rainfall rate and time interval of length `dt`.

    Parameters
    ----------
    rain : scalar
        The constant rainfall rate during the current time interval of
        length `dt`
    psi : scalar
        Soil suction head at wetting front.
    eff_theta : scalar
        Effective porosity.
    eff_sat : scalar
        Effective saturation.
    K : scalar
        Saturated hydraulic conductivity.
    dt : scalar
        Length of the time-step
    F_t0 : scalar
        The cumulative infiltration depth at the beginning of the
        current interval. Default value is zero.

    Returns
    -------
    F_t1 : scalar
        Cumulative infiltration at the end of the current time-step.

    Raises
    ------
    ValueError - If no solution can be found.

    """
    if rain == 0:
        # shortcut: cum infiltration during interval must be zero
        F_t1 = F_t0
    else:
        # compute potential infiltration rate at the start of the
        # interval
        if F_t0 == 0:
            # infinite rate - just make f_t0 > rainfall rate
            f_t0 = rain + 100.0
        else:
            f_t0 = _green_ampt_infiltration_rate(F_t0, psi,
                                                 eff_theta, eff_sat, K)

        if f_t0 <= rain:
            # ponding occurs throughout interval
            dtheta = (1 - eff_sat)*eff_theta

            F = K*dt # initial guess
            soln, infodict, ierr, mesg = fsolve(_green_ampt_cum_eq, F,
                                                args=(F_t0, psi,
                                                      dtheta, K, dt),
                                                full_output=True)
            F_t1 = soln[0]

            if ierr != 1:
                raise ValueError(mesg)
        else:
            # check whether ponding occurs during interval
            Fprime = F_t0 + rain*dt
            fprime = _green_ampt_infiltration_rate(Fprime, psi,
                                                   eff_theta, eff_sat, K)

            if fprime > rain:
                # no ponding infiltrate at rain rate
                F_t1 = Fprime
            else:
                # ponding has occurred
                dtheta = (1 - eff_sat)*eff_theta
                Fp = (K*psi*dtheta)/(rain - K)

                dtprime = (Fp - F_t0)/rain

                F = K*dt # initial guess
                dtp = dt - dtprime
                soln, infodict, ierr, mesg = fsolve(_green_ampt_cum_eq, F,
                                                    args=(Fp, psi,
                                                          dtheta, K, dtp),
                                                    full_output=True)
                F_t1 = soln[0]

                if ierr != 1:
                    raise ValueError(mesg)

    return F_t1

def test_basic_green_ampt():
    """Test the Green-Ampt infiltration solution

    This test solves Example 4.3.1 pg 116 from Chow et al. 1988,
    'Applied Hydrology'

    """
    psi = 16.7
    eff_theta = 0.486
    eff_sat = 0.3
    K = 0.65
    dt = 1

    dtheta = (1 - eff_sat)*eff_theta

    F = K*dt # initial guess
    soln, infodict, ierr, mesg = fsolve(_green_ampt_cum_eq, F,
                                        args=(0, psi,
                                              dtheta, K, dt),
                                        full_output=True)
    cum_infil = soln[0]

    assert np.allclose(cum_infil, [3.16721373])
{'fraction_non_alphanumeric': 0.08540783898305085, 'fraction_numerical': 0.01615466101694915, 'mean_word_length': 1.5991052993805919, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1748791', 'n_tokens_mistral': 2456, 'n_tokens_neox': 2193, 'n_words': 684}