calculate surface angles using unit normal

Author: cwhanse

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
@@ -210,6 +210,50 @@ 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. 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.
@@ -223,8 +267,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,7 +277,9 @@ 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``.
 
     References
     ----------
@@ -245,16 +290,17 @@ def calc_surface_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0):
     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
+    # 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,
@@ -378,15 +424,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
     -------

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': 0, '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': 0,
               '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': 0,
               '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([0, 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., 0.,  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.)
+    np.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': 90,
                            '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., 180.,   0.],
          [ nan,  nan,  nan,  nan]]),
         columns=['tracker_theta', 'aoi', 'surface_azimuth', 'surface_tilt'])
     assert_frame_equal(out, expected)
@@ -408,7 +474,7 @@ 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, 0, 90], dtype=float)
     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 +511,7 @@ 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, 0, 0, 90, 90]
     out = tracking.calc_surface_orientation(rotations)
     np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
     np.testing.assert_allclose(out['surface_azimuth'], expected_azimuths)
@@ -461,7 +527,7 @@ def test_calc_surface_orientation_special():
     # 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, 0, 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,