<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
 r""" Estimate DNI and DHI from GHI using the Erbs model. The Erbs model [1]_ estimates the diffuse fraction DF from global horizontal irradiance through an empirical relationship between DF and the ratio of GHI to extraterrestrial irradiance, Kt. The function uses the diffuse fraction to compute DHI as .. math:: DHI = DF \times GHI DNI is then estimated as .. math:: DNI = (GHI - DHI)/\cos(Z) where Z is the zenith angle. Parameters ghi: numeric Global horizontal irradiance in W/m^2. zenith: numeric True (not refraction-corrected) zenith angles in decimal degrees. datetime_or_doy : int, float, array, pd.DatetimeIndex Day of year or array of days of year e.g. pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex. min_cos_zenith : numeric, default 0.065 Minimum value of cos(zenith) to allow when calculating global clearness index `kt`. Equivalent to zenith = 86.273 degrees. max_zenith : numeric, default 87 Maximum value of zenith to allow in DNI calculation. DNI will be set to 0 for times with zenith values greater than `max_zenith`. Returns ------- data : OrderedDict or DataFrame Contains the following keys/columns: * ``dni``: the modeled direct normal irradiance in W/m^2. * ``dhi``: the modeled diffuse horizontal irradiance in W/m^2. * ``kt``: Ratio of global to extraterrestrial irradiance on a horizontal plane. References .. [1] D. G. Erbs, S. A. Klein and J. A. Duffie, Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation, Solar Energy 28(4), pp 293-302, 1982. Eq. 1 See also -------- dirint disc """


    dni_extra = get_extra_radiation(datetime_or_doy)

    kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
                         max_clearness_index=1)

    # For Kt <= 0.22, set the diffuse fraction
    df = 1 - 0.09*kt

    # For Kt > 0.22 and Kt <= 0.8, set the diffuse fraction
    df = np.where((kt > 0.22) & (kt <= 0.8),
                  0.9511 - 0.1604*kt + 4.388*kt**2 -
                  16.638*kt**3 + 12.336*kt**4,
                  df)

    # For Kt > 0.8, set the diffuse fraction
    df = np.where(kt > 0.8, 0.165, df)

    dhi = df * ghi

    dni = (ghi - dhi) / tools.cosd(zenith)
    bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
    dni = np.where(bad_values, 0, dni)
    # ensure that closure relationship remains valid
    dhi = np.where(bad_values, ghi, dhi)

    data = OrderedDict()
    data['dni'] = dni
    data['dhi'] = dhi
    data['kt'] = kt

    if isinstance(datetime_or_doy, pd.DatetimeIndex):
        data = pd.DataFrame(data, index=datetime_or_doy)

    return data