iam.py

r"""
The ``iam`` module contains functions that implement models for the incidence
angle modifier (IAM). The IAM quantifies the fraction of direct irradiance on
a module's front surface that is transmitted through the module materials to
the cells. Stated differently, the quantity 1 - IAM is the fraction of direct
irradiance that is reflected away or absorbed by the module's front materials.
IAM is typically a function of the angle of incidence (AOI) of the direct
irradiance to the module's surface.
"""

import numpy as np
import pandas as pd
import functools
from pvlib.tools import cosd, sind

# a dict of required parameter names for each IAM model
# keys are the function names for the IAM models
_IAM_MODEL_PARAMS = {
    'ashrae': {'b'},
    'physical': {'n', 'K', 'L'},
    'martin_ruiz': {'a_r'},
    'sapm': {'B0', 'B1', 'B2', 'B3', 'B4', 'B5'},
    'interp': set()
}


def ashrae(aoi, b=0.05):
    r"""
    Determine the incidence angle modifier using the ASHRAE transmission
    model.

    The ASHRAE (American Society of Heating, Refrigeration, and Air
    Conditioning Engineers) transmission model is developed in
    [1]_, and in [2]_. The model has been used in software such as PVSyst [3]_.

    Parameters
    ----------
    aoi : numeric
        The angle of incidence (AOI) between the module normal vector and the
        sun-beam vector in degrees. Angles of nan will result in nan.

    b : float, default 0.05
        A parameter to adjust the incidence angle modifier as a function of
        angle of incidence. Typical values are on the order of 0.05 [3].

    Returns
    -------
    iam : numeric
        The incident angle modifier (IAM). Returns zero for all abs(aoi) >= 90
        and for all ``iam`` values that would be less than 0.

    Notes
    -----
    The incidence angle modifier is calculated as

    .. math::

        IAM = 1 - b (\sec(aoi) - 1)

    As AOI approaches 90 degrees, the model yields negative values for IAM;
    negative IAM values are set to zero in this implementation.

    References
    ----------
    .. [1] Souka A.F., Safwat H.H., "Determination of the optimum
       orientations for the double exposure flat-plate collector and its
       reflections". Solar Energy vol .10, pp 170-174. 1966.

    .. [2] ASHRAE standard 93-77

    .. [3] PVsyst Contextual Help.
       https://files.pvsyst.com/help/index.html?iam_loss.htm retrieved on
       October 14, 2019

    See Also
    --------
    pvlib.iam.physical
    pvlib.iam.martin_ruiz
    pvlib.iam.interp
    """

    iam = 1 - b * (1 / np.cos(np.radians(aoi)) - 1)
    aoi_gte_90 = np.full_like(aoi, False, dtype='bool')
    np.greater_equal(np.abs(aoi), 90, where=~np.isnan(aoi), out=aoi_gte_90)
    iam = np.where(aoi_gte_90, 0, iam)
    iam = np.maximum(0, iam)

    if isinstance(aoi, pd.Series):
        iam = pd.Series(iam, index=aoi.index)

    return iam


def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None):
    r"""
    Determine the incidence angle modifier using refractive index ``n``,
    extinction coefficient ``K``, glazing thickness ``L`` and refractive
    index ``n_ar`` of an optional anti-reflective coating.

    ``iam.physical`` calculates the incidence angle modifier as described in
    [1]_, Section 3, with additional support of an anti-reflective coating.
    The calculation is based on a physical model of reflections, absorption,
    and transmission through a transparent cover.

    Parameters
    ----------
    aoi : numeric
        The angle of incidence between the module normal vector and the
        sun-beam vector in degrees. Angles of nan will result in nan.

    n : numeric, default 1.526
        The effective index of refraction (unitless). Reference [1]_
        indicates that a value of 1.526 is acceptable for glass.

    K : numeric, default 4.0
        The glazing extinction coefficient in units of 1/meters.
        Reference [1] indicates that a value of 4 is reasonable for
        "water white" glass.

    L : numeric, default 0.002
        The glazing thickness in units of meters. Reference [1]_
        indicates that 0.002 meters (2 mm) is reasonable for most
        glass-covered PV panels.

    n_ar : numeric, optional
        The effective index of refraction of the anti-reflective (AR) coating
        (unitless). If n_ar is None (default), no AR coating is applied.
        A typical value for the effective index of an AR coating is 1.29.

    Returns
    -------
    iam : numeric
        The incident angle modifier

    Notes
    -----
    The pvlib python authors believe that Eqn. 14 in [1]_ is
    incorrect, which presents :math:`\theta_{r} = \arcsin(n \sin(AOI))`.
    Here, :math:`\theta_{r} = \arcsin(1/n \times \sin(AOI))`

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

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

    See Also
    --------
    pvlib.iam.martin_ruiz
    pvlib.iam.ashrae
    pvlib.iam.interp
    pvlib.iam.sapm
    """
    n1, n3 = 1, n
    if n_ar is None or np.allclose(n_ar, n1):
        # no AR coating
        n2 = n
    else:
        n2 = n_ar

    # incidence angle
    costheta = np.maximum(0, cosd(aoi))  # always >= 0
    sintheta = np.sqrt(1 - costheta**2)  # always >= 0
    n1costheta1 = n1 * costheta
    n2costheta1 = n2 * costheta

    # refraction angle of first interface
    sintheta = n1 / n2 * sintheta
    costheta = np.sqrt(1 - sintheta**2)
    n1costheta2 = n1 * costheta
    n2costheta2 = n2 * costheta

    # reflectance of s-, p-polarized, and normal light by the first interface
    with np.errstate(divide='ignore', invalid='ignore'):
        rho12_s = \
            ((n1costheta1 - n2costheta2) / (n1costheta1 + n2costheta2)) ** 2
        rho12_p = \
            ((n1costheta2 - n2costheta1) / (n1costheta2 + n2costheta1)) ** 2

    rho12_0 = ((n1 - n2) / (n1 + n2)) ** 2

    # transmittance through the first interface
    tau_s = 1 - rho12_s
    tau_p = 1 - rho12_p
    tau_0 = 1 - rho12_0

    if not np.allclose(n3, n2):  # AR coated glass
        n3costheta2 = n3 * costheta
        # refraction angle of second interface
        sintheta = n2 / n3 * sintheta
        costheta = np.sqrt(1 - sintheta**2)
        n2costheta3 = n2 * costheta
        n3costheta3 = n3 * costheta

        # reflectance by the second interface
        rho23_s = (
            (n2costheta2 - n3costheta3) / (n2costheta2 + n3costheta3)
        ) ** 2
        rho23_p = (
            (n2costheta3 - n3costheta2) / (n2costheta3 + n3costheta2)
        ) ** 2
        rho23_0 = ((n2 - n3) / (n2 + n3)) ** 2

        # transmittance through the coating, including internal reflections
        # 1 + rho23*rho12 + (rho23*rho12)^2 + ... = 1/(1 - rho23*rho12)
        tau_s *= (1 - rho23_s) / (1 - rho23_s * rho12_s)
        tau_p *= (1 - rho23_p) / (1 - rho23_p * rho12_p)
        tau_0 *= (1 - rho23_0) / (1 - rho23_0 * rho12_0)

    # transmittance after absorption in the glass
    with np.errstate(divide='ignore', invalid='ignore'):
        tau_s *= np.exp(-K * L / costheta)
        tau_p *= np.exp(-K * L / costheta)

    tau_0 *= np.exp(-K * L)

    # incidence angle modifier
    iam = (tau_s + tau_p) / 2 / tau_0

    # for light coming from behind the plane, none can enter the module
    # when n2 > 1, this is already the case
    if np.isclose(n2, 1).any():
        iam = np.where(aoi >= 90, 0, iam)
        if isinstance(aoi, pd.Series):
            iam = pd.Series(iam, index=aoi.index)

    return iam


def martin_ruiz(aoi, a_r=0.16):
    r'''
    Determine the incidence angle modifier (IAM) using the Martin
    and Ruiz incident angle model.

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

    a_r : numeric
        The angular losses coefficient described in equation 3 of [1]_.
        This is an empirical dimensionless parameter. Values of ``a_r`` are
        generally on the order of 0.08 to 0.25 for flat-plate PV modules.

    Returns
    -------
    iam : numeric
        The incident angle modifier(s)

    Notes
    -----
    `martin_ruiz` calculates the incidence angle modifier (IAM) as described in
    [1]_. The information required is the incident angle (AOI) and the angular
    losses coefficient (a_r). Note that [1]_ has a corrigendum [2]_ which
    clarifies a mix-up of 'alpha's and 'a's in the former.

    The incident angle modifier is defined as

    .. math::

       IAM = \frac{1 - \exp(-\frac{\cos(aoi)}{a_r})}
       {1 - \exp(\frac{-1}{a_r})}

    which is presented as :math:`AL(\alpha) = 1 - IAM` in equation 4 of [1]_,
    with :math:`\alpha` representing the angle of incidence AOI. Thus IAM = 1
    at AOI = 0, and IAM = 0 at AOI = 90.  This equation is only valid for
    -90 <= aoi <= 90, therefore `iam` is constrained to 0.0 outside this
    interval.

    References
    ----------
    .. [1] N. Martin and J. M. Ruiz, "Calculation of the PV modules angular
       losses under field conditions by means of an analytical model", Solar
       Energy Materials & Solar Cells, vol. 70, pp. 25-38, 2001.

    .. [2] N. Martin and J. M. Ruiz, "Corrigendum to 'Calculation of the PV
       modules angular losses under field conditions by means of an
       analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
       pp. 154, 2013.

    See Also
    --------
    pvlib.iam.martin_ruiz_diffuse
    pvlib.iam.physical
    pvlib.iam.ashrae
    pvlib.iam.interp
    pvlib.iam.sapm
    '''
    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019

    aoi_input = aoi

    aoi = np.asanyarray(aoi)
    a_r = np.asanyarray(a_r)

    if np.any(np.less_equal(a_r, 0)):
        raise ValueError("The parameter 'a_r' cannot be zero or negative.")

    with np.errstate(invalid='ignore'):
        iam = (1 - np.exp(-cosd(aoi) / a_r)) / (1 - np.exp(-1 / a_r))
        iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)

    if isinstance(aoi_input, pd.Series):
        iam = pd.Series(iam, index=aoi_input.index)

    return iam


def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
    '''
    Determine the incidence angle modifiers (iam) for diffuse sky and
    ground-reflected irradiance using the Martin and Ruiz incident angle model.

    Parameters
    ----------
    surface_tilt : float or array-like, default 0
        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_tilt must be in the range [0, 180]

    a_r : numeric
        The angular losses coefficient described in equation 3 of [1]_.
        This is an empirical dimensionless parameter. Values of a_r are
        generally on the order of 0.08 to 0.25 for flat-plate PV modules.
        a_r must be greater than zero.

    c1 : float
        First fitting parameter for the expressions that approximate the
        integral of diffuse irradiance coming from different directions.
        c1 is given as the constant 4 / 3 / pi (0.4244) in [1]_.

    c2 : float
        Second fitting parameter for the expressions that approximate the
        integral of diffuse irradiance coming from different directions.
        If c2 is None, it will be calculated according to the linear
        relationship given in [3]_.

    Returns
    -------
    iam_sky : numeric
        The incident angle modifier for sky diffuse

    iam_ground : numeric
        The incident angle modifier for ground-reflected diffuse

    Notes
    -----
    Sky and ground modifiers are complementary: iam_sky for tilt = 30 is
    equal to iam_ground for tilt = 180 - 30.  For vertical surfaces,
    tilt = 90, the two factors are equal.

    References
    ----------
    .. [1] N. Martin and J. M. Ruiz, "Calculation of the PV modules angular
       losses under field conditions by means of an analytical model", Solar
       Energy Materials & Solar Cells, vol. 70, pp. 25-38, 2001.

    .. [2] N. Martin and J. M. Ruiz, "Corrigendum to 'Calculation of the PV
       modules angular losses under field conditions by means of an
       analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
       pp. 154, 2013.

    .. [3] "IEC 61853-3 Photovoltaic (PV) module performance testing and energy
       rating - Part 3: Energy rating of PV modules". IEC, Geneva, 2018.

    See Also
    --------
    pvlib.iam.martin_ruiz
    pvlib.iam.physical
    pvlib.iam.ashrae
    pvlib.iam.interp
    pvlib.iam.sapm
    '''
    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019

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

    surface_tilt = np.asanyarray(surface_tilt)

    # avoid undefined results for horizontal or upside-down surfaces
    zeroang = 1e-06

    surface_tilt = np.where(surface_tilt == 0, zeroang, surface_tilt)
    surface_tilt = np.where(surface_tilt == 180, 180 - zeroang, surface_tilt)

    if c2 is None:
        # This equation is from [3] Sect. 7.2
        c2 = 0.5 * a_r - 0.154

    beta = np.radians(surface_tilt)
    sin = np.sin
    pi = np.pi
    cos = np.cos

    # avoid RuntimeWarnings for <, sin, and cos with nan
    with np.errstate(invalid='ignore'):
        # because sin(pi) isn't exactly zero
        sin_beta = np.where(surface_tilt < 90, sin(beta), sin(pi - beta))

        trig_term_sky = sin_beta + (pi - beta - sin_beta) / (1 + cos(beta))
        trig_term_gnd = sin_beta +      (beta - sin_beta) / (1 - cos(beta))  # noqa: E222 E261 E501

    iam_sky = 1 - np.exp(-(c1 + c2 * trig_term_sky) * trig_term_sky / a_r)
    iam_gnd = 1 - np.exp(-(c1 + c2 * trig_term_gnd) * trig_term_gnd / a_r)

    if out_index is not None:
        iam_sky = pd.Series(iam_sky, index=out_index, name='iam_sky')
        iam_gnd = pd.Series(iam_gnd, index=out_index, name='iam_ground')

    return iam_sky, iam_gnd


def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
    r'''
    Determine the incidence angle modifier (IAM) by interpolating a set of
    reference values, which are usually measured values.

    Parameters
    ----------
    aoi : numeric
        The angle of incidence between the module normal vector and the
        sun-beam vector [degrees].

    theta_ref : numeric
        Vector of angles at which the IAM is known [degrees].

    iam_ref : numeric
        IAM values for each angle in ``theta_ref`` [unitless].

    method : str, default 'linear'
        Specifies the interpolation method.
        Useful options are: 'linear', 'quadratic', 'cubic'.
        See scipy.interpolate.interp1d for more options.

    normalize : boolean, default True
        When true, the interpolated values are divided by the interpolated
        value at zero degrees.  This ensures that ``iam=1.0`` at normal
        incidence.

    Returns
    -------
    iam : numeric
        The incident angle modifier(s) [unitless]

    Notes
    -----
    ``theta_ref`` must have two or more points and may span any range of
    angles. Typically there will be a dozen or more points in the range 0-90
    degrees. Beyond the range of ``theta_ref``, IAM values are extrapolated,
    but constrained to be non-negative.

    The sign of ``aoi`` is ignored; only the magnitude is used.

    See Also
    --------
    pvlib.iam.physical
    pvlib.iam.ashrae
    pvlib.iam.martin_ruiz
    pvlib.iam.sapm
    '''
    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019

    from scipy.interpolate import interp1d

    # Scipy doesn't give the clearest feedback, so check number of points here.
    MIN_REF_VALS = {'linear': 2, 'quadratic': 3, 'cubic': 4, 1: 2, 2: 3, 3: 4}

    if len(theta_ref) < MIN_REF_VALS.get(method, 2):
        raise ValueError("Too few reference points defined "
                         "for interpolation method '%s'." % method)

    if np.any(np.less(iam_ref, 0)):
        raise ValueError("Negative value(s) found in 'iam_ref'. "
                         "This is not physically possible.")

    interpolator = interp1d(theta_ref, iam_ref, kind=method,
                            fill_value='extrapolate')
    aoi_input = aoi

    aoi = np.asanyarray(aoi)
    aoi = np.abs(aoi)
    iam = interpolator(aoi)
    iam = np.clip(iam, 0, None)

    if normalize:
        iam /= interpolator(0)

    if isinstance(aoi_input, pd.Series):
        iam = pd.Series(iam, index=aoi_input.index)

    return iam


def sapm(aoi, module, upper=None):
    r"""
    Determine the incidence angle modifier (IAM) using the SAPM model.

    Parameters
    ----------
    aoi : numeric
        Angle of incidence in degrees. Negative input angles will return
        zeros.

    module : dict-like
        A dict or Series with the SAPM IAM model parameters.
        See the :py:func:`sapm` notes section for more details.

    upper : float, optional
        Upper limit on the results.

    Returns
    -------
    iam : numeric
        The SAPM angle of incidence loss coefficient, termed F2 in [1]_.

    Notes
    -----
    The SAPM [1]_ traditionally does not define an upper limit on the AOI
    loss function and values slightly exceeding 1 may exist for moderate
    angles of incidence (15-40 degrees). However, users may consider
    imposing an upper limit of 1.

    References
    ----------
    .. [1] King, D. et al, 2004, "Sandia Photovoltaic Array Performance
       Model", SAND Report 3535, Sandia National Laboratories, Albuquerque,
       NM.

    .. [2] B.H. King et al, "Procedure to Determine Coefficients for the
       Sandia Array Performance Model (SAPM)," SAND2016-5284, Sandia
       National Laboratories (2016).

    .. [3] B.H. King et al, "Recent Advancements in Outdoor Measurement
       Techniques for Angle of Incidence Effects," 42nd IEEE PVSC (2015).
       DOI: 10.1109/PVSC.2015.7355849

    See Also
    --------
    pvlib.iam.physical
    pvlib.iam.ashrae
    pvlib.iam.martin_ruiz
    pvlib.iam.interp
    """

    aoi_coeff = [module['B5'], module['B4'], module['B3'], module['B2'],
                 module['B1'], module['B0']]

    iam = np.polyval(aoi_coeff, aoi)
    iam = np.clip(iam, 0, upper)
    # nan tolerant masking
    aoi_lt_0 = np.full_like(aoi, False, dtype='bool')
    np.less(aoi, 0, where=~np.isnan(aoi), out=aoi_lt_0)
    iam = np.where(aoi_lt_0, 0, iam)

    if isinstance(aoi, pd.Series):
        iam = pd.Series(iam, aoi.index)

    return iam


def marion_diffuse(model, surface_tilt, **kwargs):
    """
    Determine diffuse irradiance incidence angle modifiers using Marion's
    method of integrating over solid angle.

    Parameters
    ----------
    model : str
        The IAM function to evaluate across solid angle. Must be one of
        `'ashrae', 'physical', 'martin_ruiz', 'sapm', 'schlick'`.

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

    **kwargs
        Extra parameters passed to the IAM function.

    Returns
    -------
    iam : dict
        IAM values for each type of diffuse irradiance:

            * 'sky': radiation from the sky dome (zenith <= 90)
            * 'horizon': radiation from the region of the sky near the horizon
              (89.5 <= zenith <= 90)
            * 'ground': radiation reflected from the ground (zenith >= 90)

        See [1]_ for a detailed description of each class.

    See Also
    --------
    pvlib.iam.marion_integrate

    References
    ----------
    .. [1] B. Marion "Numerical method for angle-of-incidence correction
       factors for diffuse radiation incident photovoltaic modules",
       Solar Energy, Volume 147, Pages 344-348. 2017.
       DOI: 10.1016/j.solener.2017.03.027

    Examples
    --------
    >>> marion_diffuse('physical', surface_tilt=20)
    {'sky': 0.9539178294437575,
     'horizon': 0.7652650139134007,
     'ground': 0.6387140117795903}

    >>> marion_diffuse('ashrae', [20, 30], b=0.04)
    {'sky': array([0.96748999, 0.96938408]),
     'horizon': array([0.86478428, 0.91825792]),
     'ground': array([0.77004435, 0.8522436 ])}
    """

    models = {
        'physical': physical,
        'ashrae': ashrae,
        'sapm': sapm,
        'martin_ruiz': martin_ruiz,
        'schlick': schlick,
    }

    try:
        iam_model = models[model]
    except KeyError:
        raise ValueError('model must be one of: ' + str(list(models.keys())))

    iam_function = functools.partial(iam_model, **kwargs)
    iam = {}
    for region in ['sky', 'horizon', 'ground']:
        iam[region] = marion_integrate(iam_function, surface_tilt, region)

    return iam


def marion_integrate(function, surface_tilt, region, num=None):
    """
    Integrate an incidence angle modifier (IAM) function over solid angle
    to determine a diffuse irradiance correction factor using Marion's method.

    This lower-level function actually performs the IAM integration for the
    specified solid angle region.

    Parameters
    ----------
    function : callable(aoi)
        The IAM function to evaluate across solid angle. The function must
        be vectorized and take only one parameter, the angle of incidence in
        degrees.

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

    region : {'sky', 'horizon', 'ground'}
        The region to integrate over. Must be one of:

            * 'sky': radiation from the sky dome (zenith <= 90)
            * 'horizon': radiation from the region of the sky near the horizon
              (89.5 <= zenith <= 90)
            * 'ground': radiation reflected from the ground (zenith >= 90)

        See [1]_ for a detailed description of each class.

    num : int, optional
        The number of increments in the zenith integration.
        If not specified, N will follow the values used in [1]_:

            * 'sky' or 'ground': num = 180
            * 'horizon': num = 1800

    Returns
    -------
    iam : numeric
        AOI diffuse correction factor for the specified region.

    See Also
    --------
    pvlib.iam.marion_diffuse

    References
    ----------
    .. [1] B. Marion "Numerical method for angle-of-incidence correction
       factors for diffuse radiation incident photovoltaic modules",
       Solar Energy, Volume 147, Pages 344-348. 2017.
       DOI: 10.1016/j.solener.2017.03.027

    Examples
    --------
    >>> marion_integrate(pvlib.iam.ashrae, 20, 'sky')
    0.9596085829811408

    >>> from functools import partial
    >>> f = partial(pvlib.iam.physical, n=1.3)
    >>> marion_integrate(f, [20, 30], 'sky')
    array([0.96225034, 0.9653219 ])
    """

    if num is None:
        if region in ['sky', 'ground']:
            num = 180
        elif region == 'horizon':
            num = 1800
        else:
            raise ValueError(f'Invalid region: {region}')

    beta = np.radians(surface_tilt)
    if isinstance(beta, pd.Series):
        # convert Series to np array for broadcasting later
        beta = beta.values
    ai = np.pi/num  # angular increment

    phi_range = np.linspace(0, np.pi, num, endpoint=False)
    psi_range = np.linspace(0, 2*np.pi, 2*num, endpoint=False)

    # the pseudocode in [1] do these checks at the end, but it's
    # faster to do this criteria check up front instead of later.
    if region == 'sky':
        mask = phi_range + ai <= np.pi/2
    elif region == 'horizon':
        lo = 89.5 * np.pi/180
        hi = np.pi/2
        mask = (lo <= phi_range) & (phi_range + ai <= hi)
    elif region == 'ground':
        mask = (phi_range >= np.pi/2)
    else:
        raise ValueError(f'Invalid region: {region}')
    phi_range = phi_range[mask]

    # fast Cartesian product of phi and psi
    angles = np.array(np.meshgrid(phi_range, psi_range)).T.reshape(-1, 2)
    # index with single-element lists to maintain 2nd dimension so that
    # these angle arrays broadcast across the beta array
    phi_1 = angles[:, [0]]
    psi_1 = angles[:, [1]]
    phi_2 = phi_1 + ai
    # psi_2 = psi_1 + ai  # not needed
    phi_avg = phi_1 + 0.5*ai
    psi_avg = psi_1 + 0.5*ai
    term_1 = np.cos(beta) * np.cos(phi_avg)
    # The AOI formula includes a term based on the difference between
    # panel azimuth and the photon azimuth, but because we assume each class
    # of diffuse irradiance is isotropic and we are integrating over all
    # angles, it doesn't matter what panel azimuth we choose (i.e., the
    # system is rotationally invariant).  So we choose gamma to be zero so
    # that we can omit it from the cos(psi_avg) term.
    # Marion's paper mentions this in the Section 3 pseudocode:
    # "set gamma to pi (or any value between 0 and 2pi)"
    term_2 = np.sin(beta) * np.sin(phi_avg) * np.cos(psi_avg)
    cosaoi = term_1 + term_2
    aoi = np.arccos(cosaoi)
    # simplify Eq 8, (psi_2 - psi_1) is always ai
    dAs = ai * (np.cos(phi_1) - np.cos(phi_2))
    cosaoi_dAs = cosaoi * dAs
    # apply the final AOI check, zeroing out non-passing points
    mask = aoi < np.pi/2
    cosaoi_dAs = np.where(mask, cosaoi_dAs, 0)
    numerator = np.sum(function(np.degrees(aoi)) * cosaoi_dAs, axis=0)
    denominator = np.sum(cosaoi_dAs, axis=0)

    with np.errstate(invalid='ignore'):
        # in some cases, no points pass the criteria
        # (e.g. region='ground', surface_tilt=0), so we override the division
        # by zero to set Fd=0.  Also, preserve nans in beta.
        Fd = np.where((denominator != 0) | ~np.isfinite(beta),
                      numerator / denominator,
                      0)

    # preserve input type
    if np.isscalar(surface_tilt):
        Fd = Fd.item()
    elif isinstance(surface_tilt, pd.Series):
        Fd = pd.Series(Fd, surface_tilt.index)

    return Fd


def schlick(aoi):
    """
    Determine incidence angle modifier (IAM) for direct irradiance using the
    Schlick approximation to the Fresnel equations.

    The Schlick approximation was proposed in [1]_ as a computationally
    efficient alternative to computing the Fresnel factor in computer
    graphics contexts.  This implementation is a normalized form of the
    equation in [1]_ so that it can be used as a PV IAM model.
    Unlike other IAM models, this model has no ability to describe
    different reflection profiles.

    In PV contexts, the Schlick approximation has been used as an analytically
    integrable alternative to the Fresnel equations for estimating IAM
    for diffuse irradiance [2]_.

    Parameters
    ----------
    aoi : numeric
        The angle of incidence (AOI) between the module normal vector and the
        sun-beam vector. Angles of nan will result in nan. [degrees]

    Returns
    -------
    iam : numeric
        The incident angle modifier.

    References
    ----------
    .. [1] Schlick, C. An inexpensive BRDF model for physically-based
       rendering. Computer graphics forum 13 (1994).

    .. [2] Xie, Y., M. Sengupta, A. Habte, A. Andreas, "The 'Fresnel Equations'
       for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
       Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
       :doi:`10.1016/j.rser.2022.112362`

    See Also
    --------
    pvlib.iam.schlick_diffuse
    """
    iam = 1 - (1 - cosd(aoi)) ** 5
    iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)

    # preserve input type
    if np.isscalar(aoi):
        iam = iam.item()
    elif isinstance(aoi, pd.Series):
        iam = pd.Series(iam, aoi.index)

    return iam


def schlick_diffuse(surface_tilt):
    """
    Determine the incidence angle modifiers (IAM) for diffuse sky and
    ground-reflected irradiance on a tilted surface using the Schlick
    incident angle model.

    The diffuse iam values are calculated using an analytical integration
    of the Schlick equation [1]_ over the portion of an isotropic sky and
    isotropic foreground that is visible from the tilted surface [2]_.

    Parameters
    ----------
    surface_tilt : numeric
        Surface tilt angle measured from horizontal (e.g. surface facing
        up = 0, surface facing horizon = 90). [degrees]

    Returns
    -------
    iam_sky : numeric
        The incident angle modifier for sky diffuse.

    iam_ground : numeric
        The incident angle modifier for ground-reflected diffuse.

    References
    ----------
    .. [1] Schlick, C. An inexpensive BRDF model for physically-based
       rendering. Computer graphics forum 13 (1994).

    .. [2] Xie, Y., M. Sengupta, A. Habte, A. Andreas, "The 'Fresnel Equations'
       for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
       Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
       :doi:`10.1016/j.rser.2022.112362`

    See Also
    --------
    pvlib.iam.schlick
    """
    # these calculations are as in [2]_, but with the refractive index
    # weighting coefficient w set to 1.0 (so it is omitted)

    # relative transmittance of sky diffuse radiation by PV cover:
    cosB = cosd(surface_tilt)
    sinB = sind(surface_tilt)
    cuk = (2 / (np.pi * (1 + cosB))) * (
        (30/7)*np.pi - (160/21)*np.radians(surface_tilt) - (10/3)*np.pi*cosB
        + (160/21)*cosB*sinB - (5/3)*np.pi*cosB*sinB**2 + (20/7)*cosB*sinB**3
        - (5/16)*np.pi*cosB*sinB**4 + (16/105)*cosB*sinB**5
    )  # Eq 4 in [2]

    # relative transmittance of ground-reflected radiation by PV cover:
    with np.errstate(divide='ignore', invalid='ignore'):  # Eq 6 in [2]
        cug = 40 / (21 * (1 - cosB)) - (1 + cosB) / (1 - cosB) * cuk

    cug = np.where(surface_tilt < 1e-6, 0, cug)

    # respect input types:
    if np.isscalar(surface_tilt):
        cuk = cuk.item()
        cug = cug.item()
    elif isinstance(surface_tilt, pd.Series):
        cuk = pd.Series(cuk, surface_tilt.index)
        cug = pd.Series(cug, surface_tilt.index)

    return cuk, cug