Use Hottel's method for bifacial.utils.*_integ functions

pvlib/bifacial/infinite_sheds.py

@@ -182,7 +182,7 @@ def _shaded_fraction(solar_zenith, solar_azimuth, surface_tilt,
 def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
                        solar_azimuth, gcr, height, pitch, ghi, dhi, dni,
                        albedo, model='isotropic', dni_extra=None, iam=1.0,
-                       npoints=100, vectorize=False):
+                       npoints=None, vectorize=None):
     r"""
     Calculate plane-of-array (POA) irradiance on one side of a row of modules.
 
@@ -250,12 +250,19 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
         on the surface that is not reflected away. [unitless]
 
     npoints : int, default 100
-        Number of discretization points for calculating integrated view
-        factors.
+
+        .. deprecated:: v0.11.2
+
+           This parameter has no effect; integrated view factors are now
+           calculated exactly instead of with discretized approximations.
 
     vectorize : bool, default False
-        If True, vectorize the view factor calculation across ``surface_tilt``.
-        This increases speed with the cost of increased memory usage.
+
+        .. deprecated:: v0.11.2
+
+           This parameter has no effect; calculations are now vectorized
+           with no memory usage penality.
+
 
     Returns
     -------
@@ -381,7 +388,7 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                    gcr, height, pitch, ghi, dhi, dni,
                    albedo, model='isotropic', dni_extra=None, iam_front=1.0,
                    iam_back=1.0, bifaciality=0.8, shade_factor=-0.02,
-                   transmission_factor=0, npoints=100, vectorize=False):
+                   transmission_factor=0, npoints=None, vectorize=None):
     """
     Get front and rear irradiance using the infinite sheds model.
 
@@ -472,12 +479,18 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         etc. A negative value is a reduction in back irradiance. [unitless]
 
     npoints : int, default 100
-        Number of discretization points for calculating integrated view
-        factors.
+
+        .. deprecated:: v0.11.2
+
+           This parameter has no effect; integrated view factors are now
+           calculated exactly instead of with discretized approximations.
 
     vectorize : bool, default False
-        If True, vectorize the view factor calculation across ``surface_tilt``.
-        This increases speed with the cost of increased memory usage.
+
+        .. deprecated:: v0.11.2
+
+           This parameter has no effect; calculations are now vectorized
+           with no memory usage penality.
 
     Returns
     -------

pvlib/bifacial/utils.py

@@ -4,7 +4,8 @@
 """
 import numpy as np
 from pvlib.tools import sind, cosd, tand
-from scipy.integrate import trapezoid
+import warnings
+from pvlib._deprecation import pvlibDeprecationWarning
 
 
 def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth):
@@ -173,8 +174,34 @@ def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10):
     return vf
 
 
+def _dist(p1, p2):
+    return ((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)**0.5
+
+
+def _angle(p1, p2):
+    return np.arctan2(p2[1] - p1[1], p2[0] - p1[0])
+
+
+def _obstructed_string_length(p1, p2, ob_left, ob_right):
+    # string length calculations for Hottel's crossed strings method,
+    # considering view obstructions from the left and right.
+    # all inputs are (x, y) points.
+
+    # unobstructed length
+    d = _dist(p1, p2)
+    # obstructed on the left
+    d = np.where(_angle(p1, p2) > _angle(p1, ob_left),
+                 _dist(p1, ob_left) + _dist(ob_left, p2),
+                 d)
+    # obstructed on the right
+    d = np.where(_angle(p1, p2) < _angle(p1, ob_right),
+                 _dist(p1, ob_right) + _dist(ob_right, p2),
+                 d)
+    return d
+
+
 def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10,
-                           npoints=100, vectorize=False):
+                           npoints=None, vectorize=None):
     """
     Integrated view factor to the sky from the ground underneath
     interior rows of the array.
@@ -205,23 +232,50 @@ def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10,
         Integration of view factor over the length between adjacent, interior
         rows.  Shape matches that of ``surface_tilt``. [unitless]
     """
-    # Abuse vf_ground_sky_2d by supplying surface_tilt in place
-    # of a signed rotation. This is OK because
-    # 1) z span the full distance between 2 rows, and
-    # 2) max_rows is set to be large upstream, and
-    # 3) _vf_ground_sky_2d considers [-max_rows, +max_rows]
-    # The VFs to the sky will thus be symmetric around z=0.5
-    z = np.linspace(0, 1, npoints)
-    rotation = np.atleast_1d(surface_tilt)
-    if vectorize:
-        fz_sky = vf_ground_sky_2d(rotation, gcr, z, pitch, height, max_rows)
-    else:
-        fz_sky = np.zeros((npoints, len(rotation)))
-        for k, r in enumerate(rotation):
-            vf = vf_ground_sky_2d(r, gcr, z, pitch, height, max_rows)
-            fz_sky[:, k] = vf[:, 0]  # remove spurious rotation dimension
-    # calculate the integrated view factor for all of the ground between rows
-    return trapezoid(fz_sky, z, axis=0)
+    if npoints is not None or vectorize is not None:
+        msg = (
+            "The `npoints` and `vectorize` parameters have no effect and will "
+            "be removed in a future version."  # TODO make this better
+        )
+        warnings.warn(msg, pvlibDeprecationWarning)
+
+    input_is_scalar = np.isscalar(surface_tilt)
+
+    collector_width = pitch * gcr
+    surface_tilt = np.atleast_2d(np.abs(surface_tilt))
+
+    # TODO seems like this should be np.arange(-max_rows, max_rows+1)?
+    # see GH #1867
+    k = np.arange(-max_rows, max_rows)[:, np.newaxis]
+
+    # primary crossed string points:
+    # a, b: boundaries of ground segment
+    # c, d: upper module edges
+    a = (0, 0)
+    b = (pitch, 0)
+    c = ((k+1)*pitch - 0.5 * collector_width * cosd(surface_tilt),
+         height + 0.5 * collector_width * sind(surface_tilt))
+    d = (c[0] - pitch, c[1])
+
+    # view obstruction points (lower module edges)
+    obs_left = (d[0] + collector_width * cosd(surface_tilt),
+                d[1] - collector_width * sind(surface_tilt))
+    obs_right = (obs_left[0] + pitch, obs_left[1])
+
+    # hottel string lengths, considering obstructions
+    ac = _obstructed_string_length(a, c, obs_left, obs_right)
+    ad = _obstructed_string_length(a, d, obs_left, obs_right)
+    bc = _obstructed_string_length(b, c, obs_left, obs_right)
+    bd = _obstructed_string_length(b, d, obs_left, obs_right)
+
+    # crossed string formula for VF
+    vf_per_slat = np.maximum(0.5 * (1/pitch) * ((ac + bd) - (bc + ad)), 0)
+    vf_total = np.sum(vf_per_slat, axis=0)
+
+    if input_is_scalar:
+        vf_total = vf_total.item()
+
+    return vf_total
 
 
 def _vf_poly(surface_tilt, gcr, x, delta):

pvlib/tests/bifacial/test_infinite_sheds.py

@@ -7,6 +7,8 @@
 from pvlib.bifacial import infinite_sheds
 from ..conftest import assert_series_equal
 
+from pvlib._deprecation import pvlibDeprecationWarning
+
 import pytest
 
 
@@ -92,11 +94,10 @@ def test_get_irradiance_poa():
     dni = 700
     albedo = 0
     iam = 1.0
-    npoints = 100
     res = infinite_sheds.get_irradiance_poa(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni,
-        albedo, iam=iam, npoints=npoints)
+        albedo, iam=iam)
     expected_diffuse = np.array([300.])
     expected_direct = np.array([700.])
     expected_global = expected_diffuse + expected_direct
@@ -122,7 +123,7 @@ def test_get_irradiance_poa():
     res = infinite_sheds.get_irradiance_poa(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni,
-        albedo, iam=iam, npoints=npoints)
+        albedo, iam=iam)
     assert np.allclose(res['poa_global'], expected_global)
     assert np.allclose(res['poa_diffuse'], expected_diffuse)
     assert np.allclose(res['poa_direct'], expected_direct)
@@ -144,7 +145,7 @@ def test_get_irradiance_poa():
     res = infinite_sheds.get_irradiance_poa(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni,
-        albedo, iam=iam, npoints=npoints)
+        albedo, iam=iam)
     assert isinstance(res, pd.DataFrame)
     assert_series_equal(res['poa_global'], expected_global)
     assert_series_equal(res['shaded_fraction'], expected_shaded_fraction)
@@ -161,8 +162,7 @@ def test__backside_tilt():
     assert np.allclose(back_az, np.array([0., 330., 90., 180.]))
 
 
-@pytest.mark.parametrize("vectorize", [True, False])
-def test_get_irradiance(vectorize):
+def test_get_irradiance():
     # singleton inputs
     solar_zenith = 0.
     solar_azimuth = 180.
@@ -177,12 +177,10 @@ def test_get_irradiance(vectorize):
     albedo = 0.
     iam_front = 1.0
     iam_back = 1.0
-    npoints = 100
     result = infinite_sheds.get_irradiance(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back,
-        bifaciality=0.8, shade_factor=-0.02, transmission_factor=0,
-        npoints=npoints, vectorize=vectorize)
+        bifaciality=0.8, shade_factor=-0.02, transmission_factor=0)
     expected_front_diffuse = np.array([300.])
     expected_front_direct = np.array([700.])
     expected_front_global = expected_front_diffuse + expected_front_direct
@@ -205,12 +203,11 @@ def test_get_irradiance(vectorize):
     result = infinite_sheds.get_irradiance(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back,
-        bifaciality=0.8, shade_factor=-0.02, transmission_factor=0,
-        npoints=npoints, vectorize=vectorize)
+        bifaciality=0.8, shade_factor=-0.02, transmission_factor=0)
     result_front = infinite_sheds.get_irradiance_poa(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni,
-        albedo, iam=iam_front, vectorize=vectorize)
+        albedo, iam=iam_front)
     assert isinstance(result, pd.DataFrame)
     expected_poa_global = pd.Series(
         [1000., 500., result_front['poa_global'][2] * (1 + 0.8 * 0.98),
@@ -223,6 +220,30 @@ def test_get_irradiance(vectorize):
                         expected_shaded_fraction)
 
 
+def test_get_irradiance_deprecated():
+    kwargs = {"surface_tilt": 0, "surface_azimuth": 0, "solar_zenith": 0,
+              "solar_azimuth": 0, "gcr": 0.5, "height": 1, "pitch": 1,
+              "ghi": 1000, "dhi": 200, "dni": 800, "albedo": 0.2}
+
+    with pytest.warns(pvlibDeprecationWarning, match="have no effect"):
+        _ = infinite_sheds.get_irradiance(**kwargs, npoints=10)
+
+    with pytest.warns(pvlibDeprecationWarning, match="have no effect"):
+        _ = infinite_sheds.get_irradiance(**kwargs, vectorize=True)
+
+
+def test_get_irradiance_poa_deprecated():
+    kwargs = {"surface_tilt": 0, "surface_azimuth": 0, "solar_zenith": 0,
+              "solar_azimuth": 0, "gcr": 0.5, "height": 1, "pitch": 1,
+              "ghi": 1000, "dhi": 200, "dni": 800, "albedo": 0.2}
+
+    with pytest.warns(pvlibDeprecationWarning, match="have no effect"):
+        _ = infinite_sheds.get_irradiance_poa(**kwargs, npoints=10)
+
+    with pytest.warns(pvlibDeprecationWarning, match="have no effect"):
+        _ = infinite_sheds.get_irradiance_poa(**kwargs, vectorize=True)
+
+
 def test_get_irradiance_limiting_gcr():
     # test confirms that irradiance on widely spaced rows is approximately
     # the same as for a single row array
@@ -239,12 +260,10 @@ def test_get_irradiance_limiting_gcr():
     albedo = 1.
     iam_front = 1.0
     iam_back = 1.0
-    npoints = 100
     result = infinite_sheds.get_irradiance(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back,
-        bifaciality=1., shade_factor=-0.00, transmission_factor=0.,
-        npoints=npoints)
+        bifaciality=1., shade_factor=-0.00, transmission_factor=0.)
     expected_ground_diffuse = np.array([500.])
     expected_sky_diffuse = np.array([150.])
     expected_direct = np.array([0.])
@@ -289,12 +308,11 @@ def test_get_irradiance_with_haydavies():
     model = 'haydavies'
     iam_front = 1.0
     iam_back = 1.0
-    npoints = 100
     result = infinite_sheds.get_irradiance(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni, albedo, model, dni_extra,
         iam_front, iam_back, bifaciality=0.8, shade_factor=-0.02,
-        transmission_factor=0, npoints=npoints)
+        transmission_factor=0)
     expected_front_diffuse = np.array([151.38])
     expected_front_direct = np.array([848.62])
     expected_front_global = expected_front_diffuse + expected_front_direct
@@ -314,4 +332,4 @@ def test_get_irradiance_with_haydavies():
             surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
             gcr, height, pitch, ghi, dhi, dni, albedo, model, None,
             iam_front, iam_back, bifaciality=0.8, shade_factor=-0.02,
-            transmission_factor=0, npoints=npoints)
+            transmission_factor=0)

pvlib/tests/bifacial/test_utils.py

@@ -8,6 +8,8 @@
 from pvlib.tools import cosd
 from scipy.integrate import trapezoid
 
+from pvlib._deprecation import pvlibDeprecationWarning
+
 
 @pytest.fixture
 def test_system_fixed_tilt():
@@ -91,17 +93,23 @@ def test__vf_ground_sky_2d(test_system_fixed_tilt):
     assert np.isclose(vf, vfs_gnd_sky[0])
 
 
-@pytest.mark.parametrize("vectorize", [True, False])
-def test_vf_ground_sky_2d_integ(test_system_fixed_tilt, vectorize):
-    ts, pts, vfs_gnd_sky = test_system_fixed_tilt
-    # pass rotation here since max_rows=1 for the hand-solved case in
-    # the fixture test_system, which means the ground-to-sky view factor
-    # isn't summed over enough rows for symmetry to hold.
-    vf_integ = utils.vf_ground_sky_2d_integ(
-        ts['rotation'], ts['gcr'], ts['height'], ts['pitch'],
-        max_rows=1, npoints=3, vectorize=vectorize)
-    expected_vf_integ = trapezoid(vfs_gnd_sky, pts, axis=0)
-    assert np.isclose(vf_integ, expected_vf_integ, rtol=0.1)
+def test_vf_ground_sky_2d_integ():
+    # test against numerical integration with vf_ground_sky_2d
+    x = np.linspace(0, 1, num=1000)
+    kwargs = dict(gcr=0.4, pitch=5.0, height=1.5)
+    vf_x = utils.vf_ground_sky_2d(rotation=-60, x=x, **kwargs)
+    vf_expected = trapezoid(vf_x, x, axis=0)
+
+    vf_actual = utils.vf_ground_sky_2d_integ(surface_tilt=60, **kwargs)
+    assert np.isclose(vf_expected, vf_actual)
+
+
+def test_vf_ground_sky_2d_integ_deprecated():
+    with pytest.warns(pvlibDeprecationWarning, match="have no effect"):
+        _ = utils.vf_ground_sky_2d_integ(0, 0.5, 1, 1, npoints=10)
+
+    with pytest.warns(pvlibDeprecationWarning, match="have no effect"):
+        _ = utils.vf_ground_sky_2d_integ(0, 0.5, 1, 1, vectorize=True)
 
 
 def test_vf_row_sky_2d(test_system_fixed_tilt):