<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:

def airmass(func, angle, H_max=86400.0, R_planet=6.371229E6, RI=1.000276):
    r'''Calculates mass of air per square meter in the atmosphere using a 
    provided atmospheric model. The lowest air mass is calculated straight up;
    as the angle is lowered to nearer and nearer the horizon, the air mass
    increases, and can approach 40x or more the minimum airmass.

    .. math::
        m(\gamma) = \int_0^\infty \rho \left\{1 - \left[1 + 2(\text{RI}-1)
        (1-\rho/\rho_0)\right]
        \left[\frac{\cos \gamma}{(1+h/R)}\right]^2\right\}^{-1/2} dH

    Parameters
    ----------
    func : float
        Function which returns the density of the atmosphere as a function of
        elevation
    angle : float
        Degrees above the horizon (90 = straight up), [degrees]
    H_max : float, optional
        Maximum height to compute the integration up to before the contribution
        of density becomes negligible, [m]
    R_planet : float, optional
        The radius of the planet for which the integration is being performed,
        [m]
    RI : float, optional
        The refractive index of the atmosphere (air on earth at 0.7 um as 
        default) assumed a constant, [-]

    Returns
    -------
    m : float
        Mass of air per square meter in the atmosphere, [kg/m^2]

    Notes
    -----
    Numerical integration via SciPy's `quad` is used to perform the 
    calculation.

    Examples
    --------
    >>> airmass(lambda Z : ATMOSPHERE_1976(Z).rho, 90)
    10356.127665863998
    
    References
    ----------
    .. [1] Kasten, Fritz, and Andrew T. Young. "Revised Optical Air Mass Tables
       and Approximation Formula." Applied Optics 28, no. 22 (November 15, 
       1989): 4735-38. https://doi.org/10.1364/AO.28.004735.
    '''
    delta0 = RI - 1.0
    rho0 = func(0.0)
    angle_term = cos(radians(angle)) 
    
    def to_int(Z):
        rho = func(Z)
        t1 = (1.0 + 2.0*delta0*(1.0 - rho/rho0))
        t2 = (angle_term/(1.0 + Z/R_planet))**2
        t3 = (1.0 - t1*t2)**-0.5
        return rho*t3

    from scipy.integrate import quad
    return float(quad(to_int, 0, 86400.0)[0])