Calculate Z array and account for tracker height in shading (#7)

* Calculate Z array and account for tracker height in shading

* Correct y_0 calculation and implement relative_slope

* Update documentation and layout plot, fix code style

* Separate original PR code into modules

* Delete two_axis_tracker_shading.py

* Update doc

* Add horizon line for tilted ground

Author: AdamRJensen

twoaxistracking/layout.py

@@ -14,7 +14,8 @@ def _rotate_origin(x, y, rotation_deg):
 
 def generate_field_layout(gcr, collector_area, L_min, neighbor_order,
                           aspect_ratio=None, offset=None, rotation=None,
-                          layout_type=None, plot=False):
+                          layout_type=None, slope_azimuth=0,
+                          slope_tilt=0, plot=False):
     """
     Generate a regularly-spaced collector field layout.
 
@@ -51,19 +52,31 @@ def generate_field_layout(gcr, collector_area, L_min, neighbor_order,
         Counterclockwise rotation of the field in degrees. 0 <= rotation < 180
     layout_type: {square, square_rotated, hexagon_e_w, hexagon_n_s}, optional
         Specification of the special layout type (only depend on gcr).
-    plot: bool, default: True
+    slope_azimuth : float, optional
+        Direction of normal to slope on horizontal [degrees]
+    slope_tilt : float, optional
+        Tilt of slope relative to horizontal [degrees]
+    plot: bool, default: False
         Whether to plot the field layout.
 
     Returns
     -------
     X: array of floats
-        x coordinates of neighboring trackers.
+        Distance of neighboring trackers to the reference tracker in the east-
+        west direction. East is positive.
     Y: array of floats
-        y coordinates of neighboring trackers.
+        Distance of neighboring trackers to the reference tracker in the north-
+        south direction. North is positive.
+    Z: array of floats
+        Relative heights of neighboring trackers.
     tracker_distance: array of floats
-        Distances between neighboring trackers and reference tracker.
+        Distances between neighboring trackers and the reference tracker.
     relative_azimuth: array of floats
-        Relative azimuth between neigboring trackers and reference tracker.
+        Relative azimuth of neighboring trackers - measured clockwise from
+        north [degrees].
+    relative_slope: array of floats
+        Slope between neighboring trackers and reference tracker. A positive
+        slope means neighboring collector is higher than reference collector.
 
     References
     ----------
@@ -128,6 +141,11 @@ def generate_field_layout(gcr, collector_area, L_min, neighbor_order,
     X = X * aspect_ratio
     # Apply field rotation
     X, Y = _rotate_origin(X, Y, rotation)
+    # Calculate relative tracker height based on surface slope
+    Z = - X * np.sin(np.deg2rad(slope_azimuth)) * \
+        np.tan(np.deg2rad(slope_tilt)) \
+        - Y * np.cos(np.deg2rad(slope_azimuth)) * \
+        np.tan(np.deg2rad(slope_tilt))
     # Calculate and apply the scaling factor based on GCR
     scaling = np.sqrt(collector_area / (gcr * aspect_ratio))
     X, Y = X*scaling, Y*scaling
@@ -136,9 +154,12 @@ def generate_field_layout(gcr, collector_area, L_min, neighbor_order,
     tracker_distance = np.sqrt(X**2 + Y**2)
     # The relative azimuth is defined clockwise eastwards from north
     relative_azimuth = np.mod(450-np.rad2deg(np.arctan2(Y, X)), 360)
+    # Relative slope of collectors
+    # positive means collector is higher than reference collector
+    relative_slope = -np.cos(np.deg2rad(slope_azimuth - relative_azimuth)) * slope_tilt  # noqa: E501
 
     # Visualize layout
     if plot:
-        plotting._plot_field_layout(X, Y, L_min)
+        plotting._plot_field_layout(X, Y, Z, L_min)
 
-    return X, Y, tracker_distance, relative_azimuth
+    return X, Y, Z, tracker_distance, relative_azimuth, relative_slope

twoaxistracking/plotting.py

@@ -2,30 +2,40 @@
 from matplotlib.collections import EllipseCollection
 from matplotlib.collections import PatchCollection
 from matplotlib.patches import Polygon
+import matplotlib.colors as mcolors
+from matplotlib import cm
 
 
-def _plot_field_layout(X, Y, L_min):
+def _plot_field_layout(X, Y, Z, L_min):
     """Plot field layout."""
-    fig, ax = plt.subplots(figsize=(6, 6))
+    # 0.000001 is added/subtracted for the limits in order for the colormap
+    # to correctly display the middle color when all tracker Z coords are zero
+    norm = mcolors.Normalize(vmin=min(Z)-0.000001, vmax=max(Z)+0.000001)
+    cmap = cm.viridis_r
+    colors = cmap(norm(Z))
+    fig, ax = plt.subplots(figsize=(6, 6), subplot_kw={'aspect': 'equal'})
     # Plot a circle with a diameter equal to L_min
     ax.add_collection(EllipseCollection(widths=L_min, heights=L_min,
                                         angles=0, units='xy',
-                                        facecolors='white',
+                                        facecolors=colors,
                                         edgecolors=("black",),
                                         linewidths=(1,),
                                         offsets=list(zip(X, Y)),
                                         transOffset=ax.transData))
-    # Add a circle for the origin
+    # Similarly, add a circle for the origin
     ax.add_collection(EllipseCollection(widths=L_min, heights=L_min,
                                         angles=0, units='xy',
                                         facecolors='red',
                                         edgecolors=("black",),
                                         linewidths=(1,), offsets=[0, 0],
                                         transOffset=ax.transData))
     plt.axis('equal')
-    ax.grid()
-    ax.set_ylim(min(Y)-L_min, max(Y)+L_min)
-    ax.set_xlim(min(X)-L_min, max(X)+L_min)
+    fig.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax, shrink=0.8,
+                 label='Relative tracker height (vertical)')
+    lower_lim = min(min(X), min(Y)) - L_min
+    upper_lim = max(max(X), max(Y)) + L_min
+    ax.set_ylim(lower_lim, upper_lim)
+    ax.set_xlim(lower_lim, upper_lim)
 
 
 def _plot_shading(collector_geometry, unshaded_geomtry, shade_geometries):

twoaxistracking/shading.py

@@ -12,17 +12,18 @@ def _rotate_origin(x, y, rotation_deg):
     return xx, yy
 
 
-def shaded_fraction(solar_azimuth, solar_elevation,
+def shaded_fraction(solar_elevation, solar_azimuth,
                     collector_geometry, L_min, tracker_distance,
-                    relative_azimuth, plot=False):
+                    relative_azimuth, relative_slope,
+                    slope_azimuth=0, slope_tilt=0, plot=False):
     """Calculate the shaded fraction for any layout of two-axis tracking collectors.
 
     Parameters
     ----------
-    solar_azimuth: float
-        Solar azimuth angle in degrees.
     solar_elevation: float
         Solar elevation angle in degrees.
+    solar_azimuth: float
+        Solar azimuth angle in degrees.
     collector_geometry: Shapely geometry object
         The collector aperture geometry.
     L_min: float
@@ -32,6 +33,13 @@ def shaded_fraction(solar_azimuth, solar_elevation,
         Distances between neighboring trackers and reference tracker.
     relative_azimuth: array of floats
         Relative azimuth between neigboring trackers and reference tracker.
+    relative_slope: array of floats
+        Slope between neighboring trackers and reference tracker. A positive
+        slope means neighboring collector is higher than reference collector.
+    slope_azimuth : float
+        Direction of normal to slope on horizontal [degrees].
+    slope_tilt : float
+        Tilt of slope relative to horizontal [degrees].
     plot: bool, default: True
         Whether to plot the projected shadows.
 
@@ -43,16 +51,21 @@ def shaded_fraction(solar_azimuth, solar_elevation,
     # If the sun is below the horizon, set the shaded fraction to nan
     if solar_elevation < 0:
         return np.nan
+    # Set shading fraction to 1 (fully shaded) if the solar elevation is below
+    # the horizon line caused by the tilted ground
+    elif solar_elevation < - np.cos(np.deg2rad(slope_azimuth-solar_azimuth)) * slope_tilt:
+        return 1
 
     azimuth_difference = solar_azimuth - relative_azimuth
 
     # Create mask to only calculate shading for collectors within +/-90° view
     mask = np.where(np.cos(np.deg2rad(azimuth_difference)) > 0)
 
     xoff = tracker_distance[mask]*np.sin(np.deg2rad(azimuth_difference[mask]))
-    yoff = -tracker_distance[mask]\
-        * np.cos(np.deg2rad(azimuth_difference[mask]))\
-        * np.sin(np.deg2rad(solar_elevation))
+    yoff = - tracker_distance[mask] *\
+        np.cos(np.deg2rad(azimuth_difference[mask])) * \
+        np.sin(np.deg2rad(solar_elevation-relative_slope[mask])) / \
+        np.cos(np.deg2rad(relative_slope[mask]))
 
     # Initialize the unshaded area as the collector area
     unshaded_geomtry = collector_geometry