ENH: BUG: calculate angles and VF from ground between panels to sky

- calculate ground-sky angles to previous and next rows, assuming height is nonzero - calculate limits on ground where it can see the sky - calculate the view factor as a function of z on the ground to the sky - fix places where it still says degrees, bad, no! - add fixme for pv-sky view factor, still has wrong formula - add tests for angles from point z on the ground to tops of current row, and limits of previous and next rows - add a script to make the plot of ground-sky view factor versus z
pvlib · mikofski · Feb 22, 2022 · May 10, 2019 · May 16, 2019 · May 16, 2019
commit 135e074813ca1eac369eca08e3910323c741e471
diff --git a/pvlib/infinite_sheds.py b/pvlib/infinite_sheds.py
@@ -139,9 +139,9 @@ def _gcr_prime(gcr, height, tilt, pitch):
     #  :    \\                      \\
     #  :.....\\......................\\........ module lower edge
     #  :       \                       \    :
-    #  :        \                      /\   h = height above ground
-    #  :         \                  tilt \  :
-    #  +----------\<---------P--------|-->\---- ground
+    #  :        \                       \   h = height above ground
+    #  :         \                 tilt  \  :
+    #  +----------\<---------P----------->\---- ground
 
     return gcr + height / np.sin(tilt) / pitch
 
@@ -184,13 +184,159 @@ def ground_sky_angles(f_z, gcr, height, tilt, pitch):
     return psi_0, psi_1
 
 
-def ground_diffuse_view_factor():
+def ground_sky_angles_prev(f_z, gcr, height, tilt, pitch):
+    """
+    Angles from point z on ground to top and bottom of previous rows beyond.
+
+    .. math::
+
+        \\tan{\\psi_0} = \\frac{\\sin{\\beta^\\prime}}{\\frac{F_z}
+        {\\text{GCR}^\\prime} + \\cos{\\beta^\\prime}}
+
+        0 < F_z < F_{z0,limit}
+
+        \\tan \\psi_1 = \\frac{h}{\\frac{h}{\\tan\\beta} - z}
+
+    Parameters
+    ----------
+    f_z : numeric
+        fraction of ground from previous to next row
+    gcr : numeric
+        ground coverage ratio
+    height : numeric
+        height of module lower edge above the ground
+    tilt : numeric
+        module tilt in radians, between 0 and 180-degrees
+    pitch : numeric
+        row spacing
+    """
+    gcr_prime = _gcr_prime(gcr, height, tilt, pitch)
+    tilt_prime = np.pi - tilt
+    # angle to bottom of panel looking at sky between rows beyond
+    psi_0 = np.arctan2(
+        np.sin(tilt_prime), (1+f_z)/gcr_prime + np.cos(tilt_prime))
+    # angle to front edge of row beyond
+    z = f_z*pitch
+    # other forms raise division by zero errors
+    # avoid division by zero errors
+    psi_1 = np.arctan2(height, height/np.tan(tilt) - z)
+    return psi_0, psi_1
+
+
+def f_z0_limit(gcr, height, tilt, pitch):
+    """
+    Limit from :math:`z_0` where sky is visible between previous rows.
+
+    .. math::
+        F_{z0,limit} = \\frac{h}{P} \\left(
+        \\frac{1}{\\tan \\beta} + \\frac{1}{\\tan \\psi_t}\\right)
+
+    <
8000
span class='blob-code-inner blob-code-marker ' data-code-marker="+">    Parameters
+    ----------
+    gcr : numeric
+        ground coverage ratio
+    height : numeric
+        height of module lower edge above the ground
+    tilt : numeric
+        module tilt in radians, between 0 and 180-degrees
+    pitch : numeric
+        row spacing
+    """
+    tan_psi_t_x0 = sky_angle_0_tangent(gcr, tilt)
+    # tan_psi_t_x0 = gcr * np.sin(tilt) / (1.0 - gcr * np.cos(tilt))
+    return height/pitch * (1/np.tan(tilt) + 1/tan_psi_t_x0)
+
+
+def ground_sky_angles_next(f_z, gcr, height, tilt, pitch):
+    """
+    Angles from point z on the ground to top and bottom of next row beyond.
+
+    .. math::
+        \\tan \\psi_0 = \\frac{h}{\\frac{h}{\\tan\\beta^\\prime}
+        - \\left(P-z\\right)}
+
+        \\tan{\\psi_1} = \\frac{\\sin{\\beta}}
+        {\\frac{F_z^\\prime}{\\text{GCR}^\\prime} + \\cos{\\beta}}
+    """
+    gcr_prime = _gcr_prime(gcr, height, tilt, pitch)
+    tilt_prime = np.pi - tilt
+    # angle to bottom of panel looking at sky between rows beyond
+    fzprime = 1-f_z
+    zprime = fzprime*pitch
+    # other forms raise division by zero errors
+    # avoid division by zero errors
+    psi_0 = np.arctan2(height, height/np.tan(tilt_prime) - zprime)
+    # angle to front edge of row beyond
+    psi_1 = np.arctan2(np.sin(tilt), (1+fzprime)/gcr_prime + np.cos(tilt))
+    return psi_0, psi_1
+
+
+def f_z1_limit(gcr, height, tilt, pitch):
+    """
+    Limit from :math:`z_1^\\prime` where sky is visible between next rows.
+
+    .. math::
+        F_{z1,limit} = \\frac{h}{P} \\left(
+        \\frac{1}{\\tan \\psi_t} - \\frac{1}{\\tan \\beta}\\right)
+
+    Parameters
+    ----------
+    gcr : numeric
+        ground coverage ratio
+    height : numeric
+        height of module lower edge above the ground
+    tilt : numeric
+        module tilt in radians, between 0 and 180-degrees
+    pitch : numeric
+        row spacing
+    """
+   
8000
 tan_psi_t_x1 = sky_angle_0_tangent(gcr, np.pi-tilt)
+    # tan_psi_t_x1 = gcr * np.sin(pi-tilt) / (1.0 - gcr * np.cos(pi-tilt))
+    return height/pitch * (1/tan_psi_t_x1 - 1/np.tan(tilt))
+
+
+def calc_fx_sky(psi_0, psi_1):
+    return (np.cos(psi_0) + np.cos(psi_1))/2
+
+
+def ground_sky_diffuse_view_factor(gcr, height, tilt, pitch, npoints=100):
     """
     Calculate the fraction of diffuse irradiance from the sky incident on the
     ground.
     :return:
     """
-    pass
+    args = gcr, height, tilt, pitch
+    gcr_prime = _gcr_prime(*args)
+    fz0_limit = f_z0_limit(*args)
+    fz1_limit = f_z1_limit(*args)
+    # divide the space between rows into N points, but include extra
+    # space to account for sky visible from adjacent rows
+    anyx = np.linspace(
+        0.0 if (1-fz1_limit) > 0 else (1-fz1_limit),
+        1.0 if fz0_limit < 1 else fz0_limit,
+        npoints)
+    # calculate the angles psi_0 and psi_1 that subtend the sky visible
+    # from between rows
+    psi_x = ground_sky_angles(anyx, *args)
+    # front edge
+    psi_x0 = ground_sky_angles_prev(anyx, *args)
+    fx0_sky_next = []
+    next_row = 0.0
+    while (fz0_limit - next_row) > 0:
+        fx0_sky_next.append(
+            np.interp(anyx + next_row, anyx, calc_fx_sky(*psi_x0)))
+        next_row += 1.0
+    # back edge
+    psi_x1 = ground_sky_angles_next(anyx, *args)
+    fx1_sky_prev = []
+    prev_row = 0.0
+    while (fz1_limit - prev_row) > 0:
+        fx1_sky_prev.append(
+            np.interp(anyx - prev_row, anyx, calc_fx_sky(*psi_x1)))
+        prev_row += 1.0
+    fx_sky = calc_fx_sky(*psi_x) + np.sum(fx0_sky_next, axis=0) + np.sum(fx1_sky_prev, axis=0)
+    x = np.linspace(0, 1, npoints)
+    return x, np.interp(x, anyx, fx_sky)
 
 
 def diffuse_fraction(ghi, dhi):
@@ -305,7 +451,7 @@ def sky_angle_tangent(gcr, tilt, f_x):
     gcr : numeric
         ratio of module length versus row spacing
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
     f_x : numeric
         fraction of module shaded from bottom
 
@@ -333,7 +479,7 @@ def sky_angle_0_tangent(gcr, tilt):
     gcr : numeric
         ratio of module length to row spacing
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
 
     Returns
     -------
@@ -352,7 +498,7 @@ def f_sky_diffuse_pv(tilt, tan_psi_top, tan_psi_top_0):
     Parameters
     ----------
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
     tan_psi_top : numeric
         tangent of angle from shade line to top of next row
     tan_psi_top_0 : numeric
@@ -385,6 +531,7 @@ def f_sky_diffuse_pv(tilt, tan_psi_top, tan_psi_top_0):
         \\large{F_{sky \\rightarrow no\\ shade} = \\frac{1 + \\frac{1 +
         \\cos \\left(\\psi_t + \\beta \\right)}{1 + \\cos \\beta} }{2}}
     """
+    # FIXME: according to Marion, should be difference of cosines
     psi_top = np.arctan(tan_psi_top)
     psi_top_0 = np.arctan(tan_psi_top_0)
     f_sky_pv_shade = (
@@ -429,7 +576,7 @@ def ground_angle(gcr, tilt, f_x):
     gcr : numeric
         ratio of module length to row spacing
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
     f_x : numeric
         fraction of module shaded from bottom, ``f_x = 0`` if shade line at
         bottom and no shade, ``f_x = 1`` if shade line at top and all shade
@@ -461,7 +608,7 @@ def ground_angle_tangent(gcr, tilt, f_x):
     gcr : numeric
         ratio of module length to row spacing
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
     f_x : numeric
         fraction of module shaded from bottom, ``f_x = 0`` if shade line at
         bottom and no shade, ``f_x = 1`` if shade line at top and all shade
@@ -489,7 +636,7 @@ def ground_angle_1_tangent(gcr, tilt):
     gcr : numeric
         ratio of module length to row spacing
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
 
     Returns
     -------
@@ -507,7 +654,7 @@ def f_ground_pv(tilt, tan_psi_bottom, tan_psi_bottom_1):
     Parameters
     ----------
     tilt : numeric
-        angle of surface normal from vertical in degrees
+        angle of surface normal from vertical in radians
     tan_psi_bottom : numeric
         tangent of angle from shade line to bottom of next row
     tan_psi_bottom_1 : numeric

diff --git a/pvlib/test/test_infinite_sheds.py b/pvlib/test/test_infinite_sheds.py
@@ -20,6 +20,8 @@
 #   system-azimuth: positive counter-clockwise from north
 TILT, SYSAZ = 20.0, 250.0
 GCR = 0.5  # ground coverage ratio
+HEIGHT = 1  # height above ground
+PITCH = 4  # row spacing
 
 # IAM parameters
 B0 = 0.05
@@ -97,6 +99,87 @@ def test_unshaded_ground_fraction():
     assert np.allclose(f_sky_gnd, F_GND_SKY)
 
 
+ARGS = (GCR, HEIGHT, TILT_RAD, PITCH)
+GCR_PRIME = pvlib.infinite_sheds._gcr_prime(*ARGS)
+
+
+def calc_ground_sky_angles_at_edges(tilt_rad=TILT_RAD, gcr_prime=GCR_PRIME):
+    back_tilt_rad = np.pi - tilt_rad
+    psi_0_x0 = back_tilt_rad
+    opposite_side = gcr_prime * np.sin(back_tilt_rad)
+    adjacent_side = 1 - gcr_prime * np.cos(back_tilt_rad)
+    # tan_psi_1_x0 = opposite_side / adjacent_side
+    psi_1_x0 = np.arctan2(opposite_side, adjacent_side)
+    opposite_side = gcr_prime * np.sin(tilt_rad)
+    adjacent_side = 1 - gcr_prime * np.cos(tilt_rad)
+    # tan_psi_0_x1 = opposite_side / adjacent_side
+    psi_0_x1 = np.arctan2(opposite_side, adjacent_side)
+    psi_1_x1 = tilt_rad
+    return psi_0_x0, psi_1_x0, psi_0_x1, psi_1_x1
+
+
+PSI_0_X0, PSI_1_X0, PSI_0_X1, PSI_1_X1 = calc_ground_sky_angles_at_edges()
+
+
+def test_ground_sky_angles():
+    # check limit at x=0, these are the same as the back edge of the row beyond
+    assert np.allclose(
+        pvlib.infinite_sheds.ground_sky_angles(0, *ARGS), (PSI_0_X0, PSI_1_X0))
+    assert np.allclose(
+        pvlib.infinite_sheds.ground_sky_angles(1, *ARGS), (PSI_0_X1, PSI_1_X1))
+
+
+FZ0_LIMIT = pvlib.infinite_sheds.f_z0_limit(*ARGS)
+PSI_TOP = np.arctan2(GCR * np.sin(TILT_RAD), (1.0 - GCR * np.cos(TILT_RAD)))
+
+
+def test_ground_sky_angles_prev():
+    if HEIGHT > 0:
+        # check limit at z=0, these are the same as z=1 of the previous row
+        assert np.allclose(
+            pvlib.infinite_sheds.ground_sky_angles_prev(0, *ARGS),
+            (PSI_0_X1, PSI_1_X1))
+        # check limit at z=z0_limit, angles must sum to 180
+        assert np.isclose(
+            sum(pvlib.infinite_sheds.ground_sky_angles_prev(FZ0_LIMIT, *ARGS)),
+            np.pi)
+        # directly under panel, angle should be 90 straight upward!
+        z_panel = HEIGHT / PITCH / np.tan(TILT_RAD)
+        assert np.isclose(
+            pvlib.infinite_sheds.ground_sky_angles_prev(z_panel, *ARGS)[1],
+            np.pi / 2.)
+    # angles must be the same as psi_top
+    assert np.isclose(
+        pvlib.infinite_sheds.ground_sky_angles_prev(FZ0_LIMIT, *ARGS)[0],
+        PSI_TOP)
+
+
+FZ1_LIMIT = pvlib.infinite_sheds.f_z1_limit(*ARGS)
+PSI_TOP_BACK = np.arctan2(
+    GCR * np.sin(BACK_TILT_RAD), (1.0 - GCR * np.cos(BACK_TILT_RAD)))
+
+
+def test_ground_sky_angles_next():
+    if HEIGHT > 0:
+        # check limit at z=1, these are the same as z=0 of the next row beyond
+        assert np.allclose(
+            pvlib.infinite_sheds.ground_sky_angles_next(1, *ARGS),
+            (PSI_0_X0, PSI_1_X0))
+        # check limit at zprime=z1_limit, angles must sum to 180
+        assert np.isclose(
+            sum(pvlib.infinite_sheds.ground_sky_angles_next(1 - FZ1_LIMIT, *ARGS)),
+            np.pi)
+        # directly under panel, angle should be 90 straight upward!
+        z_panel = 1 + HEIGHT / PITCH / np.tan(TILT_RAD)
+        assert np.isclose(
+            pvlib.infinite_sheds.ground_sky_angles_next(z_panel, *ARGS)[0],
+            np.pi / 2.)
+    # angles must be the same as psi_top
+    assert np.isclose(
+        pvlib.infinite_sheds.ground_sky_angles_next(1 - FZ1_LIMIT, *ARGS)[1],
+        PSI_TOP_BACK)
+
+
 def test_diffuse_fraction():
     df = pvlib.infinite_sheds.diffuse_fraction(GHI, DHI)
     assert np.allclose(df, DF, equal_nan=True)
@@ -156,3 +239,14 @@ def test_sky_angle_0_tangent():
     # backside
     tan_psi_top = pvlib.infinite_sheds.sky_angle_0_tangent(GCR, BACK_TILT_RAD)
     assert np.allclose(tan_psi_top, TAN_PSI_TOP0_B)
+
+
+if __name__ == '__main__':
+    from matplotlib import pyplot as plt
+    plt.ion()
+    plt.plot(*pvlib.infinite_sheds.ground_sky_diffuse_view_factor(*ARGS))
+    plt.title(
+        'combined sky view factor, not including horizon and first/last row')
+    plt.xlabel('fraction of pitch from front to back')
+    plt.ylabel('view factor')
+    plt.grid()