def volumetric_heat_capacity(se_top=1.0, porosity=0.4):
    """
    Computes the volumetric heat capacity of the soil

    .. math ::
        \\rho c_{p}=10e^{6}\\left[\\left(1-\\phi\\right)^{2}+
        2.5\\phi+4.2\\phi S_{e,top}\\right]

    Parameters
    ----------
    se_top : float
        effective saturation of the topsoil
        :math:`S_{e,top}`
        [-]
    porosity : float
        porosity of the soil
        :math:`\\phi`
        [-]

    Returns
    -------
    vhc : float
        volumetric heat capacity
        :math:`\\rho c_{p}`
        [J m-3 K-1]

    Examples
    --------
    >>> import ETLook.radiation as rad
    >>> rad.volumetric_heat_capacity(se_top=0.4, porosity = 0.5)
    23400000.0
    """
    return ((1 - porosity) ** 2 + 2.5 * porosity + 4.2 * porosity * se_top
        ) * 10 ** 6
