diff --git a/docs/examples/plot_partial_module_shading_simple.py b/docs/examples/plot_partial_module_shading_simple.py
new file mode 100644
index 0000000000..ec96d9ca45
--- /dev/null
+++ b/docs/examples/plot_partial_module_shading_simple.py
@@ -0,0 +1,344 @@
+"""
+Calculating power loss from partial module shading
+==================================================
+
+Example of modeling cell-to-cell mismatch loss from partial module shading.
+"""
+
+# %%
+# Even though the PV cell is the primary power generation unit, PV modeling is
+# often done at the module level for simplicity because module-level parameters
+# are much more available and it significantly reduces the computational scope
+# of the simulation.  However, module-level simulations are too coarse to be
+# able to model effects like cell to cell mismatch or partial shading.  This
+# example calculates cell-level IV curves and combines them to reconstruct
+# the module-level IV curve.  It uses this approach to find the maximum power
+# under various shading and irradiance conditions.
+#
+# The primary functions used here are:
+#
+# - :py:meth:`pvlib.pvsystem.calcparams_desoto` to estimate the single
+#   diode equation parameters at some specified operating conditions.
+# - :py:meth:`pvlib.singlediode.bishop88` to calculate the full cell IV curve,
+#   including the reverse bias region.
+#
+# .. note::
+#
+#     This example requires the reverse bias functionality added in pvlib 0.7.2
+#
+# .. warning::
+#
+#     Modeling partial module shading is complicated and depends significantly
+#     on the module's electrical topology.  This example makes some simplifying
+#     assumptions that are not generally applicable.  For instance, it assumes
+#     that shading only applies to beam irradiance (*i.e.* all cells receive
+#     the same amount of diffuse irradiance) and cell temperature is uniform
+#     and not affected by cell-level irradiance variation.
+
+from pvlib import pvsystem, singlediode
+import pandas as pd
+import numpy as np
+from scipy.interpolate import interp1d
+import matplotlib.pyplot as plt
+
+from scipy.constants import e as qe, k as kB
+
+# For simplicity, use cell temperature of 25C for all calculations.
+# kB is J/K, qe is C=J/V
+# kB * T / qe -> V
+Vth = kB * (273.15+25) / qe
+
+cell_parameters = {
+    'I_L_ref': 8.24,
+    'I_o_ref': 2.36e-9,
+    'a_ref': 1.3*Vth,
+    'R_sh_ref': 1000,
+    'R_s': 0.00181,
+    'alpha_sc': 0.0042,
+    'breakdown_factor': 2e-3,
+    'breakdown_exp': 3,
+    'breakdown_voltage': -15,
+}
+
+# %%
+# Simulating a cell IV curve
+# --------------------------
+#
+# First, calculate IV curves for individual cells.  The process is as follows:
+#
+# 1) Given a set of cell parameters at reference conditions and the operating
+#    conditions of interest (irradiance and temperature), use a single-diode
+#    model to calculate the single diode equation parameters for the cell at
+#    the operating conditions.  Here we use the De Soto model via
+#    :py:func:`pvlib.pvsystem.calcparams_desoto`.
+# 2) The single diode equation cannot be solved analytically, so pvlib has
+#    implemented a couple methods of solving it for us.  However, currently
+#    only the Bishop '88 method (:py:func:`pvlib.singlediode.bishop88`) has
+#    the ability to model the reverse bias characteristic in addition to the
+#    forward characteristic.  Depending on the nature of the shadow, it is
+#    sometimes necessary to model the reverse bias portion of the IV curve,
+#    so we use the Bishop '88 method here.  This gives us a set of (V, I)
+#    points on the cell's IV curve.
+
+
+def simulate_full_curve(parameters, Geff, Tcell, ivcurve_pnts=1000):
+    """
+    Use De Soto and Bishop to simulate a full IV curve with both
+    forward and reverse bias regions.
+    """
+    # adjust the reference parameters according to the operating
+    # conditions using the De Soto model:
+    sde_args = pvsystem.calcparams_desoto(
+        Geff,
+        Tcell,
+        alpha_sc=parameters['alpha_sc'],
+        a_ref=parameters['a_ref'],
+        I_L_ref=parameters['I_L_ref'],
+        I_o_ref=parameters['I_o_ref'],
+        R_sh_ref=parameters['R_sh_ref'],
+        R_s=parameters['R_s'],
+    )
+    # sde_args has values:
+    # (photocurrent, saturation_current, resistance_series,
+    # resistance_shunt, nNsVth)
+
+    # Use Bishop's method to calculate points on the IV curve with V ranging
+    # from the reverse breakdown voltage to open circuit
+    kwargs = {
+        'breakdown_factor': parameters['breakdown_factor'],
+        'breakdown_exp': parameters['breakdown_exp'],
+        'breakdown_voltage': parameters['breakdown_voltage'],
+    }
+    v_oc = singlediode.bishop88_v_from_i(
+        0.0, *sde_args, **kwargs
+    )
+    # ideally would use some intelligent log-spacing to concentrate points
+    # around the forward- and reverse-bias knees, but this is good enough:
+    vd = np.linspace(0.99*kwargs['breakdown_voltage'], v_oc, ivcurve_pnts)
+
+    ivcurve_i, ivcurve_v, _ = singlediode.bishop88(vd, *sde_args, **kwargs)
+    return pd.DataFrame({
+        'i': ivcurve_i,
+        'v': ivcurve_v,
+    })
+
+
+# %%
+# Now that we can calculate cell-level IV curves, let's compare a
+# fully-illuminated cell's curve to a shaded cell's curve.  Note that shading
+# typically does not reduce a cell's illumination to zero -- tree shading and
+# row-to-row shading block the beam portion of irradiance but leave the diffuse
+# portion largely intact.  In this example plot, we choose :math:`200 W/m^2`
+# as the amount of irradiance received by a shaded cell.
+
+def plot_curves(dfs, labels, title):
+    """plot the forward- and reverse-bias portions of an IV curve"""
+    fig, axes = plt.subplots(1, 2, sharey=True, figsize=(5, 3))
+    for df, label in zip(dfs, labels):
+        df.plot('v', 'i', label=label, ax=axes[0])
+        df.plot('v', 'i', label=label, ax=axes[1])
+        axes[0].set_xlim(right=0)
+        axes[0].set_ylim([0, 25])
+        axes[1].set_xlim([0, df['v'].max()*1.5])
+    axes[0].set_ylabel('current [A]')
+    axes[0].set_xlabel('voltage [V]')
+    axes[1].set_xlabel('voltage [V]')
+    fig.suptitle(title)
+    fig.tight_layout()
+    return axes
+
+
+cell_curve_full_sun = simulate_full_curve(cell_parameters, Geff=1000, Tcell=25)
+cell_curve_shaded = simulate_full_curve(cell_parameters, Geff=200, Tcell=25)
+ax = plot_curves([cell_curve_full_sun, cell_curve_shaded],
+                 labels=['Full Sun', 'Shaded'],
+                 title='Cell-level reverse- and forward-biased IV curves')
+
+# %%
+# This figure shows how a cell's current decreases roughly in proportion to
+# the irradiance reduction from shading, but voltage changes much less.
+# At the cell level, the effect of shading is essentially to shift the I-V
+# curve down to lower currents rather than change the curve's shape.
+#
+# Note that the forward and reverse curves are plotted separately to
+# accommodate the different voltage scales involved -- a normal crystalline
+# silicon cell reaches only ~0.6V in forward bias, but can get to -10 to -20V
+# in reverse bias.
+#
+# Combining cell IV curves to create a module IV curve
+# ----------------------------------------------------
+#
+# To combine the individual cell IV curves and form a module's IV curve,
+# the cells in each substring must be added in series.  The substrings are
+# in series as well, but with parallel bypass diodes to protect from reverse
+# bias voltages.  To add in series, the voltages for a given current are
+# added.  However, because each cell's curve is discretized and the currents
+# might not line up, we align each curve to a common set of current values
+# with interpolation.
+
+
+def interpolate(df, i):
+    """convenience wrapper around scipy.interpolate.interp1d"""
+    f_interp = interp1d(np.flipud(df['i']), np.flipud(df['v']), kind='linear',
+                        fill_value='extrapolate')
+    return f_interp(i)
+
+
+def combine_series(dfs):
+    """
+    Combine IV curves in series by aligning currents and summing voltages.
+    The current range is based on the first curve's current range.
+    """
+    df1 = dfs[0]
+    imin = df1['i'].min()
+    imax = df1['i'].max()
+    i = np.linspace(imin, imax, 1000)
+    v = 0
+    for df2 in dfs:
+        v_cell = interpolate(df2, i)
+        v += v_cell
+    return pd.DataFrame({'i': i, 'v': v})
+
+
+# %%
+# Rather than simulate all 72 cells in the module, we'll assume that there
+# are only three types of cells (fully illuminated, fully shaded, and
+# partially shaded), and within each type all cells behave identically.  This
+# means that simulating one cell from each type (for three cell simulations
+# total) is sufficient to model the module as a whole.
+#
+# This function also models the effect of bypass diodes in parallel with each
+# substring.  Bypass diodes are normally inactive but conduct when substring
+# voltage becomes sufficiently negative, presumably due to the substring
+# entering reverse bias from mismatch between substrings.  In that case the
+# substring's voltage is clamped to the diode's trigger voltage (assumed to
+# be 0.5V here).
+
+def simulate_module(cell_parameters, poa_direct, poa_diffuse, Tcell,
+                    shaded_fraction, cells_per_string=24, strings=3):
+    """
+    Simulate the IV curve for a partially shaded module.
+    The shade is assumed to be coming up from the bottom of the module when in
+    portrait orientation, so it affects all substrings equally.
+    For simplicity, cell temperature is assumed to be uniform across the
+    module, regardless of variation in cell-level illumination.
+    Substrings are assumed to be "down and back", so the number of cells per
+    string is divided between two columns of cells.
+    """
+    # find the number of cells per column that are in full shadow
+    nrow = cells_per_string // 2
+    nrow_full_shade = int(shaded_fraction * nrow)
+    # find the fraction of shade in the border row
+    partial_shade_fraction = 1 - (shaded_fraction * nrow - nrow_full_shade)
+
+    df_lit = simulate_full_curve(
+        cell_parameters,
+        poa_diffuse + poa_direct,
+        Tcell)
+    df_partial = simulate_full_curve(
+        cell_parameters,
+        poa_diffuse + partial_shade_fraction * poa_direct,
+        Tcell)
+    df_shaded = simulate_full_curve(
+        cell_parameters,
+        poa_diffuse,
+        Tcell)
+    # build a list of IV curves for a single column of cells (half a substring)
+    include_partial_cell = (shaded_fraction < 1)
+    half_substring_curves = (
+        [df_lit] * (nrow - nrow_full_shade - 1)
+        + ([df_partial] if include_partial_cell else [])  # noqa: W503
+        + [df_shaded] * nrow_full_shade  # noqa: W503
+    )
+    substring_curve = combine_series(half_substring_curves)
+    substring_curve['v'] *= 2  # turn half strings into whole strings
+    # bypass diode:
+    substring_curve['v'] = substring_curve['v'].clip(lower=-0.5)
+    # no need to interpolate since we're just scaling voltage directly:
+    substring_curve['v'] *= strings
+    return substring_curve
+
+# %%
+# Now let's see how shade affects the IV curves at the module level.  For this
+# example, the bottom 10% of the module is shaded.  Assuming 12 cells per
+# column, that means one row of cells is fully shaded and another row is
+# partially shaded.  Even though only 10% of the module is shaded, the
+# maximum power is decreased by roughly 80%!
+#
+# Note the effect of the bypass diodes.  Without bypass diodes, operating the
+# shaded module at the same current as the fully illuminated module would
+# create a reverse-bias voltage of several hundred volts!  However, the diodes
+# prevent the reverse voltage from exceeding 1.5V (three diodes at 0.5V each).
+
+
+kwargs = {
+    'cell_parameters': cell_parameters,
+    'poa_direct': 800,
+    'poa_diffuse': 200,
+    'Tcell': 25
+}
+module_curve_full_sun = simulate_module(shaded_fraction=0, **kwargs)
+module_curve_shaded = simulate_module(shaded_fraction=0.1, **kwargs)
+ax = plot_curves([module_curve_full_sun, module_curve_shaded],
+                 labels=['Full Sun', 'Shaded'],
+                 title='Module-level reverse- and forward-biased IV curves')
+
+# %%
+# Calculating shading loss across shading scenarios
+# -------------------------------------------------
+#
+# Clearly the module-level IV-curve is strongly affected by partial shading.
+# This heatmap shows the module maximum power under a range of partial shade
+# conditions, where "diffuse fraction" refers to the ratio
+# :math:`poa_{diffuse} / poa_{global}` and "shaded fraction" refers to the
+# fraction of the module that receives only diffuse irradiance.
+
+
+def find_pmp(df):
+    """simple function to find Pmp on an IV curve"""
+    return df.product(axis=1).max()
+
+
+# find Pmp under different shading conditions
+data = []
+for diffuse_fraction in np.linspace(0, 1, 11):
+    for shaded_fraction in np.linspace(0, 1, 51):
+
+        df = simulate_module(cell_parameters,
+                             poa_direct=(1-diffuse_fraction)*1000,
+                             poa_diffuse=diffuse_fraction*1000,
+                             Tcell=25,
+                             shaded_fraction=shaded_fraction)
+        data.append({
+            'fd': diffuse_fraction,
+            'fs': shaded_fraction,
+            'pmp': find_pmp(df)
+        })
+
+results = pd.DataFrame(data)
+results['pmp'] /= results['pmp'].max()  # normalize power to 0-1
+results_pivot = results.pivot('fd', 'fs', 'pmp')
+plt.figure()
+plt.imshow(results_pivot, origin='lower', aspect='auto')
+plt.xlabel('shaded fraction')
+plt.ylabel('diffuse fraction')
+xlabels = ["{:0.02f}".format(fs) for fs in results_pivot.columns[::5]]
+ylabels = ["{:0.02f}".format(fd) for fd in results_pivot.index]
+plt.xticks(range(0, 5*len(xlabels), 5), xlabels)
+plt.yticks(range(0, len(ylabels)), ylabels)
+plt.title('Module P_mp across shading conditions')
+plt.colorbar()
+plt.show()
+# use this figure as the thumbnail:
+# sphinx_gallery_thumbnail_number = 3
+
+# %%
+# The heatmap makes a few things evident:
+#
+# - When diffuse fraction is equal to 1, there is no beam irradiance to lose,
+#   so shading has no effect on production.
+# - When shaded fraction is equal to 0, no irradiance is blocked, so module
+#   output does not change with the diffuse fraction.
+# - Under sunny conditions (diffuse fraction < 0.5), module output is
+#   significantly reduced after just the first cell is shaded
+#   (1/12 = ~8% shaded fraction).