diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 50f02426de..3efb1ea206 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -873,7 +873,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.
@@ -912,10 +912,35 @@ 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
+        Flag used to decide whether to return the calculated diffuse components
+        or not.
+
     Returns
     -------
-    poa_sky_diffuse : numeric
-        The sky diffuse component of the solar radiation. [Wm⁻²]
+<<<<<<< HEAD
+    numeric, Dict, or DataFrame
+=======
+    numeric, OrderedDict, or DataFrame
+>>>>>>> 3cbc4adc (Update docstrings for Reindl)
+        Return type controlled by ``return_components`` argument.
+        If ``return_components=False``, `sky_diffuse` is returned.
+        If ``return_components=True``, `diffuse_components` is returned.
+
+    sky_diffuse : numeric
+        The sky diffuse component of the solar radiation on a tilted
+        surface.
+
+<<<<<<< HEAD
+    diffuse_components : Dict (array input) or DataFrame (Series input)
+=======
+    diffuse_components : OrderedDict (array input) or DataFrame (Series input)
+>>>>>>> 3cbc4adc (Update docstrings for Reindl)
+        Keys/columns are:
+            * poa_sky_diffuse: Total sky diffuse
+            * poa_isotropic
+            * poa_circumsolar
+            * poa_horizon
 
     Notes
     -----
@@ -939,8 +964,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
     ----------
@@ -973,16 +1002,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 = (1 - AI) * SVF * h
+
+    sky_diffuse = dhi * (term1 + term2 + term3)
+
+    if return_components:
+        diffuse_components = {
+            'poa_sky_diffuse': sky_diffuse,
+            'poa_isotropic': dhi * term1,
+            'poa_circumsolar': dhi * term2,
+            'poa_horizon': dhi * term3
+        }
+
+        if isinstance(sky_diffuse, pd.Series):
+            diffuse_components = pd.DataFrame(diffuse_components)
+        return diffuse_components
+    else:
+        return sky_diffuse
 
 
 def king(surface_tilt, dhi, ghi, solar_zenith):
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index a416636ae9..efebd3caa2 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -251,6 +251,43 @@ 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'].values, irrad_data['dni'].values,
+        irrad_data['ghi'].values, dni_et, ephem_data['apparent_zenith'].values,
+        ephem_data['azimuth'].values, 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'].values[-1], irrad_data['dni'].values[-1],
+        irrad_data['ghi'].values[-1], dni_et[-1],
+        ephem_data['apparent_zenith'].values[-1],
+        ephem_data['azimuth'].values[-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'])