def aoi_projection(surf_tilt, surf_az, sun_zen, sun_az):
    """
    Calculates the dot product of the solar vector and the surface normal.

    Input all angles in degrees.

    Parameters
    ==========

    surf_tilt : float or Series.
        Panel tilt from horizontal.
    surf_az : float or Series.
        Panel azimuth from north.
    sun_zen : float or Series.
        Solar zenith angle.
    sun_az : float or Series.
        Solar azimuth angle.

    Returns
    =======
    float or Series. Dot product of panel normal and solar angle.
    """

    projection = (
        tools.cosd(surf_tilt) * tools.cosd(sun_zen) + tools.sind(surf_tilt) *
        tools.sind(sun_zen) * tools.cosd(sun_az - surf_az))

    try:
        projection.name = 'aoi_projection'
    except AttributeError:
        pass

    return projection

def king(surf_tilt, DHI, GHI, sun_zen):
    '''
    Determine diffuse irradiance from the sky on a tilted surface using the
    King model.

    King's model determines the diffuse irradiance from the sky
    (ground reflected irradiance is not included in this algorithm) on a
    tilted surface using the surface tilt angle, diffuse horizontal
    irradiance, global horizontal irradiance, and sun zenith angle. Note
    that this model is not well documented and has not been published in
    any fashion (as of January 2012).

    Parameters
    ----------

    surf_tilt : float or Series
          Surface tilt angles in decimal degrees.
          The tilt angle is defined as
          degrees from horizontal (e.g. surface facing up = 0, surface facing
          horizon = 90)

    DHI : float or Series
          diffuse horizontal irradiance in W/m^2.

    GHI : float or Series
          global horizontal irradiance in W/m^2.

    sun_zen : float or Series
          apparent (refraction-corrected) zenith
          angles in decimal degrees.

    Returns
    --------

    SkyDiffuse : float or Series

            the diffuse component of the solar radiation  on an
            arbitrarily tilted surface as given by a model developed by
            David L. King at Sandia National Laboratories.
    '''

    pvl_logger.debug('diffuse_sky.king()')

    sky_diffuse = (DHI * ((1 + tools.cosd(surf_tilt))) / 2 + GHI *
                   ((0.012 * sun_zen - 0.04)) *
                   ((1 - tools.cosd(surf_tilt))) / 2)
    sky_diffuse[sky_diffuse < 0] = 0

    return sky_diffuse

def isotropic(surf_tilt, DHI):
    r'''
    Determine diffuse irradiance from the sky on a
    tilted surface using the isotropic sky model.

    .. math::

       I_{d} = DHI \frac{1 + \cos\beta}{2}

    Hottel and Woertz's model treats the sky as a uniform source of diffuse
    irradiance. Thus the diffuse irradiance from the sky (ground reflected
    irradiance is not included in this algorithm) on a tilted surface can
    be found from the diffuse horizontal irradiance and the tilt angle of
    the surface.

    Parameters
    ----------

    surf_tilt : float or Series
        Surface tilt angle in decimal degrees.
        surf_tilt must be >=0 and <=180. The tilt angle is defined as
        degrees from horizontal (e.g. surface facing up = 0, surface facing
        horizon = 90)

    DHI : float or Series
        Diffuse horizontal irradiance in W/m^2.
        DHI must be >=0.


    Returns
    -------
    float or Series

    The diffuse component of the solar radiation  on an
    arbitrarily tilted surface defined by the isotropic sky model as
    given in Loutzenhiser et. al (2007) equation 3.
    SkyDiffuse is the diffuse component ONLY and does not include the ground
    reflected irradiance or the irradiance due to the beam.
    SkyDiffuse is a column vector vector with a number of elements equal to
    the input vector(s).


    References
    ----------

    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
    solar irradiance on inclined surfaces for building energy simulation"
    2007, Solar Energy vol. 81. pp. 254-267

    [2] Hottel, H.C., Woertz, B.B., 1942. Evaluation of flat-plate solar heat
    collector. Trans. ASME 64, 91.
    '''

    pvl_logger.debug('diffuse_sky.isotropic()')

    sky_diffuse = DHI * (1 + tools.cosd(surf_tilt)) * 0.5

    return sky_diffuse

def haurwitz(ApparentZenith):
    """
    Determine clear sky GHI from Haurwitz model
   
    Implements the Haurwitz clear sky model for global horizontal
    irradiance (GHI) as presented in [1, 2]. A report on clear
    sky models found the Haurwitz model to have the best performance of
    models which require only zenith angle [3]. Extreme care should
    be taken in the interpretation of this result! 

    Parameters
    ----------
    ApparentZenith : DataFrame
        The apparent (refraction corrected) sun zenith angle
        in degrees.

    Returns
    -------        
    pd.Series. The modeled global horizonal irradiance in W/m^2 provided
    by the Haurwitz clear-sky model.

    Initial implementation of this algorithm by Matthew Reno.

    References
    ----------

    [1] B. Haurwitz, "Insolation in Relation to Cloudiness and Cloud 
     Density," Journal of Meteorology, vol. 2, pp. 154-166, 1945.

    [2] B. Haurwitz, "Insolation in Relation to Cloud Type," Journal of 
     Meteorology, vol. 3, pp. 123-124, 1946.

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

    See Also
    ---------
    maketimestruct    
    makelocationstruct   
    ephemeris   
    spa
    ineichen
    """

    cos_zenith = tools.cosd(ApparentZenith)

    clearsky_GHI = 1098.0 * cos_zenith * np.exp(-0.059 / cos_zenith)

    clearsky_GHI[clearsky_GHI < 0] = 0

    df_out = pd.DataFrame({"GHI": clearsky_GHI})

    return df_out

def haurwitz(apparent_zenith):
    '''
    Determine clear sky GHI from Haurwitz model.

    Implements the Haurwitz clear sky model for global horizontal
    irradiance (GHI) as presented in [1, 2]. A report on clear
    sky models found the Haurwitz model to have the best performance
    in terms of average monthly error among models which require only
    zenith angle [3].

    Parameters
    ----------
    apparent_zenith : Series
        The apparent (refraction corrected) sun zenith angle
        in degrees.

    Returns
    -------
    pd.DataFrame
    The modeled global horizonal irradiance in W/m^2 provided
    by the Haurwitz clear-sky model.

    Initial implementation of this algorithm by Matthew Reno.

    References
    ----------

    [1] B. Haurwitz, "Insolation in Relation to Cloudiness and Cloud
     Density," Journal of Meteorology, vol. 2, pp. 154-166, 1945.

    [2] B. Haurwitz, "Insolation in Relation to Cloud Type," Journal of
     Meteorology, vol. 3, pp. 123-124, 1946.

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

    cos_zenith = tools.cosd(apparent_zenith.values)
    clearsky_ghi = np.zeros_like(apparent_zenith.values)
    cos_zen_gte_0 = cos_zenith > 0
    clearsky_ghi[cos_zen_gte_0] = (1098.0 * cos_zenith[cos_zen_gte_0] *
                                   np.exp(-0.059/cos_zenith[cos_zen_gte_0]))

    df_out = pd.DataFrame(index=apparent_zenith.index,
                          data=clearsky_ghi,
                          columns=['ghi'])

    return df_out

def poa_horizontal_ratio(surf_tilt, surf_az, sun_zen, sun_az):
    """
    Calculates the ratio of the beam components of the
    plane of array irradiance and the horizontal irradiance.

    Input all angles in degrees.

    Parameters
    ==========

    surf_tilt : float or Series.
        Panel tilt from horizontal.
    surf_az : float or Series.
        Panel azimuth from north.
    sun_zen : float or Series.
        Solar zenith angle.
    sun_az : float or Series.
        Solar azimuth angle.

    Returns
    =======
    float or Series. Ratio of the plane of array irradiance to the
    horizontal plane irradiance
    """

    cos_poa_zen = aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)

    cos_sun_zen = tools.cosd(sun_zen)

    # ratio of titled and horizontal beam irradiance
    ratio = cos_poa_zen / cos_sun_zen

    try:
        ratio.name = 'poa_ratio'
    except AttributeError:
        pass

    return ratio

def physicaliam(K, L, n, aoi):
    '''
    Determine the incidence angle modifier using refractive 
    index, glazing thickness, and extinction coefficient

    physicaliam calculates the incidence angle modifier as described in
    De Soto et al. "Improvement and validation of a model for photovoltaic
    array performance", section 3. The calculation is based upon a physical
    model of absorbtion and transmission through a cover. Required
    information includes, incident angle, cover extinction coefficient,
    cover thickness

    Note: The authors of this function believe that eqn. 14 in [1] is
    incorrect. This function uses the following equation in its place:
    theta_r = arcsin(1/n * sin(theta))

    Parameters
    ----------
    K : float
        The glazing extinction coefficient in units of 1/meters. Reference
        [1] indicates that a value of  4 is reasonable for "water white"
        glass. K must be a numeric scalar or vector with all values >=0. If K
        is a vector, it must be the same size as all other input vectors.

    L : float
        The glazing thickness in units of meters. Reference [1] indicates
        that 0.002 meters (2 mm) is reasonable for most glass-covered
        PV panels. L must be a numeric scalar or vector with all values >=0. 
        If L is a vector, it must be the same size as all other input vectors.

    n : float
        The effective index of refraction (unitless). Reference [1]
        indicates that a value of 1.526 is acceptable for glass. n must be a 
        numeric scalar or vector with all values >=0. If n is a vector, it 
        must be the same size as all other input vectors.

    aoi : Series
        The angle of incidence between the module normal vector and the
        sun-beam vector in degrees. 

    Returns
    -------
    IAM : float or Series
        The incident angle modifier as specified in eqns. 14-16 of [1].
        IAM is a column vector with the same number of elements as the
        largest input vector.
        
        Theta must be a numeric scalar or vector.
        For any values of theta where abs(aoi)>90, IAM is set to 0. For any
        values of aoi where -90 < aoi < 0, theta is set to abs(aoi) and
        evaluated.

    References
    ----------
    [1] W. De Soto et al., "Improvement and validation of a model for
    photovoltaic array performance", Solar Energy, vol 80, pp. 78-88,
    2006.

    [2] Duffie, John A. & Beckman, William A.. (2006). Solar Engineering 
    of Thermal Processes, third edition. [Books24x7 version] Available 
    from http://common.books24x7.com/toc.aspx?bookid=17160. 

    See Also 
    --------
    getaoi   
    ephemeris   
    spa    
    ashraeiam
    '''
    thetar_deg = tools.asind(1.0 / n*(tools.sind(aoi)))

    tau = ( np.exp(- 1.0 * (K*L / tools.cosd(thetar_deg))) *
            ((1 - 0.5*((((tools.sind(thetar_deg - aoi)) ** 2) /
            ((tools.sind(thetar_deg + aoi)) ** 2) +
            ((tools.tand(thetar_deg - aoi)) ** 2) /
            ((tools.tand(thetar_deg + aoi)) ** 2))))) )
    
    zeroang = 1e-06
    
    thetar_deg0 = tools.asind(1.0 / n*(tools.sind(zeroang)))
    
    tau0 = ( np.exp(- 1.0 * (K*L / tools.cosd(thetar_deg0))) *
             ((1 - 0.5*((((tools.sind(thetar_deg0 - zeroang)) ** 2) /
             ((tools.sind(thetar_deg0 + zeroang)) ** 2) +
             ((tools.tand(thetar_deg0 - zeroang)) ** 2) /
             ((tools.tand(thetar_deg0 + zeroang)) ** 2))))) )
    
    IAM = tau / tau0
    
    IAM[abs(aoi) >= 90] = np.nan
    IAM[IAM < 0] = np.nan
    
    return IAM

def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
             altitude=0, dni_extra=1364., perez_enhancement=False):
    '''
    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 monthly Linke turbidity provided by SoDa [4, 5].

    Parameters
    -----------
    apparent_zenith : numeric
        Refraction corrected solar zenith angle in degrees.

    airmass_absolute : numeric
        Pressure corrected airmass.

    linke_turbidity : numeric
        Linke Turbidity.

    altitude : numeric, default 0
        Altitude above sea level in meters.

    dni_extra : numeric, default 1364
        Extraterrestrial irradiance. The units of ``dni_extra``
        determine the units of the output.

    perez_enhancement : bool, default False
        Controls if the Perez enhancement factor should be applied.
        Setting to True may produce spurious results for times when
        the Sun is near the horizon and the airmass is high.
        See https://github.com/pvlib/pvlib-python/issues/435

    Returns
    -------
    clearsky : DataFrame (if Series input) or OrderedDict of arrays
        DataFrame/OrderedDict contains the columns/keys
        ``'dhi', 'dni', 'ghi'``.

    See also
    --------
    lookup_linke_turbidity
    pvlib.location.Location.get_clearsky

    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.
    '''

    # ghi is calculated using either the equations in [1] by setting
    # perez_enhancement=False (default behavior) or using the model
    # in [2] by setting perez_enhancement=True.

    # The NaN handling is a little subtle. The AM input is likely to
    # have NaNs that we'll want to map to 0s in the output. However, we
    # want NaNs in other inputs to propagate through to the output. This
    # is accomplished by judicious use and placement of np.maximum,
    # np.minimum, and np.fmax

    # use max so that nighttime values will result in 0s instead of
    # negatives. propagates nans.
    cos_zenith = np.maximum(tools.cosd(apparent_zenith), 0)

    tl = linke_turbidity

    fh1 = np.exp(-altitude/8000.)
    fh2 = np.exp(-altitude/1250.)
    cg1 = 5.09e-05 * altitude + 0.868
    cg2 = 3.92e-05 * altitude + 0.0387

    ghi = np.exp(-cg2*airmass_absolute*(fh1 + fh2*(tl - 1)))

    # https://github.com/pvlib/pvlib-python/issues/435
    if perez_enhancement:
        ghi *= np.exp(0.01*airmass_absolute**1.8)

    # use fmax to map airmass nans to 0s. multiply and divide by tl to
    # reinsert tl nans
    ghi = cg1 * dni_extra * cos_zenith * tl / tl * np.fmax(ghi, 0)
#.........这里部分代码省略.........

#.........这里部分代码省略.........
    # MATLAB to Python conversion by
    # Will Holmgren (@wholmgren), U. Arizona. March, 2015.

    if isinstance(apparent_zenith, pd.Series):
        index = apparent_zenith.index
    else:
        index = None

    # convert scalars to arrays
    apparent_azimuth = np.atleast_1d(apparent_azimuth)
    apparent_zenith = np.atleast_1d(apparent_zenith)

    if apparent_azimuth.ndim > 1 or apparent_zenith.ndim > 1:
        raise ValueError('Input dimensions must not exceed 1')

    # Calculate sun position x, y, z using coordinate system as in [1], Eq 2.

    # Positive y axis is oriented parallel to earth surface along tracking axis
    # (for the purpose of illustration, assume y is oriented to the south);
    # positive x axis is orthogonal, 90 deg clockwise from y-axis, and parallel
    # to the earth's surface (if y axis is south, x axis is west);
    # positive z axis is normal to x, y axes, pointed upward.

    # Equations in [1] assume solar azimuth is relative to reference vector
    # pointed south, with clockwise positive.
    # Here, the input solar azimuth is degrees East of North,
    # i.e., relative to a reference vector pointed
    # north with clockwise positive.
    # Rotate sun azimuth to coordinate system as in [1]
    # to calculate sun position.

    az = apparent_azimuth - 180
    apparent_elevation = 90 - apparent_zenith
    x = cosd(apparent_elevation) * sind(az)
    y = cosd(apparent_elevation) * cosd(az)
    z = sind(apparent_elevation)

    # translate array azimuth from compass bearing to [1] coord system
    # wholmgren: strange to see axis_azimuth calculated differently from az,
    # (not that it matters, or at least it shouldn't...).
    axis_azimuth_south = axis_azimuth - 180

    # translate input array tilt angle axis_tilt to [1] coordinate system.

    # In [1] coordinates, axis_tilt is a rotation about the x-axis.
    # For a system with array azimuth (y-axis) oriented south,
    # the x-axis is oriented west, and a positive axis_tilt is a
    # counterclockwise rotation, i.e, lifting the north edge of the panel.
    # Thus, in [1] coordinate system, in the northern hemisphere a positive
    # axis_tilt indicates a rotation toward the equator,
    # whereas in the southern hemisphere rotation toward the equator is
    # indicated by axis_tilt<0.  Here, the input axis_tilt is
    # always positive and is a rotation toward the equator.

    # Calculate sun position (xp, yp, zp) in panel-oriented coordinate system:
    # positive y-axis is oriented along tracking axis at panel tilt;
    # positive x-axis is orthogonal, clockwise, parallel to earth surface;
    # positive z-axis is normal to x-y axes, pointed upward.
    # Calculate sun position (xp,yp,zp) in panel coordinates using [1] Eq 11
    # note that equation for yp (y' in Eq. 11 of Lorenzo et al 2011) is
    # corrected, after conversation with paper's authors.

    xp = x*cosd(axis_azimuth_south) - y*sind(axis_azimuth_south)
    yp = (x*cosd(axis_tilt)*sind(axis_azimuth_south) +
          y*cosd(axis_tilt)*cosd(axis_azimuth_south) -
          z*sind(axis_tilt))

#.........这里部分代码省略.........
        # Note that the numbers within the matrix are 20 * Linke Turbidity,
        # so divide the number from the file by 20 to get the
        # turbidity.

        try:
            import scipy.io
        except ImportError:
            raise ImportError(
                "The Linke turbidity lookup table requires scipy. "
                + "You can still use clearsky.ineichen if you "
                + "supply your own turbidities."
            )

        # consider putting this code at module level
        this_path = os.path.dirname(os.path.abspath(__file__))
        logger.debug("this_path={}".format(this_path))

        mat = scipy.io.loadmat(os.path.join(this_path, "data", "LinkeTurbidities.mat"))
        linke_turbidity = mat["LinkeTurbidity"]
        LatitudeIndex = np.round_(_linearly_scale(location.latitude, 90, -90, 1, 2160))
        LongitudeIndex = np.round_(_linearly_scale(location.longitude, -180, 180, 1, 4320))
        g = linke_turbidity[LatitudeIndex][LongitudeIndex]
        if interp_turbidity:
            logger.info("interpolating turbidity to the day")
            g2 = np.concatenate([[g[-1]], g, [g[0]]])  # wrap ends around
            days = np.linspace(-15, 380, num=14)  # map day of year onto month (approximate)
            LT = pd.Series(np.interp(time.dayofyear, days, g2), index=time)
        else:
            logger.info("using monthly turbidity")
            ApplyMonth = lambda x: g[x[0] - 1]
            LT = pd.DataFrame(time.month, index=time)
            LT = LT.apply(ApplyMonth, axis=1)
        TL = LT / 20.0
        logger.info("using TL=\n{}".format(TL))
    else:
        TL = linke_turbidity

    # Get the absolute airmass assuming standard local pressure (per
    # alt2pres) using Kasten and Young's 1989 formula for airmass.

    if airmass_data is None:
        AMabsolute = atmosphere.absoluteairmass(
            AMrelative=atmosphere.relativeairmass(ApparentZenith, airmass_model),
            pressure=atmosphere.alt2pres(location.altitude),
        )
    else:
        AMabsolute = airmass_data

    fh1 = np.exp(-location.altitude / 8000.0)
    fh2 = np.exp(-location.altitude / 1250.0)
    cg1 = 5.09e-05 * location.altitude + 0.868
    cg2 = 3.92e-05 * location.altitude + 0.0387
    logger.debug("fh1={}, fh2={}, cg1={}, cg2={}".format(fh1, fh2, cg1, cg2))

    #  Dan's note on the TL correction: By my reading of the publication on
    #  pages 151-157, Ineichen and Perez introduce (among other things) three
    #  things. 1) Beam model in eqn. 8, 2) new turbidity factor in eqn 9 and
    #  appendix A, and 3) Global horizontal model in eqn. 11. They do NOT appear
    #  to use the new turbidity factor (item 2 above) in either the beam or GHI
    #  models. The phrasing of appendix A seems as if there are two separate
    #  corrections, the first correction is used to correct the beam/GHI models,
    #  and the second correction is used to correct the revised turibidity
    #  factor. In my estimation, there is no need to correct the turbidity
    #  factor used in the beam/GHI models.

    #  Create the corrected TL for TL < 2
    #  TLcorr = TL;
    #  TLcorr(TL < 2) = TLcorr(TL < 2) - 0.25 .* (2-TLcorr(TL < 2)) .^ (0.5);

    #  This equation is found in Solar Energy 73, pg 311.
    #  Full ref: Perez et. al., Vol. 73, pp. 307-317 (2002).
    #  It is slightly different than the equation given in Solar Energy 73, pg 156.
    #  We used the equation from pg 311 because of the existence of known typos
    #  in the pg 156 publication (notably the fh2-(TL-1) should be fh2 * (TL-1)).

    cos_zenith = tools.cosd(ApparentZenith)

    clearsky_GHI = (
        cg1 * I0 * cos_zenith * np.exp(-cg2 * AMabsolute * (fh1 + fh2 * (TL - 1))) * np.exp(0.01 * AMabsolute ** 1.8)
    )
    clearsky_GHI[clearsky_GHI < 0] = 0

    # BncI == "normal beam clear sky radiation"
    b = 0.664 + 0.163 / fh1
    BncI = b * I0 * np.exp(-0.09 * AMabsolute * (TL - 1))
    logger.debug("b={}".format(b))

    # "empirical correction" SE 73, 157 & SE 73, 312.
    BncI_2 = clearsky_GHI * (1 - (0.1 - 0.2 * np.exp(-TL)) / (0.1 + 0.882 / fh1)) / cos_zenith
    # return BncI, BncI_2
    clearsky_DNI = np.minimum(BncI, BncI_2)  # Will H: use np.minimum explicitly

    clearsky_DHI = clearsky_GHI - clearsky_DNI * cos_zenith

    df_out = pd.DataFrame({"GHI": clearsky_GHI, "DNI": clearsky_DNI, "DHI": clearsky_DHI})
    df_out.fillna(0, inplace=True)
    # df_out['BncI'] = BncI
    # df_out['BncI_2'] = BncI

    return df_out

#.........这里部分代码省略.........
    float or Series

          the diffuse component of the solar radiation  on an
          arbitrarily tilted surface defined by the Perez model as given in
          reference [3].
          SkyDiffuse is the diffuse component ONLY and does not include the
          ground reflected irradiance or the irradiance due to the beam.


    References
    ----------

    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
    solar irradiance on inclined surfaces for building energy simulation"
    2007, Solar Energy vol. 81. pp. 254-267

    [2] Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D., 1987.
    A new simplified version of the Perez diffuse irradiance model for tilted
    surfaces. Solar Energy 39(3), 221-232.

    [3] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R., 1990.
    Modeling daylight availability and irradiance components from direct
    and global irradiance. Solar Energy 44 (5), 271-289.

    [4] Perez, R. et. al 1988. "The Development and Verification of the
    Perez Diffuse Radiation Model". SAND88-7030
    '''

    pvl_logger.debug('diffuse_sky.perez()')

    kappa = 1.041  # for sun_zen in radians
    z = np.radians(sun_zen)  # convert to radians

    # epsilon is the sky's "clearness"
    eps = ((DHI + DNI) / DHI + kappa * (z ** 3)) / (1 + kappa * (z ** 3))

    # Perez et al define clearness bins according to the following rules.
    # 1 = overcast ... 8 = clear
    # (these names really only make sense for small zenith angles, but...)
    # these values will eventually be used as indicies for coeffecient look ups
    ebin = eps.copy()
    ebin[(eps < 1.065)] = 1
    ebin[(eps >= 1.065) & (eps < 1.23)] = 2
    ebin[(eps >= 1.23) & (eps < 1.5)] = 3
    ebin[(eps >= 1.5) & (eps < 1.95)] = 4
    ebin[(eps >= 1.95) & (eps < 2.8)] = 5
    ebin[(eps >= 2.8) & (eps < 4.5)] = 6
    ebin[(eps >= 4.5) & (eps < 6.2)] = 7
    ebin[eps >= 6.2] = 8

    ebin = ebin - 1  # correct for 0 indexing in coeffecient lookup

    # remove night time values
    ebin = ebin.dropna().astype(int)

    # This is added because in cases where the sun is below the horizon
    # (var.sun_zen > 90) but there is still diffuse horizontal light
    # (var.DHI>0), it is possible that the airmass (var.AM) could be NaN, which
    # messes up later calculations. Instead, if the sun is down, and there is
    # still var.DHI, we set the airmass to the airmass value on the horizon
    # (approximately 37-38).
    # var.AM(var.sun_zen >=90 & var.DHI >0) = 37;

    # var.DNI_ET[var.DNI_ET==0] = .00000001 #very hacky, fix this

    # delta is the sky's "brightness"
    delta = DHI * AM / DNI_ET

    # keep only valid times
    delta = delta[ebin.index]
    z = z[ebin.index]

    # The various possible sets of Perez coefficients are contained
    # in a subfunction to clean up the code.
    F1c, F2c = _get_perez_coefficients(modelt)

    F1 = F1c[ebin, 0] + F1c[ebin, 1] * delta + F1c[ebin, 2] * z
    F1[F1 < 0] = 0
    F1 = F1.astype(float)

    F2 = F2c[ebin, 0] + F2c[ebin, 1] * delta + F2c[ebin, 2] * z
    F2[F2 < 0] = 0
    F2 = F2.astype(float)

    A = aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)
    A[A < 0] = 0

    B = tools.cosd(sun_zen)
    B[B < tools.cosd(85)] = tools.cosd(85)

    # Calculate Diffuse POA from sky dome

    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surf_tilt))
    term2 = F1 * A[ebin.index] / B[ebin.index]
    term3 = F2 * tools.sind(surf_tilt)

    sky_diffuse = DHI[ebin.index] * (term1 + term2 + term3)
    sky_diffuse[sky_diffuse < 0] = 0

    return sky_diffuse

    def complete_irradiance(self, times=None, weather=None):
        """
        Determine the missing irradiation columns. Only two of the
        following data columns (dni, ghi, dhi) are needed to calculate
        the missing data.

        This function is not safe at the moment. Results can be too high
        or negative. Please contribute and help to improve this function
        on https://github.com/pvlib/pvlib-python

        Parameters
        ----------
        times : None or DatetimeIndex, default None
            Times at which to evaluate the model. Can be None if
            attribute `times` is already set.
        weather : None or pandas.DataFrame, default None
            Table with at least two columns containing one of the
            following data sets: dni, dhi, ghi. Can be None if attribute
            `weather` is already set.

        Returns
        -------
        self

        Assigns attributes: times, weather

        Examples
        --------
        This example does not work until the parameters `my_system`,
        `my_location`, `my_datetime` and `my_weather` are not defined
        properly but shows the basic idea how this method can be used.

        >>> from pvlib.modelchain import ModelChain

        >>> # my_weather containing 'dhi' and 'ghi'.
        >>> mc = ModelChain(my_system, my_location)  # doctest: +SKIP
        >>> mc.complete_irradiance(my_datetime, my_weather)  # doctest: +SKIP
        >>> mc.run_model()  # doctest: +SKIP

        >>> # my_weather containing 'dhi', 'ghi' and 'dni'.
        >>> mc = ModelChain(my_system, my_location)  # doctest: +SKIP
        >>> mc.run_model(my_datetime, my_weather)  # doctest: +SKIP
        """
        if weather is not None:
            self.weather = weather
        if times is not None:
            self.times = times
        self.solar_position = self.location.get_solarposition(
            self.times, method=self.solar_position_method)
        icolumns = set(self.weather.columns)
        wrn_txt = ("This function is not safe at the moment.\n" +
                   "Results can be too high or negative.\n" +
                   "Help to improve this function on github:\n" +
                   "https://github.com/pvlib/pvlib-python \n")

        if {'ghi', 'dhi'} <= icolumns and 'dni' not in icolumns:
            clearsky = self.location.get_clearsky(
                times, solar_position=self.solar_position)
            self.weather.loc[:, 'dni'] = pvlib.irradiance.dni(
                self.weather.loc[:, 'ghi'], self.weather.loc[:, 'dhi'],
                self.solar_position.zenith,
                clearsky_dni=clearsky['dni'],
                clearsky_tolerance=1.1)
        elif {'dni', 'dhi'} <= icolumns and 'ghi' not in icolumns:
            warnings.warn(wrn_txt, UserWarning)
            self.weather.loc[:, 'ghi'] = (
                self.weather.dni * tools.cosd(self.solar_position.zenith) +
                self.weather.dhi)
        elif {'dni', 'ghi'} <= icolumns and 'dhi' not in icolumns:
            warnings.warn(wrn_txt, UserWarning)
            self.weather.loc[:, 'dhi'] = (
                self.weather.ghi - self.weather.dni *
                tools.cosd(self.solar_position.zenith))

        return self

def reindl(surf_tilt, surf_az, DHI, DNI, GHI, DNI_ET, sun_zen, sun_az):
    r'''
    Determine diffuse irradiance from the sky on a
    tilted surface using Reindl's 1990 model

    .. math::

       I_{d} = DHI (A R_b + (1 - A) (\frac{1 + \cos\beta}{2})
       (1 + \sqrt{\frac{I_{hb}}{I_h}} \sin^3(\beta/2)) )

    Reindl's 1990 model determines the diffuse irradiance from the sky
    (ground reflected irradiance is not included in this algorithm) on a
    tilted surface using the surface tilt angle, surface azimuth angle,
    diffuse horizontal irradiance, direct normal irradiance, global
    horizontal irradiance, extraterrestrial irradiance, sun zenith angle,
    and sun azimuth angle.

    Parameters
    ----------

    surf_tilt : float or Series.
        Surface tilt angles in decimal degrees.
        The tilt angle is defined as
        degrees from horizontal (e.g. surface facing up = 0, surface facing
        horizon = 90)

    surf_az : float or Series.
        Surface azimuth angles in decimal degrees.
        The Azimuth convention is defined
        as degrees east of north (e.g. North = 0, South=180 East = 90,
        West = 270).

    DHI : float or Series.
        diffuse horizontal irradiance in W/m^2.

    DNI : float or Series.
        direct normal irradiance in W/m^2.

    GHI: float or Series.
        Global irradiance in W/m^2.

    DNI_ET : float or Series.
        extraterrestrial normal irradiance in W/m^2.

    sun_zen : float or Series.
        apparent (refraction-corrected) zenith
        angles in decimal degrees.

    sun_az : float or Series.
        Sun azimuth angles in decimal degrees.
        The Azimuth convention is defined
        as degrees east of north (e.g. North = 0, East = 90, West = 270).

    Returns
    -------

    SkyDiffuse : float or Series.

        The diffuse component of the solar radiation  on an
        arbitrarily tilted surface defined by the Reindl model as given in
        Loutzenhiser et. al (2007) equation 8.
        SkyDiffuse is the diffuse component ONLY and does not include the
        ground reflected irradiance or the irradiance due to the beam.
        SkyDiffuse is a column vector vector with a number of elements equal to
        the input vector(s).


    Notes
    -----

    The POAskydiffuse calculation is generated from the Loutzenhiser et al.
    (2007) paper, equation 8. Note that I have removed the beam and ground
    reflectance portion of the equation and this generates ONLY the diffuse
    radiation from the sky and circumsolar, so the form of the equation
    varies slightly from equation 8.

    References
    ----------

    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
    solar irradiance on inclined surfaces for building energy simulation"
    2007, Solar Energy vol. 81. pp. 254-267

    [2] Reindl, D.T., Beckmann, W.A., Duffie, J.A., 1990a. Diffuse fraction
    correlations. Solar Energy 45(1), 1-7.

    [3] Reindl, D.T., Beckmann, W.A., Duffie, J.A., 1990b. Evaluation of hourly
    tilted surface radiation models. Solar Energy 45(1), 9-17.
    '''

    pvl_logger.debug('diffuse_sky.reindl()')

    cos_tt = aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)

    cos_sun_zen = tools.cosd(sun_zen)

    # ratio of titled and horizontal beam irradiance
    Rb = cos_tt / cos_sun_zen

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

def haydavies(surf_tilt, surf_az, DHI, DNI, DNI_ET, sun_zen, sun_az):
    r'''
    Determine diffuse irradiance from the sky on a
    tilted surface using Hay & Davies' 1980 model

    .. math::
        I_{d} = DHI ( A R_b + (1 - A) (\frac{1 + \cos\beta}{2}) )

    Hay and Davies' 1980 model determines the diffuse irradiance from the sky
    (ground reflected irradiance is not included in this algorithm) on a
    tilted surface using the surface tilt angle, surface azimuth angle,
    diffuse horizontal irradiance, direct normal irradiance,
    extraterrestrial irradiance, sun zenith angle, and sun azimuth angle.

    Parameters
    ----------

    surf_tilt : float or Series
        Surface tilt angles in decimal degrees.
        The tilt angle is defined as
        degrees from horizontal (e.g. surface facing up = 0, surface facing
        horizon = 90)

    surf_az : float or Series
          Surface azimuth angles in decimal degrees.
          The Azimuth convention is defined
          as degrees east of north (e.g. North = 0, South=180 East = 90,
          West = 270).

    DHI : float or Series
          diffuse horizontal irradiance in W/m^2.

    DNI : float or Series
          direct normal irradiance in W/m^2.

    DNI_ET : float or Series
          extraterrestrial normal irradiance in W/m^2.

    sun_zen : float or Series
          apparent (refraction-corrected) zenith
          angles in decimal degrees.

    sun_az : float or Series
          Sun azimuth angles in decimal degrees.
          The Azimuth convention is defined
          as degrees east of north (e.g. North = 0, East = 90, West = 270).

    Returns
    --------

    SkyDiffuse : float or Series

          the diffuse component of the solar radiation  on an
          arbitrarily tilted surface defined by the Perez model as given in
          reference [3].
          SkyDiffuse is the diffuse component ONLY and does not include the
          ground reflected irradiance or the irradiance due to the beam.

    References
    -----------
    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
    solar irradiance on inclined surfaces for building energy simulation"
    2007, Solar Energy vol. 81. pp. 254-267

    [2] Hay, J.E., Davies, J.A., 1980. Calculations of the solar radiation
    incident on an inclined surface. In: Hay, J.E., Won, T.K. (Eds.), Proc. of
    First Canadian Solar Radiation Data Workshop, 59. Ministry of Supply
    and Services, Canada.
    '''

    pvl_logger.debug('diffuse_sky.haydavies()')

    cos_tt = aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)

    cos_sun_zen = tools.cosd(sun_zen)

    # ratio of titled and horizontal beam irradiance
    Rb = cos_tt / cos_sun_zen

    # Anisotropy Index
    AI = DNI / DNI_ET

    # these are actually the () and [] sub-terms of the second term of eqn 7
    term1 = 1 - AI
    term2 = 0.5 * (1 + tools.cosd(surf_tilt))

    sky_diffuse = DHI * (AI * Rb + term1 * term2)
    sky_diffuse[sky_diffuse < 0] = 0

    return sky_diffuse

def klucher(surf_tilt, surf_az, DHI, GHI, sun_zen, sun_az):
    r'''
    Determine diffuse irradiance from the sky on a tilted surface
    using Klucher's 1979 model

    .. math::

       I_{d} = DHI \frac{1 + \cos\beta}{2} (1 + F' \sin^3(\beta/2))
       (1 + F' \cos^2\theta\sin^3\theta_z)

    where

    .. math::

       F' = 1 - (I_{d0} / GHI)

    Klucher's 1979 model determines the diffuse irradiance from the sky
    (ground reflected irradiance is not included in this algorithm) on a
    tilted surface using the surface tilt angle, surface azimuth angle,
    diffuse horizontal irradiance, direct normal irradiance, global
    horizontal irradiance, extraterrestrial irradiance, sun zenith angle,
    and sun azimuth angle.

    Parameters
    ----------

    surf_tilt : float or Series
        Surface tilt angles in decimal degrees.
        surf_tilt must be >=0 and <=180. The tilt angle is defined as
        degrees from horizontal (e.g. surface facing up = 0, surface facing
        horizon = 90)

    surf_az : float or Series
        Surface azimuth angles in decimal degrees.
        surf_az must be >=0 and <=360. The Azimuth convention is defined
        as degrees east of north (e.g. North = 0, South=180 East = 90,
        West = 270).

    DHI : float or Series
        diffuse horizontal irradiance in W/m^2.
        DHI must be >=0.

    GHI : float or Series
        Global  irradiance in W/m^2.
        DNI must be >=0.

    sun_zen : float or Series
        apparent (refraction-corrected) zenith
        angles in decimal degrees.
        sun_zen must be >=0 and <=180.

    sun_az : float or Series
        Sun azimuth angles in decimal degrees.
        sun_az must be >=0 and <=360. The Azimuth convention is defined
        as degrees east of north (e.g. North = 0, East = 90, West = 270).

    Returns
    -------
    float or Series.

    The diffuse component of the solar radiation on an
    arbitrarily tilted surface defined by the Klucher model as given in
    Loutzenhiser et. al (2007) equation 4.
    SkyDiffuse is the diffuse component ONLY and does not include the ground
    reflected irradiance or the irradiance due to the beam.
    SkyDiffuse is a column vector vector with a number of elements equal to
    the input vector(s).

    References
    ----------
    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
    solar irradiance on inclined surfaces for building energy simulation"
    2007, Solar Energy vol. 81. pp. 254-267

    [2] Klucher, T.M., 1979. Evaluation of models to predict insolation on
    tilted surfaces. Solar Energy 23 (2), 111-114.
    '''

    pvl_logger.debug('diffuse_sky.klucher()')

    # zenith angle with respect to panel normal.
    cos_tt = aoi_projection(surf_tilt, surf_az, sun_zen, sun_az)

    F = 1 - ((DHI / GHI) ** 2)
    try:
        # fails with single point input
        F.fillna(0, inplace=True)
    except AttributeError:
        F = 0

    term1 = 0.5 * (1 + tools.cosd(surf_tilt))
    term2 = 1 + F * (tools.sind(0.5 * surf_tilt) ** 3)
    term3 = 1 + F * (cos_tt ** 2) * (tools.sind(sun_zen) ** 3)

    sky_diffuse = DHI * term1 * term2 * term3

    return sky_diffuse

def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
              solar_zenith=None, solar_azimuth=None, projection_ratio=None):
    r'''
    Determine diffuse irradiance from the sky on a
    tilted surface using Hay & Davies' 1980 model

    .. math::
        I_{d} = DHI ( A R_b + (1 - A) (\frac{1 + \cos\beta}{2}) )

    Hay and Davies' 1980 model determines the diffuse irradiance from the sky
    (ground reflected irradiance is not included in this algorithm) on a
    tilted surface using the surface tilt angle, surface azimuth angle,
    diffuse horizontal irradiance, direct normal irradiance,
    extraterrestrial irradiance, sun zenith angle, and sun azimuth angle.

    Parameters
    ----------

    surface_tilt : float or Series
        Surface tilt angles in decimal degrees.
        The tilt angle is defined as
        degrees from horizontal (e.g. surface facing up = 0, surface facing
        horizon = 90)

    surface_azimuth : float or Series
        Surface azimuth angles in decimal degrees.
        The azimuth convention is defined
        as degrees east of north (e.g. North=0, South=180, East=90,
        West=270).

    dhi : float