Extend infinite_sheds to optionally use haydavies transposition

docs/sphinx/source/user_guide/bifacial.rst

@@ -91,8 +91,16 @@ considered to be towards the "front" of the array. Array height differs in this
 code from the description in [1], where array height is described at the row's
 lower edge.
 
+If ``model='isotropic'`` (the default), ``dhi`` is assumed to be isotropically
+distributed across the sky dome as in [1]_. This implementation provides an
+optional extension to [1]_ to model sky anisotropy: if ``model='haydavies'``,
+the input ``dhi`` is decomposed into circumsolar and isotropic components using
+:py:func:`~pvlib.irradiance.haydavies`, with the circumsolar component treated
+as additional ``dni`` for transposition and shading purposes.
+
 This model is influenced by the 2D model published by Marion, *et al.* in [2].
 
+
 References
 ----------
 .. [1] Mikofski, M., Darawali, R., Hamer, M., Neubert, A., and Newmiller,

docs/sphinx/source/whatsnew/v0.9.5.rst

@@ -19,6 +19,12 @@ Enhancements
 * Added the optional `string_factor` parameter to
   :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
 
+* Added an optional ``model`` parameter to
+  :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
+  :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa`
+  to enable use of the hay-davies sky diffuse irradiance model
+  instead of the default isotropic model. (:pull:`1668`)
+
 Bug fixes
 ~~~~~~~~~
 * Added a limit to :py:func:`pvlib.snow.loss_townsend` to guard against
@@ -69,3 +75,5 @@ Contributors
 * Anton Driesse (:ghuser:`adriesse`)
 * Adam R. Jensen (:ghuser:`AdamRJensen`)
 * Michael Deceglie (:ghuser:`mdeceglie`)
+* Saurabh Aneja (:ghuser:`spaneja`)
+* John Moseley (:ghuser:`johnMoseleyArray`)

pvlib/bifacial/infinite_sheds.py

@@ -7,8 +7,7 @@
 from pvlib.tools import cosd, sind, tand
 from pvlib.bifacial import utils
 from pvlib.shading import masking_angle
-from pvlib.irradiance import beam_component, aoi
-
+from pvlib.irradiance import beam_component, aoi, haydavies
 
 def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
                          pitch, max_rows=10, npoints=100):
@@ -401,7 +400,8 @@ 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, iam=1.0, npoints=100):
+                       albedo, model='isotropic', dni_extra=None, iam=1.0,
+                       npoints=100):
     r"""
     Calculate plane-of-array (POA) irradiance on one side of a row of modules.
 
@@ -457,6 +457,13 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
     albedo : numeric
         Surface albedo. [unitless]
 
+    model : str, default 'isotropic'
+        Irradiance model - can be one of 'isotropic' or 'haydavies'.
+
+    dni_extra : numeric, optional
+        Extraterrestrial direct normal irradiance. Required when
+        ``model='haydavies'``. [W/m2]
+
     iam : numeric, default 1.0
         Incidence angle modifier, the fraction of direct irradiance incident
         on the surface that is not reflected away. [unitless]
@@ -495,6 +502,27 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
     --------
     get_irradiance
     """
+    if model == 'haydavies':
+        if dni_extra is None:
+            raise ValueError(f'must supply dni_extra for {model} model')
+        # Call haydavies first time within the horizontal plane - to subtract
+        # circumsolar_horizontal from DHI
+        sky_diffuse_comps_horizontal = haydavies(0, 180, dhi, dni, dni_extra,
+                                                 solar_zenith, solar_azimuth,
+                                                 return_components=True)
+        circumsolar_horizontal = sky_diffuse_comps_horizontal['circumsolar']
+
+        # Call haydavies a second time where circumsolar_normal is facing
+        # directly towards sun, and can be added to DNI
+        sky_diffuse_comps_normal = haydavies(solar_zenith, solar_azimuth, dhi,
+                                             dni, dni_extra, solar_zenith,
+                                             solar_azimuth,
+                                             return_components=True)
+        circumsolar_normal = sky_diffuse_comps_normal['circumsolar']
+
+        dhi = dhi - circumsolar_horizontal
+        dni = dni + circumsolar_normal
+
     # Calculate some geometric quantities
     # rows to consider in front and behind current row
     # ensures that view factors to the sky are computed to within 5 degrees
@@ -580,8 +608,8 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
 
 def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                    gcr, height, pitch, ghi, dhi, dni,
-                   albedo, iam_front=1.0, iam_back=1.0,
-                   bifaciality=0.8, shade_factor=-0.02,
+                   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):
     """
     Get front and rear irradiance using the infinite sheds model.
@@ -643,6 +671,13 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
     albedo : numeric
         Surface albedo. [unitless]
 
+    model : str, default 'isotropic'
+        Irradiance model - can be one of 'isotropic' or 'haydavies'.
+
+    dni_extra : numeric, optional
+        Extraterrestrial direct normal irradiance. Required when
+        ``model='haydavies'``. [W/m2]
+
     iam_front : numeric, default 1.0
         Incidence angle modifier, the fraction of direct irradiance incident
         on the front surface that is not reflected away. [unitless]
@@ -720,13 +755,15 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         surface_tilt=surface_tilt, surface_azimuth=surface_azimuth,
         solar_zenith=solar_zenith, solar_azimuth=solar_azimuth,
         gcr=gcr, height=height, pitch=pitch, ghi=ghi, dhi=dhi, dni=dni,
-        albedo=albedo, iam=iam_front, npoints=npoints)
+        albedo=albedo, model=model, dni_extra=dni_extra, iam=iam_front,
+        npoints=npoints)
     # back side POA irradiance
     irrad_back = get_irradiance_poa(
         surface_tilt=backside_tilt, surface_azimuth=backside_sysaz,
         solar_zenith=solar_zenith, solar_azimuth=solar_azimuth,
         gcr=gcr, height=height, pitch=pitch, ghi=ghi, dhi=dhi, dni=dni,
-        albedo=albedo, iam=iam_back, npoints=npoints)
+        albedo=albedo, model=model, dni_extra=dni_extra, iam=iam_back,
+        npoints=npoints)
 
     colmap_front = {
         'poa_global': 'poa_front',

pvlib/tests/bifacial/test_infinite_sheds.py

@@ -353,3 +353,42 @@ def test_get_irradiance_limiting_gcr():
                       result['poa_back_sky_diffuse'])
     assert np.isclose(result['poa_front_ground_diffuse'],
                       result['poa_back_ground_diffuse'])
+
+
+def test_get_irradiance_with_haydavies():
+    # singleton inputs
+    solar_zenith = 0.
+    solar_azimuth = 180.
+    surface_tilt = 0.
+    surface_azimuth = 180.
+    gcr = 0.5
+    height = 1.
+    pitch = 1.
+    ghi = 1000.
+    dhi = 300.
+    dni = 700.
+    albedo = 0.
+    dni_extra = 1413.
+    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)
+    expected_front_diffuse = np.array([151.38])
+    expected_front_direct = np.array([848.62])
+    expected_front_global = expected_front_diffuse + expected_front_direct
+    assert np.isclose(result['poa_front'], expected_front_global)
+    assert np.isclose(result['poa_front_diffuse'], expected_front_diffuse)
+    assert np.isclose(result['poa_front_direct'], expected_front_direct)
+    assert np.isclose(result['poa_global'], result['poa_front'])
+    # test for when dni_extra is not supplied
+    with pytest.raises(ValueError, match='supply dni_extra for haydavies'):
+        result = infinite_sheds.get_irradiance(
+            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)