infinite_sheds.py

r"""
Functions for the infinite sheds bifacial irradiance model.
"""

import numpy as np
import pandas as pd
from pvlib.tools import cosd, sind, tand
from pvlib.bifacial import utils
from pvlib.irradiance import beam_component, aoi, haydavies


def _poa_ground_shadows(ghi, dhi, albedo, f_gnd_beam, vf_gnd_sky):
    """
    Reduce ground-reflected irradiance to the tilted plane (poa_ground) to
    account for shadows on the ground.

    Parameters
    ----------
    ghi : numeric
        Global horizontal irradiance, with no adjustments. [W/m^2]
    dhi : numeric
        Diffuse horizontal irradiance, with no adjustments. [W/m^2]
    albedo : numeric
        Ground albedo, the ratio of reflected to incident irradiance of the
        ground surface. [W/m^2]
    f_gnd_beam : numeric
        Fraction of the distance between rows that is illuminated (unshaded).
        [unitless]
    vf_gnd_sky : numeric
        View factor from the ground to the sky, integrated along the distance
        between rows. [unitless]

    Returns
    -------
    poa_gnd_sky : numeric
        Adjusted ground-reflected irradiance accounting for shadows on the
        ground. [W/m^2]

    """
    return albedo * (f_gnd_beam * (ghi - dhi) + vf_gnd_sky * dhi)


def _poa_sky_diffuse_pv(dhi, gcr, surface_tilt):
    r"""
    Integrated view factors from the shaded and unshaded parts of
    the row slant height to the sky.

    Parameters
    ----------
    f_x : numeric
        Fraction of row slant height from the bottom that is shaded from
        direct irradiance. [unitless]
    surface_tilt : numeric
        Surface tilt angle in degrees from horizontal, e.g., surface facing up
        = 0, surface facing horizon = 90. [degree]
    gcr : float
        Ratio of row slant length to row spacing (pitch). [unitless]
    npoints : int, default 100
        Number of points for integration. [unitless]

    A detailed calculation would be

        dhi * (f_x * vf_shade_sky_integ + (1 - f_x) * vf_noshade_sky_integ)

    where vf_shade_sky_integ is the average view factor between 0 and f_x
    (the shaded portion). But the average view factor is

        1/(f_x - 0) Integral_0^f_x vf(x) dx

    so the detailed calculation is equivalent to

        dhi * 1/(1 - 0) Integral_0^1 vf(x) dx

    Parameters
    ----------
    f_x : numeric
        Fraction of row slant height from the bottom that is shaded from
        direct irradiance. [unitless]
    dhi : numeric
        Diffuse horizontal irradiance (DHI). [W/m^2]
    gcr : float
        ground coverage ratio, ratio of row slant length to row spacing.
        [unitless]
    surface_tilt : numeric
        Surface tilt angle in degrees from horizontal, e.g., surface facing up
        = 0, surface facing horizon = 90. [degree]

    Returns
    -------
    poa_sky_diffuse_pv : numeric
        Total sky diffuse irradiance incident on the PV surface. [W/m^2]
    """
    vf_integ = utils.vf_row_sky_2d_integ(surface_tilt, gcr, 0., 1.)
    return dhi * vf_integ


def _poa_ground_pv(poa_ground, gcr, surface_tilt):
    """
    Reduce ground-reflected irradiance to account for limited view of the
    ground from the row surface.

    Parameters
    ----------
    poa_ground : numeric
        Ground-reflected irradiance that would reach the row surface if the
        full ground was visible. poa_gnd_sky accounts for limited view of the
        sky from the ground. [W/m^2]
    gcr : float
        ground coverage ratio, ratio of row slant length to row spacing.
        [unitless]
    surface_tilt : numeric
        Surface tilt angle in degrees from horizontal, e.g., surface facing up
        = 0, surface facing horizon = 90. [degree]

    Returns
    -------
    numeric
        Ground diffuse irradiance on the row plane. [W/m^2]
    """
    vf_integ = utils.vf_row_ground_2d_integ(surface_tilt, gcr, 0., 1.)
    return poa_ground * vf_integ


def _shaded_fraction(solar_zenith, solar_azimuth, surface_tilt,
                     surface_azimuth, gcr):
    """
    Calculate fraction (from the bottom) of row slant height that is shaded
    from direct irradiance by the row in front toward the sun.

    See [1], Eq. 14 and also [2], Eq. 32.

    .. math::
        F_x = \\max \\left( 0, \\min \\left(\\frac{\\text{GCR} \\cos \\theta
        + \\left( \\text{GCR} \\sin \\theta - \\tan \\beta_{c} \\right)
        \\tan Z - 1}
        {\\text{GCR} \\left( \\cos \\theta + \\sin \\theta \\tan Z \\right)},
        1 \\right) \\right)

    Parameters
    ----------
    solar_zenith : numeric
        Apparent (refraction-corrected) solar zenith. [degrees]
    solar_azimuth : numeric
        Solar azimuth. [degrees]
    surface_tilt : numeric
        Row tilt from horizontal, e.g. surface facing up = 0, surface facing
        horizon = 90. [degrees]
    surface_azimuth : numeric
        Azimuth angle of the row surface. North=0, East=90, South=180,
        West=270. [degrees]
    gcr : numeric
        Ground coverage ratio, which is the ratio of row slant length to row
        spacing (pitch). [unitless]

    Returns
    -------
    f_x : numeric
        Fraction of row slant height from the bottom that is shaded from
        direct irradiance.

    References
    ----------
    .. [1] Mikofski, M., Darawali, R., Hamer, M., Neubert, A., and Newmiller,
       J. "Bifacial Performance Modeling in Large Arrays". 2019 IEEE 46th
       Photovoltaic Specialists Conference (PVSC), 2019, pp. 1282-1287.
       :doi:`10.1109/PVSC40753.2019.8980572`.
    .. [2] Kevin Anderson and Mark Mikofski, "Slope-Aware Backtracking for
       Single-Axis Trackers", Technical Report NREL/TP-5K00-76626, July 2020.
       :doi:`10.2172/1660126`
    """
    tan_phi = utils._solar_projection_tangent(
        solar_zenith, solar_azimuth, surface_azimuth)
    # length of shadow behind a row as a fraction of pitch
    x = gcr * (sind(surface_tilt) * tan_phi + cosd(surface_tilt))
    f_x = 1 - 1. / x
    # set f_x to be 1 when sun is behind the array
    ao = aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)
    f_x = np.where(ao < 90, f_x, 1.)
    # when x < 1, the shadow is not long enough to fall on the row surface
    f_x = np.where(x > 1., f_x, 0.)
    return f_x


def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
                       solar_azimuth, gcr, height, pitch, ghi, dhi, dni,
                       albedo, model='isotropic', dni_extra=None, iam=1.0,
                       npoints=100, vectorize=False):
    r"""
    Calculate plane-of-array (POA) irradiance on one side of a row of modules.

    The infinite sheds model [1] assumes the PV system comprises parallel,
    evenly spaced rows on a level, horizontal surface. Rows can be on fixed
    racking or single axis trackers. The model calculates irradiance at a
    location far from the ends of any rows, in effect, assuming that the
    rows (sheds) are infinitely long.

    POA irradiance components include direct, diffuse and global (total).
    Irradiance values are reduced to account for reflection of direct light,
    but are not adjusted for solar spectrum or reduced by a module's
    bifaciality factor.

    Parameters
    ----------
    surface_tilt : numeric
        Tilt of the surface from horizontal. Must be between 0 and 180. For
        example, for a fixed tilt module mounted at 30 degrees from
        horizontal, use ``surface_tilt=30`` to get front-side irradiance and
        ``surface_tilt=150`` to get rear-side irradiance. [degree]

    surface_azimuth : numeric
        Surface azimuth in decimal degrees east of north
        (e.g. North = 0, South = 180, East = 90, West = 270). [degree]

    solar_zenith : numeric
        Refraction-corrected solar zenith. [degree]

    solar_azimuth : numeric
        Solar azimuth. [degree]

    gcr : float
        Ground coverage ratio, ratio of row slant length to row spacing.
        [unitless]

    height : float
        Height of the center point of the row above the ground; must be in the
        same units as ``pitch``.

    pitch : float
        Distance between two rows; must be in the same units as ``height``.

    ghi : numeric
        Global horizontal irradiance. [W/m2]

    dhi : numeric
        Diffuse horizontal irradiance. [W/m2]

    dni : numeric
        Direct normal irradiance. [W/m2]

    albedo : numeric
        Surface albedo. [unitless]

    model : str, default 'isotropic'
        Irradiance model - can be one of 'isotropic' or 'haydavies'.

    dni_extra : numeric, optional
        Extraterrestrial direct normal irradiance. Required when
        ``model='haydavies'``. [W/m2]

    iam : numeric, default 1.0
        Incidence angle modifier, the fraction of direct irradiance incident
        on the surface that is not reflected away. [unitless]

    npoints : int, default 100
        Number of discretization points for calculating integrated view
        factors.

    vectorize : bool, default False
        If True, vectorize the view factor calculation across ``surface_tilt``.
        This increases speed with the cost of increased memory usage.

    Returns
    -------
    output : dict or DataFrame
        Output is a ``pandas.DataFrame`` when ``ghi`` is a Series.
        Otherwise it is a dict of ``numpy.ndarray``
        See Notes for descriptions of content.

    Notes
    -----
    Input parameters ``height`` and ``pitch`` must have the same unit.

    ``output`` always includes:

    - ``poa_global`` : total POA irradiance. [W/m^2]
    - ``poa_diffuse`` : total diffuse POA irradiance from all sources. [W/m^2]
    - ``poa_direct`` : total direct POA irradiance. [W/m^2]
    - ``poa_sky_diffuse`` : total sky diffuse irradiance on the plane of array.
      [W/m^2]
    - ``poa_ground_diffuse`` : total ground-reflected diffuse irradiance on the
      plane of array. [W/m^2]
    - ``shaded_fraction`` : fraction of row slant height from the bottom that
      is shaded from direct irradiance by adjacent rows. [unitless]

    References
    ----------
    .. [1] Mikofski, M., Darawali, R., Hamer, M., Neubert, A., and Newmiller,
       J. "Bifacial Performance Modeling in Large Arrays". 2019 IEEE 46th
       Photovoltaic Specialists Conference (PVSC), 2019, pp. 1282-1287.
       :doi:`10.1109/PVSC40753.2019.8980572`.

    See also
    --------
    get_irradiance
    """
    if model == 'haydavies':
        if dni_extra is None:
            raise ValueError(f'must supply dni_extra for {model} model')
        # Call haydavies first time within the horizontal plane - to subtract
        # circumsolar_horizontal from DHI
        sky_diffuse_comps_horizontal = haydavies(0, 180, dhi, dni, dni_extra,
                                                 solar_zenith, solar_azimuth,
                                                 return_components=True)
        circumsolar_horizontal = \
            sky_diffuse_comps_horizontal['poa_circumsolar']

        # Call haydavies a second time where circumsolar_normal is facing
        # directly towards sun, and can be added to DNI
        sky_diffuse_comps_normal = haydavies(solar_zenith, solar_azimuth, dhi,
                                             dni, dni_extra, solar_zenith,
                                             solar_azimuth,
                                             return_components=True)
        circumsolar_normal = sky_diffuse_comps_normal['poa_circumsolar']

        dhi = dhi - circumsolar_horizontal
        dni = dni + circumsolar_normal

    # Calculate some geometric quantities
    # rows to consider in front and behind current row
    # ensures that view factors to the sky are computed to within 5 degrees
    # of the horizon
    max_rows = np.ceil(height / (pitch * tand(5)))
    # fraction of ground between rows that is illuminated accounting for
    # shade from panels. [1], Eq. 4
    f_gnd_beam = utils._unshaded_ground_fraction(
        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, gcr)
    # integrated view factor from the ground to the sky, integrated between
    # adjacent rows interior to the array
    # method differs from [1], Eq. 7 and Eq. 8; height is defined at row
    # center rather than at row lower edge as in [1].
    vf_gnd_sky = utils.vf_ground_sky_2d_integ(
        surface_tilt, gcr, height, pitch, max_rows, npoints,
        vectorize)
    # fraction of row slant height that is shaded from direct irradiance
    f_x = _shaded_fraction(solar_zenith, solar_azimuth, surface_tilt,
                           surface_azimuth, gcr)

    # Total sky diffuse received by both shaded and unshaded portions
    poa_sky_pv = _poa_sky_diffuse_pv(dhi, gcr, surface_tilt)

    # irradiance reflected from the ground before accounting for shadows
    # and restricted views
    # this is a deviation from [1], because the row to ground view factor
    # is accounted for in a different manner
    # Reduce ground-reflected irradiance because other rows in the array
    # block irradiance from reaching the ground.
    # [2], Eq. 9
    ground_diffuse = _poa_ground_shadows(ghi, dhi, albedo, f_gnd_beam,
                                         vf_gnd_sky)

    # Ground-reflected irradiance on the row surface accounting for
    # the view to the ground. This deviates from [1], Eq. 10, 11 and
    # subsequent. Here, the row to ground view factor is computed. In [1],
    # the usual ground-reflected irradiance includes the single row to ground
    # view factor (1 - cos(tilt))/2, and Eq. 10, 11 and later multiply
    # this quantity by a ratio of view factors.
    poa_gnd_pv = _poa_ground_pv(ground_diffuse, gcr, surface_tilt)

    # add sky and ground-reflected irradiance on the row by irradiance
    # component
    poa_diffuse = poa_gnd_pv + poa_sky_pv
    # beam on plane, make an array for consistency with poa_diffuse
    poa_beam = np.atleast_1d(beam_component(
        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni))
    poa_direct = poa_beam * (1 - f_x) * iam  # direct only on the unshaded part
    poa_global = poa_direct + poa_diffuse

    output = {
        'poa_global': poa_global, 'poa_direct': poa_direct,
        'poa_diffuse': poa_diffuse, 'poa_ground_diffuse': poa_gnd_pv,
        'poa_sky_diffuse': poa_sky_pv, 'shaded_fraction': f_x}
    if isinstance(ghi, pd.Series):
        output = pd.DataFrame(output)
    return output


def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                   gcr, height, pitch, ghi, dhi, dni,
                   albedo, model='isotropic', dni_extra=None, iam_front=1.0,
                   iam_back=1.0, bifaciality=0.8, shade_factor=-0.02,
                   transmission_factor=0, npoints=100, vectorize=False):
    """
    Get front and rear irradiance using the infinite sheds model.

    The infinite sheds model [1] assumes the PV system comprises parallel,
    evenly spaced rows on a level, horizontal surface. Rows can be on fixed
    racking or single axis trackers. The model calculates irradiance at a
    location far from the ends of any rows, in effect, assuming that the
    rows (sheds) are infinitely long.

    The model accounts for the following effects:

    - restricted view of the sky from module surfaces due to the nearby rows.
    - restricted view of the ground from module surfaces due to nearby rows.
    - restricted view of the sky from the ground due to rows.
    - shading of module surfaces by nearby rows.
    - shading of rear cells of a module by mounting structure and by
      module features.

    The model implicitly assumes that diffuse irradiance from the sky is
    isotropic, and that module surfaces do not allow irradiance to transmit
    through the module to the ground through gaps between cells.

    Parameters
    ----------
    surface_tilt : numeric
        Tilt from horizontal of the front-side surface. [degree]

    surface_azimuth : numeric
        Surface azimuth in decimal degrees east of north
        (e.g. North = 0, South = 180, East = 90, West = 270). [degree]

    solar_zenith : numeric
        Refraction-corrected solar zenith. [degree]

    solar_azimuth : numeric
        Solar azimuth. [degree]

    gcr : float
        Ground coverage ratio, ratio of row slant length to row spacing.
        [unitless]

    height : float
        Height of the center point of the row above the ground; must be in the
        same units as ``pitch``.

    pitch : float
        Distance between two rows; must be in the same units as ``height``.

    ghi : numeric
        Global horizontal irradiance. [W/m2]

    dhi : numeric
        Diffuse horizontal irradiance. [W/m2]

    dni : numeric
        Direct normal irradiance. [W/m2]

    albedo : numeric
        Surface albedo. [unitless]

    model : str, default 'isotropic'
        Irradiance model - can be one of 'isotropic' or 'haydavies'.

    dni_extra : numeric, optional
        Extraterrestrial direct normal irradiance. Required when
        ``model='haydavies'``. [W/m2]

    iam_front : numeric, default 1.0
        Incidence angle modifier, the fraction of direct irradiance incident
        on the front surface that is not reflected away. [unitless]

    iam_back : numeric, default 1.0
        Incidence angle modifier, the fraction of direct irradiance incident
        on the back surface that is not reflected away. [unitless]

    bifaciality : numeric, default 0.8
        Ratio of the efficiency of the module's rear surface to the efficiency
        of the front surface. [unitless]

    shade_factor : numeric, default -0.02
        Fraction of back surface irradiance that is blocked by array mounting
        structures. Negative value is a reduction in back irradiance.
        [unitless]

    transmission_factor : numeric, default 0.0
        Fraction of irradiance on the back surface that does not reach the
        module's cells due to module features such as busbars, junction box,
        etc. A negative value is a reduction in back irradiance. [unitless]

    npoints : int, default 100
        Number of discretization points for calculating integrated view
        factors.

    vectorize : bool, default False
        If True, vectorize the view factor calculation across ``surface_tilt``.
        This increases speed with the cost of increased memory usage.

    Returns
    -------
    output : dict or DataFrame
        Output is a DataFrame when input ghi is a Series. See Notes for
        descriptions of content.

    Notes
    -----

    ``output`` includes:

    - ``poa_global`` : total irradiance reaching the module cells from both
      front and back surfaces. [W/m^2]
    - ``poa_front`` : total irradiance reaching the module cells from the front
      surface. [W/m^2]
    - ``poa_back`` : total irradiance reaching the module cells from the back
      surface. [W/m^2]
    - ``poa_front_direct`` : direct irradiance reaching the module cells from
      the front surface. [W/m^2]
    - ``poa_front_diffuse`` : total diffuse irradiance reaching the module
      cells from the front surface. [W/m^2]
    - ``poa_front_sky_diffuse`` : sky diffuse irradiance reaching the module
      cells from the front surface. [W/m^2]
    - ``poa_front_ground_diffuse`` : ground-reflected diffuse irradiance
      reaching the module cells from the front surface. [W/m^2]
    - ``shaded_fraction_front`` : fraction of row slant height from the bottom
      that is shaded from direct irradiance on the front surface by adjacent
      rows. [unitless]
    - ``poa_back_direct`` : direct irradiance reaching the module cells from
      the back surface. [W/m^2]
    - ``poa_back_diffuse`` : total diffuse irradiance reaching the module
      cells from the back surface. [W/m^2]
    - ``poa_back_sky_diffuse`` : sky diffuse irradiance reaching the module
      cells from the back surface. [W/m^2]
    - ``poa_back_ground_diffuse`` : ground-reflected diffuse irradiance
      reaching the module cells from the back surface. [W/m^2]
    - ``shaded_fraction_back`` : fraction of row slant height from the bottom
      that is shaded from direct irradiance on the back surface by adjacent
      rows. [unitless]

    References
    ----------
    .. [1] Mikofski, M., Darawali, R., Hamer, M., Neubert, A., and Newmiller,
       J. "Bifacial Performance Modeling in Large Arrays". 2019 IEEE 46th
       Photovoltaic Specialists Conference (PVSC), 2019, pp. 1282-1287.
       :doi:`10.1109/PVSC40753.2019.8980572`.

    See also
    --------
    get_irradiance_poa
    """
    # backside is rotated and flipped relative to front
    backside_tilt, backside_sysaz = _backside(surface_tilt, surface_azimuth)
    # front side POA irradiance
    irrad_front = get_irradiance_poa(
        surface_tilt=surface_tilt, surface_azimuth=surface_azimuth,
        solar_zenith=solar_zenith, solar_azimuth=solar_azimuth,
        gcr=gcr, height=height, pitch=pitch, ghi=ghi, dhi=dhi, dni=dni,
        albedo=albedo, model=model, dni_extra=dni_extra, iam=iam_front,
        npoints=npoints, vectorize=vectorize)
    # back side POA irradiance
    irrad_back = get_irradiance_poa(
        surface_tilt=backside_tilt, surface_azimuth=backside_sysaz,
        solar_zenith=solar_zenith, solar_azimuth=solar_azimuth,
        gcr=gcr, height=height, pitch=pitch, ghi=ghi, dhi=dhi, dni=dni,
        albedo=albedo, model=model, dni_extra=dni_extra, iam=iam_back,
        npoints=npoints, vectorize=vectorize)

    colmap_front = {
        'poa_global': 'poa_front',
        'poa_direct': 'poa_front_direct',
        'poa_diffuse': 'poa_front_diffuse',
        'poa_sky_diffuse': 'poa_front_sky_diffuse',
        'poa_ground_diffuse': 'poa_front_ground_diffuse',
        'shaded_fraction': 'shaded_fraction_front',
    }
    colmap_back = {
        'poa_global': 'poa_back',
        'poa_direct': 'poa_back_direct',
        'poa_diffuse': 'poa_back_diffuse',
        'poa_sky_diffuse': 'poa_back_sky_diffuse',
        'poa_ground_diffuse': 'poa_back_ground_diffuse',
        'shaded_fraction': 'shaded_fraction_back',
    }

    if isinstance(ghi, pd.Series):
        irrad_front = irrad_front.rename(columns=colmap_front)
        irrad_back = irrad_back.rename(columns=colmap_back)
        output = pd.concat([irrad_front, irrad_back], axis=1)
    else:
        for old_key, new_key in colmap_front.items():
            irrad_front[new_key] = irrad_front.pop(old_key)
        for old_key, new_key in colmap_back.items():
            irrad_back[new_key] = irrad_back.pop(old_key)
        irrad_front.update(irrad_back)
        output = irrad_front

    effects = (1 + shade_factor) * (1 + transmission_factor)
    output['poa_global'] = output['poa_front'] + \
        output['poa_back'] * bifaciality * effects
    return output


def _backside(tilt, surface_azimuth):
    backside_tilt = 180. - tilt
    backside_sysaz = (180. + surface_azimuth) % 360.
    return backside_tilt, backside_sysaz