def compute_Rn(albedo_c, albedo_s, T_air, T_c, T_s, e_atm, Rs_c, Rs_s, F):
    """Compute Soil and Canopy Net Radiation

    Parameters
    ----------
    albedo_c : ee.Image
    albedo_s : ee.Image
    T_air : ee.Image
        Air temperature (Kelvin).
    T_c : ee.Image
        Canopy temperature (Kelvin).
    T_s : ee.Image
        Soil temperature (Kelvin).
    e_atm : ee.Image
    Rs_c : ee.Image
    Rs_s : ee.Image
    F : ee.Image

    Returns
    -------
    Rn_s : ee.Image
        Soil net radiation (W m-2)
    Rn_c : ee.Image
        Canopy net radiation (W m-2)
    Rn : ee.Image
        Net radiation (W m-2)

    """
    # Long-wave extinction coefficient [-]
    kL = 0.95
    # Soil Emissivity [-]
    eps_s = 0.94
    # Canopy emissivity [-]
    eps_c = 0.99
      
    L_c = T_c.expression(
        'eps_c * 0.0000000567 * (T_c ** 4)', {'eps_c': eps_c, 'T_c': T_c})
    L_s = T_s.expression(
        'eps_s * 0.0000000567 * (T_s ** 4)', {'eps_s': eps_s, 'T_s': T_s})
    Rle = T_air.expression(
        'e_atm * 0.0000000567 * (T_air ** 4)',
        {'e_atm': e_atm, 'T_air': T_air})
    Rn_c = albedo_c.expression(
        '((1 - albedo_c) * Rs_c) + '
        '((1 - exp(-kL * F)) * (Rle + Ls - 2 * L_c))',
        {'albedo_c': albedo_c, 'F': F, 'kL': kL, 'L_c': L_c, 'L_s': L_s,
         'Rle': Rle, 'Rs_c': Rs_c})
    Rn_s = albedo_s.expression(
        '((1 - albedo_s) * Rs_s) + '
        '((exp(-kL * F)) * Rle) + ((1 - exp(-kL * F)) * L_c) - L_s',
        {'albedo_s': albedo_s, 'F': F, 'kL': kL, 'Lc': L_c, 'Ls': L_s,
         'Rle': Rle, 'Rs_s': Rs_s})
    Rn = Rn_s.expression(
        'Rn_s + Rn_c', {'Rn_s': Rn_s, 'Rn_c': Rn_c})

    return Rn_s, Rn_c, Rn