Add support for return_components to Reindl

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

@@ -23,6 +23,9 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
+* Add  ``return_components`` kwarg to :py:func:`pvlib.irradiance.reindl` to
+  support returning the components of sky diffuse irradiance.
+  (:issue:`2750`, :pull:`2775`)
 * Add iotools functions to retrieve irradiance and weather data from NSRDB PSM4 Polar,
   which provides satellite-derived irradiance data above 60 degree latitude.
   :py:func:`~pvlib.iotools.get_nsrdb_psm4_polar` and

pvlib/irradiance.py

@@ -872,7 +872,7 @@ def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
 
 
 def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
-           solar_zenith, solar_azimuth):
+           solar_zenith, solar_azimuth, return_components=False):
     r'''
     Determine the diffuse irradiance from the sky on a tilted surface using
     the Reindl (1990) model.
@@ -911,10 +911,30 @@ def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
     solar_azimuth : numeric
         Solar azimuth angles. See :term:`solar_azimuth`. [°]
 
+    return_components : bool, default ``False``
+        If ``False``, ``poa_sky_diffuse`` is returned.
+        If ``True``, ``diffuse_components`` is returned.
+
     Returns
     -------
+    numeric, dict, or DataFrame
+        Return type controlled by ``return_components`` argument.
+        If ``return_components=False``, ``poa_sky_diffuse`` is returned.
+        If ``return_components=True``, ``diffuse_components`` is returned.
+
     poa_sky_diffuse : numeric
-        The sky diffuse component of the solar radiation. [Wm⁻²]
+        The sky diffuse component of irradiance on a tilted plane. [Wm⁻²]
+
+    diffuse_components : dict (array input) or DataFrame (Series input)
+        Keys/columns are:
+            * poa_sky_diffuse: The sky diffuse component of irradiance on a
+              tilted plane. [Wm⁻²]
+            * poa_isotropic: The portion of sky diffuse irradiance on a tilted
+              plane from the isotropic sky dome. [Wm⁻²]
+            * poa_circumsolar: The portion of sky diffuse irradiance on a
+              tilted plane from the circumsolar region. [Wm⁻²]
+            * poa_horizon: The portion of sky diffuse irradiance on a tilted
+              plane from the horizon. [Wm⁻²]
 
     Notes
     -----
@@ -938,8 +958,12 @@ def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
     Implementation is based on Loutzenhiser et al.
     (2007) [3]_, Equation 8. The beam and ground reflectance portion of the
     equation have been removed, therefore the model described here generates
-    ONLY the diffuse radiation from the sky and circumsolar, so the form of the
-    equation varies slightly from Equation 8 in [3]_.
+    ONLY the diffuse radiation from the sky, circumsolar, and horizon
+    brightening, so the form of the equation varies slightly from Equation 8
+    in [3]_.
+
+    For clarity, the horizon component in ``reindl`` corresponds to the term
+    added on top of the ``haydavies`` formulation, on which ``reindl`` builds.
 
     References
     ----------
@@ -975,16 +999,31 @@ def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
     HB = dni * cos_solar_zenith
     HB = np.maximum(HB, 0)
 
-    # these are the () and [] sub-terms of the second term of eqn 8
-    term1 = 1 - AI
-    term2 = 0.5 * (1 + tools.cosd(surface_tilt))
+    SVF = (1 + tools.cosd(surface_tilt)) / 2
+
     with np.errstate(invalid='ignore', divide='ignore'):
         hb_to_ghi = np.where(ghi == 0, 0, np.divide(HB, ghi))
-    term3 = 1 + np.sqrt(hb_to_ghi) * (tools.sind(0.5 * surface_tilt)**3)
-    sky_diffuse = dhi * (AI * Rb + term1 * term2 * term3)
-    sky_diffuse = np.maximum(sky_diffuse, 0)
+    h = np.sqrt(hb_to_ghi) * (tools.sind(surface_tilt / 2) ** 3)
 
-    return sky_diffuse
+    term1 = (1 - AI) * SVF
+    term2 = AI * Rb
+    term3 = term1 * h
+
+    poa_sky_diffuse = dhi * (term1 + term2 + term3)
+
+    if return_components:
+        diffuse_components = {
+            'poa_sky_diffuse': poa_sky_diffuse,
+            'poa_isotropic': dhi * term1,
+            'poa_circumsolar': dhi * term2,
+            'poa_horizon': dhi * term3
+        }
+
+        if isinstance(poa_sky_diffuse, pd.Series):
+            diffuse_components = pd.DataFrame(diffuse_components)
+        return diffuse_components
+    else:
+        return poa_sky_diffuse
 
 
 def king(surface_tilt, dhi, ghi, solar_zenith):

tests/test_irradiance.py

@@ -249,6 +249,45 @@ def test_reindl(irrad_data, ephem_data, dni_et):
     assert_allclose(result, [0., 27.9412, 104.1317, 34.1663], atol=1e-4)
 
 
+def test_reindl_components(irrad_data, ephem_data, dni_et):
+    keys = ['poa_sky_diffuse', 'poa_isotropic', 'poa_circumsolar',
+            'poa_horizon']
+    expected = pd.DataFrame(np.array(
+        [[0, 27.941170, 104.131724, 34.166258],
+         [0, 27.177514, 30.181807, 27.983728],
+         [0, 0, 72.813055, 5.207138],
+         [0, 0.763656, 1.136862, 0.975393]]).T,
+        columns=keys,
+        index=irrad_data.index
+    )
+    # pandas
+    result = irradiance.reindl(
+        40, 180, irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'],
+        dni_et, ephem_data['apparent_zenith'], ephem_data['azimuth'],
+        return_components=True)
+    assert_frame_equal(result, expected, check_less_precise=4)
+    # numpy
+    result = irradiance.reindl(
+        40, 180, irrad_data['dhi'].to_numpy(), irrad_data['dni'].to_numpy(),
+        irrad_data['ghi'].to_numpy(), dni_et,
+        ephem_data['apparent_zenith'].to_numpy(),
+        ephem_data['azimuth'].to_numpy(), return_components=True)
+    for key in keys:
+        assert_allclose(result[key], expected[key], atol=1e-4)
+    assert isinstance(result, dict)
+    # scalar
+    result = irradiance.reindl(
+        40, 180, irrad_data['dhi'].to_numpy()[-1],
+        irrad_data['dni'].to_numpy()[-1],
+        irrad_data['ghi'].to_numpy()[-1], dni_et[-1],
+        ephem_data['apparent_zenith'].to_numpy()[-1],
+        ephem_data['azimuth'].to_numpy()[-1],
+        return_components=True)
+    for key in keys:
+        assert_allclose(result[key], expected[key].iloc[-1], atol=1e-4)
+    assert isinstance(result, dict)
+
+
 def test_king(irrad_data, ephem_data):
     result = irradiance.king(40, irrad_data['dhi'], irrad_data['ghi'],
                              ephem_data['apparent_zenith'])