Projected zenith convenience function

docs/sphinx/source/reference/effects_on_pv_system_output/shading.rst

@@ -9,4 +9,5 @@ Shading
    shading.ground_angle
    shading.masking_angle
    shading.masking_angle_passias
-   shading.sky_diffuse_passias
+   shading.sky_diffuse_passias
+   shading.projected_solar_zenith_angle

docs/sphinx/source/whatsnew/v0.10.4.rst

@@ -7,6 +7,8 @@ v0.10.4 (Anticipated March, 2024)
 
 Enhancements
 ~~~~~~~~~~~~
+* Added function :py:func:`pvlib.shading.projected_solar_zenith_angle`,
+  a common calculation in shading and tracking. (:issue:`1734`, :pull:`1904`)
 
 
 Bug fixes

pvlib/shading.py

@@ -232,3 +232,76 @@ def sky_diffuse_passias(masking_angle):
        Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
     """
     return 1 - cosd(masking_angle/2)**2
+
+
+def projected_solar_zenith_angle(axis_tilt, axis_azimuth,
+                                 solar_zenith, solar_azimuth):
+    r"""
+    Calculate projected solar zenith angle in degrees.
+
+    This solar zenith angle is projected onto the plane whose normal vector is
+    defined by ``axis_tilt`` and ``axis_azimuth``. The normal vector is in the
+    direction of ``axis_azimuth`` (clockwise from north) and tilted from
+    horizontal by ``axis_tilt``. See Figure 5 in [1]_:
+
+    .. figure:: ../../_images/Anderson_Mikofski_2020_Fig5.jpg
+       :alt: Wire diagram of coordinates systems to obtain the projected angle.
+       :align: center
+       :scale: 75 %
+
+       Fig. 5, [1]_: Solar coordinates projection onto tracker rotation plane.
+
+    Parameters
+    ----------
+    axis_tilt : numeric
+        Axis tilt angle in degrees. From horizontal plane to array plane.
+    axis_azimuth : numeric
+        Axis azimuth angle in degrees.
+        North = 0°; East = 90°; South = 180°; West = 270°
+    solar_zenith : numeric
+        Sun's apparent zenith in degrees.
+    solar_azimuth : numeric
+        Sun's azimuth in degrees.
+
+    Returns
+    -------
+    Projected_solar_zenith : numeric
+        In degrees.
+
+    References
+    ----------
+    .. [1] K. Anderson and M. Mikofski, 'Slope-Aware Backtracking for
+       Single-Axis Trackers', National Renewable Energy Lab. (NREL), Golden,
+       CO (United States);
+       NREL/TP-5K00-76626, Jul. 2020. :doi:`10.2172/1660126`.
+    """
+    # Assume the tracker reference frame is right-handed. Positive y-axis is
+    # oriented along tracking axis; from north, the y-axis is rotated clockwise
+    # by the axis azimuth and tilted from horizontal by the axis tilt. The
+    # positive x-axis is 90 deg clockwise from the y-axis and parallel to
+    # horizontal (e.g., if the y-axis is south, the x-axis is west); the
+    # positive z-axis is normal to the x and y axes, pointed upward.
+
+    # Since elevation = 90 - zenith, sin(90-x) = cos(x) & cos(90-x) = sin(x):
+    # Notation from [1], modified to use zenith instead of elevation
+    # cos(elevation) = sin(zenith) and sin(elevation) = cos(zenith)
+    # Avoid recalculating these values
+    sind_solar_zenith = sind(solar_zenith)
+    cosd_axis_azimuth = cosd(axis_azimuth)
+    sind_axis_azimuth = sind(axis_azimuth)
+    sind_axis_tilt = sind(axis_tilt)
+
+    # Sun's x, y, z coords
+    sx = sind_solar_zenith * sind(solar_azimuth)
+    sy = sind_solar_zenith * cosd(solar_azimuth)
+    sz = cosd(solar_zenith)
+    # Eq. (4); sx', sz' values from sun coordinates projected onto surface
+    sx_prime = sx * cosd_axis_azimuth - sy * sind_axis_azimuth
+    sz_prime = (
+        sx * sind_axis_azimuth * sind_axis_tilt
+        + sy * sind_axis_tilt * cosd_axis_azimuth
+        + sz * cosd(axis_tilt)
+    )
+    # Eq. (5); angle between sun's beam and surface
+    theta_T = np.degrees(np.arctan2(sx_prime, sz_prime))
+    return theta_T

pvlib/tests/test_shading.py

@@ -2,8 +2,11 @@
 import pandas as pd
 
 from pandas.testing import assert_series_equal
+from numpy.testing import assert_allclose
 import pytest
+from datetime import timezone, timedelta
 
+import pvlib
 from pvlib import shading
 
 
@@ -37,7 +40,7 @@ def test__ground_angle_zero_gcr():
 
 @pytest.fixture
 def surface_tilt():
-    idx = pd.date_range('2019-01-01', freq='h', periods=3)
+    idx = pd.date_range("2019-01-01", freq="h", periods=3)
     return pd.Series([0, 20, 90], index=idx)
 
 
@@ -104,3 +107,155 @@ def test_sky_diffuse_passias_scalar(average_masking_angle, shading_loss):
     for angle, loss in zip(average_masking_angle, shading_loss):
         actual_loss = shading.sky_diffuse_passias(angle)
         assert np.isclose(loss, actual_loss)
+
+
+@pytest.fixture
+def true_tracking_angle_and_inputs_NREL():
+    # data from NREL 'Slope-Aware Backtracking for Single-Axis Trackers'
+    # doi.org/10.2172/1660126 ; Accessed on 2023-11-06.
+    tzinfo = timezone(timedelta(hours=-5))
+    axis_tilt_angle = 9.666  # deg
+    axis_azimuth_angle = 195.0  # deg
+    timedata = pd.DataFrame(
+        columns=("Apparent Elevation", "Solar Azimuth", "True-Tracking"),
+        data=(
+            (2.404287, 122.791770, -84.440),
+            (11.263058, 133.288729, -72.604),
+            (18.733558, 145.285552, -59.861),
+            (24.109076, 158.939435, -45.578),
+            (26.810735, 173.931802, -28.764),
+            (26.482495, 189.371536, -8.475),
+            (23.170447, 204.136810, 15.120),
+            (17.296785, 217.446538, 39.562),
+            (9.461862, 229.102218, 61.587),
+            (0.524817, 239.330401, 79.530),
+        ),
+    )
+    timedata.index = pd.date_range(
+        "2019-01-01T08", "2019-01-01T17", freq="1H", tz=tzinfo
+    )
+    timedata["Apparent Zenith"] = 90.0 - timedata["Apparent Elevation"]
+    return (axis_tilt_angle, axis_azimuth_angle, timedata)
+
+
+@pytest.fixture
+def singleaxis_psz_implementation_port_data():
+    # data generated with the PSZ angle implementation in tracking.singleaxis
+    # See GitHub issue #1734 & PR #1904
+    axis_tilt_angle = 12.224
+    axis_azimuth_angle = 187.2
+
+    singleaxis_result = pd.DataFrame(
+        columns=[
+            "Apparent Zenith",
+            "Solar Azimuth",
+            "tracker_theta",
+            "surface_azimuth",
+            "surface_tilt",
+        ],
+        data=[
+            [88.86131915, 116.14911543, -84.67346, 98.330924, 84.794565],
+            [85.67558254, 119.46577753, -80.544188, 99.219659, 80.760477],
+            [82.4784391, 122.90558458, -76.226064, 100.171259, 76.5443],
+            [79.37555806, 126.48822166, -71.79054, 101.184411, 72.217365],
+            [76.40491865, 130.23239671, -67.237442, 102.276947, 67.781439],
+            [73.59273783, 134.15525777, -62.55178, 103.476096, 63.224495],
+            [70.96318968, 138.2715258, -57.713941, 104.819827, 58.53107],
+            [68.54068323, 142.59233032, -52.702658, 106.361922, 53.685798],
+            [66.35031258, 147.12377575, -47.496592, 108.18131, 48.676053],
+            [64.41759166, 151.8653323, -42.07579, 110.39903, 43.495367],
+            [62.76775062, 156.80824414, -36.423404, 113.210504, 38.148938],
+            [61.42469841, 161.9342438, -30.527799, 116.950922, 32.663696],
+            [60.40974474, 167.21493901, -24.385012, 122.236817, 27.108957],
+            [59.74022062, 172.61222482, -18.001341, 130.288224, 21.645102],
+            [59.42818646, 178.07994717, -11.395651, 143.610698, 16.652493],
+            [59.47944177, 183.56677914, -4.600779, 166.390187, 13.048796],
+            [59.89302187, 189.01995634, 2.336615, 198.108, 12.441979],
+            [60.66128258, 194.38926277, 9.358232, 225.094855, 15.351466],
+            [61.77055542, 199.63057627, 16.398369, 241.465486, 20.352345],
+            [63.20224386, 204.70842576, 23.389598, 251.116742, 26.231294],
+            [64.93416116, 209.59729217, 30.268795, 257.259578, 32.425598],
+            [66.94189859, 214.28170196, 36.982274, 261.49605, 38.674352],
+            [69.20004673, 218.75538494, 43.489104, 264.617474, 44.841832],
+            [71.68314725, 223.01986867, 49.762279, 267.042188, 50.852813],
+            [74.36628597, 227.08285659, 55.787916, 269.007999, 56.666604],
+            [77.22520074, 230.95665462, 61.562937, 270.658956, 62.264111],
+            [80.23550305, 234.65680797, 67.091395, 272.086933, 67.639267],
+            [83.3693091, 238.20102038, 72.378024, 273.352342, 72.790188],
+            [86.57992299, 241.60837123, 77.408775, 274.492262, 77.698775],
+            [89.70940444, 244.89880789, 82.045935, 275.505443, 82.227402],
+        ],
+    )
+    singleaxis_result.index = pd.date_range(
+        "2024-01-25 08:40",
+        "2024-01-25 18:20",
+        freq="20min",
+        tz=timezone(timedelta(hours=1)),
+    )
+    return (axis_tilt_angle, axis_azimuth_angle, singleaxis_result)
+
+
+def test_projected_solar_zenith_angle_numeric(
+    true_tracking_angle_and_inputs_NREL,
+    singleaxis_psz_implementation_port_data
+):
+    psz_func = shading.projected_solar_zenith_angle
+    axis_tilt, axis_azimuth, timedata = true_tracking_angle_and_inputs_NREL
+    # test against data provided by NREL
+    psz = psz_func(
+        axis_tilt,
+        axis_azimuth,
+        timedata["Apparent Zenith"],
+        timedata["Solar Azimuth"],
+    )
+    assert_allclose(psz, timedata["True-Tracking"], atol=1e-3)
+    # test by changing axis azimuth and tilt
+    psz = psz_func(
+        -axis_tilt,
+        axis_azimuth - 180,
+        timedata["Apparent Zenith"],
+        timedata["Solar Azimuth"],
+    )
+    assert_allclose(psz, -timedata["True-Tracking"], atol=1e-3)
+
+    # test implementation port from tracking.singleaxis
+    axis_tilt, axis_azimuth, singleaxis = \
+        singleaxis_psz_implementation_port_data
+    psz = pvlib.tracking.singleaxis(
+        singleaxis["Apparent Zenith"],
+        singleaxis["Solar Azimuth"],
+        axis_tilt,
+        axis_azimuth,
+        backtrack=False,
+    )
+    assert_allclose(
+        psz["tracker_theta"],
+        singleaxis["tracker_theta"],
+        atol=1e-6
+    )
+
+
+@pytest.mark.parametrize(
+    "cast_type, cast_func",
+    [
+        (float, lambda x: float(x)),
+        (np.ndarray, lambda x: np.array([x])),
+        (pd.Series, lambda x: pd.Series(data=[x])),
+    ],
+)
+def test_projected_solar_zenith_angle_datatypes(
+    cast_type, cast_func, true_tracking_angle_and_inputs_NREL
+):
+    psz_func = shading.projected_solar_zenith_angle
+    axis_tilt, axis_azimuth, timedata = true_tracking_angle_and_inputs_NREL
+    sun_apparent_zenith = timedata["Apparent Zenith"].iloc[0]
+    sun_azimuth = timedata["Solar Azimuth"].iloc[0]
+
+    axis_tilt, axis_azimuth, sun_apparent_zenith, sun_azimuth = (
+        cast_func(axis_tilt),
+        cast_func(axis_azimuth),
+        cast_func(sun_apparent_zenith),
+        cast_func(sun_azimuth),
+    )
+    psz = psz_func(axis_tilt, axis_azimuth, sun_apparent_zenith, axis_azimuth)
+    assert isinstance(psz, cast_type)

pvlib/tracking.py

@@ -3,6 +3,7 @@
 
 from pvlib.tools import cosd, sind, tand, acosd, asind
 from pvlib import irradiance
+from pvlib import shading
 
 
 def singleaxis(apparent_zenith, apparent_azimuth,
@@ -126,51 +127,20 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     if apparent_azimuth.ndim > 1 or apparent_zenith.ndim > 1:
         raise ValueError('Input dimensions must not exceed 1')
 
-    # Calculate sun position x, y, z using coordinate system as in [1], Eq 1.
-
-    # NOTE: solar elevation = 90 - solar zenith, then use trig identities:
-    # sin(90-x) = cos(x) & cos(90-x) = sin(x)
-    sin_zenith = sind(apparent_zenith)
-    x = sin_zenith * sind(apparent_azimuth)
-    y = sin_zenith * cosd(apparent_azimuth)
-    z = cosd(apparent_zenith)
-
-    # Assume the tracker reference frame is right-handed. Positive y-axis is
-    # oriented along tracking axis; from north, the y-axis is rotated clockwise
-    # by the axis azimuth and tilted from horizontal by the axis tilt. The
-    # positive x-axis is 90 deg clockwise from the y-axis and parallel to
-    # horizontal (e.g., if the y-axis is south, the x-axis is west); the
-    # positive z-axis is normal to the x and y axes, pointed upward.
-
-    # Calculate sun position (xp, yp, zp) in tracker coordinate system using
-    # [1] Eq 4.
-
-    cos_axis_azimuth = cosd(axis_azimuth)
-    sin_axis_azimuth = sind(axis_azimuth)
-    cos_axis_tilt = cosd(axis_tilt)
-    sin_axis_tilt = sind(axis_tilt)
-    xp = x*cos_axis_azimuth - y*sin_axis_azimuth
-    # not necessary to calculate y'
-    # yp = (x*cos_axis_tilt*sin_axis_azimuth
-    #       + y*cos_axis_tilt*cos_axis_azimuth
-    #       - z*sin_axis_tilt)
-    zp = (x*sin_axis_tilt*sin_axis_azimuth
-          + y*sin_axis_tilt*cos_axis_azimuth
-          + z*cos_axis_tilt)
-
     # The ideal tracking angle wid is the rotation to place the sun position
-    # vector (xp, yp, zp) in the (y, z) plane, which is normal to the panel and
+    # vector (xp, yp, zp) in the (x, z) plane, which is normal to the panel and
     # contains the axis of rotation.  wid = 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.
-
-    # Calculate angle from x-y plane to projection of sun vector onto x-z plane
-    # using [1] Eq. 5.
-
-    wid = np.degrees(np.arctan2(xp, zp))
+    wid = shading.projected_solar_zenith_angle(
+        axis_tilt=axis_tilt,
+        axis_azimuth=axis_azimuth,
+        solar_zenith=apparent_zenith,
+        solar_azimuth=apparent_azimuth,
+    )
 
     # filter for sun above panel horizon
     zen_gt_90 = apparent_zenith > 90