irradiance.py

"""
The ``irradiance`` module contains functions for modeling global
horizontal irradiance, direct normal irradiance, diffuse horizontal
irradiance, and total irradiance under various conditions.
"""

import datetime
from collections import OrderedDict
from functools import partial

import numpy as np
import pandas as pd

from pvlib import atmosphere, solarposition, tools

# see References section of grounddiffuse function
SURFACE_ALBEDOS = {'urban': 0.18,
                   'grass': 0.20,
                   'fresh grass': 0.26,
                   'soil': 0.17,
                   'sand': 0.40,
                   'snow': 0.65,
                   'fresh snow': 0.75,
                   'asphalt': 0.12,
                   'concrete': 0.30,
                   'aluminum': 0.85,
                   'copper': 0.74,
                   'fresh steel': 0.35,
                   'dirty steel': 0.08,
                   'sea': 0.06}


def get_extra_radiation(datetime_or_doy, solar_constant=1366.1,
                        method='spencer', epoch_year=2014, **kwargs):
    """
    Determine extraterrestrial radiation from day of year.

    Parameters
    ----------
    datetime_or_doy : numeric, array, date, datetime, Timestamp, DatetimeIndex
        Day of year, array of days of year, or datetime-like object

    solar_constant : float, default 1366.1
        The solar constant.

    method : string, default 'spencer'
        The method by which the ET radiation should be calculated.
        Options include ``'pyephem', 'spencer', 'asce', 'nrel'``.

    epoch_year : int, default 2014
        The year in which a day of year input will be calculated. Only
        applies to day of year input used with the pyephem or nrel
        methods.

    kwargs :
        Passed to solarposition.nrel_earthsun_distance

    Returns
    -------
    dni_extra : float, array, or Series
        The extraterrestrial radiation present in watts per square meter
        on a surface which is normal to the sun. Pandas Timestamp and
        DatetimeIndex inputs will yield a Pandas TimeSeries. All other
        inputs will yield a float or an array of floats.

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

    .. [2] <http://solardat.uoregon.edu/SolarRadiationBasics.html>, Eqs.
       SR1 and SR2

    .. [3] Partridge, G. W. and Platt, C. M. R. 1976. Radiative Processes
       in Meteorology and Climatology.

    .. [4] Duffie, J. A. and Beckman, W. A. 1991. Solar Engineering of
       Thermal Processes, 2nd edn. J. Wiley and Sons, New York.

    .. [5] ASCE, 2005. The ASCE Standardized Reference Evapotranspiration
       Equation, Environmental and Water Resources Institute of the American
       Civil Engineers, Ed. R. G. Allen et al.
    """

    to_doy, to_datetimeindex, to_output = \
        _handle_extra_radiation_types(datetime_or_doy, epoch_year)

    # consider putting asce and spencer methods in their own functions
    method = method.lower()
    if method == 'asce':
        B = solarposition._calculate_simple_day_angle(to_doy(datetime_or_doy),
                                                      offset=0)
        RoverR0sqrd = 1 + 0.033 * np.cos(B)
    elif method == 'spencer':
        B = solarposition._calculate_simple_day_angle(to_doy(datetime_or_doy))
        RoverR0sqrd = (1.00011 + 0.034221 * np.cos(B) + 0.00128 * np.sin(B) +
                       0.000719 * np.cos(2 * B) + 7.7e-05 * np.sin(2 * B))
    elif method == 'pyephem':
        times = to_datetimeindex(datetime_or_doy)
        RoverR0sqrd = solarposition.pyephem_earthsun_distance(times) ** (-2)
    elif method == 'nrel':
        times = to_datetimeindex(datetime_or_doy)
        RoverR0sqrd = \
            solarposition.nrel_earthsun_distance(times, **kwargs) ** (-2)
    else:
        raise ValueError('Invalid method: %s', method)

    Ea = solar_constant * RoverR0sqrd

    Ea = to_output(Ea)

    return Ea


def _handle_extra_radiation_types(datetime_or_doy, epoch_year):
    # This block will set the functions that can be used to convert the
    # inputs to either day of year or pandas DatetimeIndex, and the
    # functions that will yield the appropriate output type. It's
    # complicated because there are many day-of-year-like input types,
    # and the different algorithms need different types. Maybe you have
    # a better way to do it.
    if isinstance(datetime_or_doy, pd.DatetimeIndex):
        to_doy = tools._pandas_to_doy  # won't be evaluated unless necessary
        def to_datetimeindex(x): return x                       # noqa: E306
        to_output = partial(pd.Series, index=datetime_or_doy)
    elif isinstance(datetime_or_doy, pd.Timestamp):
        to_doy = tools._pandas_to_doy
        to_datetimeindex = \
            tools._datetimelike_scalar_to_datetimeindex
        to_output = tools._scalar_out
    elif isinstance(datetime_or_doy,
                    (datetime.date, datetime.datetime, np.datetime64)):
        to_doy = tools._datetimelike_scalar_to_doy
        to_datetimeindex = \
            tools._datetimelike_scalar_to_datetimeindex
        to_output = tools._scalar_out
    elif np.isscalar(datetime_or_doy):  # ints and floats of various types
        def to_doy(x): return x                                 # noqa: E306
        to_datetimeindex = partial(tools._doy_to_datetimeindex,
                                   epoch_year=epoch_year)
        to_output = tools._scalar_out
    else:  # assume that we have an array-like object of doy
        def to_doy(x): return x                                 # noqa: E306
        to_datetimeindex = partial(tools._doy_to_datetimeindex,
                                   epoch_year=epoch_year)
        to_output = tools._array_out

    return to_doy, to_datetimeindex, to_output


def aoi_projection(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    """
    Calculates the dot product of the sun position unit vector and the surface
    normal unit vector; in other words, the cosine of the angle of incidence.

    Usage note: When the sun is behind the surface the value returned is
    negative.  For many uses negative values must be set to zero.

    Input all angles in degrees.

    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.

    Returns
    -------
    projection : numeric
        Dot product of panel normal and solar angle.
    """

    projection = (
        tools.cosd(surface_tilt) * tools.cosd(solar_zenith) +
        tools.sind(surface_tilt) * tools.sind(solar_zenith) *
        tools.cosd(solar_azimuth - surface_azimuth))

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

    return projection


def aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    """
    Calculates the angle of incidence of the solar vector on a surface.
    This is the angle between the solar vector and the surface normal.

    Input all angles in degrees.

    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.

    Returns
    -------
    aoi : numeric
        Angle of incidence in degrees.
    """

    projection = aoi_projection(surface_tilt, surface_azimuth,
                                solar_zenith, solar_azimuth)
    aoi_value = np.rad2deg(np.arccos(projection))

    try:
        aoi_value.name = 'aoi'
    except AttributeError:
        pass

    return aoi_value


def poa_horizontal_ratio(surface_tilt, surface_azimuth,
                         solar_zenith, solar_azimuth):
    """
    Calculates the ratio of the beam components of the plane of array
    irradiance and the horizontal irradiance.

    Input all angles in degrees.

    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.

    Returns
    -------
    ratio : numeric
        Ratio of the plane of array irradiance to the horizontal plane
        irradiance
    """

    cos_poa_zen = aoi_projection(surface_tilt, surface_azimuth,
                                 solar_zenith, solar_azimuth)

    cos_solar_zenith = tools.cosd(solar_zenith)

    # ratio of tilted and horizontal beam irradiance
    ratio = cos_poa_zen / cos_solar_zenith

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

    return ratio


def beam_component(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                   dni):
    """
    Calculates the beam component of the plane of array irradiance.

    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 : numeric
        Direct Normal Irradiance

    Returns
    -------
    beam : numeric
        Beam component
    """
    beam = dni * aoi_projection(surface_tilt, surface_azimuth,
                                solar_zenith, solar_azimuth)
    beam = np.maximum(beam, 0)

    return beam


def get_total_irradiance(surface_tilt, surface_azimuth,
                         solar_zenith, solar_azimuth,
                         dni, ghi, dhi, dni_extra=None, airmass=None,
                         albedo=.25, surface_type=None,
                         model='isotropic',
                         model_perez='allsitescomposite1990', **kwargs):
    r"""
    Determine total in-plane irradiance and its beam, sky diffuse and ground
    reflected components, using the specified sky diffuse irradiance model.

    .. math::

       I_{tot} = I_{beam} + I_{sky diffuse} + I_{ground}

    Sky diffuse models include:
        * isotropic (default)
        * klucher
        * haydavies
        * reindl
        * king
        * perez

    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 : numeric
        Direct Normal Irradiance
    ghi : numeric
        Global horizontal irradiance
    dhi : numeric
        Diffuse horizontal irradiance
    dni_extra : None or numeric, default None
        Extraterrestrial direct normal irradiance
    airmass : None or numeric, default None
        Airmass
    albedo : numeric, default 0.25
        Surface albedo
    surface_type : None or String, default None
        Surface type. See grounddiffuse.
    model : String, default 'isotropic'
        Irradiance model.
    model_perez : String, default 'allsitescomposite1990'
        Used only if model='perez'. See :py:func:`perez`.

    Returns
    -------
    total_irrad : OrderedDict or DataFrame
        Contains keys/columns ``'poa_global', 'poa_direct', 'poa_diffuse',
        'poa_sky_diffuse', 'poa_ground_diffuse'``.
    """
    poa_sky_diffuse = get_sky_diffuse(
        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
        dni, ghi, dhi, dni_extra=dni_extra, airmass=airmass, model=model,
        model_perez=model_perez)

    poa_ground_diffuse = get_ground_diffuse(surface_tilt, ghi, albedo,
                                            surface_type)
    aoi_ = aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)
    irrads = poa_components(aoi_, dni, poa_sky_diffuse, poa_ground_diffuse)
    return irrads


def get_sky_diffuse(surface_tilt, surface_azimuth,
                    solar_zenith, solar_azimuth,
                    dni, ghi, dhi, dni_extra=None, airmass=None,
                    model='isotropic',
                    model_perez='allsitescomposite1990'):
    r"""
    Determine in-plane sky diffuse irradiance component
    using the specified sky diffuse irradiance model.

    Sky diffuse models include:
        * isotropic (default)
        * klucher
        * haydavies
        * reindl
        * king
        * perez

    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 : numeric
        Direct Normal Irradiance
    ghi : numeric
        Global horizontal irradiance
    dhi : numeric
        Diffuse horizontal irradiance
    dni_extra : None or numeric, default None
        Extraterrestrial direct normal irradiance
    airmass : None or numeric, default None
        Airmass
    model : String, default 'isotropic'
        Irradiance model.
    model_perez : String, default 'allsitescomposite1990'
        See perez.

    Returns
    -------
    poa_sky_diffuse : numeric
    """

    model = model.lower()
    if model == 'isotropic':
        sky = isotropic(surface_tilt, dhi)
    elif model == 'klucher':
        sky = klucher(surface_tilt, surface_azimuth, dhi, ghi,
                      solar_zenith, solar_azimuth)
    elif model == 'haydavies':
        sky = haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
                        solar_zenith, solar_azimuth)
    elif model == 'reindl':
        sky = reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
                     solar_zenith, solar_azimuth)
    elif model == 'king':
        sky = king(surface_tilt, dhi, ghi, solar_zenith)
    elif model == 'perez':
        sky = perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
                    solar_zenith, solar_azimuth, airmass,
                    model=model_perez)
    else:
        raise ValueError('invalid model selection {}'.format(model))

    return sky


def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
    r'''
    Determine in-plane irradiance components.

    Combines DNI with sky diffuse and ground-reflected irradiance to calculate
    total, direct and diffuse irradiance components in the plane of array.

    Parameters
    ----------
    aoi : numeric
        Angle of incidence of solar rays with respect to the module
        surface, from :func:`aoi`.

    dni : numeric
        Direct normal irradiance (W/m^2), as measured from a TMY file or
        calculated with a clearsky model.

    poa_sky_diffuse : numeric
        Diffuse irradiance (W/m^2) in the plane of the modules, as
        calculated by a diffuse irradiance translation function

    poa_ground_diffuse : numeric
        Ground reflected irradiance (W/m^2) in the plane of the modules,
        as calculated by an albedo model (eg. :func:`grounddiffuse`)

    Returns
    -------
    irrads : OrderedDict or DataFrame
        Contains the following keys:

        * ``poa_global`` : Total in-plane irradiance (W/m^2)
        * ``poa_direct`` : Total in-plane beam irradiance (W/m^2)
        * ``poa_diffuse`` : Total in-plane diffuse irradiance (W/m^2)
        * ``poa_sky_diffuse`` : In-plane diffuse irradiance from sky (W/m^2)
        * ``poa_ground_diffuse`` : In-plane diffuse irradiance from ground
          (W/m^2)

    Notes
    ------
    Negative beam irradiation due to aoi :math:`> 90^{\circ}` or AOI
    :math:`< 0^{\circ}` is set to zero.
    '''

    poa_direct = np.maximum(dni * np.cos(np.radians(aoi)), 0)
    poa_diffuse = poa_sky_diffuse + poa_ground_diffuse
    poa_global = poa_direct + poa_diffuse

    irrads = OrderedDict()
    irrads['poa_global'] = poa_global
    irrads['poa_direct'] = poa_direct
    irrads['poa_diffuse'] = poa_diffuse
    irrads['poa_sky_diffuse'] = poa_sky_diffuse
    irrads['poa_ground_diffuse'] = poa_ground_diffuse

    if isinstance(poa_direct, pd.Series):
        irrads = pd.DataFrame(irrads)

    return irrads


def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None):
    '''
    Estimate diffuse irradiance from ground reflections given
    irradiance, albedo, and surface tilt

    Function to determine the portion of irradiance on a tilted surface
    due to ground reflections. Any of the inputs may be DataFrames or
    scalars.

    Parameters
    ----------
    surface_tilt : numeric
        Surface tilt angles in decimal degrees. 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).

    ghi : numeric
        Global horizontal irradiance in W/m^2.

    albedo : numeric, default 0.25
        Ground reflectance, typically 0.1-0.4 for surfaces on Earth
        (land), may increase over snow, ice, etc. May also be known as
        the reflection coefficient. Must be >=0 and <=1. Will be
        overridden if surface_type is supplied.

    surface_type: None or string, default None
        If not None, overrides albedo. String can be one of 'urban',
        'grass', 'fresh grass', 'snow', 'fresh snow', 'asphalt', 'concrete',
        'aluminum', 'copper', 'fresh steel', 'dirty steel', 'sea'.

    Returns
    -------
    grounddiffuse : numeric
        Ground reflected irradiances in W/m^2.


    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.

    The calculation is the last term of equations 3, 4, 7, 8, 10, 11, and 12.

    .. [2] albedos from:
       http://files.pvsyst.com/help/albedo.htm
       and
       http://en.wikipedia.org/wiki/Albedo
       and
       https://doi.org/10.1175/1520-0469(1972)029<0959:AOTSS>2.0.CO;2
    '''

    if surface_type is not None:
        albedo = SURFACE_ALBEDOS[surface_type]

    diffuse_irrad = ghi * albedo * (1 - np.cos(np.radians(surface_tilt))) * 0.5

    try:
        diffuse_irrad.name = 'diffuse_ground'
    except AttributeError:
        pass

    return diffuse_irrad


def isotropic(surface_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
    ----------
    surface_tilt : numeric
        Surface tilt angle in decimal degrees. 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 : numeric
        Diffuse horizontal irradiance in W/m^2. DHI must be >=0.

    Returns
    -------
    diffuse : numeric
        The sky diffuse component of the solar radiation.

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

    sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5

    return sky_diffuse


def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,
            solar_azimuth):
    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)^2

    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
    ----------
    surface_tilt : numeric
        Surface tilt angles in decimal degrees. surface_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)

    surface_azimuth : numeric
        Surface azimuth angles in decimal degrees. surface_azimuth 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 : numeric
        Diffuse horizontal irradiance in W/m^2. DHI must be >=0.

    ghi : numeric
        Global irradiance in W/m^2. DNI must be >=0.

    solar_zenith : numeric
        Apparent (refraction-corrected) zenith angles in decimal
        degrees. solar_zenith must be >=0 and <=180.

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

    Returns
    -------
    diffuse : numeric
        The sky diffuse component of the solar radiation.

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

    # zenith angle with respect to panel normal.
    cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                            solar_zenith, solar_azimuth)
    cos_tt = np.maximum(cos_tt, 0)  # GH 526

    # silence warning from 0 / 0
    with np.errstate(invalid='ignore'):
        F = 1 - ((dhi / ghi) ** 2)

    try:
        # fails with single point input
        F.fillna(0, inplace=True)
    except AttributeError:
        F = np.where(np.isnan(F), 0, F)

    term1 = 0.5 * (1 + tools.cosd(surface_tilt))
    term2 = 1 + F * (tools.sind(0.5 * surface_tilt) ** 3)
    term3 = 1 + F * (cos_tt ** 2) * (tools.sind(solar_zenith) ** 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 : numeric
        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 : numeric
        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 : numeric
        Diffuse horizontal irradiance in W/m^2.

    dni : numeric
        Direct normal irradiance in W/m^2.

    dni_extra : numeric
        Extraterrestrial normal irradiance in W/m^2.

    solar_zenith : None or numeric, default None
        Solar apparent (refraction-corrected) zenith angles in decimal
        degrees. Must supply ``solar_zenith`` and ``solar_azimuth`` or
        supply ``projection_ratio``.

    solar_azimuth : None or numeric, default None
        Solar azimuth angles in decimal degrees. Must supply
        ``solar_zenith`` and ``solar_azimuth`` or supply
        ``projection_ratio``.

    projection_ratio : None or numeric, default None
        Ratio of angle of incidence projection to solar zenith angle
        projection. Must supply ``solar_zenith`` and ``solar_azimuth``
        or supply ``projection_ratio``.

    Returns
    --------
    sky_diffuse : numeric
        The sky diffuse component of the solar radiation.

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

    # if necessary, calculate ratio of titled and horizontal beam irradiance
    if projection_ratio is None:
        cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                                solar_zenith, solar_azimuth)
        cos_tt = np.maximum(cos_tt, 0)  # GH 526
        cos_solar_zenith = tools.cosd(solar_zenith)
        Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745)  # GH 432
    else:
        Rb = projection_ratio

    # Anisotropy Index
    AI = dni / dni_extra

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

    sky_diffuse = dhi * (AI * Rb + term1 * term2)
    sky_diffuse = np.maximum(sky_diffuse, 0)

    return sky_diffuse


def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
           solar_zenith, solar_azimuth):
    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
    ----------
    surface_tilt : numeric
        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 : numeric
        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 : numeric
        diffuse horizontal irradiance in W/m^2.

    dni : numeric
        direct normal irradiance in W/m^2.

    ghi: numeric
        Global irradiance in W/m^2.

    dni_extra : numeric
        Extraterrestrial normal irradiance in W/m^2.

    solar_zenith : numeric
        Apparent (refraction-corrected) zenith angles in decimal degrees.

    solar_azimuth : numeric
        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
    -------
    poa_sky_diffuse : numeric
        The sky diffuse component of the solar radiation.

    Notes
    -----
    The poa_sky_diffuse 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.
    '''

    cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                            solar_zenith, solar_azimuth)
    cos_tt = np.maximum(cos_tt, 0)  # GH 526

    # do not apply cos(zen) limit here (needed for HB below)
    cos_solar_zenith = tools.cosd(solar_zenith)

    # ratio of titled and horizontal beam irradiance
    Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745)  # GH 432

    # Anisotropy Index
    AI = dni / dni_extra

    # DNI projected onto horizontal
    HB = dni * cos_solar_zenith
    HB = np.maximum(HB, 0)

    # these are the () and [] sub-terms of the second term of eqn 8
    term1 = 1 - AI
    term2 = 0.5 * (1 + tools.cosd(surface_tilt))
    term3 = 1 + np.sqrt(HB / ghi) * (tools.sind(0.5 * surface_tilt) ** 3)

    sky_diffuse = dhi * (AI * Rb + term1 * term2 * term3)
    sky_diffuse = np.maximum(sky_diffuse, 0)

    return sky_diffuse


def king(surface_tilt, dhi, ghi, solar_zenith):
    '''
    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
    ----------
    surface_tilt : numeric
        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 : numeric
        Diffuse horizontal irradiance in W/m^2.

    ghi : numeric
        Global horizontal irradiance in W/m^2.

    solar_zenith : numeric
        Apparent (refraction-corrected) zenith angles in decimal degrees.

    Returns
    --------
    poa_sky_diffuse : numeric
        The diffuse component of the solar radiation.
    '''

    sky_diffuse = (dhi * ((1 + tools.cosd(surface_tilt))) / 2 + ghi *
                   ((0.012 * solar_zenith - 0.04)) *
                   ((1 - tools.cosd(surface_tilt))) / 2)
    sky_diffuse = np.maximum(sky_diffuse, 0)

    return sky_diffuse


def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
          solar_zenith, solar_azimuth, airmass,
          model='allsitescomposite1990', return_components=False):
    '''
    Determine diffuse irradiance from the sky on a tilted surface using
    one of the Perez models.

    Perez models determine 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, sun azimuth angle, and relative (not
    pressure-corrected) airmass. Optionally a selector may be used to
    use any of Perez's model coefficient sets.

    Parameters
    ----------
    surface_tilt : numeric
        Surface tilt angles in decimal degrees. surface_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)

    surface_azimuth : numeric
        Surface azimuth angles in decimal degrees. surface_azimuth 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 : numeric
        Diffuse horizontal irradiance in W/m^2. DHI must be >=0.

    dni : numeric
        Direct normal irradiance in W/m^2. DNI must be >=0.

    dni_extra : numeric
        Extraterrestrial normal irradiance in W/m^2.

    solar_zenith : numeric
        apparent (refraction-corrected) zenith angles in decimal
        degrees. solar_zenith must be >=0 and <=180.

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

    airmass : numeric
        Relative (not pressure-corrected) airmass values. If AM is a
        DataFrame it must be of the same size as all other DataFrame
        inputs. AM must be >=0 (careful using the 1/sec(z) model of AM
        generation)

    model : string (optional, default='allsitescomposite1990')
        A string which selects the desired set of Perez coefficients. If
        model is not provided as an input, the default, '1990' will be
        used. All possible model selections are:

        * '1990'
        * 'allsitescomposite1990' (same as '1990')
        * 'allsitescomposite1988'
        * 'sandiacomposite1988'
        * 'usacomposite1988'
        * 'france1988'
        * 'phoenix1988'
        * 'elmonte1988'
        * 'osage1988'
        * 'albuquerque1988'
        * 'capecanaveral1988'
        * 'albany1988'

    return_components: bool (optional, default=False)
        Flag used to decide whether to return the calculated diffuse components
        or not.

    Returns
    --------
    numeric, OrderedDict, or DataFrame
        Return type controlled by `return_components` argument.
        If ``return_components=False``, `sky_diffuse` is returned.
        If ``return_components=True``, `diffuse_components` is returned.

    sky_diffuse : numeric
        The sky diffuse component of the solar radiation on a tilted
        surface.

    diffuse_components : OrderedDict (array input) or DataFrame (Series input)
        Keys/columns are:
            * sky_diffuse: Total sky diffuse
            * isotropic
            * circumsolar
            * horizon


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

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

    # delta is the sky's "brightness"
    delta = dhi * airmass / dni_extra

    # epsilon is the sky's "clearness"
    with np.errstate(invalid='ignore'):
        eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))

    # numpy indexing below will not work with a Series
    if isinstance(eps, pd.Series):
        eps = eps.values

    # 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 = np.digitize(eps, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
    ebin = np.array(ebin)  # GH 642
    ebin[np.isnan(eps)] = 0

    # correct for 0 indexing in coeffecient lookup
    # later, ebin = -1 will yield nan coefficients
    ebin -= 1

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

    # results in invalid eps (ebin = -1) being mapped to nans
    nans = np.array([np.nan, np.nan, np.nan])
    F1c = np.vstack((F1c, nans))
    F2c = np.vstack((F2c, nans))

    F1 = (F1c[ebin, 0] + F1c[ebin, 1] * delta + F1c[ebin, 2] * z)
    F1 = np.maximum(F1, 0)

    F2 = (F2c[ebin, 0] + F2c[ebin, 1] * delta + F2c[ebin, 2] * z)
    F2 = np.maximum(F2, 0)

    A = aoi_projection(surface_tilt, surface_azimuth,
                       solar_zenith, solar_azimuth)
    A = np.maximum(A, 0)

    B = tools.cosd(solar_zenith)
    B = np.maximum(B, tools.cosd(85))

    # Calculate Diffuse POA from sky dome
    term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
    term2 = F1 * A / B
    term3 = F2 * tools.sind(surface_tilt)

    sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)

    # we've preserved the input type until now, so don't ruin it!
    if isinstance(sky_diffuse, pd.Series):
        sky_diffuse[np.isnan(airmass)] = 0
    else:
        sky_diffuse = np.where(np.isnan(airmass), 0, sky_diffuse)

    if return_components:
        diffuse_components = OrderedDict()
        diffuse_components['sky_diffuse'] = sky_diffuse

        # Calculate the different components
        diffuse_components['isotropic'] = dhi * term1
        diffuse_components['circumsolar'] = dhi * term2
        diffuse_components['horizon'] = dhi * term3

        # Set values of components to 0 when sky_diffuse is 0
        mask = sky_diffuse == 0
        if isinstance(sky_diffuse, pd.Series):
            diffuse_components = pd.DataFrame(diffuse_components)
            diffuse_components.loc[mask] = 0
        else:
            diffuse_components = {k: np.where(mask, 0, v) for k, v in
                                  diffuse_components.items()}
        return diffuse_components
    else:
        return sky_diffuse


def clearsky_index(ghi, clearsky_ghi, max_clearsky_index=2.0):
    """
    Calculate the clearsky index.

    The clearsky index is the ratio of global to clearsky global irradiance.
    Negative and non-finite clearsky index values will be truncated to zero.

    Parameters
    ----------
    ghi : numeric
        Global horizontal irradiance in W/m^2.

    clearsky_ghi : numeric
        Modeled clearsky GHI

    max_clearsky_index : numeric, default 2.0
        Maximum value of the clearsky index. The default, 2.0, allows
        for over-irradiance events typically seen in sub-hourly data.

    Returns
    -------
    clearsky_index : numeric
        Clearsky index
    """
    clearsky_index = ghi / clearsky_ghi
    # set +inf, -inf, and nans to zero
    clearsky_index = np.where(~np.isfinite(clearsky_index), 0,
                              clearsky_index)
    # but preserve nans in the input arrays
    input_is_nan = ~np.isfinite(ghi) | ~np.isfinite(clearsky_ghi)
    clearsky_index = np.where(input_is_nan, np.nan, clearsky_index)

    clearsky_index = np.maximum(clearsky_index, 0)
    clearsky_index = np.minimum(clearsky_index, max_clearsky_index)

    # preserve input type
    if isinstance(ghi, pd.Series):
        clearsky_index = pd.Series(clearsky_index, index=ghi.index)

    return clearsky_index


def clearness_index(ghi, solar_zenith, extra_radiation, min_cos_zenith=0.065,
                    max_clearness_index=2.0):
    """
    Calculate the clearness index.

    The clearness index is the ratio of global to extraterrestrial
    irradiance on a horizontal plane [1]_.

    Parameters
    ----------
    ghi : numeric
        Global horizontal irradiance in W/m^2.

    solar_zenith : numeric
        True (not refraction-corrected) solar zenith angle in decimal
        degrees.

    extra_radiation : numeric
        Irradiance incident at the top of the atmosphere

    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_clearness_index : numeric, default 2.0
        Maximum value of the clearness index. The default, 2.0, allows
        for over-irradiance events typically seen in sub-hourly data.
        NREL's SRRL Fortran code used 0.82 for hourly data.

    Returns
    -------
    kt : numeric
        Clearness index

    References
    ----------
    .. [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
           Global Horizontal to Direct Normal Insolation", Technical
           Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
           Institute, 1987.
    """
    cos_zenith = tools.cosd(solar_zenith)
    I0h = extra_radiation * np.maximum(cos_zenith, min_cos_zenith)
    # consider adding
    # with np.errstate(invalid='ignore', divide='ignore'):
    # to kt calculation, but perhaps it's good to allow these
    # warnings to the users that override min_cos_zenith
    kt = ghi / I0h
    kt = np.maximum(kt, 0)
    kt = np.minimum(kt, max_clearness_index)
    return kt


def clearness_index_zenith_independent(clearness_index, airmass,
                                       max_clearness_index=2.0):
    """
    Calculate the zenith angle independent clearness index.

    See [1]_ for details.

    Parameters
    ----------
    clearness_index : numeric
        Ratio of global to extraterrestrial irradiance on a horizontal
        plane

    airmass : numeric
        Airmass

    max_clearness_index : numeric, default 2.0
        Maximum value of the clearness index. The default, 2.0, allows
        for over-irradiance events typically seen in sub-hourly data.
        NREL's SRRL Fortran code used 0.82 for hourly data.

    Returns
    -------
    kt_prime : numeric
        Zenith independent clearness index

    References
    ----------
    .. [1] Perez, R., P. Ineichen, E. Maxwell, R. Seals and A. Zelenka,
           (1992). "Dynamic Global-to-Direct Irradiance Conversion Models".
           ASHRAE Transactions-Research Series, pp. 354-369
    """
    # Perez eqn 1
    kt_prime = clearness_index / _kt_kt_prime_factor(airmass)
    kt_prime = np.maximum(kt_prime, 0)
    kt_prime = np.minimum(kt_prime, max_clearness_index)
    return kt_prime


def _kt_kt_prime_factor(airmass):
    """
    Calculate the conversion factor between kt and kt prime.
    Function is useful because DIRINT and GTI-DIRINT both use this.
    """
    # consider adding
    # airmass = np.maximum(airmass, 12)  # GH 450
    return 1.031 * np.exp(-1.4 / (0.9 + 9.4 / airmass)) + 0.1


def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
         min_cos_zenith=0.065, max_zenith=87, max_airmass=12):
    """
    Estimate Direct Normal Irradiance from Global Horizontal Irradiance
    using the DISC model.

    The DISC algorithm converts global horizontal irradiance to direct
    normal irradiance through empirical relationships between the global
    and direct clearness indices.

    The pvlib implementation limits the clearness index to 1.

    The original report describing the DISC model [1]_ uses the
    relative airmass rather than the absolute (pressure-corrected)
    airmass. However, the NREL implementation of the DISC model [2]_
    uses absolute airmass. PVLib Matlab also uses the absolute airmass.
    pvlib python defaults to absolute airmass, but the relative airmass
    can be used by supplying `pressure=None`.

    Parameters
    ----------
    ghi : numeric
        Global horizontal irradiance in W/m^2.

    solar_zenith : numeric
        True (not refraction-corrected) solar 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.

    pressure : None or numeric, default 101325
        Site pressure in Pascal. If None, relative airmass is used
        instead of absolute (pressure-corrected) airmass.

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

    max_airmass : numeric, default 12
        Maximum value of the airmass to allow in Kn calculation.
        Default value (12) comes from range over which Kn was fit
        to airmass in the original paper.

    Returns
    -------
    output : OrderedDict or DataFrame
        Contains the following keys:

        * ``dni``: The modeled direct normal irradiance
          in W/m^2 provided by the
          Direct Insolation Simulation Code (DISC) model.
        * ``kt``: Ratio of global to extraterrestrial
          irradiance on a horizontal plane.
        * ``airmass``: Airmass

    References
    ----------
    .. [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
       Global Horizontal to Direct Normal Insolation", Technical
       Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
       Institute, 1987.

    .. [2] Maxwell, E. "DISC Model", Excel Worksheet.
       https://www.nrel.gov/grid/solar-resource/disc.html

    See Also
    --------
    dirint
    """

    # this is the I0 calculation from the reference
    # SSC uses solar constant = 1367.0 (checked 2018 08 15)
    I0 = get_extra_radiation(datetime_or_doy, 1370., 'spencer')

    kt = clearness_index(ghi, solar_zenith, I0, min_cos_zenith=min_cos_zenith,
                         max_clearness_index=1)

    am = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
    if pressure is not None:
        am = atmosphere.get_absolute_airmass(am, pressure)

    Kn, am = _disc_kn(kt, am, max_airmass=max_airmass)
    dni = Kn * I0

    bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
    dni = np.where(bad_values, 0, dni)

    output = OrderedDict()
    output['dni'] = dni
    output['kt'] = kt
    output['airmass'] = am

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

    return output


def _disc_kn(clearness_index, airmass, max_airmass=12):
    """
    Calculate Kn for `disc`

    Parameters
    ----------
    clearness_index : numeric
    airmass : numeric
    max_airmass : float
        airmass > max_airmass is set to max_airmass before being used
        in calculating Kn.

    Returns
    -------
    Kn : numeric
    am : numeric
        airmass used in the calculation of Kn. am <= max_airmass.
    """
    # short names for equations
    kt = clearness_index
    am = airmass

    am = np.minimum(am, max_airmass)  # GH 450

    # powers of kt will be used repeatedly, so compute only once
    kt2 = kt * kt  # about the same as kt ** 2
    kt3 = kt2 * kt  # 5-10x faster than kt ** 3

    bools = (kt <= 0.6)
    a = np.where(bools,
                 0.512 - 1.56*kt + 2.286*kt2 - 2.222*kt3,
                 -5.743 + 21.77*kt - 27.49*kt2 + 11.56*kt3)
    b = np.where(bools,
                 0.37 + 0.962*kt,
                 41.4 - 118.5*kt + 66.05*kt2 + 31.9*kt3)
    c = np.where(bools,
                 -0.28 + 0.932*kt - 2.048*kt2,
                 -47.01 + 184.2*kt - 222.0*kt2 + 73.81*kt3)

    delta_kn = a + b * np.exp(c*am)

    Knc = 0.866 - 0.122*am + 0.0121*am**2 - 0.000653*am**3 + 1.4e-05*am**4
    Kn = Knc - delta_kn
    return Kn, am


def dirint(ghi, solar_zenith, times, pressure=101325., use_delta_kt_prime=True,
           temp_dew=None, min_cos_zenith=0.065, max_zenith=87):
    """
    Determine DNI from GHI using the DIRINT modification of the DISC
    model.

    Implements the modified DISC model known as "DIRINT" introduced in
    [1]. DIRINT predicts direct normal irradiance (DNI) from measured
    global horizontal irradiance (GHI). DIRINT improves upon the DISC
    model by using time-series GHI data and dew point temperature
    information. The effectiveness of the DIRINT model improves with
    each piece of information provided.

    The pvlib implementation limits the clearness index to 1.

    Parameters
    ----------
    ghi : array-like
        Global horizontal irradiance in W/m^2.

    solar_zenith : array-like
        True (not refraction-corrected) solar_zenith angles in decimal
        degrees.

    times : DatetimeIndex

    pressure : float or array-like, default 101325.0
        The site pressure in Pascal. Pressure may be measured or an
        average pressure may be calculated from site altitude.

    use_delta_kt_prime : bool, default True
        If True, indicates that the stability index delta_kt_prime is
        included in the model. The stability index adjusts the estimated
        DNI in response to dynamics in the time series of GHI. It is
        recommended that delta_kt_prime is not used if the time between
        GHI points is 1.5 hours or greater. If use_delta_kt_prime=True,
        input data must be Series.

    temp_dew : None, float, or array-like, default None
        Surface dew point temperatures, in degrees C. Values of temp_dew
        may be numeric or NaN. Any single time period point with a
        temp_dew=NaN does not have dew point improvements applied. If
        temp_dew is not provided, then dew point improvements are not
        applied.

    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
    -------
    dni : array-like
        The modeled direct normal irradiance in W/m^2 provided by the
        DIRINT model.

    Notes
    -----
    DIRINT model requires time series data (ie. one of the inputs must
    be a vector of length > 2).

    References
    ----------
    .. [1] Perez, R., P. Ineichen, E. Maxwell, R. Seals and A. Zelenka,
       (1992). "Dynamic Global-to-Direct Irradiance Conversion Models".
       ASHRAE Transactions-Research Series, pp. 354-369

    .. [2] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
       Global Horizontal to Direct Normal Insolation", Technical Report No.
       SERI/TR-215-3087, Golden, CO: Solar Energy Research Institute, 1987.
    """

    disc_out = disc(ghi, solar_zenith, times, pressure=pressure,
                    min_cos_zenith=min_cos_zenith, max_zenith=max_zenith)
    airmass = disc_out['airmass']
    kt = disc_out['kt']

    kt_prime = clearness_index_zenith_independent(
        kt, airmass, max_clearness_index=1)
    delta_kt_prime = _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime,
                                            times)
    w = _temp_dew_dirint(temp_dew, times)

    dirint_coeffs = _dirint_coeffs(times, kt_prime, solar_zenith, w,
                                   delta_kt_prime)

    # Perez eqn 5
    dni = disc_out['dni'] * dirint_coeffs

    return dni


def _dirint_from_dni_ktprime(dni, kt_prime, solar_zenith, use_delta_kt_prime,
                             temp_dew):
    """
    Calculate DIRINT DNI from supplied DISC DNI and Kt'.

    Supports :py:func:`gti_dirint`
    """
    times = dni.index
    delta_kt_prime = _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime,
                                            times)
    w = _temp_dew_dirint(temp_dew, times)
    dirint_coeffs = _dirint_coeffs(times, kt_prime, solar_zenith, w,
                                   delta_kt_prime)
    dni_dirint = dni * dirint_coeffs
    return dni_dirint


def _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime, times):
    """
    Calculate delta_kt_prime (Perez eqn 2 and eqn 3), or return a default value
    for use with :py:func:`_dirint_bins`.
    """
    if use_delta_kt_prime:
        # Perez eqn 2
        kt_next = kt_prime.shift(-1)
        kt_previous = kt_prime.shift(1)
        # replace nan with values that implement Perez Eq 3 for first and last
        # positions. Use kt_previous and kt_next to handle series of length 1
        kt_next.iloc[-1] = kt_previous.iloc[-1]
        kt_previous.iloc[0] = kt_next.iloc[0]
        delta_kt_prime = 0.5 * ((kt_prime - kt_next).abs().add(
                                (kt_prime - kt_previous).abs(),
                                fill_value=0))
    else:
        # do not change unless also modifying _dirint_bins
        delta_kt_prime = pd.Series(-1, index=times)
    return delta_kt_prime


def _temp_dew_dirint(temp_dew, times):
    """
    Calculate precipitable water from surface dew point temp (Perez eqn 4),
    or return a default value for use with :py:func:`_dirint_bins`.
    """
    if temp_dew is not None:
        # Perez eqn 4
        w = pd.Series(np.exp(0.07 * temp_dew - 0.075), index=times)
    else:
        # do not change unless also modifying _dirint_bins
        w = pd.Series(-1, index=times)
    return w


def _dirint_coeffs(times, kt_prime, solar_zenith, w, delta_kt_prime):
    """
    Determine the DISC to DIRINT multiplier `dirint_coeffs`.

    dni = disc_out['dni'] * dirint_coeffs

    Parameters
    ----------
    times : pd.DatetimeIndex
    kt_prime : Zenith-independent clearness index
    solar_zenith : Solar zenith angle
    w : precipitable water estimated from surface dew-point temperature
    delta_kt_prime : stability index

    Returns
    -------
    dirint_coeffs : array-like
    """
    kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin = \
        _dirint_bins(times, kt_prime, solar_zenith, w, delta_kt_prime)

    # get the coefficients
    coeffs = _get_dirint_coeffs()

    # subtract 1 to account for difference between MATLAB-style bin
    # assignment and Python-style array lookup.
    dirint_coeffs = coeffs[kt_prime_bin-1, zenith_bin-1,
                           delta_kt_prime_bin-1, w_bin-1]

    # convert unassigned bins to nan
    dirint_coeffs = np.where((kt_prime_bin == 0) | (zenith_bin == 0) |
                             (w_bin == 0) | (delta_kt_prime_bin == 0),
                             np.nan, dirint_coeffs)
    return dirint_coeffs


def _dirint_bins(times, kt_prime, zenith, w, delta_kt_prime):
    """
    Determine the bins for the DIRINT coefficients.

    Parameters
    ----------
    times : pd.DatetimeIndex
    kt_prime : Zenith-independent clearness index
    zenith : Solar zenith angle
    w : precipitable water estimated from surface dew-point temperature
    delta_kt_prime : stability index

    Returns
    -------
    tuple of kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin
    """
    # @wholmgren: the following bin assignments use MATLAB's 1-indexing.
    # Later, we'll subtract 1 to conform to Python's 0-indexing.

    # Create kt_prime bins
    kt_prime_bin = pd.Series(0, index=times, dtype=np.int64)
    kt_prime_bin[(kt_prime >= 0) & (kt_prime < 0.24)] = 1
    kt_prime_bin[(kt_prime >= 0.24) & (kt_prime < 0.4)] = 2
    kt_prime_bin[(kt_prime >= 0.4) & (kt_prime < 0.56)] = 3
    kt_prime_bin[(kt_prime >= 0.56) & (kt_prime < 0.7)] = 4
    kt_prime_bin[(kt_prime >= 0.7) & (kt_prime < 0.8)] = 5
    kt_prime_bin[(kt_prime >= 0.8) & (kt_prime <= 1)] = 6

    # Create zenith angle bins
    zenith_bin = pd.Series(0, index=times, dtype=np.int64)
    zenith_bin[(zenith >= 0) & (zenith < 25)] = 1
    zenith_bin[(zenith >= 25) & (zenith < 40)] = 2
    zenith_bin[(zenith >= 40) & (zenith < 55)] = 3
    zenith_bin[(zenith >= 55) & (zenith < 70)] = 4
    zenith_bin[(zenith >= 70) & (zenith < 80)] = 5
    zenith_bin[(zenith >= 80)] = 6

    # Create the bins for w based on dew point temperature
    w_bin = pd.Series(0, index=times, dtype=np.int64)
    w_bin[(w >= 0) & (w < 1)] = 1
    w_bin[(w >= 1) & (w < 2)] = 2
    w_bin[(w >= 2) & (w < 3)] = 3
    w_bin[(w >= 3)] = 4
    w_bin[(w == -1)] = 5

    # Create delta_kt_prime binning.
    delta_kt_prime_bin = pd.Series(0, index=times, dtype=np.int64)
    delta_kt_prime_bin[(delta_kt_prime >= 0) & (delta_kt_prime < 0.015)] = 1
    delta_kt_prime_bin[(delta_kt_prime >= 0.015) &
                       (delta_kt_prime < 0.035)] = 2
    delta_kt_prime_bin[(delta_kt_prime >= 0.035) & (delta_kt_prime < 0.07)] = 3
    delta_kt_prime_bin[(delta_kt_prime >= 0.07) & (delta_kt_prime < 0.15)] = 4
    delta_kt_prime_bin[(delta_kt_prime >= 0.15) & (delta_kt_prime < 0.3)] = 5
    delta_kt_prime_bin[(delta_kt_prime >= 0.3) & (delta_kt_prime <= 1)] = 6
    delta_kt_prime_bin[delta_kt_prime == -1] = 7

    return kt_prime_bin, zenith_bin, w_bin, delta_kt_prime_bin


def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
             use_delta_kt_prime=True, temp_dew=None, min_cos_zenith=0.065,
             max_zenith=87):
    """
    Determine DNI from GHI using the DIRINDEX model.

    The DIRINDEX model [1] modifies the DIRINT model implemented in
    :py:func:``pvlib.irradiance.dirint`` by taking into account information
    from a clear sky model. It is recommended that ``ghi_clearsky`` be
    calculated using the Ineichen clear sky model
    :py:func:``pvlib.clearsky.ineichen`` with ``perez_enhancement=True``.

    The pvlib implementation limits the clearness index to 1.

    Parameters
    ----------
    ghi : array-like
        Global horizontal irradiance in W/m^2.

    ghi_clearsky : array-like
        Global horizontal irradiance from clear sky model, in W/m^2.

    dni_clearsky : array-like
        Direct normal irradiance from clear sky model, in W/m^2.

    zenith : array-like
        True (not refraction-corrected) zenith angles in decimal
        degrees. If Z is a vector it must be of the same size as all
        other vector inputs. Z must be >=0 and <=180.

    times : DatetimeIndex

    pressure : float or array-like, default 101325.0
        The site pressure in Pascal. Pressure may be measured or an
        average pressure may be calculated from site altitude.

    use_delta_kt_prime : bool, default True
        If True, indicates that the stability index delta_kt_prime is
        included in the model. The stability index adjusts the estimated
        DNI in response to dynamics in the time series of GHI. It is
        recommended that delta_kt_prime is not used if the time between
        GHI points is 1.5 hours or greater. If use_delta_kt_prime=True,
        input data must be Series.

    temp_dew : None, float, or array-like, default None
        Surface dew point temperatures, in degrees C. Values of temp_dew
        may be numeric or NaN. Any single time period point with a
        temp_dew=NaN does not have dew point improvements applied. If
        temp_dew is not provided, then dew point improvements are not
        applied.

    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
    -------
    dni : array-like
        The modeled direct normal irradiance in W/m^2.

    Notes
    -----
    DIRINDEX model requires time series data (ie. one of the inputs must
    be a vector of length > 2).

    References
    ----------
    .. [1] Perez, R., Ineichen, P., Moore, K., Kmiecik, M., Chain, C., George,
       R., & Vignola, F. (2002). A new operational model for satellite-derived
       irradiances: description and validation. Solar Energy, 73(5), 307-317.
    """

    dni_dirint = dirint(ghi, zenith, times, pressure=pressure,
                        use_delta_kt_prime=use_delta_kt_prime,
                        temp_dew=temp_dew, min_cos_zenith=min_cos_zenith,
                        max_zenith=max_zenith)

    dni_dirint_clearsky = dirint(ghi_clearsky, zenith, times,
                                 pressure=pressure,
                                 use_delta_kt_prime=use_delta_kt_prime,
                                 temp_dew=temp_dew,
                                 min_cos_zenith=min_cos_zenith,
                                 max_zenith=max_zenith)

    dni_dirindex = dni_clearsky * dni_dirint / dni_dirint_clearsky

    dni_dirindex[dni_dirindex < 0] = 0.

    return dni_dirindex


def gti_dirint(poa_global, aoi, solar_zenith, solar_azimuth, times,
               surface_tilt, surface_azimuth, pressure=101325.,
               use_delta_kt_prime=True, temp_dew=None, albedo=.25,
               model='perez', model_perez='allsitescomposite1990',
               calculate_gt_90=True, max_iterations=30):
    """
    Determine GHI, DNI, DHI from POA global using the GTI DIRINT model.

    The GTI DIRINT model is described in [1]_.

    .. warning::

        Model performance is poor for AOI greater than approximately
        80 degrees `and` plane of array irradiance greater than
        approximately 200 W/m^2.

    Parameters
    ----------
    poa_global : array-like
        Plane of array global irradiance in W/m^2.

    aoi : array-like
        Angle of incidence of solar rays with respect to the module
        surface normal.

    solar_zenith : array-like
        True (not refraction-corrected) solar zenith angles in decimal
        degrees.

    solar_azimuth : array-like
        Solar azimuth angles in decimal degrees.

    times : DatetimeIndex
        Time indices for the input array-like data.

    surface_tilt : numeric
        Surface tilt angles in decimal degrees. 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).

    surface_azimuth : numeric
        Surface azimuth angles in decimal degrees. surface_azimuth 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).

    pressure : numeric, default 101325.0
        The site pressure in Pascal. Pressure may be measured or an
        average pressure may be calculated from site altitude.

    use_delta_kt_prime : bool, default True
        If True, indicates that the stability index delta_kt_prime is
        included in the model. The stability index adjusts the estimated
        DNI in response to dynamics in the time series of GHI. It is
        recommended that delta_kt_prime is not used if the time between
        GHI points is 1.5 hours or greater. If use_delta_kt_prime=True,
        input data must be Series.

    temp_dew : None, float, or array-like, default None
        Surface dew point temperatures, in degrees C. Values of temp_dew
        may be numeric or NaN. Any single time period point with a
        temp_dew=NaN does not have dew point improvements applied. If
        temp_dew is not provided, then dew point improvements are not
        applied.

    albedo : numeric, default 0.25
        Surface albedo

    model : String, default 'isotropic'
        Irradiance model.

    model_perez : String, default 'allsitescomposite1990'
        Used only if model='perez'. See :py:func:`perez`.

    calculate_gt_90 : bool, default True
        Controls if the algorithm evaluates inputs with AOI >= 90 degrees.
        If False, returns nan for AOI >= 90 degrees. Significant speed ups
        can be achieved by setting this parameter to False.

    max_iterations : int, default 30
        Maximum number of iterations for the aoi < 90 deg algorithm.

    Returns
    -------
    data : OrderedDict or DataFrame
        Contains the following keys/columns:

            * ``ghi``: the modeled global horizontal irradiance in W/m^2.
            * ``dni``: the modeled direct normal irradiance in W/m^2.
            * ``dhi``: the modeled diffuse horizontal irradiance in
              W/m^2.

    References
    ----------
    .. [1] B. Marion, A model for deriving the direct normal and
           diffuse horizontal irradiance from the global tilted
           irradiance, Solar Energy 122, 1037-1046.
           :doi:`10.1016/j.solener.2015.10.024`
    """

    aoi_lt_90 = aoi < 90

    # for AOI less than 90 degrees
    ghi, dni, dhi, kt_prime = _gti_dirint_lt_90(
        poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth, times,
        surface_tilt, surface_azimuth, pressure=pressure,
        use_delta_kt_prime=use_delta_kt_prime, temp_dew=temp_dew,
        albedo=albedo, model=model, model_perez=model_perez,
        max_iterations=max_iterations)

    # for AOI greater than or equal to 90 degrees
    if calculate_gt_90:
        ghi_gte_90, dni_gte_90, dhi_gte_90 = _gti_dirint_gte_90(
            poa_global, aoi, solar_zenith, solar_azimuth,
            surface_tilt, times, kt_prime,
            pressure=pressure, temp_dew=temp_dew, albedo=albedo)
    else:
        ghi_gte_90, dni_gte_90, dhi_gte_90 = np.nan, np.nan, np.nan

    # put the AOI < 90 and AOI >= 90 conditions together
    output = OrderedDict()
    output['ghi'] = ghi.where(aoi_lt_90, ghi_gte_90)
    output['dni'] = dni.where(aoi_lt_90, dni_gte_90)
    output['dhi'] = dhi.where(aoi_lt_90, dhi_gte_90)

    output = pd.DataFrame(output, index=times)

    return output


def _gti_dirint_lt_90(poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth,
                      times, surface_tilt, surface_azimuth, pressure=101325.,
                      use_delta_kt_prime=True, temp_dew=None, albedo=.25,
                      model='perez', model_perez='allsitescomposite1990',
                      max_iterations=30):
    """
    GTI-DIRINT model for AOI < 90 degrees. See Marion 2015 Section 2.1.

    See gti_dirint signature for parameter details.
    """
    I0 = get_extra_radiation(times, 1370, 'spencer')
    cos_zenith = tools.cosd(solar_zenith)
    # I0h as in Marion 2015 eqns 1, 3
    I0h = I0 * np.maximum(0.065, cos_zenith)

    airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
    airmass = atmosphere.get_absolute_airmass(airmass, pressure)

    # these coeffs and diff variables and the loop below
    # implement figure 1 of Marion 2015

    # make coeffs that is at least 30 elements long so that all
    # coeffs can be assigned as specified in Marion 2015.
    # slice below will limit iterations if necessary
    coeffs = np.empty(max(30, max_iterations))
    coeffs[0:3] = 1
    coeffs[3:10] = 0.5
    coeffs[10:20] = 0.25
    coeffs[20:] = 0.125
    coeffs = coeffs[:max_iterations]  # covers case where max_iterations < 30

    # initialize diff
    diff = pd.Series(9999, index=times)
    best_diff = diff

    # initialize poa_global_i
    poa_global_i = poa_global

    for iteration, coeff in enumerate(coeffs):

        # test if difference between modeled GTI and
        # measured GTI (poa_global) is less than 1 W/m^2
        # only test for aoi less than 90 deg
        best_diff_lte_1 = best_diff <= 1
        best_diff_lte_1_lt_90 = best_diff_lte_1[aoi_lt_90]
        if best_diff_lte_1_lt_90.all():
            # all aoi < 90 points have a difference <= 1, so break loop
            break

        # calculate kt and DNI from GTI
        kt = clearness_index(poa_global_i, aoi, I0)  # kt from Marion eqn 2
        disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0)
        kt_prime = clearness_index_zenith_independent(kt, airmass)
        # dirint DNI in Marion eqn 3
        dni = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
                                       use_delta_kt_prime, temp_dew)

        # calculate DHI using Marion eqn 3 (identify 1st term on RHS as GHI)
        # I0h has a minimum zenith projection, but multiplier of DNI does not
        ghi = kt * I0h                  # Kt * I0 * max(0.065, cos(zen))
        dhi = ghi - dni * cos_zenith    # no cos(zen) restriction here

        # following SSC code
        dni = np.maximum(dni, 0)
        ghi = np.maximum(ghi, 0)
        dhi = np.maximum(dhi, 0)

        # use DNI and DHI to model GTI
        # GTI-DIRINT uses perez transposition model, but we allow for
        # any model here
        all_irrad = get_total_irradiance(
            surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
            dni, ghi, dhi, dni_extra=I0, airmass=airmass,
            albedo=albedo, model=model, model_perez=model_perez)

        gti_model = all_irrad['poa_global']

        # calculate new diff
        diff = gti_model - poa_global

        # determine if the new diff is smaller in magnitude
        # than the old diff
        diff_abs = diff.abs()
        smallest_diff = diff_abs < best_diff

        # save the best differences
        best_diff = diff_abs.where(smallest_diff, best_diff)

        # on first iteration, the best values are the only values
        if iteration == 0:
            best_ghi = ghi
            best_dni = dni
            best_dhi = dhi
            best_kt_prime = kt_prime
        else:
            # save new DNI, DHI, DHI if they provide the best consistency
            # otherwise use the older values.
            best_ghi = ghi.where(smallest_diff, best_ghi)
            best_dni = dni.where(smallest_diff, best_dni)
            best_dhi = dhi.where(smallest_diff, best_dhi)
            best_kt_prime = kt_prime.where(smallest_diff, best_kt_prime)

        # calculate adjusted inputs for next iteration. Marion eqn 4
        poa_global_i = np.maximum(1.0, poa_global_i - coeff * diff)
    else:
        # we are here because we ran out of coeffs to loop over and
        # therefore we have exceeded max_iterations
        import warnings
        failed_points = best_diff[aoi_lt_90][~best_diff_lte_1_lt_90]
        warnings.warn(
            ('%s points failed to converge after %s iterations. best_diff:\n%s'
             % (len(failed_points), max_iterations, failed_points)),
            RuntimeWarning)

    # return the best data, whether or not the solution converged
    return best_ghi, best_dni, best_dhi, best_kt_prime


def _gti_dirint_gte_90(poa_global, aoi, solar_zenith, solar_azimuth,
                       surface_tilt, times, kt_prime,
                       pressure=101325., temp_dew=None, albedo=.25):
    """
    GTI-DIRINT model for AOI >= 90 degrees. See Marion 2015 Section 2.2.

    See gti_dirint signature for parameter details.
    """
    kt_prime_gte_90 = _gti_dirint_gte_90_kt_prime(aoi, solar_zenith,
                                                  solar_azimuth, times,
                                                  kt_prime)

    I0 = get_extra_radiation(times, 1370, 'spencer')
    airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
    airmass = atmosphere.get_absolute_airmass(airmass, pressure)
    kt = kt_prime_gte_90 * _kt_kt_prime_factor(airmass)
    disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0)

    dni_gte_90 = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
                                          False, temp_dew)

    dni_gte_90_proj = dni_gte_90 * tools.cosd(solar_zenith)
    cos_surface_tilt = tools.cosd(surface_tilt)

    # isotropic sky plus ground diffuse
    dhi_gte_90 = (
        (2 * poa_global - dni_gte_90_proj * albedo * (1 - cos_surface_tilt)) /
        (1 + cos_surface_tilt + albedo * (1 - cos_surface_tilt)))

    ghi_gte_90 = dni_gte_90_proj + dhi_gte_90

    return ghi_gte_90, dni_gte_90, dhi_gte_90


def _gti_dirint_gte_90_kt_prime(aoi, solar_zenith, solar_azimuth, times,
                                kt_prime):
    """
    Determine kt' values to be used in GTI-DIRINT AOI >= 90 deg case.
    See Marion 2015 Section 2.2.

    For AOI >= 90 deg: average of the kt_prime values for 65 < AOI < 80
    in each day's morning and afternoon. Morning and afternoon are treated
    separately.

    For AOI < 90 deg: NaN.

    See gti_dirint signature for parameter details.

    Returns
    -------
    kt_prime_gte_90 : Series
        Index is `times`.
    """
    # kt_prime values from DIRINT calculation for AOI < 90 case
    # set the kt_prime from sunrise to AOI=90 to be equal to
    # the kt_prime for 65 < AOI < 80 during the morning.
    # similar for the afternoon. repeat for every day.
    aoi_gte_90 = aoi >= 90
    aoi_65_80 = (aoi > 65) & (aoi < 80)
    zenith_lt_90 = solar_zenith < 90
    morning = solar_azimuth < 180
    afternoon = solar_azimuth > 180
    aoi_65_80_morning = aoi_65_80 & morning
    aoi_65_80_afternoon = aoi_65_80 & afternoon
    zenith_lt_90_aoi_gte_90_morning = zenith_lt_90 & aoi_gte_90 & morning
    zenith_lt_90_aoi_gte_90_afternoon = zenith_lt_90 & aoi_gte_90 & afternoon

    kt_prime_gte_90 = []
    for date, data in kt_prime.groupby(times.date):
        kt_prime_am_avg = data[aoi_65_80_morning].mean()
        kt_prime_pm_avg = data[aoi_65_80_afternoon].mean()

        kt_prime_by_date = pd.Series(np.nan, index=data.index)
        kt_prime_by_date[zenith_lt_90_aoi_gte_90_morning] = kt_prime_am_avg
        kt_prime_by_date[zenith_lt_90_aoi_gte_90_afternoon] = kt_prime_pm_avg
        kt_prime_gte_90.append(kt_prime_by_date)
    kt_prime_gte_90 = pd.concat(kt_prime_gte_90)

    return kt_prime_gte_90


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


def liujordan(zenith, transmittance, airmass, dni_extra=1367.0):
    '''
    Determine DNI, DHI, GHI from extraterrestrial flux, transmittance,
    and optical air mass number.

    Liu and Jordan, 1960, developed a simplified direct radiation model.
    DHI is from an empirical equation for diffuse radiation from Liu and
    Jordan, 1960.

    Parameters
    ----------
    zenith: pd.Series
        True (not refraction-corrected) zenith angles in decimal
        degrees. If Z is a vector it must be of the same size as all
        other vector inputs. Z must be >=0 and <=180.

    transmittance: float
        Atmospheric transmittance between 0 and 1.

    pressure: float, default 101325.0
        Air pressure

    dni_extra: float, default 1367.0
        Direct irradiance incident at the top of the atmosphere.

    Returns
    -------
    irradiance: DataFrame
        Modeled direct normal irradiance, direct horizontal irradiance,
        and global horizontal irradiance in W/m^2

    References
    ----------
    .. [1] Campbell, G. S., J. M. Norman (1998) An Introduction to
       Environmental Biophysics. 2nd Ed. New York: Springer.

    .. [2] Liu, B. Y., R. C. Jordan, (1960). "The interrelationship and
       characteristic distribution of direct, diffuse, and total solar
       radiation".  Solar Energy 4:1-19
    '''

    tau = transmittance

    dni = dni_extra*tau**airmass
    dhi = 0.3 * (1.0 - tau**airmass) * dni_extra * np.cos(np.radians(zenith))
    ghi = dhi + dni * np.cos(np.radians(zenith))

    irrads = OrderedDict()
    irrads['ghi'] = ghi
    irrads['dni'] = dni
    irrads['dhi'] = dhi

    if isinstance(ghi, pd.Series):
        irrads = pd.DataFrame(irrads)

    return irrads


def _get_perez_coefficients(perezmodel):
    '''
    Find coefficients for the Perez model

    Parameters
    ----------

    perezmodel : string (optional, default='allsitescomposite1990')

          a character string which selects the desired set of Perez
          coefficients. If model is not provided as an input, the default,
          '1990' will be used.

    All possible model selections are:

          * '1990'
          * 'allsitescomposite1990' (same as '1990')
          * 'allsitescomposite1988'
          * 'sandiacomposite1988'
          * 'usacomposite1988'
          * 'france1988'
          * 'phoenix1988'
          * 'elmonte1988'
          * 'osage1988'
          * 'albuquerque1988'
          * 'capecanaveral1988'
          * 'albany1988'

    Returns
    --------
    F1coeffs, F2coeffs : (array, array)
          F1 and F2 coefficients for the Perez model

    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

    '''
    coeffdict = {
        'allsitescomposite1990': [
            [-0.0080,    0.5880,   -0.0620,   -0.0600,    0.0720,   -0.0220],
            [0.1300,    0.6830,   -0.1510,   -0.0190,    0.0660,   -0.0290],
            [0.3300,    0.4870,   -0.2210,    0.0550,   -0.0640,   -0.0260],
            [0.5680,    0.1870,   -0.2950,    0.1090,   -0.1520,   -0.0140],
            [0.8730,   -0.3920,   -0.3620,    0.2260,   -0.4620,    0.0010],
            [1.1320,   -1.2370,   -0.4120,    0.2880,   -0.8230,    0.0560],
            [1.0600,   -1.6000,   -0.3590,    0.2640,   -1.1270,    0.1310],
            [0.6780,   -0.3270,   -0.2500,    0.1560,   -1.3770,    0.2510]],
        'allsitescomposite1988': [
            [-0.0180,    0.7050,   -0.071,   -0.0580,    0.1020,   -0.0260],
            [0.1910,    0.6450,   -0.1710,    0.0120,    0.0090,   -0.0270],
            [0.4400,    0.3780,   -0.2560,    0.0870,   -0.1040,   -0.0250],
            [0.7560,   -0.1210,   -0.3460,    0.1790,   -0.3210,   -0.0080],
            [0.9960,   -0.6450,   -0.4050,    0.2600,   -0.5900,    0.0170],
            [1.0980,   -1.2900,   -0.3930,    0.2690,   -0.8320,    0.0750],
            [0.9730,   -1.1350,   -0.3780,    0.1240,   -0.2580,    0.1490],
            [0.6890,   -0.4120,   -0.2730,    0.1990,   -1.6750,    0.2370]],
        'sandiacomposite1988': [
            [-0.1960,    1.0840,   -0.0060,   -0.1140,    0.1800,   -0.0190],
            [0.2360,    0.5190,   -0.1800,   -0.0110,    0.0200,   -0.0380],
            [0.4540,    0.3210,   -0.2550,    0.0720,   -0.0980,   -0.0460],
            [0.8660,   -0.3810,   -0.3750,    0.2030,   -0.4030,   -0.0490],
            [1.0260,   -0.7110,   -0.4260,    0.2730,   -0.6020,   -0.0610],
            [0.9780,   -0.9860,   -0.3500,    0.2800,   -0.9150,   -0.0240],
            [0.7480,   -0.9130,   -0.2360,    0.1730,   -1.0450,    0.0650],
            [0.3180,   -0.7570,    0.1030,    0.0620,   -1.6980,    0.2360]],
        'usacomposite1988': [
            [-0.0340,    0.6710,   -0.0590,   -0.0590,    0.0860,   -0.0280],
            [0.2550,    0.4740,   -0.1910,    0.0180,   -0.0140,   -0.0330],
            [0.4270,    0.3490,   -0.2450,    0.0930,   -0.1210,   -0.0390],
            [0.7560,   -0.2130,   -0.3280,    0.1750,   -0.3040,   -0.0270],
            [1.0200,   -0.8570,   -0.3850,    0.2800,   -0.6380,   -0.0190],
            [1.0500,   -1.3440,   -0.3480,    0.2800,   -0.8930,    0.0370],
            [0.9740,   -1.5070,   -0.3700,    0.1540,   -0.5680,    0.1090],
            [0.7440,   -1.8170,   -0.2560,    0.2460,   -2.6180,    0.2300]],
        'france1988': [
            [0.0130,    0.7640,   -0.1000,   -0.0580,    0.1270,   -0.0230],
            [0.0950,    0.9200,   -0.1520,         0,    0.0510,   -0.0200],
            [0.4640,    0.4210,   -0.2800,    0.0640,   -0.0510,   -0.0020],
            [0.7590,   -0.0090,   -0.3730,    0.2010,   -0.3820,    0.0100],
            [0.9760,   -0.4000,   -0.4360,    0.2710,   -0.6380,    0.0510],
            [1.1760,   -1.2540,   -0.4620,    0.2950,   -0.9750,    0.1290],
            [1.1060,   -1.5630,   -0.3980,    0.3010,   -1.4420,    0.2120],
            [0.9340,   -1.5010,   -0.2710,    0.4200,   -2.9170,    0.2490]],
        'phoenix1988': [
            [-0.0030,    0.7280,   -0.0970,   -0.0750,    0.1420,   -0.0430],
            [0.2790,    0.3540,   -0.1760,    0.0300,   -0.0550,   -0.0540],
            [0.4690,    0.1680,   -0.2460,    0.0480,   -0.0420,   -0.0570],
            [0.8560,   -0.5190,   -0.3400,    0.1760,   -0.3800,   -0.0310],
            [0.9410,   -0.6250,   -0.3910,    0.1880,   -0.3600,   -0.0490],
            [1.0560,   -1.1340,   -0.4100,    0.2810,   -0.7940,   -0.0650],
            [0.9010,   -2.1390,   -0.2690,    0.1180,   -0.6650,    0.0460],
            [0.1070,    0.4810,    0.1430,   -0.1110,   -0.1370,    0.2340]],
        'elmonte1988': [
            [0.0270,    0.7010,   -0.1190,   -0.0580,    0.1070,  -0.0600],
            [0.1810,    0.6710,   -0.1780,   -0.0790,    0.1940,  -0.0350],
            [0.4760,    0.4070,   -0.2880,    0.0540,   -0.0320,  -0.0550],
            [0.8750,   -0.2180,   -0.4030,    0.1870,   -0.3090,  -0.0610],
            [1.1660,   -1.0140,   -0.4540,    0.2110,   -0.4100,  -0.0440],
            [1.1430,   -2.0640,   -0.2910,    0.0970,   -0.3190,   0.0530],
            [1.0940,   -2.6320,   -0.2590,    0.0290,   -0.4220,   0.1470],
            [0.1550,    1.7230,    0.1630,   -0.1310,   -0.0190,   0.2770]],
        'osage1988': [
            [-0.3530,    1.4740,   0.0570,   -0.1750,    0.3120,   0.0090],
            [0.3630,    0.2180,  -0.2120,    0.0190,   -0.0340,  -0.0590],
            [-0.0310,    1.2620,  -0.0840,   -0.0820,    0.2310,  -0.0170],
            [0.6910,    0.0390,  -0.2950,    0.0910,   -0.1310,  -0.0350],
            [1.1820,   -1.3500,  -0.3210,    0.4080,   -0.9850,  -0.0880],
            [0.7640,    0.0190,  -0.2030,    0.2170,   -0.2940,  -0.1030],
            [0.2190,    1.4120,   0.2440,    0.4710,   -2.9880,   0.0340],
            [3.5780,   22.2310, -10.7450,    2.4260,    4.8920,  -5.6870]],
        'albuquerque1988': [
            [0.0340,    0.5010,  -0.0940,   -0.0630,    0.1060,  -0.0440],
            [0.2290,    0.4670,  -0.1560,   -0.0050,   -0.0190,  -0.0230],
            [0.4860,    0.2410,  -0.2530,    0.0530,   -0.0640,  -0.0220],
            [0.8740,   -0.3930,  -0.3970,    0.1810,   -0.3270,  -0.0370],
            [1.1930,   -1.2960,  -0.5010,    0.2810,   -0.6560,  -0.0450],
            [1.0560,   -1.7580,  -0.3740,    0.2260,   -0.7590,   0.0340],
            [0.9010,   -4.7830,  -0.1090,    0.0630,   -0.9700,   0.1960],
            [0.8510,   -7.0550,  -0.0530,    0.0600,   -2.8330,   0.3300]],
        'capecanaveral1988': [
            [0.0750,    0.5330,   -0.1240,  -0.0670,   0.0420,  -0.0200],
            [0.2950,    0.4970,   -0.2180,  -0.0080,   0.0030,  -0.0290],
            [0.5140,    0.0810,   -0.2610,   0.0750,  -0.1600,  -0.0290],
            [0.7470,   -0.3290,   -0.3250,   0.1810,  -0.4160,  -0.0300],
            [0.9010,   -0.8830,   -0.2970,   0.1780,  -0.4890,   0.0080],
            [0.5910,   -0.0440,   -0.1160,   0.2350,  -0.9990,   0.0980],
            [0.5370,   -2.4020,    0.3200,   0.1690,  -1.9710,   0.3100],
            [-0.8050,    4.5460,    1.0720,  -0.2580,  -0.9500,    0.7530]],
        'albany1988': [
            [0.0120,    0.5540,   -0.0760, -0.0520,   0.0840,  -0.0290],
            [0.2670,    0.4370,   -0.1940,  0.0160,   0.0220,  -0.0360],
            [0.4200,    0.3360,   -0.2370,  0.0740,  -0.0520,  -0.0320],
            [0.6380,   -0.0010,   -0.2810,  0.1380,  -0.1890,  -0.0120],
            [1.0190,   -1.0270,   -0.3420,  0.2710,  -0.6280,   0.0140],
            [1.1490,   -1.9400,   -0.3310,  0.3220,  -1.0970,   0.0800],
            [1.4340,   -3.9940,   -0.4920,  0.4530,  -2.3760,   0.1170],
            [1.0070,   -2.2920,   -0.4820,  0.3900,  -3.3680,   0.2290]], }

    array = np.array(coeffdict[perezmodel])

    F1coeffs = array[:, 0:3]
    F2coeffs = array[:, 3:7]

    return F1coeffs, F2coeffs


def _get_dirint_coeffs():
    """
    A place to stash the dirint coefficients.

    Returns
    -------
    np.array with shape ``(6, 6, 7, 5)``.
    Ordering is ``[kt_prime_bin, zenith_bin, delta_kt_prime_bin, w_bin]``
    """

    # To allow for maximum copy/paste from the MATLAB 1-indexed code,
    # we create and assign values to an oversized array.
    # Then, we return the [1:, 1:, :, :] slice.

    coeffs = np.zeros((7, 7, 7, 5))

    coeffs[1, 1, :, :] = [
        [0.385230, 0.385230, 0.385230, 0.462880, 0.317440],
        [0.338390, 0.338390, 0.221270, 0.316730, 0.503650],
        [0.235680, 0.235680, 0.241280, 0.157830, 0.269440],
        [0.830130, 0.830130, 0.171970, 0.841070, 0.457370],
        [0.548010, 0.548010, 0.478000, 0.966880, 1.036370],
        [0.548010, 0.548010, 1.000000, 3.012370, 1.976540],
        [0.582690, 0.582690, 0.229720, 0.892710, 0.569950]]

    coeffs[1, 2, :, :] = [
        [0.131280, 0.131280, 0.385460, 0.511070, 0.127940],
        [0.223710, 0.223710, 0.193560, 0.304560, 0.193940],
        [0.229970, 0.229970, 0.275020, 0.312730, 0.244610],
        [0.090100, 0.184580, 0.260500, 0.687480, 0.579440],
        [0.131530, 0.131530, 0.370190, 1.380350, 1.052270],
        [1.116250, 1.116250, 0.928030, 3.525490, 2.316920],
        [0.090100, 0.237000, 0.300040, 0.812470, 0.664970]]

    coeffs[1, 3, :, :] = [
        [0.587510, 0.130000, 0.400000, 0.537210, 0.832490],
        [0.306210, 0.129830, 0.204460, 0.500000, 0.681640],
        [0.224020, 0.260620, 0.334080, 0.501040, 0.350470],
        [0.421540, 0.753970, 0.750660, 3.706840, 0.983790],
        [0.706680, 0.373530, 1.245670, 0.864860, 1.992630],
        [4.864400, 0.117390, 0.265180, 0.359180, 3.310820],
        [0.392080, 0.493290, 0.651560, 1.932780, 0.898730]]

    coeffs[1, 4, :, :] = [
        [0.126970, 0.126970, 0.126970, 0.126970, 0.126970],
        [0.810820, 0.810820, 0.810820, 0.810820, 0.810820],
        [3.241680, 2.500000, 2.291440, 2.291440, 2.291440],
        [4.000000, 3.000000, 2.000000, 0.975430, 1.965570],
        [12.494170, 12.494170, 8.000000, 5.083520, 8.792390],
        [21.744240, 21.744240, 21.744240, 21.744240, 21.744240],
        [3.241680, 12.494170, 1.620760, 1.375250, 2.331620]]

    coeffs[1, 5, :, :] = [
        [0.126970, 0.126970, 0.126970, 0.126970, 0.126970],
        [0.810820, 0.810820, 0.810820, 0.810820, 0.810820],
        [3.241680, 2.500000, 2.291440, 2.291440, 2.291440],
        [4.000000, 3.000000, 2.000000, 0.975430, 1.965570],
        [12.494170, 12.494170, 8.000000, 5.083520, 8.792390],
        [21.744240, 21.744240, 21.744240, 21.744240, 21.744240],
        [3.241680, 12.494170, 1.620760, 1.375250, 2.331620]]

    coeffs[1, 6, :, :] = [
        [0.126970, 0.126970, 0.126970, 0.126970, 0.126970],
        [0.810820, 0.810820, 0.810820, 0.810820, 0.810820],
        [3.241680, 2.500000, 2.291440, 2.291440, 2.291440],
        [4.000000, 3.000000, 2.000000, 0.975430, 1.965570],
        [12.494170, 12.494170, 8.000000, 5.083520, 8.792390],
        [21.744240, 21.744240, 21.744240, 21.744240, 21.744240],
        [3.241680, 12.494170, 1.620760, 1.375250, 2.331620]]

    coeffs[2, 1, :, :] = [
        [0.337440, 0.337440, 0.969110, 1.097190, 1.116080],
        [0.337440, 0.337440, 0.969110, 1.116030, 0.623900],
        [0.337440, 0.337440, 1.530590, 1.024420, 0.908480],
        [0.584040, 0.584040, 0.847250, 0.914940, 1.289300],
        [0.337440, 0.337440, 0.310240, 1.435020, 1.852830],
        [0.337440, 0.337440, 1.015010, 1.097190, 2.117230],
        [0.337440, 0.337440, 0.969110, 1.145730, 1.476400]]

    coeffs[2, 2, :, :] = [
        [0.300000, 0.300000, 0.700000, 1.100000, 0.796940],
        [0.219870, 0.219870, 0.526530, 0.809610, 0.649300],
        [0.386650, 0.386650, 0.119320, 0.576120, 0.685460],
        [0.746730, 0.399830, 0.470970, 0.986530, 0.785370],
        [0.575420, 0.936700, 1.649200, 1.495840, 1.335590],
        [1.319670, 4.002570, 1.276390, 2.644550, 2.518670],
        [0.665190, 0.678910, 1.012360, 1.199940, 0.986580]]

    coeffs[2, 3, :, :] = [
        [0.378870, 0.974060, 0.500000, 0.491880, 0.665290],
        [0.105210, 0.263470, 0.407040, 0.553460, 0.582590],
        [0.312900, 0.345240, 1.144180, 0.854790, 0.612280],
        [0.119070, 0.365120, 0.560520, 0.793720, 0.802600],
        [0.781610, 0.837390, 1.270420, 1.537980, 1.292950],
        [1.152290, 1.152290, 1.492080, 1.245370, 2.177100],
        [0.424660, 0.529550, 0.966910, 1.033460, 0.958730]]

    coeffs[2, 4, :, :] = [
        [0.310590, 0.714410, 0.252450, 0.500000, 0.607600],
        [0.975190, 0.363420, 0.500000, 0.400000, 0.502800],
        [0.175580, 0.196250, 0.476360, 1.072470, 0.490510],
        [0.719280, 0.698620, 0.657770, 1.190840, 0.681110],
        [0.426240, 1.464840, 0.678550, 1.157730, 0.978430],
        [2.501120, 1.789130, 1.387090, 2.394180, 2.394180],
        [0.491640, 0.677610, 0.685610, 1.082400, 0.735410]]

    coeffs[2, 5, :, :] = [
        [0.597000, 0.500000, 0.300000, 0.310050, 0.413510],
        [0.314790, 0.336310, 0.400000, 0.400000, 0.442460],
        [0.166510, 0.460440, 0.552570, 1.000000, 0.461610],
        [0.401020, 0.559110, 0.403630, 1.016710, 0.671490],
        [0.400360, 0.750830, 0.842640, 1.802600, 1.023830],
        [3.315300, 1.510380, 2.443650, 1.638820, 2.133990],
        [0.530790, 0.745850, 0.693050, 1.458040, 0.804500]]

    coeffs[2, 6, :, :] = [
        [0.597000, 0.500000, 0.300000, 0.310050, 0.800920],
        [0.314790, 0.336310, 0.400000, 0.400000, 0.237040],
        [0.166510, 0.460440, 0.552570, 1.000000, 0.581990],
        [0.401020, 0.559110, 0.403630, 1.016710, 0.898570],
        [0.400360, 0.750830, 0.842640, 1.802600, 3.400390],
        [3.315300, 1.510380, 2.443650, 1.638820, 2.508780],
        [0.204340, 1.157740, 2.003080, 2.622080, 1.409380]]

    coeffs[3, 1, :, :] = [
        [1.242210, 1.242210, 1.242210, 1.242210, 1.242210],
        [0.056980, 0.056980, 0.656990, 0.656990, 0.925160],
        [0.089090, 0.089090, 1.040430, 1.232480, 1.205300],
        [1.053850, 1.053850, 1.399690, 1.084640, 1.233340],
        [1.151540, 1.151540, 1.118290, 1.531640, 1.411840],
        [1.494980, 1.494980, 1.700000, 1.800810, 1.671600],
        [1.018450, 1.018450, 1.153600, 1.321890, 1.294670]]

    coeffs[3, 2, :, :] = [
        [0.700000, 0.700000, 1.023460, 0.700000, 0.945830],
        [0.886300, 0.886300, 1.333620, 0.800000, 1.066620],
        [0.902180, 0.902180, 0.954330, 1.126690, 1.097310],
        [1.095300, 1.075060, 1.176490, 1.139470, 1.096110],
        [1.201660, 1.201660, 1.438200, 1.256280, 1.198060],
        [1.525850, 1.525850, 1.869160, 1.985410, 1.911590],
        [1.288220, 1.082810, 1.286370, 1.166170, 1.119330]]

    coeffs[3, 3, :, :] = [
        [0.600000, 1.029910, 0.859890, 0.550000, 0.813600],
        [0.604450, 1.029910, 0.859890, 0.656700, 0.928840],
        [0.455850, 0.750580, 0.804930, 0.823000, 0.911000],
        [0.526580, 0.932310, 0.908620, 0.983520, 0.988090],
        [1.036110, 1.100690, 0.848380, 1.035270, 1.042380],
        [1.048440, 1.652720, 0.900000, 2.350410, 1.082950],
        [0.817410, 0.976160, 0.861300, 0.974780, 1.004580]]

    coeffs[3, 4, :, :] = [
        [0.782110, 0.564280, 0.600000, 0.600000, 0.665740],
        [0.894480, 0.680730, 0.541990, 0.800000, 0.669140],
        [0.487460, 0.818950, 0.841830, 0.872540, 0.709040],
        [0.709310, 0.872780, 0.908480, 0.953290, 0.844350],
        [0.863920, 0.947770, 0.876220, 1.078750, 0.936910],
        [1.280350, 0.866720, 0.769790, 1.078750, 0.975130],
        [0.725420, 0.869970, 0.868810, 0.951190, 0.829220]]

    coeffs[3, 5, :, :] = [
        [0.791750, 0.654040, 0.483170, 0.409000, 0.597180],
        [0.566140, 0.948990, 0.971820, 0.653570, 0.718550],
        [0.648710, 0.637730, 0.870510, 0.860600, 0.694300],
        [0.637630, 0.767610, 0.925670, 0.990310, 0.847670],
        [0.736380, 0.946060, 1.117590, 1.029340, 0.947020],
        [1.180970, 0.850000, 1.050000, 0.950000, 0.888580],
        [0.700560, 0.801440, 0.961970, 0.906140, 0.823880]]

    coeffs[3, 6, :, :] = [
        [0.500000, 0.500000, 0.586770, 0.470550, 0.629790],
        [0.500000, 0.500000, 1.056220, 1.260140, 0.658140],
        [0.500000, 0.500000, 0.631830, 0.842620, 0.582780],
        [0.554710, 0.734730, 0.985820, 0.915640, 0.898260],
        [0.712510, 1.205990, 0.909510, 1.078260, 0.885610],
        [1.899260, 1.559710, 1.000000, 1.150000, 1.120390],
        [0.653880, 0.793120, 0.903320, 0.944070, 0.796130]]

    coeffs[4, 1, :, :] = [
        [1.000000, 1.000000, 1.050000, 1.170380, 1.178090],
        [0.960580, 0.960580, 1.059530, 1.179030, 1.131690],
        [0.871470, 0.871470, 0.995860, 1.141910, 1.114600],
        [1.201590, 1.201590, 0.993610, 1.109380, 1.126320],
        [1.065010, 1.065010, 0.828660, 0.939970, 1.017930],
        [1.065010, 1.065010, 0.623690, 1.119620, 1.132260],
        [1.071570, 1.071570, 0.958070, 1.114130, 1.127110]]

    coeffs[4, 2, :, :] = [
        [0.950000, 0.973390, 0.852520, 1.092200, 1.096590],
        [0.804120, 0.913870, 0.980990, 1.094580, 1.042420],
        [0.737540, 0.935970, 0.999940, 1.056490, 1.050060],
        [1.032980, 1.034540, 0.968460, 1.032080, 1.015780],
        [0.900000, 0.977210, 0.945960, 1.008840, 0.969960],
        [0.600000, 0.750000, 0.750000, 0.844710, 0.899100],
        [0.926800, 0.965030, 0.968520, 1.044910, 1.032310]]

    coeffs[4, 3, :, :] = [
        [0.850000, 1.029710, 0.961100, 1.055670, 1.009700],
        [0.818530, 0.960010, 0.996450, 1.081970, 1.036470],
        [0.765380, 0.953500, 0.948260, 1.052110, 1.000140],
        [0.775610, 0.909610, 0.927800, 0.987800, 0.952100],
        [1.000990, 0.881880, 0.875950, 0.949100, 0.893690],
        [0.902370, 0.875960, 0.807990, 0.942410, 0.917920],
        [0.856580, 0.928270, 0.946820, 1.032260, 0.972990]]

    coeffs[4, 4, :, :] = [
        [0.750000, 0.857930, 0.983800, 1.056540, 0.980240],
        [0.750000, 0.987010, 1.013730, 1.133780, 1.038250],
        [0.800000, 0.947380, 1.012380, 1.091270, 0.999840],
        [0.800000, 0.914550, 0.908570, 0.999190, 0.915230],
        [0.778540, 0.800590, 0.799070, 0.902180, 0.851560],
        [0.680190, 0.317410, 0.507680, 0.388910, 0.646710],
        [0.794920, 0.912780, 0.960830, 1.057110, 0.947950]]

    coeffs[4, 5, :, :] = [
        [0.750000, 0.833890, 0.867530, 1.059890, 0.932840],
        [0.979700, 0.971470, 0.995510, 1.068490, 1.030150],
        [0.858850, 0.987920, 1.043220, 1.108700, 1.044900],
        [0.802400, 0.955110, 0.911660, 1.045070, 0.944470],
        [0.884890, 0.766210, 0.885390, 0.859070, 0.818190],
        [0.615680, 0.700000, 0.850000, 0.624620, 0.669300],
        [0.835570, 0.946150, 0.977090, 1.049350, 0.979970]]

    coeffs[4, 6, :, :] = [
        [0.689220, 0.809600, 0.900000, 0.789500, 0.853990],
        [0.854660, 0.852840, 0.938200, 0.923110, 0.955010],
        [0.938600, 0.932980, 1.010390, 1.043950, 1.041640],
        [0.843620, 0.981300, 0.951590, 0.946100, 0.966330],
        [0.694740, 0.814690, 0.572650, 0.400000, 0.726830],
        [0.211370, 0.671780, 0.416340, 0.297290, 0.498050],
        [0.843540, 0.882330, 0.911760, 0.898420, 0.960210]]

    coeffs[5, 1, :, :] = [
        [1.054880, 1.075210, 1.068460, 1.153370, 1.069220],
        [1.000000, 1.062220, 1.013470, 1.088170, 1.046200],
        [0.885090, 0.993530, 0.942590, 1.054990, 1.012740],
        [0.920000, 0.950000, 0.978720, 1.020280, 0.984440],
        [0.850000, 0.908500, 0.839940, 0.985570, 0.962180],
        [0.800000, 0.800000, 0.810080, 0.950000, 0.961550],
        [1.038590, 1.063200, 1.034440, 1.112780, 1.037800]]

    coeffs[5, 2, :, :] = [
        [1.017610, 1.028360, 1.058960, 1.133180, 1.045620],
        [0.920000, 0.998970, 1.033590, 1.089030, 1.022060],
        [0.912370, 0.949930, 0.979770, 1.020420, 0.981770],
        [0.847160, 0.935300, 0.930540, 0.955050, 0.946560],
        [0.880260, 0.867110, 0.874130, 0.972650, 0.883420],
        [0.627150, 0.627150, 0.700000, 0.774070, 0.845130],
        [0.973700, 1.006240, 1.026190, 1.071960, 1.017240]]

    coeffs[5, 3, :, :] = [
        [1.028710, 1.017570, 1.025900, 1.081790, 1.024240],
        [0.924980, 0.985500, 1.014100, 1.092210, 0.999610],
        [0.828570, 0.934920, 0.994950, 1.024590, 0.949710],
        [0.900810, 0.901330, 0.928830, 0.979570, 0.913100],
        [0.761030, 0.845150, 0.805360, 0.936790, 0.853460],
        [0.626400, 0.546750, 0.730500, 0.850000, 0.689050],
        [0.957630, 0.985480, 0.991790, 1.050220, 0.987900]]

    coeffs[5, 4, :, :] = [
        [0.992730, 0.993880, 1.017150, 1.059120, 1.017450],
        [0.975610, 0.987160, 1.026820, 1.075440, 1.007250],
        [0.871090, 0.933190, 0.974690, 0.979840, 0.952730],
        [0.828750, 0.868090, 0.834920, 0.905510, 0.871530],
        [0.781540, 0.782470, 0.767910, 0.764140, 0.795890],
        [0.743460, 0.693390, 0.514870, 0.630150, 0.715660],
        [0.934760, 0.957870, 0.959640, 0.972510, 0.981640]]

    coeffs[5, 5, :, :] = [
        [0.965840, 0.941240, 0.987100, 1.022540, 1.011160],
        [0.988630, 0.994770, 0.976590, 0.950000, 1.034840],
        [0.958200, 1.018080, 0.974480, 0.920000, 0.989870],
        [0.811720, 0.869090, 0.812020, 0.850000, 0.821050],
        [0.682030, 0.679480, 0.632450, 0.746580, 0.738550],
        [0.668290, 0.445860, 0.500000, 0.678920, 0.696510],
        [0.926940, 0.953350, 0.959050, 0.876210, 0.991490]]

    coeffs[5, 6, :, :] = [
        [0.948940, 0.997760, 0.850000, 0.826520, 0.998470],
        [1.017860, 0.970000, 0.850000, 0.700000, 0.988560],
        [1.000000, 0.950000, 0.850000, 0.606240, 0.947260],
        [1.000000, 0.746140, 0.751740, 0.598390, 0.725230],
        [0.922210, 0.500000, 0.376800, 0.517110, 0.548630],
        [0.500000, 0.450000, 0.429970, 0.404490, 0.539940],
        [0.960430, 0.881630, 0.775640, 0.596350, 0.937680]]

    coeffs[6, 1, :, :] = [
        [1.030000, 1.040000, 1.000000, 1.000000, 1.049510],
        [1.050000, 0.990000, 0.990000, 0.950000, 0.996530],
        [1.050000, 0.990000, 0.990000, 0.820000, 0.971940],
        [1.050000, 0.790000, 0.880000, 0.820000, 0.951840],
        [1.000000, 0.530000, 0.440000, 0.710000, 0.928730],
        [0.540000, 0.470000, 0.500000, 0.550000, 0.773950],
        [1.038270, 0.920180, 0.910930, 0.821140, 1.034560]]

    coeffs[6, 2, :, :] = [
        [1.041020, 0.997520, 0.961600, 1.000000, 1.035780],
        [0.948030, 0.980000, 0.900000, 0.950360, 0.977460],
        [0.950000, 0.977250, 0.869270, 0.800000, 0.951680],
        [0.951870, 0.850000, 0.748770, 0.700000, 0.883850],
        [0.900000, 0.823190, 0.727450, 0.600000, 0.839870],
        [0.850000, 0.805020, 0.692310, 0.500000, 0.788410],
        [1.010090, 0.895270, 0.773030, 0.816280, 1.011680]]

    coeffs[6, 3, :, :] = [
        [1.022450, 1.004600, 0.983650, 1.000000, 1.032940],
        [0.943960, 0.999240, 0.983920, 0.905990, 0.978150],
        [0.936240, 0.946480, 0.850000, 0.850000, 0.930320],
        [0.816420, 0.885000, 0.644950, 0.817650, 0.865310],
        [0.742960, 0.765690, 0.561520, 0.700000, 0.827140],
        [0.643870, 0.596710, 0.474460, 0.600000, 0.651200],
        [0.971740, 0.940560, 0.714880, 0.864380, 1.001650]]

    coeffs[6, 4, :, :] = [
        [0.995260, 0.977010, 1.000000, 1.000000, 1.035250],
        [0.939810, 0.975250, 0.939980, 0.950000, 0.982550],
        [0.876870, 0.879440, 0.850000, 0.900000, 0.917810],
        [0.873480, 0.873450, 0.751470, 0.850000, 0.863040],
        [0.761470, 0.702360, 0.638770, 0.750000, 0.783120],
        [0.734080, 0.650000, 0.600000, 0.650000, 0.715660],
        [0.942160, 0.919100, 0.770340, 0.731170, 0.995180]]

    coeffs[6, 5, :, :] = [
        [0.952560, 0.916780, 0.920000, 0.900000, 1.005880],
        [0.928620, 0.994420, 0.900000, 0.900000, 0.983720],
        [0.913070, 0.850000, 0.850000, 0.800000, 0.924280],
        [0.868090, 0.807170, 0.823550, 0.600000, 0.844520],
        [0.769570, 0.719870, 0.650000, 0.550000, 0.733500],
        [0.580250, 0.650000, 0.600000, 0.500000, 0.628850],
        [0.904770, 0.852650, 0.708370, 0.493730, 0.949030]]

    coeffs[6, 6, :, :] = [
        [0.911970, 0.800000, 0.800000, 0.800000, 0.956320],
        [0.912620, 0.682610, 0.750000, 0.700000, 0.950110],
        [0.653450, 0.659330, 0.700000, 0.600000, 0.856110],
        [0.648440, 0.600000, 0.641120, 0.500000, 0.695780],
        [0.570000, 0.550000, 0.598800, 0.400000, 0.560150],
        [0.475230, 0.500000, 0.518640, 0.339970, 0.520230],
        [0.743440, 0.592190, 0.603060, 0.316930, 0.794390]]

    return coeffs[1:, 1:, :, :]


def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1,
        zenith_threshold_for_zero_dni=88.0,
        zenith_threshold_for_clearsky_limit=80.0):
    """
    Determine DNI from GHI and DHI.

    When calculating the DNI from GHI and DHI the calculated DNI may be
    unreasonably high or negative for zenith angles close to 90 degrees
    (sunrise/sunset transitions). This function identifies unreasonable DNI
    values and sets them to NaN. If the clearsky DNI is given unreasonably high
    values are cut off.

    Parameters
    ----------
    ghi : Series
        Global horizontal irradiance.

    dhi : Series
        Diffuse horizontal irradiance.

    zenith : Series
        True (not refraction-corrected) zenith angles in decimal
        degrees. Angles must be >=0 and <=180.

    clearsky_dni : None or Series, default None
        Clearsky direct normal irradiance.

    clearsky_tolerance : float, default 1.1
        If 'clearsky_dni' is given this parameter can be used to allow a
        tolerance by how much the calculated DNI value can be greater than
        the clearsky value before it is identified as an unreasonable value.

    zenith_threshold_for_zero_dni : float, default 88.0
        Non-zero DNI values for zenith angles greater than or equal to
        'zenith_threshold_for_zero_dni' will be set to NaN.

    zenith_threshold_for_clearsky_limit : float, default 80.0
        DNI values for zenith angles greater than or equal to
        'zenith_threshold_for_clearsky_limit' and smaller the
        'zenith_threshold_for_zero_dni' that are greater than the clearsky DNI
        (times allowed tolerance) will be corrected. Only applies if
        'clearsky_dni' is not None.

    Returns
    -------
    dni : Series
        The modeled direct normal irradiance.
    """

    # calculate DNI
    dni = (ghi - dhi) / tools.cosd(zenith)

    # cutoff negative values
    dni[dni < 0] = float('nan')

    # set non-zero DNI values for zenith angles >=
    # zenith_threshold_for_zero_dni to NaN
    dni[(zenith >= zenith_threshold_for_zero_dni) & (dni != 0)] = float('nan')

    # correct DNI values for zenith angles greater or equal to the
    # zenith_threshold_for_clearsky_limit and smaller than the
    # upper_cutoff_zenith that are greater than the clearsky DNI (times
    # clearsky_tolerance)
    if clearsky_dni is not None:
        max_dni = clearsky_dni * clearsky_tolerance
        dni[(zenith >= zenith_threshold_for_clearsky_limit) &
            (zenith < zenith_threshold_for_zero_dni) &
            (dni > max_dni)] = max_dni
    return dni