tracking.py

import numpy as np
import pandas as pd

from pvlib.tools import cosd, sind, tand, acosd, asind
from pvlib import irradiance
from pvlib import shading
from pvlib._deprecation import renamed_kwarg_warning


@renamed_kwarg_warning(
    since='0.13.1',
    old_param_name='apparent_azimuth',
    new_param_name='solar_azimuth')
def singleaxis(apparent_zenith, solar_azimuth,
               axis_tilt=0, axis_azimuth=0, max_angle=90,
               backtrack=True, gcr=2.0/7.0, cross_axis_tilt=0):
    """
    Determine the rotation angle of a single-axis tracker when given particular
    solar zenith and azimuth angles.

    See [1]_ and [2]_ for details about the equations. Backtracking may be
    specified, in which case a ground coverage ratio is required.

    Rotation angle is determined in a right-handed coordinate system. The
    tracker ``axis_azimuth`` defines the positive y-axis, the positive x-axis
    is 90 degrees clockwise from the y-axis and parallel to the Earth's
    surface, and the positive z-axis is normal to both x and y-axes and
    oriented skyward. Rotation angle ``tracker_theta`` is a right-handed
    rotation around the y-axis in the x, y, z coordinate system and indicates
    tracker position relative to horizontal. For example, if tracker
    ``axis_azimuth`` is 180 (oriented south) and ``axis_tilt`` is zero, then a
    ``tracker_theta`` of zero is horizontal, a ``tracker_theta`` of 30 degrees
    is a rotation of 30 degrees towards the west, and a ``tracker_theta`` of
    -90 degrees is a rotation to the vertical plane facing east.

    Parameters
    ----------
    apparent_zenith : float, 1d array, or Series
        Solar apparent zenith angles in decimal degrees.

    solar_azimuth : float, 1d array, or Series
        Solar apparent azimuth angles in decimal degrees.

    axis_tilt : float, default 0
        The tilt of the axis of rotation (i.e, the y-axis defined by
        ``axis_azimuth``) with respect to horizontal (degrees). Positive
        ``axis_tilt`` is *downward* in the direction of ``axis_azimuth``. For
        example, for a tracker with ``axis_azimuth``=180 and ``axis_tilt``=10,
        the north end is higher than the south end of the axis.

    axis_azimuth : float, default 0
        A value denoting the compass direction along which the axis of
        rotation lies. Measured in decimal degrees east of north.

    max_angle : float or tuple, default 90
        A value denoting the maximum rotation angle, in decimal degrees,
        of the one-axis tracker from its horizontal position (horizontal
        if axis_tilt = 0). If a float is provided, it represents the maximum
        rotation angle, and the minimum rotation angle is assumed to be the
        opposite of the maximum angle. If a tuple of (min_angle, max_angle) is
        provided, it represents both the minimum and maximum rotation angles.

        A rotation to ``max_angle`` is a counter-clockwise rotation about the
        y-axis of the tracker coordinate system. For example, for a tracker
        with ``axis_azimuth`` oriented to the south, a rotation to
        ``max_angle`` is towards the west, and a rotation toward ``-max_angle``
        is in the opposite direction, toward the east.

    backtrack : bool, default True
        Controls whether the tracker has the capability to "backtrack"
        to avoid row-to-row shading. False denotes no backtrack
        capability. True denotes backtrack capability.

    gcr : float, default 2.0/7.0
        A value denoting the ground coverage ratio of a tracker system that
        utilizes backtracking; i.e. the ratio between the PV array surface area
        to the total ground area. A tracker system with modules 2 meters wide,
        centered on the tracking axis, with 6 meters between the tracking axes
        has a ``gcr`` of 2/6=0.333. If ``gcr`` is not provided, a ``gcr`` of
        2/7 is default. ``gcr`` must be <=1.

    cross_axis_tilt : float, default 0.0
        The angle, relative to horizontal, of the line formed by the
        intersection between the slope containing the tracker axes and a plane
        perpendicular to the tracker axes. The cross-axis tilt should be
        specified using a right-handed convention. For example, trackers with
        ``axis_azimuth`` of 180 degrees (heading south) will have a negative
        cross-axis tilt if the tracker axes plane slopes down to the east and
        positive cross-axis tilt if the tracker axes plane slopes down to the
        west. Use :func:`~pvlib.tracking.calc_cross_axis_tilt` to calculate
        ``cross_axis_tilt``. [degrees]

    Returns
    -------
    dict or DataFrame with the following columns:
        * `tracker_theta`: The rotation angle of the tracker is a right-handed
          rotation defined by `axis_azimuth`.
          tracker_theta = 0 is horizontal. [degrees]
        * `aoi`: The angle-of-incidence of direct irradiance onto the
          rotated panel surface. [degrees]
        * `surface_tilt`: The angle between the panel surface and the earth
          surface, accounting for panel rotation. [degrees]
        * `surface_azimuth`: The azimuth of the rotated panel, determined by
          projecting the vector normal to the panel's surface to the earth's
          surface. [degrees]

    See also
    --------
    pvlib.tracking.calc_axis_tilt
    pvlib.tracking.calc_cross_axis_tilt
    pvlib.tracking.calc_surface_orientation

    References
    ----------
    .. [1] Anderson, K., and Mikofski, M., "Slope-Aware Backtracking for
       Single-Axis Trackers", Technical Report NREL/TP-5K00-76626, July 2020.
       https://www.nrel.gov/docs/fy20osti/76626.pdf
    .. [2] Lorenzo, E., Narvarte, L., and Muñoz, J. (2011). Tracking and
       back-tracking 19(6), 747–753. :doi:`10.1002/pip.1085`
    """

    # MATLAB to Python conversion by
    # Will Holmgren (@wholmgren), U. Arizona. March, 2015.

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

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

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

    # The ideal tracking angle, omega_ideal, is the rotation to place the sun
    # position vector (xp, yp, zp) in the (x, z) plane, which is normal to
    # the panel and contains the axis of rotation. omega_ideal=0 indicates
    # that the panel is horizontal. Here, our convention is that a clockwise
    # rotation is positive, to view rotation angles in the same frame of
    # reference as azimuth. For example, for a system with tracking
    # axis oriented south, a rotation toward the east is negative, and a
    # rotation to the west is positive. This is a right-handed rotation
    # around the tracker y-axis.
    omega_ideal = shading.projected_solar_zenith_angle(
        axis_tilt=axis_tilt,
        axis_azimuth=axis_azimuth,
        solar_zenith=apparent_zenith,
        solar_azimuth=solar_azimuth,
    )

    # filter for sun above panel horizon
    zen_gt_90 = apparent_zenith > 90
    omega_ideal[zen_gt_90] = np.nan

    # Account for backtracking
    if backtrack:
        # distance between rows in terms of rack lengths relative to cross-axis
        # tilt
        axes_distance = 1/(gcr * cosd(cross_axis_tilt))

        # NOTE: account for rare angles below array, see GH 824
        temp = np.abs(axes_distance * cosd(omega_ideal - cross_axis_tilt))

        # backtrack angle using [1], Eq. 14
        with np.errstate(invalid='ignore'):
            omega_correction = np.degrees(
                -np.sign(omega_ideal)*np.arccos(temp))

        # NOTE: in the middle of the day, arccos(temp) is out of range because
        # there's no row-to-row shade to avoid, & backtracking is unnecessary
        # [1], Eqs. 15-16
        with np.errstate(invalid='ignore'):
            tracker_theta = omega_ideal + np.where(
                temp < 1, omega_correction,
                0)
    else:
        tracker_theta = omega_ideal

    # NOTE: max_angle defined relative to zero-point rotation, not the
    # system-plane normal

    # Determine minimum and maximum rotation angles based on max_angle.
    # If max_angle is a single value, assume min_angle is the negative.
    if np.isscalar(max_angle):
        min_angle = -max_angle
    else:
        min_angle, max_angle = max_angle

    # Clip tracker_theta between the minimum and maximum angles.
    tracker_theta = np.clip(tracker_theta, min_angle, max_angle)

    # Calculate auxiliary angles
    surface = calc_surface_orientation(tracker_theta, axis_tilt, axis_azimuth)
    surface_tilt = surface['surface_tilt']
    surface_azimuth = surface['surface_azimuth']
    aoi = irradiance.aoi(surface_tilt, surface_azimuth,
                         apparent_zenith, solar_azimuth)

    # Bundle DataFrame for return values and filter for sun below horizon.
    out = {'tracker_theta': tracker_theta, 'aoi': aoi,
           'surface_azimuth': surface_azimuth, 'surface_tilt': surface_tilt}
    if index is not None:
        out = pd.DataFrame(out, index=index)
        out[zen_gt_90] = np.nan
    else:
        out = {k: np.where(zen_gt_90, np.nan, v) for k, v in out.items()}

    return out


def _unit_normal(axis_azimuth, axis_tilt, theta):
    """
    Unit normal to rotated tracker surface, in global E-N-Up coordinates,
    given by R*(0, 0, 1)^T, where:

        R = Rz(-axis_azimuth) Rx(-axis_tilt) Ry(theta) *

    Rz is a rotation by -axis_azimuth about the z-axis (axis_azimuth
    is negated to convert from an azimuth angle to a rotation angle). Rx is a
    rotation by -axis_tilt about the x-axis, where axis_tilt is negated
    because pvlib's convention is that the positive y-axis is tilted
    downwards. Ry is a rotation by theta
    about the y-axis. theta is negated so that a negative.

    Parameters
    ----------
    axis_azimuth : scalar
    axis_tilt    : scalar
    theta        : scalar or array-like

    Returns
    -------
    ndarray
        Shape (3,) if theta scalar
        Shape (N,3) if theta has length N
    """

    theta = np.asarray(theta)

    cA, sA = cosd(-axis_azimuth), sind(-axis_azimuth)
    cT, sT = cosd(-axis_tilt), sind(-axis_tilt)

    cTh = cosd(theta)
    sTh = sind(theta)

    x = sA * sT * cTh + cA * sTh
    y = sA * sTh - cA * sT * cTh
    z = cT * cTh

    result = np.column_stack((x, y, z))

    return result


def calc_surface_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0):
    """
    Calculate the surface tilt and azimuth angles for a given tracker rotation.

    Parameters
    ----------
    tracker_theta : numeric
        Tracker rotation angle as a right-handed rotation around
        the axis defined by ``axis_tilt`` and ``axis_azimuth``.  For example,
        with ``axis_tilt=0`` and ``axis_azimuth=180``, ``tracker_theta > 0``
        results in ``surface_azimuth`` to the West while ``tracker_theta < 0``
        results in ``surface_azimuth`` to the East. [degree]
    axis_tilt : float, default 0
        The tilt of the axis of rotation with respect to horizontal. [degree]
    axis_azimuth : float, default 0
        A value denoting the compass direction along which the axis of
        rotation lies. Measured east of north. [degree]

    Returns
    -------
    dict or DataFrame
        Contains keys ``'surface_tilt'`` and ``'surface_azimuth'`` representing
        the module orientation accounting for tracker rotation and axis
        orientation (degree).
        Where ``surface_tilt``=0, ``surface_azimuth`` is set equal to
        ``axis_azimuth``.

    References
    ----------
    .. [1] William F. Marion and Aron P. Dobos, "Rotation Angle for the Optimum
       Tracking of One-Axis Trackers", Technical Report NREL/TP-6A20-58891,
       July 2013. :doi:`10.2172/1089596`
    """
    with np.errstate(invalid='ignore', divide='ignore'):
        surface_tilt = acosd(cosd(tracker_theta) * cosd(axis_tilt))

    # unit normal to rotated tracker surface
    unit_normal = _unit_normal(axis_azimuth, axis_tilt, tracker_theta)

    # deviate from [1] to allow for negative tilt.
    # project unit_normal to x-y plane to calculate azimuth
    surface_azimuth = np.degrees(
        np.arctan2(unit_normal[:, 0], unit_normal[:, 1]))
    # constrain angles to [0, 360)
    surface_azimuth = np.mod(surface_azimuth, 360.0)

    surface_azimuth = np.where(surface_tilt == 0., axis_azimuth,
                               surface_azimuth)

    out = {
        'surface_tilt': surface_tilt,
        'surface_azimuth': surface_azimuth,
    }
    if hasattr(tracker_theta, 'index'):
        out = pd.DataFrame(out)
    return out


def calc_axis_tilt(slope_azimuth, slope_tilt, axis_azimuth):
    """
    Calculate tracker axis tilt in the global reference frame when on a sloped
    plane. Axis tilt is the inclination of the tracker rotation axis with
    respect to horizontal, ranging from 0 degrees (horizontal axis) to 90
    degrees (vertical axis).

    Parameters
    ----------
    slope_azimuth : float
        direction of normal to slope on horizontal [degrees]
    slope_tilt : float
        tilt of normal to slope relative to vertical [degrees]
    axis_azimuth : float
        direction of tracker axes on horizontal [degrees]

    Returns
    -------
    axis_tilt : float
        tilt of tracker [degrees]

    See also
    --------
    pvlib.tracking.singleaxis
    pvlib.tracking.calc_cross_axis_tilt

    Notes
    -----
    See [1]_ for derivation of equations.

    References
    ----------
    .. [1] Kevin Anderson and Mark Mikofski, "Slope-Aware Backtracking for
       Single-Axis Trackers", Technical Report NREL/TP-5K00-76626, July 2020.
       https://www.nrel.gov/docs/fy20osti/76626.pdf
    """
    delta_gamma = axis_azimuth - slope_azimuth
    # equations 18-19
    tan_axis_tilt = cosd(delta_gamma) * tand(slope_tilt)
    return np.degrees(np.arctan(tan_axis_tilt))


def _calc_tracker_norm(ba, bg, dg):
    """
    Calculate tracker normal, v, cross product of tracker axis and unit normal,
    N, to the system slope plane.

    Parameters
    ----------
    ba : float
        axis tilt [degrees]
    bg : float
        ground tilt [degrees]
    dg : float
        delta gamma, difference between axis and ground azimuths [degrees]

    Returns
    -------
    vector : tuple
        vx, vy, vz
    """
    cos_ba = cosd(ba)
    cos_bg = cosd(bg)
    sin_bg = sind(bg)
    sin_dg = sind(dg)
    vx = sin_dg * cos_ba * cos_bg
    vy = sind(ba)*sin_bg + cosd(dg)*cos_ba*cos_bg
    vz = -sin_dg*sin_bg*cos_ba
    return vx, vy, vz


def _calc_beta_c(v, dg, ba):
    """
    Calculate the cross-axis tilt angle.

    Parameters
    ----------
    v : tuple
        tracker normal
    dg : float
        delta gamma, difference between axis and ground azimuths [degrees]
    ba : float
        axis tilt [degrees]

    Returns
    -------
    beta_c : float
        cross-axis tilt angle [radians]
    """
    vnorm = np.sqrt(np.dot(v, v))
    beta_c = np.arcsin(
        ((v[0]*cosd(dg) - v[1]*sind(dg)) * sind(ba) + v[2]*cosd(ba)) / vnorm)
    return beta_c


def calc_cross_axis_tilt(
        slope_azimuth, slope_tilt, axis_azimuth, axis_tilt):
    """
    Calculate the angle, relative to horizontal, of the line formed by the
    intersection between the slope containing the tracker axes and a plane
    perpendicular to the tracker axes.

    Use the cross-axis tilt to avoid row-to-row shade when backtracking on a
    slope not parallel with the axis azimuth. Cross-axis tilt should be
    specified using a right-handed convention. For example, trackers with axis
    azimuth of 180 degrees (heading south) will have a negative cross-axis tilt
    if the tracker axes plane slopes down to the east and positive cross-axis
    tilt if the tracker axes plane slopes down to the west.

    Parameters
    ----------
    slope_azimuth : float
        direction of the normal to the slope containing the tracker axes, when
        projected on the horizontal. [degrees]
    slope_tilt : float
        angle of the slope containing the tracker axes, relative to horizontal.
        [degrees]
    axis_azimuth : float
        direction of tracker axes projected on the horizontal. [degrees]
    axis_tilt : float
        tilt of trackers relative to horizontal. [degree]

    Returns
    -------
    cross_axis_tilt : float
        angle, relative to horizontal, of the line formed by the intersection
        between the slope containing the tracker axes and a plane perpendicular
        to the tracker axes [degrees]

    See also
    --------
    pvlib.tracking.singleaxis
    pvlib.tracking.calc_axis_tilt

    Notes
    -----
    See [1]_ for derivation of equations.

    References
    ----------
    .. [1] Kevin Anderson and Mark Mikofski, "Slope-Aware Backtracking for
       Single-Axis Trackers", Technical Report NREL/TP-5K00-76626, July 2020.
       https://www.nrel.gov/docs/fy20osti/76626.pdf
    """
    # delta-gamma, difference between axis and slope azimuths
    delta_gamma = axis_azimuth - slope_azimuth
    # equation 22
    v = _calc_tracker_norm(axis_tilt, slope_tilt, delta_gamma)
    # equation 26
    beta_c = _calc_beta_c(v, delta_gamma, axis_tilt)
    return np.degrees(beta_c)