diff --git a/docs/sphinx/source/whatsnew/v0.9.4.rst b/docs/sphinx/source/whatsnew/v0.9.4.rst
index d62586107e..ae6b710bed 100644
--- a/docs/sphinx/source/whatsnew/v0.9.4.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.4.rst
@@ -10,6 +10,8 @@ Deprecations
 Enhancements
 ~~~~~~~~~~~~
 * Multiple code style issues fixed that were reported by LGTM analysis. (:issue:`1275`, :pull:`1559`)
+* Add optional ``return_components`` parameter to :py:func:`pvlib.irradiance.haydavies` to return
+  individual diffuse irradiance components (:issue:`1553`, :pull:`1568`)
 
 Bug fixes
 ~~~~~~~~~
@@ -35,4 +37,5 @@ Requirements
 Contributors
 ~~~~~~~~~~~~
 * Christian Orner (:ghuser:`chrisorner`)
-* Marcus Boumans (:ghuser:`bowie2211`)
\ No newline at end of file
+* Saurabh Aneja (:ghuser:`spaneja`)
+* Marcus Boumans (:ghuser:`bowie2211`)
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 03ddd13f5a..ddcb23a7ff 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -739,7 +739,8 @@ def klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith,
 
 
 def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
-              solar_zenith=None, solar_azimuth=None, projection_ratio=None):
+              solar_zenith=None, solar_azimuth=None, projection_ratio=None,
+              return_components=False):
     r'''
     Determine diffuse irradiance from the sky on a tilted surface using
     Hay & Davies' 1980 model
@@ -790,10 +791,27 @@ def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
         projection. Must supply ``solar_zenith`` and ``solar_azimuth``
         or supply ``projection_ratio``.
 
+    return_components : bool, default False
+        Flag used to decide whether to return the calculated diffuse components
+        or not.
+
     Returns
     --------
+    numeric, OrderedDict, or DataFrame
+        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.
+        The sky diffuse component of the solar radiation on a tilted
+        surface.
+
+    diffuse_components : OrderedDict (array input) or DataFrame (Series input)
+        Keys/columns are:
+            * sky_diffuse: Total sky diffuse
+            * isotropic
+            * circumsolar
+            * horizon
 
     Notes
     ------
@@ -830,10 +848,25 @@ def haydavies(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     term1 = 1 - AI
     term2 = 0.5 * (1 + tools.cosd(surface_tilt))
 
-    sky_diffuse = dhi * (AI * Rb + term1 * term2)
-    sky_diffuse = np.maximum(sky_diffuse, 0)
+    poa_isotropic = np.maximum(dhi * term1 * term2, 0)
+    poa_circumsolar = np.maximum(dhi * (AI * Rb), 0)
+    sky_diffuse = poa_isotropic + poa_circumsolar
 
-    return sky_diffuse
+    if return_components:
+        diffuse_components = OrderedDict()
+        diffuse_components['sky_diffuse'] = sky_diffuse
+
+        # Calculate the individual components
+        diffuse_components['isotropic'] = poa_isotropic
+        diffuse_components['circumsolar'] = poa_circumsolar
+        diffuse_components['horizon'] = np.where(
+            np.isnan(diffuse_components['isotropic']), np.nan, 0.)
+
+        if isinstance(sky_diffuse, pd.Series):
+            diffuse_components = pd.DataFrame(diffuse_components)
+        return diffuse_components
+    else:
+        return sky_diffuse
 
 
 def reindl(surface_tilt, surface_azimuth, dhi, dni, ghi, dni_extra,
diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py
index 8f3c565154..c87e8dcb47 100644
--- a/pvlib/tests/test_irradiance.py
+++ b/pvlib/tests/test_irradiance.py
@@ -207,6 +207,45 @@ def test_haydavies(irrad_data, ephem_data, dni_et):
     assert_allclose(result, [0, 27.1775, 102.9949, 33.1909], atol=1e-4)
 
 
+def test_haydavies_components(irrad_data, ephem_data, dni_et):
+    expected = pd.DataFrame(np.array(
+        [[0, 27.1775, 102.9949, 33.1909],
+         [0, 27.1775, 30.1818, 27.9837],
+         [0, 0, 72.8130, 5.2071],
+         [0, 0, 0, 0]]).T,
+        columns=['sky_diffuse', 'isotropic', 'circumsolar', 'horizon'],
+        index=irrad_data.index
+    )
+    # pandas
+    result = irradiance.haydavies(
+        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
+        ephem_data['apparent_zenith'], ephem_data['azimuth'],
+        return_components=True)
+    assert_frame_equal(result, expected, check_less_precise=4)
+    # numpy
+    result = irradiance.haydavies(
+        40, 180, irrad_data['dhi'].values, irrad_data['dni'].values, dni_et,
+        ephem_data['apparent_zenith'].values, ephem_data['azimuth'].values,
+        return_components=True)
+    assert_allclose(result['sky_diffuse'], expected['sky_diffuse'], atol=1e-4)
+    assert_allclose(result['isotropic'], expected['isotropic'], atol=1e-4)
+    assert_allclose(result['circumsolar'], expected['circumsolar'], atol=1e-4)
+    assert_allclose(result['horizon'], expected['horizon'], atol=1e-4)
+    assert isinstance(result, dict)
+    # scalar
+    result = irradiance.haydavies(
+        40, 180, irrad_data['dhi'].values[-1], irrad_data['dni'].values[-1],
+        dni_et[-1], ephem_data['apparent_zenith'].values[-1],
+        ephem_data['azimuth'].values[-1], return_components=True)
+    assert_allclose(result['sky_diffuse'], expected['sky_diffuse'][-1],
+                    atol=1e-4)
+    assert_allclose(result['isotropic'], expected['isotropic'][-1],
+                    atol=1e-4)
+    assert_allclose(result['circumsolar'], expected['circumsolar'][-1],
+                    atol=1e-4)
+    assert_allclose(result['horizon'], expected['horizon'][-1], atol=1e-4)
+    assert isinstance(result, dict)
+
 def test_reindl(irrad_data, ephem_data, dni_et):
     result = irradiance.reindl(
         40, 180, irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'],