diff --git a/docs/examples/plot_diffuse_aoi_correction.py b/docs/examples/plot_diffuse_aoi_correction.py
new file mode 100644
index 0000000000..5cd2f5cb41
--- /dev/null
+++ b/docs/examples/plot_diffuse_aoi_correction.py
@@ -0,0 +1,119 @@
+"""
+Diffuse IAM Calculation
+=======================
+
+Integrating an IAM model across angles to determine the overall reflection
+loss for diffuse irradiance.
+"""
+
+# %%
+# The fraction of light reflected from the front of a module depends on the
+# angle of incidence (AOI) of the light compared to the panel surface.  The
+# greater the AOI, the larger the reflected fraction is.  The incident angle
+# modifier (IAM) is defined as the ratio of light transmitted at the given
+# AOI to transmitted light at normal incidence.
+# Several models exist to calculate the IAM for a given incidence
+# angle (e.g. :py:func:`pvlib.iam.ashrae`, :py:func:`pvlib.iam.martin_ruiz`,
+# :py:func:`pvlib.iam.sapm`, :py:func:`pvlib.iam.physical`).
+# However, evaluating the IAM for diffuse light is
+# not as straightforward because it comes from all directions and therefore
+# has a range of angles of incidence.  Here we show how to integrate the effect
+# of AOI reflection across this AOI range using the process described in [1]_.
+# In particular, we will recreate Figures 3, 4, and 5 in that paper.
+#
+# References
+# ----------
+#  .. [1] B. Marion "Numerical method for angle-of-incidence correction
+#     factors for diffuse radiation incident photovoltaic modules",
+#     Solar Energy, Volume 147, Pages 344-348. 2017.
+#     DOI: 10.1016/j.solener.2017.03.027
+#
+#  .. [2] Duffie, John A. & Beckman, William A. (2013). Solar Engineering
+#     of Thermal Processes.  DOI: 10.1002/9781118671603
+
+
+from pvlib.iam import marion_diffuse, physical
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# %%
+# IAM Model
+# ---------
+#
+# The IAM model used to generate the figures in [1]_ uses Snell's, Fresnel's,
+# and Beer's laws to determine the fraction of light transmitted through the
+# air-glass interface as a function of AOI.
+# The function :py:func:`pvlib.iam.physical` implements this model, except it
+# also includes an exponential term to model attenuation in the glazing layer.
+# To be faithful to Marion's implementation, we will disable this extinction
+# term by setting the attenuation coefficient ``K`` parameter to zero.
+# For more details on this IAM model, see [2]_.
+#
+# Marion generated diffuse irradiance modifiers for two cases:  a standard
+# uncoated glass with index of refraction n=1.526 and a glass with
+# anti-reflective (AR) coating with n=1.3.
+# Comparing the IAM model across AOI recreates Figure 3 in [1]_:
+
+aoi = np.arange(0, 91)
+iam_no_coating = physical(aoi, n=1.526, K=0)
+iam_ar_coating = physical(aoi, n=1.3, K=0)
+
+plt.plot(aoi, iam_ar_coating, c='b', label='$F_b$, AR coated, n=1.3')
+plt.plot(aoi, iam_no_coating, c='r', label='$F_b$, uncoated, n=1.526')
+plt.xlabel(r'Angle-of-Incidence, AOI $(\degree)$')
+plt.ylabel('Diffuse Incidence Angle Modifier')
+plt.legend()
+plt.ylim([0, 1.2])
+plt.grid()
+
+# %%
+# Diffuse sky, ground, and horizon IAM
+# ------------------------------------
+#
+# Now that we have an AOI model, we use :py:func:`pvlib.iam.marion_diffuse`
+# to integrate it across solid angle and determine diffuse irradiance IAM.
+# Marion defines three types of diffuse irradiance:
+# sky, horizon, and ground-reflected.  The diffuse IAM value is evaluated
+# independently for each type.
+
+tilts = np.arange(0, 91, 2.5)
+
+# marion_diffuse calculates all three IAM values (sky, horizon, ground)
+iam_no_coating = marion_diffuse('physical', tilts, n=1.526, K=0)
+iam_ar_coating = marion_diffuse('physical', tilts, n=1.3, K=0)
+
+# %%
+# First we recreate Figure 4 in [1]_, showing the dependence of the sky diffuse
+# incidence angle modifier on module tilt.
+
+plt.plot(tilts, iam_ar_coating['sky'], c='b', marker='^',
+         label='$F_{sky}$, AR coated, n=1.3')
+plt.plot(tilts, iam_no_coating['sky'], c='r', marker='x',
+         label='$F_{sky}$, uncoated, n=1.526')
+plt.ylim([0.9, 1.0])
+plt.xlabel(r'PV Module Tilt, $\beta (\degree)$')
+plt.ylabel('Diffuse Incidence Angle Modifier')
+plt.grid()
+plt.legend()
+plt.show()
+
+# %%
+# Now we recreate Figure 5 in [1]_, showing the dependence of the diffuse iam
+# values for horizon and ground diffuse irradiance on module tilt.  Note that
+# :py:func:`pvlib.iam.marion_diffuse` defaults to using 1800 points for the
+# horizon case (instead of 180 like the others) to match [1]_.
+
+plt.plot(tilts, iam_ar_coating['horizon'], c='b', marker='^',
+         label='$F_{hor}$, AR coated, n=1.3')
+plt.plot(tilts, iam_no_coating['horizon'], c='r', marker='x',
+         label='$F_{hor}$, uncoated, n=1.526')
+plt.plot(tilts, iam_ar_coating['ground'], c='b', marker='s',
+         label='$F_{grd}$, AR coated, n=1.3')
+plt.plot(tilts, iam_no_coating['ground'], c='r', marker='+',
+         label='$F_{grd}$, uncoated, n=1.526')
+plt.xlabel(r'PV Module Tilt, $\beta (\degree)$')
+plt.ylabel('Diffuse Incidence Angle Modifier')
+plt.grid()
+plt.legend()
+plt.show()
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index eb5691f567..a0460c5991 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -223,6 +223,8 @@ Incident angle modifiers
    iam.martin_ruiz_diffuse
    iam.sapm
    iam.interp
+   iam.marion_diffuse
+   iam.marion_integrate
 
 PV temperature models
 ---------------------
diff --git a/docs/sphinx/source/whatsnew/v0.8.0.rst b/docs/sphinx/source/whatsnew/v0.8.0.rst
index 4b6ce43ac4..a664fc7a76 100644
--- a/docs/sphinx/source/whatsnew/v0.8.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.8.0.rst
@@ -11,6 +11,9 @@ API Changes
 Enhancements
 ~~~~~~~~~~~~
 * Update :func:`~pvlib.bifacial.pvfactors_timeseries` to run with ``pvfactors`` v1.4.1 (:issue:`902`)(:pull:`934`)
+* Add :py:func:`pvlib.iam.marion_diffuse` and 
+  :py:func:`pvlib.iam.marion_integrate` to calculate IAM values for
+  diffuse irradiance. (:pull:`984`)
 
 Bug fixes
 ~~~~~~~~~
@@ -29,6 +32,8 @@ Documentation
 * Clarify units for heat loss factors in
   :py:func:`pvlib.temperature.pvsyst_cell` and
   :py:func:`pvlib.temperature.faiman`. (:pull:`960`)
+* Add a gallery example of calculating diffuse IAM using
+  :py:func:`pvlib.iam.marion_diffuse`. (:pull:`984`)
 
 Requirements
 ~~~~~~~~~~~~
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 6143439289..243e097443 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -10,6 +10,7 @@
 
 import numpy as np
 import pandas as pd
+import functools
 from pvlib.tools import cosd, sind, tand, asind
 
 # a dict of required parameter names for each IAM model
@@ -527,3 +528,221 @@ def sapm(aoi, module, upper=None):
         iam = pd.Series(iam, aoi.index)
 
     return iam
+
+
+def marion_diffuse(model, surface_tilt, **kwargs):
+    """
+    Determine diffuse irradiance incidence angle modifiers using Marion's
+    method of integrating over solid angle.
+
+    Parameters
+    ----------
+    model : str
+        The IAM function to evaluate across solid angle. Must be one of
+        `'ashrae', 'physical', 'martin_ruiz', 'sapm'`.
+
+    surface_tilt : numeric
+        Surface tilt angles in decimal degrees.
+        The tilt angle is defined as degrees from horizontal
+        (e.g. surface facing up = 0, surface facing horizon = 90).
+
+    **kwargs
+        Extra parameters passed to the IAM function.
+
+    Returns
+    -------
+    iam : dict
+        IAM values for each type of diffuse irradiance:
+
+            * 'sky': radiation from the sky dome (zenith <= 90)
+            * 'horizon': radiation from the region of the sky near the horizon
+              (89.5 <= zenith <= 90)
+            * 'ground': radiation reflected from the ground (zenith >= 90)
+
+        See [1]_ for a detailed description of each class.
+
+    See Also
+    --------
+    pvlib.iam.marion_integrate
+
+    References
+    ----------
+    .. [1] B. Marion "Numerical method for angle-of-incidence correction
+       factors for diffuse radiation incident photovoltaic modules",
+       Solar Energy, Volume 147, Pages 344-348. 2017.
+       DOI: 10.1016/j.solener.2017.03.027
+
+    Examples
+    --------
+    >>> marion_diffuse('physical', surface_tilt=20)
+    {'sky': 0.9539178294437575,
+     'horizon': 0.7652650139134007,
+     'ground': 0.6387140117795903}
+
+    >>> marion_diffuse('ashrae', [20, 30], b=0.04)
+    {'sky': array([0.96748999, 0.96938408]),
+     'horizon': array([0.86478428, 0.91825792]),
+     'ground': array([0.77004435, 0.8522436 ])}
+    """
+
+    models = {
+        'physical': physical,
+        'ashrae': ashrae,
+        'sapm': sapm,
+        'martin_ruiz': martin_ruiz,
+    }
+
+    try:
+        iam_model = models[model]
+    except KeyError:
+        raise ValueError('model must be one of: ' + str(list(models.keys())))
+
+    iam_function = functools.partial(iam_model, **kwargs)
+    iam = {}
+    for region in ['sky', 'horizon', 'ground']:
+        iam[region] = marion_integrate(iam_function, surface_tilt, region)
+
+    return iam
+
+
+def marion_integrate(function, surface_tilt, region, num=None):
+    """
+    Integrate an incidence angle modifier (IAM) function over solid angle
+    to determine a diffuse irradiance correction factor using Marion's method.
+
+    This lower-level function actually performs the IAM integration for the
+    specified solid angle region.
+
+    Parameters
+    ----------
+    function : callable(aoi)
+        The IAM function to evaluate across solid angle. The function must
+        be vectorized and take only one parameter, the angle of incidence in
+        degrees.
+
+    surface_tilt : numeric
+        Surface tilt angles in decimal degrees.
+        The tilt angle is defined as degrees from horizontal
+        (e.g. surface facing up = 0, surface facing horizon = 90).
+
+    region : {'sky', 'horizon', 'ground'}
+        The region to integrate over. Must be one of:
+
+            * 'sky': radiation from the sky dome (zenith <= 90)
+            * 'horizon': radiation from the region of the sky near the horizon
+              (89.5 <= zenith <= 90)
+            * 'ground': radiation reflected from the ground (zenith >= 90)
+
+        See [1]_ for a detailed description of each class.
+
+    num : int, optional
+        The number of increments in the zenith integration.
+        If not specified, N will follow the values used in [1]_:
+
+            * 'sky' or 'ground': num = 180
+            * 'horizon': num = 1800
+
+    Returns
+    -------
+    iam : numeric
+        AOI diffuse correction factor for the specified region.
+
+    See Also
+    --------
+    pvlib.iam.marion_diffuse
+
+    References
+    ----------
+    .. [1] B. Marion "Numerical method for angle-of-incidence correction
+       factors for diffuse radiation incident photovoltaic modules",
+       Solar Energy, Volume 147, Pages 344-348. 2017.
+       DOI: 10.1016/j.solener.2017.03.027
+
+    Examples
+    --------
+    >>> marion_integrate(pvlib.iam.ashrae, 20, 'sky')
+    0.9596085829811408
+
+    >>> from functools import partial
+    >>> f = partial(pvlib.iam.physical, n=1.3)
+    >>> marion_integrate(f, [20, 30], 'sky')
+    array([0.96225034, 0.9653219 ])
+    """
+
+    if num is None:
+        if region in ['sky', 'ground']:
+            num = 180
+        elif region == 'horizon':
+            num = 1800
+        else:
+            raise ValueError('Invalid region: {}'.format(region))
+
+    beta = np.radians(surface_tilt)
+    if isinstance(beta, pd.Series):
+        # convert Series to np array for broadcasting later
+        beta = beta.values
+    ai = np.pi/num  # angular increment
+
+    phi_range = np.linspace(0, np.pi, num, endpoint=False)
+    psi_range = np.linspace(0, 2*np.pi, 2*num, endpoint=False)
+
+    # the pseudocode in [1] do these checks at the end, but it's
+    # faster to do this criteria check up front instead of later.
+    if region == 'sky':
+        mask = phi_range + ai <= np.pi/2
+    elif region == 'horizon':
+        lo = 89.5 * np.pi/180
+        hi = np.pi/2
+        mask = (lo <= phi_range) & (phi_range + ai <= hi)
+    elif region == 'ground':
+        mask = (phi_range >= np.pi/2)
+    else:
+        raise ValueError('Invalid region: {}'.format(region))
+    phi_range = phi_range[mask]
+
+    # fast Cartesian product of phi and psi
+    angles = np.array(np.meshgrid(phi_range, psi_range)).T.reshape(-1, 2)
+    # index with single-element lists to maintain 2nd dimension so that
+    # these angle arrays broadcast across the beta array
+    phi_1 = angles[:, [0]]
+    psi_1 = angles[:, [1]]
+    phi_2 = phi_1 + ai
+    # psi_2 = psi_1 + ai  # not needed
+    phi_avg = phi_1 + 0.5*ai
+    psi_avg = psi_1 + 0.5*ai
+    term_1 = np.cos(beta) * np.cos(phi_avg)
+    # The AOI formula includes a term based on the difference between
+    # panel azimuth and the photon azimuth, but because we assume each class
+    # of diffuse irradiance is isotropic and we are integrating over all
+    # angles, it doesn't matter what panel azimuth we choose (i.e., the
+    # system is rotationally invariant).  So we choose gamma to be zero so
+    # that we can omit it from the cos(psi_avg) term.
+    # Marion's paper mentions this in the Section 3 pseudocode:
+    # "set gamma to pi (or any value between 0 and 2pi)"
+    term_2 = np.sin(beta) * np.sin(phi_avg) * np.cos(psi_avg)
+    cosaoi = term_1 + term_2
+    aoi = np.arccos(cosaoi)
+    # simplify Eq 8, (psi_2 - psi_1) is always ai
+    dAs = ai * (np.cos(phi_1) - np.cos(phi_2))
+    cosaoi_dAs = cosaoi * dAs
+    # apply the final AOI check, zeroing out non-passing points
+    mask = aoi < np.pi/2
+    cosaoi_dAs = np.where(mask, cosaoi_dAs, 0)
+    numerator = np.sum(function(np.degrees(aoi)) * cosaoi_dAs, axis=0)
+    denominator = np.sum(cosaoi_dAs, axis=0)
+
+    with np.errstate(invalid='ignore'):
+        # in some cases, no points pass the criteria
+        # (e.g. region='ground', surface_tilt=0), so we override the division
+        # by zero to set Fd=0.  Also, preserve nans in beta.
+        Fd = np.where((denominator != 0) | ~np.isfinite(beta),
+                      numerator / denominator,
+                      0)
+
+    # preserve input type
+    if np.isscalar(surface_tilt):
+        Fd = Fd.item()
+    elif isinstance(surface_tilt, pd.Series):
+        Fd = pd.Series(Fd, surface_tilt.index)
+
+    return Fd
diff --git a/pvlib/tests/test_iam.py b/pvlib/tests/test_iam.py
index 8e04ee2975..7fc204bc95 100644
--- a/pvlib/tests/test_iam.py
+++ b/pvlib/tests/test_iam.py
@@ -218,3 +218,114 @@ def test_sapm_limits():
 
     module_parameters = {'B0': -5, 'B1': 0, 'B2': 0, 'B3': 0, 'B4': 0, 'B5': 0}
     assert _iam.sapm(1, module_parameters) == 0
+
+
+def test_marion_diffuse_model(mocker):
+    # 1: return values are correct
+    # 2: the underlying models are called appropriately
+    ashrae_expected = {
+        'sky': 0.9596085829811408,
+        'horizon': 0.8329070417832541,
+        'ground': 0.719823559106309
+    }
+    physical_expected = {
+        'sky': 0.9539178294437575,
+        'horizon': 0.7652650139134007,
+        'ground': 0.6387140117795903
+    }
+    ashrae_spy = mocker.spy(_iam, 'ashrae')
+    physical_spy = mocker.spy(_iam, 'physical')
+
+    ashrae_actual = _iam.marion_diffuse('ashrae', 20)
+    assert ashrae_spy.call_count == 3  # one call for each of the 3 regions
+    assert physical_spy.call_count == 0
+    physical_actual = _iam.marion_diffuse('physical', 20)
+    assert ashrae_spy.call_count == 3
+    assert physical_spy.call_count == 3
+
+    for k, v in ashrae_expected.items():
+        np.testing.assert_allclose(ashrae_actual[k], v)
+
+    for k, v in physical_expected.items():
+        np.testing.assert_allclose(physical_actual[k], v)
+
+
+def test_marion_diffuse_kwargs():
+    # kwargs get passed to underlying model
+    expected = {
+        'sky': 0.967489994422575,
+        'horizon': 0.8647842827418412,
+        'ground': 0.7700443455928433
+    }
+    actual = _iam.marion_diffuse('ashrae', 20, b=0.04)
+
+    for k, v in expected.items():
+        np.testing.assert_allclose(actual[k], v)
+
+
+def test_marion_diffuse_invalid():
+    with pytest.raises(ValueError):
+        _iam.marion_diffuse('not_a_model', 20)
+
+
+@pytest.mark.parametrize('region,N,expected', [
+    ('sky', 180, 0.9596085829811408),
+    ('horizon', 1800, 0.8329070417832541),
+    ('ground', 180, 0.719823559106309)
+])
+def test_marion_integrate_scalar(region, N, expected):
+    actual = _iam.marion_integrate(_iam.ashrae, 20, region, N)
+    assert_allclose(actual, expected)
+
+    with np.errstate(invalid='ignore'):
+        actual = _iam.marion_integrate(_iam.ashrae, np.nan, region, N)
+    expected = np.nan
+    assert_allclose(actual, expected)
+
+
+@pytest.mark.parametrize('region,N,expected', [
+    ('sky', 180, [0.9523611991069362, 0.9596085829811408, 0.9619811198105501]),
+    ('horizon', 1800, [0.0, 0.8329070417832541, 0.8987287652347437]),
+    ('ground', 180, [0.0, 0.719823559106309, 0.8186360238536674])
+])
+def test_marion_integrate_list(region, N, expected):
+    actual = _iam.marion_integrate(_iam.ashrae, [0, 20, 30], region, N)
+    assert_allclose(actual, expected)
+
+    with np.errstate(invalid='ignore'):
+        actual = _iam.marion_integrate(_iam.ashrae, [0, 20, np.nan], region, N)
+    assert_allclose(actual, [expected[0], expected[1], np.nan])
+
+
+@pytest.mark.parametrize('region,N,expected', [
+    ('sky', 180, [0.9523611991069362, 0.9596085829811408, 0.9619811198105501]),
+    ('horizon', 1800, [0.0, 0.8329070417832541, 0.8987287652347437]),
+    ('ground', 180, [0.0, 0.719823559106309, 0.8186360238536674])
+])
+def test_marion_integrate_series(region, N, expected):
+    idx = pd.date_range('2019-01-01', periods=3, freq='h')
+    tilt = pd.Series([0, 20, 30], index=idx)
+    expected = pd.Series(expected, index=idx)
+    actual = _iam.marion_integrate(_iam.ashrae, tilt, region, N)
+    assert_series_equal(actual, expected)
+
+    tilt.iloc[1] = np.nan
+    expected.iloc[1] = np.nan
+    with np.errstate(invalid='ignore'):
+        actual = _iam.marion_integrate(_iam.ashrae, tilt, region, N)
+    assert_allclose(actual, expected)
+
+
+def test_marion_integrate_ground_flat():
+    iam = _iam.marion_integrate(_iam.ashrae, 0, 'horizon', 1800)
+    assert_allclose(iam, 0)
+
+
+def test_marion_integrate_invalid():
+    # check for invalid region string.  this actually gets checked twice,
+    # with the difference being whether `num` is specified or not.
+    with pytest.raises(ValueError):
+        _iam.marion_integrate(_iam.ashrae, 0, 'bad')
+
+    with pytest.raises(ValueError):
+        _iam.marion_integrate(_iam.ashrae, 0, 'bad', 180)