def test_perez():
    AM = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    irradiance.perez(40, 180, irrad_data['DHI'], irrad_data['DNI'],
                            dni_et,
                            ephem_data['apparent_zenith'],
                            ephem_data['apparent_azimuth'],
                            AM) 

def test_perez_components():
    am = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out, df_components = irradiance.perez(40, 180, irrad_data['dhi'], dni,
                     dni_et, ephem_data['apparent_zenith'],
                     ephem_data['azimuth'], am, return_components=True)
    expected = pd.Series(np.array(
        [   0.        ,   31.46046871,  np.nan,   45.45539877]),
        index=times)
    expected_components = pd.DataFrame(
        np.array([[  0.        ,  26.84138589,          np.nan,  31.72696071],
                 [ 0.        ,  0.        ,         np.nan,  4.47966439],
                 [ 0.        ,  4.62212181,         np.nan,  9.25316454]]).T,
        columns=['isotropic', 'circumsolar', 'horizon'],
        index=times
    )
    if pandas_0_22():
        expected_for_sum = expected.copy()
        expected_for_sum.iloc[2] = 0
    else:
        expected_for_sum = expected
    sum_components = df_components.sum(axis=1)

    assert_series_equal(out, expected, check_less_precise=2)
    assert_frame_equal(df_components, expected_components)
    assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)

def test_perez_arrays():
    am = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez(40, 180, irrad_data['dhi'].values, dni.values,
                     dni_et, ephem_data['apparent_zenith'].values,
                     ephem_data['azimuth'].values, am.values)
    expected = np.array(
        [   0.        ,   31.46046871,  np.nan,   45.45539877])
    assert_allclose(out, expected, atol=1e-2)

def test_globalinplane():
    AOI = irradiance.aoi(40, 180, ephem_data['apparent_zenith'],
                         ephem_data['apparent_azimuth'])
    AM = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    gr_sand = irradiance.grounddiffuse(40, ghi, surface_type='sand')
    diff_perez = irradiance.perez(
        40, 180, irrad_data['DHI'], irrad_data['DNI'], dni_et,
        ephem_data['apparent_zenith'], ephem_data['apparent_azimuth'], AM)
    irradiance.globalinplane(
        AOI=AOI, DNI=irrad_data['DNI'], In_Plane_SkyDiffuse=diff_perez,
        GR=gr_sand)

def test_globalinplane():
    aoi = irradiance.aoi(40, 180, ephem_data['apparent_zenith'],
                         ephem_data['apparent_azimuth'])
    airmass = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    gr_sand = irradiance.grounddiffuse(40, ghi, surface_type='sand')
    diff_perez = irradiance.perez(
        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
        ephem_data['apparent_zenith'], ephem_data['apparent_azimuth'], airmass)
    irradiance.globalinplane(
        aoi=aoi, dni=irrad_data['dni'], poa_sky_diffuse=diff_perez,
        poa_ground_diffuse=gr_sand)

def test_perez():
    am = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez(40, 180, irrad_data['dhi'], dni,
                     dni_et, ephem_data['apparent_zenith'],
                     ephem_data['azimuth'], am)
    expected = pd.Series(np.array(
        [   0.        ,   31.46046871,  np.nan,   45.45539877]),
        index=times)
    assert_series_equal(out, expected, check_less_precise=2)

def test_total_irrad():
    models = ['isotropic', 'klutcher', 'haydavies', 'reindl', 'king', 'perez']
    AM = atmosphere.relativeairmass(ephem_data['apparent_zenith'])

    for model in models:
        total = irradiance.total_irrad(
            32, 180, 
            ephem_data['apparent_zenith'], ephem_data['azimuth'],
            dni=irrad_data['dni'], ghi=irrad_data['ghi'], dhi=irrad_data['dhi'],
            dni_extra=dni_et, airmass=AM,
            model=model,
            surface_type='urban')

    def get_airmass(self, times=None, solar_position=None,
                    model='kastenyoung1989'):
        """
        Calculate the relative and absolute airmass.
        
        Automatically chooses zenith or apparant zenith
        depending on the selected model.

        Parameters
        ----------
        times : None or DatetimeIndex
            Only used if solar_position is not provided.
        solar_position : None or DataFrame
            DataFrame with with columns 'apparent_zenith', 'zenith'.
        model : str
            Relative airmass model
        
        Returns
        -------
        airmass : DataFrame
            Columns are 'airmass_relative', 'airmass_absolute'
        """

        if solar_position is None:
            solar_position = self.get_solarposition(times)

        if model in atmosphere.APPARENT_ZENITH_MODELS:
            zenith = solar_position['apparent_zenith']
        elif model in atmosphere.TRUE_ZENITH_MODELS:
            zenith = solar_position['zenith']
        else:
            raise ValueError('{} is not a valid airmass model'.format(model))

        airmass_relative = atmosphere.relativeairmass(zenith, model)

        pressure = atmosphere.alt2pres(self.altitude)
        airmass_absolute = atmosphere.absoluteairmass(airmass_relative,
                                                      pressure)

        airmass = pd.DataFrame()
        airmass['airmass_relative'] = airmass_relative
        airmass['airmass_absolute'] = airmass_absolute

        return airmass
            

def test_total_irrad():
    models = ['isotropic', 'klutcher', 'klucher',
              'haydavies', 'reindl', 'king', 'perez']
    AM = atmosphere.relativeairmass(ephem_data['apparent_zenith'])

    for model in models:
        total = irradiance.total_irrad(
            32, 180,
            ephem_data['apparent_zenith'], ephem_data['azimuth'],
            dni=irrad_data['dni'], ghi=irrad_data['ghi'],
            dhi=irrad_data['dhi'],
            dni_extra=dni_et, airmass=AM,
            model=model,
            surface_type='urban')

        assert total.columns.tolist() == ['poa_global', 'poa_direct',
                                          'poa_diffuse', 'poa_sky_diffuse',
                                          'poa_ground_diffuse']

def test_ineichen_series():
    tus = Location(32.2, -111, 'US/Arizona', 700)
    times = pd.date_range(start='2014-06-24', end='2014-06-25', freq='3h')
    times_localized = times.tz_localize(tus.tz)
    ephem_data = solarposition.get_solarposition(times_localized, tus.latitude,
                                                 tus.longitude)
    am = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
    am = atmosphere.absoluteairmass(am, atmosphere.alt2pres(tus.altitude))
    expected = pd.DataFrame(np.
        array([[    0.        ,     0.        ,     0.        ],
               [    0.        ,     0.        ,     0.        ],
               [   91.12492792,   321.16092181,    51.17628184],
               [  716.46580533,   888.90147035,    99.5050056 ],
               [ 1053.42066043,   953.24925854,   116.32868969],
               [  863.54692781,   922.06124712,   106.95536561],
               [  271.06382274,   655.44925241,    73.05968071],
               [    0.        ,     0.        ,     0.        ],
               [    0.        ,     0.        ,     0.        ]]),
                            columns=['ghi', 'dni', 'dhi'],
                            index=times_localized)

    out = clearsky.ineichen(ephem_data['apparent_zenith'], am, 3)
    assert_frame_equal(expected, out)

def test_absoluteairmass():
    relative_am = atmosphere.relativeairmass(ephem_data['zenith'], 'simple')
    atmosphere.absoluteairmass(relative_am)
    atmosphere.absoluteairmass(relative_am, pressure=100000)

def test_airmass_invalid():
    with pytest.raises(ValueError):
        atmosphere.relativeairmass(ephem_data['zenith'], 'invalid')

def test_airmass_scalar_nan():
    assert np.isnan(atmosphere.relativeairmass(100))

def test_airmass(model):
    out = atmosphere.relativeairmass(ephem_data['zenith'], model)
    assert isinstance(out, pd.Series)
    out = atmosphere.relativeairmass(ephem_data['zenith'].values, model)
    assert isinstance(out, np.ndarray)

def test_bird():
    """Test Bird/Hulstrom Clearsky Model"""
    times = pd.DatetimeIndex(start='1/1/2015 0:00', end='12/31/2015 23:00',
                             freq='H')
    tz = -7  # test timezone
    gmt_tz = pytz.timezone('Etc/GMT%+d' % -(tz))
    times = times.tz_localize(gmt_tz)  # set timezone
    # match test data from BIRD_08_16_2012.xls
    latitude = 40.
    longitude = -105.
    press_mB = 840.
    o3_cm = 0.3
    h2o_cm = 1.5
    aod_500nm = 0.1
    aod_380nm = 0.15
    b_a = 0.85
    alb = 0.2
    eot = solarposition.equation_of_time_spencer71(times.dayofyear)
    hour_angle = solarposition.hour_angle(times, longitude, eot) - 0.5 * 15.
    declination = solarposition.declination_spencer71(times.dayofyear)
    zenith = solarposition.solar_zenith_analytical(
        np.deg2rad(latitude), np.deg2rad(hour_angle), declination
    )
    zenith = np.rad2deg(zenith)
    airmass = atmosphere.relativeairmass(zenith, model='kasten1966')
    etr = irradiance.extraradiation(times)
    # test Bird with time series data
    field_names = ('dni', 'direct_horizontal', 'ghi', 'dhi')
    irrads = clearsky.bird(
        zenith, airmass, aod_380nm, aod_500nm, h2o_cm, o3_cm, press_mB * 100.,
        etr, b_a, alb
    )
    Eb, Ebh, Gh, Dh = (irrads[_] for _ in field_names)
    clearsky_path = os.path.dirname(os.path.abspath(__file__))
    pvlib_path = os.path.dirname(clearsky_path)
    data_path = os.path.join(pvlib_path, 'data', 'BIRD_08_16_2012.csv')
    testdata = pd.read_csv(data_path, usecols=range(1, 26), header=1).dropna()
    testdata.index = times[1:48]
    assert np.allclose(testdata['DEC'], np.rad2deg(declination[1:48]))
    assert np.allclose(testdata['EQT'], eot[1:48], rtol=1e-4)
    assert np.allclose(testdata['Hour Angle'], hour_angle[1:48])
    assert np.allclose(testdata['Zenith Ang'], zenith[1:48])
    dawn = zenith < 88.
    dusk = testdata['Zenith Ang'] < 88.
    am = pd.Series(np.where(dawn, airmass, 0.), index=times).fillna(0.0)
    assert np.allclose(
        testdata['Air Mass'].where(dusk, 0.), am[1:48], rtol=1e-3
    )
    direct_beam = pd.Series(np.where(dawn, Eb, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata['Direct Beam'].where(dusk, 0.), direct_beam[1:48], rtol=1e-3
    )
    direct_horz = pd.Series(np.where(dawn, Ebh, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata['Direct Hz'].where(dusk, 0.), direct_horz[1:48], rtol=1e-3
    )
    global_horz = pd.Series(np.where(dawn, Gh, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata['Global Hz'].where(dusk, 0.), global_horz[1:48], rtol=1e-3
    )
    diffuse_horz = pd.Series(np.where(dawn, Dh, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:48], rtol=1e-3
    )
    # test keyword parameters
    irrads2 = clearsky.bird(
        zenith, airmass, aod_380nm, aod_500nm, h2o_cm, dni_extra=etr
    )
    Eb2, Ebh2, Gh2, Dh2 = (irrads2[_] for _ in field_names)
    clearsky_path = os.path.dirname(os.path.abspath(__file__))
    pvlib_path = os.path.dirname(clearsky_path)
    data_path = os.path.join(pvlib_path, 'data', 'BIRD_08_16_2012_patm.csv')
    testdata2 = pd.read_csv(data_path, usecols=range(1, 26), header=1).dropna()
    testdata2.index = times[1:48]
    direct_beam2 = pd.Series(np.where(dawn, Eb2, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata2['Direct Beam'].where(dusk, 0.), direct_beam2[1:48], rtol=1e-3
    )
    direct_horz2 = pd.Series(np.where(dawn, Ebh2, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata2['Direct Hz'].where(dusk, 0.), direct_horz2[1:48], rtol=1e-3
    )
    global_horz2 = pd.Series(np.where(dawn, Gh2, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata2['Global Hz'].where(dusk, 0.), global_horz2[1:48], rtol=1e-3
    )
    diffuse_horz2 = pd.Series(np.where(dawn, Dh2, 0.), index=times).fillna(0.)
    assert np.allclose(
        testdata2['Dif Hz'].where(dusk, 0.), diffuse_horz2[1:48], rtol=1e-3
    )
    # test scalars just at noon
    # XXX: calculations start at 12am so noon is at index = 12
    irrads3 = clearsky.bird(
        zenith[12], airmass[12], aod_380nm, aod_500nm, h2o_cm, dni_extra=etr[12]
    )
    Eb3, Ebh3, Gh3, Dh3 = (irrads3[_] for _ in field_names)
    # XXX: testdata starts at 1am so noon is at index = 11
    np.allclose(
        [Eb3, Ebh3, Gh3, Dh3],
        testdata2[['Direct Beam', 'Direct Hz', 'Global Hz', 'Dif Hz']].iloc[11],
#.........这里部分代码省略.........

#.........这里部分代码省略.........

    solar_position_method : str
        Passed to location.get_solarposition.

    airmass_model : str
        Passed to location.get_airmass.

    altitude : None or float
        If None, computed from pressure. Assumed to be 0 m
        if pressure is also None.

    pressure : None or float
        If None, computed from altitude. Assumed to be 101325 Pa
        if altitude is also None.

    **kwargs
        Arbitrary keyword arguments.
        See code for details.

    Returns
    -------
    output : (dc, ac)
        Tuple of DC power (with SAPM parameters) (DataFrame) and AC
        power (Series).
    """

    # use surface_tilt and surface_azimuth if provided,
    # otherwise set them using the orientation_strategy
    if surface_tilt is not None and surface_azimuth is not None:
        pass
    elif orientation_strategy is not None:
        surface_tilt, surface_azimuth = \
            get_orientation(orientation_strategy, latitude=latitude)
    else:
        raise ValueError('orientation_strategy or surface_tilt and ' +
                         'surface_azimuth must be provided')

    times = times

    if altitude is None and pressure is None:
        altitude = 0.
        pressure = 101325.
    elif altitude is None:
        altitude = atmosphere.pres2alt(pressure)
    elif pressure is None:
        pressure = atmosphere.alt2pres(altitude)

    solar_position = solarposition.get_solarposition(times, latitude,
                                                     longitude,
                                                     altitude=altitude,
                                                     pressure=pressure,
                                                     **kwargs)

    # possible error with using apparent zenith with some models
    airmass = atmosphere.relativeairmass(solar_position['apparent_zenith'],
                                         model=airmass_model)
    airmass = atmosphere.absoluteairmass(airmass, pressure)
    dni_extra = pvlib.irradiance.extraradiation(solar_position.index)
    dni_extra = pd.Series(dni_extra, index=solar_position.index)

    aoi = pvlib.irradiance.aoi(surface_tilt, surface_azimuth,
                               solar_position['apparent_zenith'],
                               solar_position['azimuth'])

    if irradiance is None:
        irradiance = clearsky.ineichen(
            solar_position.index,
            latitude,
            longitude,
            zenith_data=solar_position['apparent_zenith'],
            airmass_data=airmass,
            altitude=altitude)

    total_irrad = pvlib.irradiance.total_irrad(
        surface_tilt,
        surface_azimuth,
        solar_position['apparent_zenith'],
        solar_position['azimuth'],
        irradiance['dni'],
        irradiance['ghi'],
        irradiance['dhi'],
        model=transposition_model,
        dni_extra=dni_extra)

    if weather is None:
        weather = {'wind_speed': 0, 'temp_air': 20}

    temps = pvsystem.sapm_celltemp(total_irrad['poa_global'],
                                   weather['wind_speed'],
                                   weather['temp_air'])

    dc = pvsystem.sapm(module_parameters, total_irrad['poa_direct'],
                       total_irrad['poa_diffuse'],
                       temps['temp_cell'],
                       airmass,
                       aoi)

    ac = pvsystem.snlinverter(inverter_parameters, dc['v_mp'], dc['p_mp'])

    return dc, ac

    def get_irradiance(self, surface_tilt, surface_azimuth,
                       solar_zenith, solar_azimuth, dni, ghi, dhi,
                       dni_extra=None, airmass=None, model='haydavies',
                       **kwargs):
        """
        Uses the :func:`irradiance.total_irrad` function to calculate
        the plane of array irradiance components on a tilted surface
        defined by the input data and ``self.albedo``.

        For a given set of solar zenith and azimuth angles, the
        surface tilt and azimuth parameters are typically determined
        by :py:method:`~SingleAxisTracker.singleaxis`.

        Parameters
        ----------
        surface_tilt : numeric
            Panel tilt from horizontal.
        surface_azimuth : numeric
            Panel azimuth from north
        solar_zenith : numeric
            Solar zenith angle.
        solar_azimuth : numeric
            Solar azimuth angle.
        dni : float or Series
            Direct Normal Irradiance
        ghi : float or Series
            Global horizontal irradiance
        dhi : float or Series
            Diffuse horizontal irradiance
        dni_extra : float or Series, default None
            Extraterrestrial direct normal irradiance
        airmass : float or Series, default None
            Airmass
        model : String, default 'haydavies'
            Irradiance model.

        **kwargs
            Passed to :func:`irradiance.total_irrad`.

        Returns
        -------
        poa_irradiance : DataFrame
            Column names are: ``total, beam, sky, ground``.
        """

        # not needed for all models, but this is easier
        if dni_extra is None:
            dni_extra = irradiance.extraradiation(solar_zenith.index)

        if airmass is None:
            airmass = atmosphere.relativeairmass(solar_zenith)

        return irradiance.total_irrad(surface_tilt,
                                      surface_azimuth,
                                      solar_zenith,
                                      solar_azimuth,
                                      dni, ghi, dhi,
                                      dni_extra=dni_extra, airmass=airmass,
                                      model=model,
                                      albedo=self.albedo,
                                      **kwargs)

def test_airmass_invalid():
    atmosphere.relativeairmass(ephem_data['zenith'], 'invalid')

def ineichen(
    time,
    location,
    linke_turbidity=None,
    solarposition_method="pyephem",
    zenith_data=None,
    airmass_model="young1994",
    airmass_data=None,
    interp_turbidity=True,
):
    """
    Determine clear sky GHI, DNI, and DHI from Ineichen/Perez model

    Implements the Ineichen and Perez clear sky model for global horizontal
    irradiance (GHI), direct normal irradiance (DNI), and calculates
    the clear-sky diffuse horizontal (DHI) component as the difference
    between GHI and DNI*cos(zenith) as presented in [1, 2]. A report on clear
    sky models found the Ineichen/Perez model to have excellent performance
    with a minimal input data set [3]. 
    
    Default values for montly Linke turbidity provided by SoDa [4, 5].

    Parameters
    -----------
    time : pandas.DatetimeIndex
    
    location : pvlib.Location
    
    linke_turbidity : None or float
        If None, uses ``LinkeTurbidities.mat`` lookup table.
    
    solarposition_method : string
        Sets the solar position algorithm. 
        See solarposition.get_solarposition()
    
    zenith_data : None or pandas.Series
        If None, ephemeris data will be calculated using ``solarposition_method``.
    
    airmass_model : string
        See pvlib.airmass.relativeairmass().
    
    airmass_data : None or pandas.Series
        If None, absolute air mass data will be calculated using 
        ``airmass_model`` and location.alitude.
    
    interp_turbidity : bool
        If ``True``, interpolates the monthly Linke turbidity values
        found in ``LinkeTurbidities.mat`` to daily values.

    Returns
    --------
    DataFrame with the following columns: ``GHI, DNI, DHI``.

    Notes
    -----
    If you are using this function
    in a loop, it may be faster to load LinkeTurbidities.mat outside of
    the loop and feed it in as a variable, rather than
    having the function open the file each time it is called.

    References
    ----------

    [1] P. Ineichen and R. Perez, "A New airmass independent formulation for
        the Linke turbidity coefficient", Solar Energy, vol 73, pp. 151-157, 2002.

    [2] R. Perez et. al., "A New Operational Model for Satellite-Derived
        Irradiances: Description and Validation", Solar Energy, vol 73, pp.
        307-317, 2002.

    [3] M. Reno, C. Hansen, and J. Stein, "Global Horizontal Irradiance Clear
        Sky Models: Implementation and Analysis", Sandia National
        Laboratories, SAND2012-2389, 2012.

    [4] http://www.soda-is.com/eng/services/climat_free_eng.php#c5 (obtained
        July 17, 2012).

    [5] J. Remund, et. al., "Worldwide Linke Turbidity Information", Proc.
        ISES Solar World Congress, June 2003. Goteborg, Sweden.
    """
    # Initial implementation of this algorithm by Matthew Reno.
    # Ported to python by Rob Andrews
    # Added functionality by Will Holmgren

    I0 = irradiance.extraradiation(time.dayofyear)

    if zenith_data is None:
        ephem_data = solarposition.get_solarposition(time, location, method=solarposition_method)
        time = ephem_data.index  # fixes issue with time possibly not being tz-aware
        try:
            ApparentZenith = ephem_data["apparent_zenith"]
        except KeyError:
            ApparentZenith = ephem_data["zenith"]
            logger.warning("could not find apparent_zenith. using zenith")
    else:
        ApparentZenith = zenith_data
    # ApparentZenith[ApparentZenith >= 90] = 90 # can cause problems in edge cases

    if linke_turbidity is None:
        # The .mat file 'LinkeTurbidities.mat' contains a single 2160 x 4320 x 12
#.........这里部分代码省略.........

def ineichen(time, latitude, longitude, altitude=0, linke_turbidity=None,
             solarposition_method='nrel_numpy', zenith_data=None,
             airmass_model='young1994', airmass_data=None,
             interp_turbidity=True):
    '''
    Determine clear sky GHI, DNI, and DHI from Ineichen/Perez model

    Implements the Ineichen and Perez clear sky model for global horizontal
    irradiance (GHI), direct normal irradiance (DNI), and calculates
    the clear-sky diffuse horizontal (DHI) component as the difference
    between GHI and DNI*cos(zenith) as presented in [1, 2]. A report on clear
    sky models found the Ineichen/Perez model to have excellent performance
    with a minimal input data set [3].

    Default values for montly Linke turbidity provided by SoDa [4, 5].

    Parameters
    -----------
    time : pandas.DatetimeIndex

    latitude : float

    longitude : float

    altitude : float

    linke_turbidity : None or float
        If None, uses ``LinkeTurbidities.mat`` lookup table.

    solarposition_method : string
        Sets the solar position algorithm.
        See solarposition.get_solarposition()

    zenith_data : None or Series
        If None, ephemeris data will be calculated using ``solarposition_method``.

    airmass_model : string
        See pvlib.airmass.relativeairmass().

    airmass_data : None or Series
        If None, absolute air mass data will be calculated using
        ``airmass_model`` and location.alitude.

    interp_turbidity : bool
        If ``True``, interpolates the monthly Linke turbidity values
        found in ``LinkeTurbidities.mat`` to daily values.

    Returns
    --------
    DataFrame with the following columns: ``ghi, dni, dhi``.

    Notes
    -----
    If you are using this function
    in a loop, it may be faster to load LinkeTurbidities.mat outside of
    the loop and feed it in as a keyword argument, rather than
    having the function open and process the file each time it is called.

    References
    ----------

    [1] P. Ineichen and R. Perez, "A New airmass independent formulation for
        the Linke turbidity coefficient", Solar Energy, vol 73, pp. 151-157, 2002.

    [2] R. Perez et. al., "A New Operational Model for Satellite-Derived
        Irradiances: Description and Validation", Solar Energy, vol 73, pp.
        307-317, 2002.

    [3] M. Reno, C. Hansen, and J. Stein, "Global Horizontal Irradiance Clear
        Sky Models: Implementation and Analysis", Sandia National
        Laboratories, SAND2012-2389, 2012.

    [4] http://www.soda-is.com/eng/services/climat_free_eng.php#c5 (obtained
        July 17, 2012).

    [5] J. Remund, et. al., "Worldwide Linke Turbidity Information", Proc.
        ISES Solar World Congress, June 2003. Goteborg, Sweden.
    '''
    # Initial implementation of this algorithm by Matthew Reno.
    # Ported to python by Rob Andrews
    # Added functionality by Will Holmgren (@wholmgren)

    I0 = irradiance.extraradiation(time.dayofyear)

    if zenith_data is None:
        ephem_data = solarposition.get_solarposition(time,
                                                     latitude=latitude,
                                                     longitude=longitude,
                                                     altitude=altitude,
                                                     method=solarposition_method)
        time = ephem_data.index # fixes issue with time possibly not being tz-aware
        try:
            ApparentZenith = ephem_data['apparent_zenith']
        except KeyError:
            ApparentZenith = ephem_data['zenith']
            logger.warning('could not find apparent_zenith. using zenith')
    else:
        ApparentZenith = zenith_data
    #ApparentZenith[ApparentZenith >= 90] = 90 # can cause problems in edge cases

#.........这里部分代码省略.........