calculate surface angles using unit normal (#2702)

cwhanse · kandersolar · web-flow · commit 36f20768b491 · 2026-04-20T09:59:21.000-04:00
* calculate surface angles using unit normal

* formatting

* more formatting

* comments from review

* change axis_azimuth default when surface_tilt=0

* subtract 90, not add 90

* clarify returned shape, fix a missing assert

* whatsnew

* line length

* Update docs/sphinx/source/whatsnew/v0.15.1.rst

Co-authored-by: Kevin Anderson &lt;kevin.anderso@gmail.com&gt;

---------

Co-authored-by: Kevin Anderson &lt;kevin.anderso@gmail.com&gt;
diff --git a/docs/sphinx/source/whatsnew/v0.15.1.rst b/docs/sphinx/source/whatsnew/v0.15.1.rst
@@ -38,6 +38,9 @@ Enhancements
 * Accelerate the internals of :py:func:`~pvlib.solarpostion.ephemeris`. (:pull:`2626`)
 * Accelerate the intervals of :py:func:`~pvlib.pvsystem.singlediode` when
   `method='lambertw'`. (:pull:`2732`, :pull:`2723`)
+* :py:func:`~pvlib.tracking.singleaxis` accepts negative values for 
+  parameter `axis_tilt`. (:pull:`2702`, :issue:`1976`)
+
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/pvlib/tracking.py b/pvlib/tracking.py
@@ -43,8 +43,10 @@ def singleaxis(apparent_zenith, solar_azimuth,
 
     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.
-        ``axis_tilt`` must be >= 0 and <= 90. [degrees]
+        ``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
@@ -62,9 +64,7 @@ def singleaxis(apparent_zenith, solar_azimuth,
         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. Hence, a ``max_angle``
-        of 180 degrees (equivalent to max_angle = (-180, 180)) allows the
-        tracker to achieve its full rotation capability.
+        is in the opposite direction, toward the east.
 
     backtrack : bool, default True
         Controls whether the tracker has the capability to "backtrack"
@@ -84,7 +84,7 @@ def singleaxis(apparent_zenith, solar_azimuth,
         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
+        ``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
@@ -100,9 +100,10 @@ def singleaxis(apparent_zenith, solar_azimuth,
           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]
+        * `surface_azimuth`: The azimuth of the rotated panel (degrees),
+          determined by projecting the vector normal to the panel's surface to
+          the earth's surface. Where ``surface_tilt``=0, ``surface_azimuth``
+          is set equal to ``axis_azimuth`` - 90.
 
     See also
     --------
@@ -210,6 +211,48 @@ def singleaxis(apparent_zenith, solar_azimuth,
     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.
+
+    Parameters
+    ----------
+    axis_azimuth : scalar
+    axis_tilt    : scalar
+    theta        : scalar or array-like
+
+    Returns
+    -------
+    ndarray
+        Shape (N,3) where 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.
@@ -223,8 +266,7 @@ def calc_surface_orientation(tracker_theta, axis_tilt=0, axis_azimuth=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.
-        ``axis_tilt`` must be >= 0 and <= 90.  [degree]
+        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]
@@ -234,27 +276,32 @@ def calc_surface_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0):
     dict or DataFrame
         Contains keys ``'surface_tilt'`` and ``'surface_azimuth'`` representing
         the module orientation accounting for tracker rotation and axis
-        orientation. [degree]
+        orientation (degree).
+        Where ``surface_tilt``=0, ``surface_azimuth`` is set equal to
+        ``axis_azimuth`` - 90.
 
     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`
     """
+    # from [1], Eq. 1
     with np.errstate(invalid='ignore', divide='ignore'):
         surface_tilt = acosd(cosd(tracker_theta) * cosd(axis_tilt))
 
-        # clip(..., -1, +1) to prevent arcsin(1 + epsilon) issues:
-        azimuth_delta = asind(np.clip(sind(tracker_theta) / sind(surface_tilt),
-                                      a_min=-1, a_max=1))
-        # Combine Eqs 2, 3, and 4:
-        azimuth_delta = np.where(abs(tracker_theta) < 90,
-                                 azimuth_delta,
-                                 -azimuth_delta + np.sign(tracker_theta) * 180)
-        # handle surface_tilt=0 case:
-        azimuth_delta = np.where(sind(surface_tilt) != 0, azimuth_delta, 90)
-        surface_azimuth = (axis_azimuth + azimuth_delta) % 360
+    # for surface azimuth deviate from [1] to allow for negative tilt.
+    # unit normal to rotated tracker surface
+    unit_normal = _unit_normal(axis_azimuth, axis_tilt, tracker_theta)
+
+    # project unit_normal to x-y plane to calculate azimuth
+    surface_azimuth = np.degrees(
+        np.arctan2(unit_normal[:, 0], unit_normal[:, 1]))
+
+    surface_azimuth = np.where(surface_tilt == 0., axis_azimuth - 90.,
+                               surface_azimuth)
+    # constrain angles to [0, 360)
+    surface_azimuth = np.mod(surface_azimuth, 360.0)
 
     out = {
         'surface_tilt': surface_tilt,
@@ -378,15 +425,14 @@ def calc_cross_axis_tilt(
     ----------
     slope_azimuth : float
         direction of the normal to the slope containing the tracker axes, when
-        projected on the horizontal [degrees]
+        projected on the horizontal. [degrees]
     slope_tilt : float
-        angle of the slope containing the tracker axes, relative to horizontal
+        angle of the slope containing the tracker axes, relative to horizontal.
         [degrees]
     axis_azimuth : float
-        direction of tracker axes projected on the horizontal [degrees]
+        direction of tracker axes projected on the horizontal. [degrees]
     axis_tilt : float
-        tilt of trackers relative to horizontal.  ``axis_tilt`` must be >= 0
-        and <= 90. [degree]
+        tilt of trackers relative to horizontal. [degree]
 
     Returns
     -------
diff --git a/tests/test_tracking.py b/tests/test_tracking.py
@@ -23,7 +23,7 @@ def test_solar_noon():
                                        gcr=2.0/7.0)
 
     expect = pd.DataFrame({'tracker_theta': 0, 'aoi': 10,
-                           'surface_azimuth': 90, 'surface_tilt': 0},
+                           'surface_azimuth': 270, 'surface_tilt': 0},
                           index=index, dtype=np.float64)
     expect = expect[SINGLEAXIS_COL_ORDER]
 
@@ -38,7 +38,7 @@ def test_scalars():
                                        max_angle=90, backtrack=True,
                                        gcr=2.0/7.0)
     assert isinstance(tracker_data, dict)
-    expect = {'tracker_theta': 0, 'aoi': 10, 'surface_azimuth': 90,
+    expect = {'tracker_theta': 0, 'aoi': 10, 'surface_azimuth': 270,
               'surface_tilt': 0}
     for k, v in expect.items():
         assert np.isclose(tracker_data[k], v)
@@ -52,7 +52,7 @@ def test_arrays():
                                        max_angle=90, backtrack=True,
                                        gcr=2.0/7.0)
     assert isinstance(tracker_data, dict)
-    expect = {'tracker_theta': 0, 'aoi': 10, 'surface_azimuth': 90,
+    expect = {'tracker_theta': 0, 'aoi': 10, 'surface_azimuth': 270,
               'surface_tilt': 0}
     for k, v in expect.items():
         assert_allclose(tracker_data[k], v, atol=1e-7)
@@ -68,7 +68,7 @@ def test_nans():
                                            gcr=2.0/7.0)
     expect = {'tracker_theta': np.array([0, nan, nan]),
               'aoi': np.array([10, nan, nan]),
-              'surface_azimuth': np.array([90, nan, nan]),
+              'surface_azimuth': np.array([270, nan, nan]),
               'surface_tilt': np.array([0, nan, nan])}
     for k, v in expect.items():
         assert_allclose(tracker_data[k], v, atol=1e-7)
@@ -82,7 +82,7 @@ def test_nans():
                                            max_angle=90, backtrack=True,
                                            gcr=2.0/7.0)
     expect = pd.DataFrame(np.array(
-        [[ 0., 10., 90.,  0.],
+        [[ 0., 10., 270.,  0.],
          [nan, nan, nan, nan],
          [nan, nan, nan, nan]]),
         columns=['tracker_theta', 'aoi', 'surface_azimuth', 'surface_tilt'])
@@ -195,6 +195,54 @@ def test_backtrack():
     assert_frame_equal(expect, tracker_data)
 
 
+def test__unit_normal():
+    # with scalar input
+    unorm = tracking._unit_normal(180., 45., 45.)
+    assert_allclose(unorm, np.array([[-np.sqrt(2)/2, -0.5, 0.5]]))
+    # with vector input
+    az = np.array([0., 90., 180., 270.,
+                   0., 90., 180., 270.,
+                   180., 180., 180, 180.,
+                   180., 180., 180., 180,
+                   0., 90., 180., 270.,
+                   ])
+    tilt = np.array([30., 30., 30., 30.,
+                     0., 0., 0., 0.,
+                     -30., -90., 90., 180.,
+                     0., 0., 0., 0.,
+                     30., 30., 30., 30,
+                     ])
+    theta = np.array([0., 0., 0., 0.,
+                      0., 0., 0., 0.,
+                      0., 0., 0., 0.,
+                      -30., 30., -90., 90.,
+                      30., 30., 30., 30.,
+                      ])
+    expected = np.array(
+        [[ 0., 0.5, 0.8660254],
+         [ 0.5, 0., 0.8660254],
+         [ 0., -0.5, 0.8660254],
+         [-0.5, -0., 0.8660254],
+         [ 0., 0., 1.],
+         [ 0., 0., 1.],
+         [ 0., -0., 1.],
+         [-0., 0., 1.],
+         [-0., 0.5, 0.8660254],
+         [-0., 1., 0.],
+         [ 0., -1., 0.],
+         [ 0., -0., -1.],
+         [ 0.5, 0., 0.8660254],
+         [-0.5, -0.,  0.8660254],
+         [ 1., 0., 0.],
+         [-1., -0., 0.],
+         [ 0.5, 0.4330127, 0.75],
+         [ 0.4330127, -0.5, 0.75],
+         [-0.5, -0.4330127, 0.75],
+         [-0.4330127, 0.5, 0.75]])
+    unorms = tracking._unit_normal(az, tilt, theta)
+    assert np.allclose(unorms, expected)
+
+
 def test_axis_tilt():
     apparent_zenith = pd.Series([30])
     apparent_azimuth = pd.Series([135])
@@ -226,6 +274,7 @@ def test_axis_tilt():
 
 
 def test_axis_azimuth():
+    # sun to the east, horizontal east-oriented tracker
     apparent_zenith = pd.Series([30])
     apparent_azimuth = pd.Series([90])
 
@@ -234,13 +283,30 @@ def test_axis_azimuth():
                                        max_angle=90, backtrack=True,
                                        gcr=2.0/7.0)
 
-    expect = pd.DataFrame({'aoi': 30, 'surface_azimuth': 180,
+    expect = pd.DataFrame({'aoi': 30, 'surface_azimuth': 0,
                            'surface_tilt': 0, 'tracker_theta': 0},
                           index=[0], dtype=np.float64)
     expect = expect[SINGLEAXIS_COL_ORDER]
 
     assert_frame_equal(expect, tracker_data)
 
+    # sun to the east, horizontal south-oriented tracker
+    apparent_zenith = pd.Series([30])
+    apparent_azimuth = pd.Series([90])
+
+    tracker_data = tracking.singleaxis(apparent_zenith, apparent_azimuth,
+                                       axis_tilt=0, axis_azimuth=180,
+                                       max_angle=90, backtrack=True,
+                                       gcr=2.0/7.0)
+
+    expect = pd.DataFrame({'aoi': 0, 'surface_azimuth': 90,
+                           'surface_tilt': 30, 'tracker_theta': -30},
+                          index=[0], dtype=np.float64)
+    expect = expect[SINGLEAXIS_COL_ORDER]
+
+    assert_frame_equal(expect, tracker_data)
+
+    # sun to the south, horizontal east-oriented tracker
     apparent_zenith = pd.Series([30])
     apparent_azimuth = pd.Series([180])
 
@@ -269,7 +335,7 @@ def test_horizon_flat():
                               axis_azimuth=180, backtrack=False, max_angle=180)
     expected = pd.DataFrame(np.array(
         [[ nan,  nan,  nan,  nan],
-         [  0.,  45., 270.,   0.],
+         [  0.,  45., 90.,   0.],
          [ nan,  nan,  nan,  nan]]),
         columns=['tracker_theta', 'aoi', 'surface_azimuth', 'surface_tilt'])
     assert_frame_equal(out, expected)
@@ -389,6 +455,23 @@ def test_slope_aware_backtracking():
                         check_less_precise=True)
 
 
+def test_singleaxis_neg_axis_tilt():
+    ''' Check equivalence of (negative tilt, axis azimuth) and
+    (positive tilt, axis azimuth + 180)
+    '''
+    params = dict(apparent_zenith=45, solar_azimuth=270)
+
+    tr_pos = pvlib.tracking.singleaxis(axis_tilt=10, axis_azimuth=0,
+                                       **params)
+    tr_neg = pvlib.tracking.singleaxis(axis_tilt=-10, axis_azimuth=180,
+                                       **params)
+
+    tr_neg['tracker_theta'] *= -1  # expect tracker_theta to be negated
+
+    for key in tr_pos:
+        assert_allclose(tr_pos[key], tr_neg[key])
+
+
 def test_singleaxis_aoi_gh1221():
     # vertical tracker
     loc = pvlib.location.Location(40.1134, -88.3695)
@@ -408,7 +491,8 @@ def test_calc_surface_orientation_types():
     # numpy arrays
     rotations = np.array([-10, 0, 10])
     expected_tilts = np.array([10, 0, 10], dtype=float)
-    expected_azimuths = np.array([270, 90, 90], dtype=float)
+    expected_azimuths = np.array([270, 270, 90], dtype=float)
+    # defaults to axis_azimuth=0
     out = tracking.calc_surface_orientation(tracker_theta=rotations)
     np.testing.assert_allclose(expected_tilts, out['surface_tilt'])
     np.testing.assert_allclose(expected_azimuths, out['surface_azimuth'])
@@ -445,7 +529,8 @@ def test_calc_surface_orientation_special():
     # special cases for rotations
     rotations = np.array([-180, -90, -0, 0, 90, 180])
     expected_tilts = np.array([180, 90, 0, 0, 90, 180], dtype=float)
-    expected_azimuths = [270, 270, 90, 90, 90, 90]
+    expected_azimuths = [270, 270, 270, 270, 90, 90]
+    # defaults to axis_azimuth=0
     out = tracking.calc_surface_orientation(rotations)
     np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
     np.testing.assert_allclose(out['surface_azimuth'], expected_azimuths)
@@ -454,14 +539,15 @@ def test_calc_surface_orientation_special():
     rotations = np.array([-10, 0, 10])
     expected_tilts = np.array([90, 90, 90], dtype=float)
     expected_azimuths = np.array([350, 0, 10], dtype=float)
+    # defaults to axis_azimuth=0
     out = tracking.calc_surface_orientation(rotations, axis_tilt=90)
     np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
     np.testing.assert_allclose(out['surface_azimuth'], expected_azimuths)
 
     # special cases for axis_azimuth
     rotations = np.array([-10, 0, 10])
     expected_tilts = np.array([10, 0, 10], dtype=float)
-    expected_azimuth_offsets = np.array([-90, 90, 90], dtype=float)
+    expected_azimuth_offsets = np.array([-90, -90, 90], dtype=float)
     for axis_azimuth in [0, 90, 180, 270, 360]:
         expected_azimuths = (axis_azimuth + expected_azimuth_offsets) % 360
         out = tracking.calc_surface_orientation(rotations,
