diff --git a/docs/examples/plot_passias_diffuse_shading.py b/docs/examples/plot_passias_diffuse_shading.py
new file mode 100644
index 0000000000..989e977fdb
--- /dev/null
+++ b/docs/examples/plot_passias_diffuse_shading.py
@@ -0,0 +1,84 @@
+"""
+Diffuse Self-Shading
+====================
+
+Modeling the reduction in diffuse irradiance caused by row-to-row diffuse
+shading.
+"""
+
+# %%
+# The term "self-shading" usually refers to adjacent rows blocking direct
+# irradiance and casting shadows on each other. However, the concept also
+# applies to diffuse irradiance because rows block a portion of the sky
+# dome even when the sun is high in the sky. The irradiance loss fraction
+# depends on how tightly the rows are packed and where on the module the
+# loss is evaluated -- a point near the top of edge of a module will see
+# more of the sky than a point near the bottom edge.
+#
+# This example uses the approach presented by Passias and Källbäck in [1]_
+# and recreates two figures from that paper using
+# :py:func:`pvlib.shading.masking_angle_passias` and
+# :py:func:`pvlib.shading.sky_diffuse_passias`.
+#
+# References
+# ----------
+#  .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
+#     panels", Solar Cells, Volume 11, Pages 281-291.  1984.
+#     DOI: 10.1016/0379-6787(84)90017-6
+
+from pvlib import shading, irradiance
+import matplotlib.pyplot as plt
+import numpy as np
+
+# %%
+# First we'll recreate Figure 4, showing how the average masking angle varies
+# with array tilt and array packing. The masking angle of a given point on a
+# module is the angle from horizontal to the next row's top edge and represents
+# the portion of the sky dome blocked by the next row. Because it changes
+# from the bottom to the top of a module, the average across the module is
+# calculated. In [1]_, ``k`` refers to the ratio of row pitch to row slant
+# height (i.e. 1 / GCR).
+
+surface_tilt = np.arange(0, 90, 0.5)
+
+plt.figure()
+for k in [1, 1.5, 2, 2.5, 3, 4, 5, 7, 10]:
+    gcr = 1/k
+    psi = shading.masking_angle_passias(surface_tilt, gcr)
+    plt.plot(surface_tilt, psi, label='k={}'.format(k))
+
+plt.xlabel('Inclination angle [degrees]')
+plt.ylabel('Average masking angle [degrees]')
+plt.legend()
+plt.show()
+
+# %%
+# So as the array is packed tighter (decreasing ``k``), the average masking
+# angle increases.
+#
+# Next we'll recreate Figure 5. Note that the y-axis here is the ratio of
+# diffuse plane of array irradiance (after accounting for shading) to diffuse
+# horizontal irradiance. This means that the deviation from 100% is due to the
+# combination of self-shading and the fact that being at a tilt blocks off
+# the portion of the sky behind the row. The first effect is modeled with
+# :py:func:`pvlib.shading.sky_diffuse_passias` and the second with
+# :py:func:`pvlib.irradiance.isotropic`.
+
+plt.figure()
+for k in [1, 1.5, 2, 10]:
+    gcr = 1/k
+    psi = shading.masking_angle_passias(surface_tilt, gcr)
+    shading_loss = shading.sky_diffuse_passias(psi)
+    transposition_ratio = irradiance.isotropic(surface_tilt, dhi=1.0)
+    relative_diffuse = transposition_ratio * (1-shading_loss) * 100  # %
+    plt.plot(surface_tilt, relative_diffuse, label='k={}'.format(k))
+
+plt.xlabel('Inclination angle [degrees]')
+plt.ylabel('Relative diffuse irradiance [%]')
+plt.ylim(0, 105)
+plt.legend()
+plt.show()
+
+# %%
+# As ``k`` decreases, GCR increases, so self-shading loss increases and
+# collected diffuse irradiance decreases.
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index dd7c618646..a7087edad0 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -354,6 +354,12 @@ Effects on PV System Output
    soiling.hsu
    soiling.kimber
 
+.. autosummary::
+   :toctree: generated/
+
+   shading.masking_angle
+   shading.masking_angle_passias
+   shading.sky_diffuse_passias
 
 
 Tracking
diff --git a/docs/sphinx/source/whatsnew/v0.8.0.rst b/docs/sphinx/source/whatsnew/v0.8.0.rst
index 2aa1ebfb67..1fd41f232f 100644
--- a/docs/sphinx/source/whatsnew/v0.8.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.8.0.rst
@@ -38,6 +38,10 @@ Enhancements
 * Add :py:func:`pvlib.iam.marion_diffuse` and
   :py:func:`pvlib.iam.marion_integrate` to calculate IAM values for
   diffuse irradiance. (:pull:`984`)
+* Add :py:func:`pvlib.shading.sky_diffuse_passias`,
+  :py:func:`pvlib.shading.masking_angle_passias`, and
+  :py:func:`pvlib.shading.masking_angle` to model diffuse shading loss.
+  (:pull:`1017`)
 * Add :py:func:`pvlib.inverter.fit_sandia` that fits the Sandia inverter model
   to a set of inverter efficiency curves. (:pull:`1011`)
 
@@ -75,6 +79,7 @@ Documentation
 * Add a transposition gain example to the gallery.  (:pull:`979`)
 * Add a gallery example of calculating diffuse IAM using
   :py:func:`pvlib.iam.marion_diffuse`. (:pull:`984`)
+* Add a gallery example of modeling diffuse shading loss. (:pull:`1017`)
 * Add minigalleries to API reference pages. (:pull:`991`)
 
 Requirements
diff --git a/pvlib/shading.py b/pvlib/shading.py
new file mode 100644
index 0000000000..9479eb1739
--- /dev/null
+++ b/pvlib/shading.py
@@ -0,0 +1,193 @@
+"""
+The ``shading`` module contains functions that model module shading and the
+associated effects on PV module output
+"""
+
+import numpy as np
+import pandas as pd
+from pvlib.tools import sind, cosd
+
+
+def masking_angle(surface_tilt, gcr, slant_height):
+    """
+    The elevation angle below which diffuse irradiance is blocked.
+
+    The ``height`` parameter determines how far up the module's surface to
+    evaluate the masking angle.  The lower the point, the steeper the masking
+    angle [1]_.  SAM uses a "worst-case" approach where the masking angle
+    is calculated for the bottom of the array (i.e. ``slant_height=0``) [2]_.
+
+    Parameters
+    ----------
+    surface_tilt : numeric
+        Panel tilt from horizontal [degrees].
+
+    gcr : float
+        The ground coverage ratio of the array [unitless].
+
+    slant_height : numeric
+        The distance up the module's slant height to evaluate the masking
+        angle, as a fraction [0-1] of the module slant height [unitless].
+
+    Returns
+    -------
+    mask_angle : numeric
+        Angle from horizontal where diffuse light is blocked by the
+        preceding row [degrees].
+
+    See Also
+    --------
+    masking_angle_passias
+    sky_diffuse_passias
+
+    References
+    ----------
+    .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
+       panels", Solar Cells, Volume 11, Pages 281-291.  1984.
+       DOI: 10.1016/0379-6787(84)90017-6
+    .. [2] Gilman, P. et al., (2018). "SAM Photovoltaic Model Technical
+       Reference Update", NREL Technical Report NREL/TP-6A20-67399.
+       Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
+    """
+    # The original equation (8 in [1]) requires pitch and collector width,
+    # but it's easy to non-dimensionalize it to make it a function of GCR
+    # by factoring out B from the argument to arctan.
+    numerator = (1 - slant_height) * sind(surface_tilt)
+    denominator = 1/gcr - (1 - slant_height) * cosd(surface_tilt)
+    phi = np.arctan(numerator / denominator)
+    return np.degrees(phi)
+
+
+def masking_angle_passias(surface_tilt, gcr):
+    r"""
+    The average masking angle over the slant height of a row.
+
+    The masking angle is the angle from horizontal where the sky dome is
+    blocked by the row in front. The masking angle is larger near the lower
+    edge of a row than near the upper edge. This function calculates the
+    average masking angle as described in [1]_.
+
+    Parameters
+    ----------
+    surface_tilt : numeric
+        Panel tilt from horizontal [degrees].
+
+    gcr : float
+        The ground coverage ratio of the array [unitless].
+
+    Returns
+    ----------
+    mask_angle : numeric
+        Average angle from horizontal where diffuse light is blocked by the
+        preceding row [degrees].
+
+    See Also
+    --------
+    masking_angle
+    sky_diffuse_passias
+
+    Notes
+    -----
+    The pvlib-python authors believe that Eqn. 9 in [1]_ is incorrect.
+    Here we use an independent equation.  First, Eqn. 8 is non-dimensionalized
+    (recasting in terms of GCR):
+
+    .. math::
+
+        \psi(z') = \arctan \left [
+            \frac{(1 - z') \sin \beta}
+                 {\mathrm{GCR}^{-1} + (z' - 1) \cos \beta}
+        \right ]
+
+    Where :math:`GCR = B/C` and :math:`z' = z/B`. The average masking angle
+    :math:`\overline{\psi} = \int_0^1 \psi(z') \mathrm{d}z'` is then
+    evaluated symbolically using Maxima (using :math:`X = 1/\mathrm{GCR}`):
+
+    .. code-block:: none
+
+        load(scifac)    /* for the gcfac function */
+        assume(X>0, cos(beta)>0, cos(beta)-X<0);   /* X is 1/GCR */
+        gcfac(integrate(atan((1-z)*sin(beta)/(X+(z-1)*cos(beta))), z, 0, 1))
+
+    This yields the equation implemented by this function:
+
+    .. math::
+
+        \overline{\psi} = \
+            &-\frac{X}{2} \sin\beta \log | 2 X \cos\beta - (X^2 + 1)| \\
+            &+ (X \cos\beta - 1) \arctan \frac{X \cos\beta - 1}{X \sin\beta} \\
+            &+ (1 - X \cos\beta) \arctan \frac{\cos\beta}{\sin\beta} \\
+            &+ X \log X \sin\beta
+
+    The pvlib-python authors have validated this equation against numerical
+    integration of :math:`\overline{\psi} = \int_0^1 \psi(z') \mathrm{d}z'`.
+
+    References
+    ----------
+    .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
+       panels", Solar Cells, Volume 11, Pages 281-291.  1984.
+       DOI: 10.1016/0379-6787(84)90017-6
+    """
+    # wrap it in an array so that division by zero is handled well
+    beta = np.radians(np.array(surface_tilt))
+    sin_b = np.sin(beta)
+    cos_b = np.cos(beta)
+    X = 1/gcr
+
+    with np.errstate(divide='ignore', invalid='ignore'):  # ignore beta=0
+        term1 = -X * sin_b * np.log(np.abs(2 * X * cos_b - (X**2 + 1))) / 2
+        term2 = (X * cos_b - 1) * np.arctan((X * cos_b - 1) / (X * sin_b))
+        term3 = (1 - X * cos_b) * np.arctan(cos_b / sin_b)
+        term4 = X * np.log(X) * sin_b
+
+    psi_avg = term1 + term2 + term3 + term4
+    # when beta=0, divide by zero makes psi_avg NaN.  replace with 0:
+    psi_avg = np.where(np.isfinite(psi_avg), psi_avg, 0)
+
+    if isinstance(surface_tilt, pd.Series):
+        psi_avg = pd.Series(psi_avg, index=surface_tilt.index)
+
+    return np.degrees(psi_avg)
+
+
+def sky_diffuse_passias(masking_angle):
+    r"""
+    The diffuse irradiance loss caused by row-to-row sky diffuse shading.
+
+    Even when the sun is high in the sky, a row's view of the sky dome will
+    be partially blocked by the row in front. This causes a reduction in the
+    diffuse irradiance incident on the module. The reduction depends on the
+    masking angle, the elevation angle from a point on the shaded module to
+    the top of the shading row. In [1]_ the masking angle is calculated as
+    the average across the module height. SAM assumes the "worst-case" loss
+    where the masking angle is calculated for the bottom of the array [2]_.
+
+    This function, as in [1]_, makes the assumption that sky diffuse
+    irradiance is isotropic.
+
+    Parameters
+    ----------
+    masking_angle : numeric
+        The elevation angle below which diffuse irradiance is blocked
+        [degrees].
+
+    Returns
+    -------
+    derate : numeric
+        The fraction [0-1] of blocked sky diffuse irradiance.
+
+    See Also
+    --------
+    masking_angle
+    masking_angle_passias
+
+    References
+    ----------
+    .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
+       panels", Solar Cells, Volume 11, Pages 281-291.  1984.
+       DOI: 10.1016/0379-6787(84)90017-6
+    .. [2] Gilman, P. et al., (2018). "SAM Photovoltaic Model Technical
+       Reference Update", NREL Technical Report NREL/TP-6A20-67399.
+       Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
+    """
+    return 1 - cosd(masking_angle/2)**2
diff --git a/pvlib/tests/test_shading.py b/pvlib/tests/test_shading.py
new file mode 100644
index 0000000000..8a9fd46a69
--- /dev/null
+++ b/pvlib/tests/test_shading.py
@@ -0,0 +1,71 @@
+import numpy as np
+import pandas as pd
+
+from pandas.testing import assert_series_equal
+import pytest
+
+from pvlib import shading
+
+
+@pytest.fixture
+def surface_tilt():
+    idx = pd.date_range('2019-01-01', freq='h', periods=3)
+    return pd.Series([0, 20, 90], index=idx)
+
+
+@pytest.fixture
+def masking_angle(surface_tilt):
+    # masking angles for the surface_tilt fixture,
+    # assuming GCR=0.5 and height=0.25
+    return pd.Series([0.0, 11.20223712, 20.55604522], index=surface_tilt.index)
+
+
+@pytest.fixture
+def average_masking_angle(surface_tilt):
+    # average masking angles for the surface_tilt fixture, assuming GCR=0.5
+    return pd.Series([0.0, 7.20980655, 13.779867461], index=surface_tilt.index)
+
+
+@pytest.fixture
+def shading_loss(surface_tilt):
+    # diffuse shading loss values for the average_masking_angle fixture
+    return pd.Series([0, 0.00395338, 0.01439098], index=surface_tilt.index)
+
+
+def test_masking_angle_series(surface_tilt, masking_angle):
+    # series inputs and outputs
+    masking_angle_actual = shading.masking_angle(surface_tilt, 0.5, 0.25)
+    assert_series_equal(masking_angle_actual, masking_angle)
+
+
+def test_masking_angle_scalar(surface_tilt, masking_angle):
+    # scalar inputs and outputs, including zero
+    for tilt, angle in zip(surface_tilt, masking_angle):
+        masking_angle_actual = shading.masking_angle(tilt, 0.5, 0.25)
+        assert np.isclose(masking_angle_actual, angle)
+
+
+def test_masking_angle_passias_series(surface_tilt, average_masking_angle):
+    # pandas series inputs and outputs
+    masking_angle_actual = shading.masking_angle_passias(surface_tilt, 0.5)
+    assert_series_equal(masking_angle_actual, average_masking_angle)
+
+
+def test_masking_angle_passias_scalar(surface_tilt, average_masking_angle):
+    # scalar inputs and outputs, including zero
+    for tilt, angle in zip(surface_tilt, average_masking_angle):
+        masking_angle_actual = shading.masking_angle_passias(tilt, 0.5)
+        assert np.isclose(masking_angle_actual, angle)
+
+
+def test_sky_diffuse_passias_series(average_masking_angle, shading_loss):
+    # pandas series inputs and outputs
+    actual_loss = shading.sky_diffuse_passias(average_masking_angle)
+    assert_series_equal(shading_loss, actual_loss)
+
+
+def test_sky_diffuse_passias_scalar(average_masking_angle, shading_loss):
+    # scalar inputs and outputs
+    for angle, loss in zip(average_masking_angle, shading_loss):
+        actual_loss = shading.sky_diffuse_passias(angle)
+        assert np.isclose(loss, actual_loss)