diff --git a/doc/docs/Python_Tutorials/Mode_Decomposition.md b/doc/docs/Python_Tutorials/Mode_Decomposition.md
index e12340c1f..7f1d86172 100644
--- a/doc/docs/Python_Tutorials/Mode_Decomposition.md
+++ b/doc/docs/Python_Tutorials/Mode_Decomposition.md
@@ -143,6 +143,264 @@ The reflectance values computed using the two methods are nearly identical. For
 
 In the reflected-flux calculation, we apply our usual trick of first performing a reference simulation with just the incident field and then subtracting that from our taper simulation with `load_minus_flux_data`, so that what is left over is the reflected fields (from which we obtain the reflected flux).  In *principle*, this trick would not be required for the mode-decomposition method, because the reflected mode is orthogonal to the forward mode and so the decomposition will separate the forward and reflected coefficients automatically.   However, this is only true in the limit of infinite resolution — for a *finite* resolution, the reflected mode used for the mode coefficient calculation (calculated via MPB) is not exactly orthogonal to the forward mode propagating in Meep (whose discretization scheme is different from that of MPB).   In consequence, if you did not subtract the fields of the reference simulation, the mode-coefficient could only calculate the reflected power down to a "noise floor" set by the discretization error.  With the subtraction, in contrast, you can compute much smaller reflections (limited by the floating-point precision).
 
+Phase of a Total Internal Reflected Mode
+----------------------------------------
+
+For certain applications, such as [physically based ray tracing](https://pbr-book.org/), the phase of the reflection/transmission coefficient of a surface interface (flat or textured) is an important parameter used for modeling coherent effects (i.e., constructive or destructive interference from a collection of intersecting rays scattered off the surface). This tutorial demonstrates how to use mode decomposition to compute the phase of the reflection coefficient of a [total internal reflected](https://en.wikipedia.org/wiki/Total_internal_reflection) (TIR) mode. The simulated results are verified using the [Fresnel equations](https://en.wikipedia.org/wiki/Total_internal_reflection#Phase_shifts).
+
+The example involves a 2d simulation of a flat interface of two lossless media with $n_1=1.5$ and $n_2=1.0$. The [critical angle](https://en.wikipedia.org/wiki/Total_internal_reflection#Critical_angle) for this interface is $\theta_c = 41.8^\circ$. The source is an incident planewave with linear polarization ($S$ or $P$). The 2d cell is periodic in $y$ with PML in the $x$ direction. The cell size does not affect the results and is therefore arbitrary. A schematic of the simulation layout is shown below.
+
+The key item to consider in these types of calculations is the *location of the mode monitor relative to the interface*. The mode monitor is a line extending the entire length of the cell in the $y$ direction. In order to compare the result with the Fresnel equations, the phase of the reflection coefficient must be computed *exactly* at the interface. However, it is problematic to measure the reflection coefficient exactly at the interface because the amplitude of the left-going wave drops discontinuously to zero across the interface. For this reason, the mode monitor must be positioned away from the interface (somewhere within the $n_1$ medium) and the measured phase must be corrected in post processing to account for this offset.
+
+![](../images/refl_coeff_flat_interface.png#center)
+
+The calculation of the reflection coefficient requires two separate runs: (1) a normalization run involving just the source medium ($n_1$) to obtain the incident fields, and (2) the main run including the interface whereby the incident fields from (1) are first subtracted from the monitor to obtain only the reflected fields. The mode coefficient in (2) divided by (1) is, by definition, the reflection coefficient.
+
+The phase of the reflection coefficient needs to be corrected for the offset of the mode monitor relative to the interface — labeled $L$ in the schematic above — using the formula $\exp(i 2\pi k_x 2L)$, where $k_x$ is the $x$ component of the planewave's wavevector (`k` in the script). The $2\pi$ factor is necessary because `k` does *not* include this factor (per convention in Meep). The factor 2 in front of the $L$ is necessary because the phase needs to be corrected for the monitors in the normalization and main runs separately. The correction factor is just the phase accumulated as the wave propagates in the surface-normal direction for a distance $L$ in a medium with index $n_1$.
+
+With this setup, we measure the phase of the reflection coefficient for two different source configurations (polarization and angle) and compare the result with the Fresnel equations. The location of the mode monitor is also varied in the two test configurations. Results are shown in the two tables below. There is good agreement between the simulation and theory. (As additional validation, we note that the magnitude of the reflection coefficient of a TIR mode must be one which is indeed the case.)
+
+| $S$ pol., $\theta$ = 54.3° | reflection coefficient | phase (radians) |
+|:--------------------------:|:----------------------:|:---------------:|
+|            Meep            |    0.23415-0.96597j    |     -1.33299    |
+|           Fresnel          |    0.22587-0.97416j    |     -1.34296    |
+
+
+| $P$ pol., $\theta$ = 48.5° | reflection coefficient | phase (radians) |
+|:--------------------------:|:----------------------:|:---------------:|
+|            Meep            |    0.11923+0.98495j    |      1.45033    |
+|           Fresnel          |    0.14645+0.98922j    |      1.42382    |
+
+The simulation script is in [examples/mode_coeff_phase.py](https://github.com/NanoComp/meep/blob/master/python/examples/mode_coeff_phase.py).
+
+```py
+import cmath
+from enum import Enum
+import math
+
+import meep as mp
+import numpy as np
+
+Polarization = Enum("Polarization", "S P")
+
+# refractive indices of materials 1 and 2
+n1 = 1.5
+n2 = 1.0
+
+
+def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
+    """Computes the reflection coefficient of a TIR mode using mode
+       decomposition.
+
+    Args:
+        pol: polarization of the incident planewave (S or P).
+        theta: angle of the incident planewave (radians).
+        L: position of the mode monitor relative to the flat interface. 0 is
+            interface position.
+    """
+    resolution = 50.0
+
+    # cell size is arbitrary
+    sx = 7.0
+    sy = 3.0
+    dpml = 2.0
+    cell_size = mp.Vector3(sx + 2 * dpml, sy, 0)
+    pml_layers = [mp.PML(dpml, direction=mp.X)]
+
+    fcen = 1.0  # center frequency
+    df = 0.1 * fcen
+
+    # k (in source medium) with correct length
+    # plane of incidence is xy
+    k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta)
+
+    # planewave amplitude function (for source)
+    def pw_amp(k, x0):
+        def _pw_amp(x):
+            return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
+
+        return _pw_amp
+
+    src_pt = mp.Vector3(-0.5 * sx, 0, 0)
+
+    if pol.name == "S":
+        src_cmpt = mp.Ez
+        eig_parity = mp.ODD_Z
+    elif pol.name == "P":
+        src_cmpt = mp.Hz
+        eig_parity = mp.EVEN_Z
+    else:
+        raise ValueError("pol must be S or P, only.")
+
+    sources = [
+        mp.Source(
+            mp.GaussianSource(fcen, fwidth=df),
+            component=src_cmpt,
+            center=src_pt,
+            size=mp.Vector3(0, cell_size.y, 0),
+            amp_func=pw_amp(k, src_pt),
+        ),
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        default_material=mp.Medium(index=n1),
+        boundary_layers=pml_layers,
+        k_point=k,
+        sources=sources,
+    )
+
+    # DFT monitor for incident fields
+    mode_mon = sim.add_mode_monitor(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(-L, 0, 0),
+            size=mp.Vector3(0, cell_size.y, 0),
+        ),
+    )
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            50,
+            src_cmpt,
+            mp.Vector3(-L, 0, 0),
+            1e-6,
+        ),
+    )
+
+    res = sim.get_eigenmode_coefficients(
+        mode_mon,
+        bands=[1],
+        eig_parity=eig_parity,
+        kpoint_func=lambda *not_used: k,
+        direction=mp.NO_DIRECTION,
+    )
+
+    input_mode_coeff = res.alpha[0, 0, 0]
+    input_flux_data = sim.get_flux_data(mode_mon)
+
+    sim.reset_meep()
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=n1),
+            center=mp.Vector3(-0.25 * (sx + 2 * dpml), 0, 0),
+            size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf),
+        ),
+        mp.Block(
+            material=mp.Medium(index=n2),
+            center=mp.Vector3(0.25 * (sx + 2 * dpml), 0, 0),
+            size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf),
+        ),
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        boundary_layers=pml_layers,
+        k_point=k,
+        sources=sources,
+        geometry=geometry,
+    )
+
+    # DFT monitor for reflected fields
+    mode_mon = sim.add_mode_monitor(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(-L, 0, 0),
+            size=mp.Vector3(0, cell_size.y, 0),
+        ),
+    )
+
+    sim.load_minus_flux_data(mode_mon, input_flux_data)
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            50,
+            mp.Ez,
+            mp.Vector3(-L, 0, 0),
+            1e-6,
+        ),
+    )
+
+    res = sim.get_eigenmode_coefficients(
+        mode_mon,
+        bands=[1],
+        eig_parity=eig_parity,
+        kpoint_func=lambda *not_used: k,
+        direction=mp.NO_DIRECTION,
+    )
+
+    # mode coefficient of reflected planewave
+    refl_mode_coeff = res.alpha[0, 0, 1]
+
+    # reflection coefficient
+    refl_coeff = refl_mode_coeff / input_mode_coeff
+
+    # apply phase correction factor
+    refl_coeff /= cmath.exp(1j * k.x * 2 * math.pi * 2 * L)
+
+    return refl_coeff
+
+
+def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex:
+    """Computes the reflection coefficient of a TIR mode using the Fresnel
+       equations.
+
+    Args:
+        pol: polarization of the incident planewave (S or P).
+        theta: angle of the incident planewave (degrees).
+    """
+    if pol.name == "S":
+        refl_coeff = (
+            math.cos(theta) - ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+        ) / (math.cos(theta) + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5)
+    else:
+        refl_coeff = (
+            -((n2 / n1) ** 2) * math.cos(theta)
+            + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+        ) / (
+            (n2 / n1) ** 2 * math.cos(theta)
+            + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+        )
+
+    return refl_coeff
+
+
+if __name__ == "__main__":
+    thetas = [54.3, 48.5] # angle of incident planewave (degrees)
+    Ls = [0.4, 1.2] # position of mode monitor relative to flat interface
+    pols = [Polarization.S, Polarization.P] # polarization of incident planewave
+
+    for pol, theta, L in zip(pols, thetas, Ls):
+        theta_rad = np.radians(theta)
+        rc_m = refl_coeff_meep(pol, theta_rad, L)
+        rc_f = refl_coeff_Fresnel(pol, theta_rad)
+
+        rc_m_str = f"{rc_m.real:.5f}{rc_m.imag:+.5f}j"
+        rc_f_str = f"{rc_f.real:.5f}{rc_f.imag:+.5f}j"
+        print(f"refl-coeff:, {pol.name}, {theta}, {rc_m_str} (Meep), "
+              f"{rc_f_str} (Fresnel)")
+
+        mag_m = abs(rc_m)
+        mag_f = abs(rc_f)
+        err_mag = abs(mag_m - mag_f) / mag_f
+        print(f"magnitude:, {mag_m:.5f} (Meep), {mag_f:.5f} (Fresnel), "
+              f"{err_mag:.5f} (error)")
+
+        phase_m = cmath.phase(rc_m)
+        phase_f = cmath.phase(rc_f)
+        err_phase = abs(phase_m - phase_f) / abs(phase_f)
+        print(f"phase:, {phase_m:.5f} (Meep), {phase_f:.5f} (Fresnel), "
+              f"{err_phase:.5f} (error)")
+```
+
+
 Diffraction Spectrum of a Binary Grating
 ----------------------------------------
 
diff --git a/doc/docs/images/refl_coeff_flat_interface.png b/doc/docs/images/refl_coeff_flat_interface.png
new file mode 100644
index 000000000..6a248f197
Binary files /dev/null and b/doc/docs/images/refl_coeff_flat_interface.png differ
diff --git a/python/examples/mode_coeff_phase.py b/python/examples/mode_coeff_phase.py
new file mode 100644
index 000000000..6f78f4eed
--- /dev/null
+++ b/python/examples/mode_coeff_phase.py
@@ -0,0 +1,239 @@
+"""
+Verifies that the magnitude and phase of the reflection coefficient of a
+total internal reflected (TIR) mode of a flat interface of two lossless
+materials given an incident planewave at oblique incidence computed using
+the mode-decomposition feature matches the Fresnel equations.
+
+ref: https://meep.readthedocs.io/en/latest/Python_Tutorials/Mode_Decomposition/#phase-of-a-total-internal-reflected-mode
+"""
+
+import cmath
+from enum import Enum
+import math
+
+import meep as mp
+import numpy as np
+
+Polarization = Enum("Polarization", "S P")
+
+# refractive indices of materials 1 and 2
+n1 = 1.5
+n2 = 1.0
+
+
+def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
+    """Computes the reflection coefficient of a TIR mode using mode
+       decomposition.
+
+    Args:
+        pol: polarization of the incident planewave (S or P).
+        theta: angle of the incident planewave (radians).
+        L: position of the mode monitor relative to the flat interface. 0 is
+            interface position.
+    """
+    resolution = 50.0
+
+    # cell size is arbitrary
+    sx = 7.0
+    sy = 3.0
+    dpml = 2.0
+    cell_size = mp.Vector3(sx + 2 * dpml, sy, 0)
+    pml_layers = [mp.PML(dpml, direction=mp.X)]
+
+    fcen = 1.0  # center frequency
+    df = 0.1 * fcen
+
+    # k (in source medium) with correct length
+    # plane of incidence is xy
+    k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta)
+
+    # planewave amplitude function (for source)
+    def pw_amp(k, x0):
+        def _pw_amp(x):
+            return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
+
+        return _pw_amp
+
+    src_pt = mp.Vector3(-0.5 * sx, 0, 0)
+
+    if pol.name == "S":
+        src_cmpt = mp.Ez
+        eig_parity = mp.ODD_Z
+    elif pol.name == "P":
+        src_cmpt = mp.Hz
+        eig_parity = mp.EVEN_Z
+    else:
+        raise ValueError("pol must be S or P, only.")
+
+    sources = [
+        mp.Source(
+            mp.GaussianSource(fcen, fwidth=df),
+            component=src_cmpt,
+            center=src_pt,
+            size=mp.Vector3(0, cell_size.y, 0),
+            amp_func=pw_amp(k, src_pt),
+        ),
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        default_material=mp.Medium(index=n1),
+        boundary_layers=pml_layers,
+        k_point=k,
+        sources=sources,
+    )
+
+    # DFT monitor for incident fields
+    mode_mon = sim.add_mode_monitor(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(-L, 0, 0),
+            size=mp.Vector3(0, cell_size.y, 0),
+        ),
+    )
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            50,
+            src_cmpt,
+            mp.Vector3(-L, 0, 0),
+            1e-6,
+        ),
+    )
+
+    res = sim.get_eigenmode_coefficients(
+        mode_mon,
+        bands=[1],
+        eig_parity=eig_parity,
+        kpoint_func=lambda *not_used: k,
+        direction=mp.NO_DIRECTION,
+    )
+
+    input_mode_coeff = res.alpha[0, 0, 0]
+    input_flux_data = sim.get_flux_data(mode_mon)
+
+    sim.reset_meep()
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=n1),
+            center=mp.Vector3(-0.25 * (sx + 2 * dpml), 0, 0),
+            size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf),
+        ),
+        mp.Block(
+            material=mp.Medium(index=n2),
+            center=mp.Vector3(0.25 * (sx + 2 * dpml), 0, 0),
+            size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf),
+        ),
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        boundary_layers=pml_layers,
+        k_point=k,
+        sources=sources,
+        geometry=geometry,
+    )
+
+    # DFT monitor for reflected fields
+    mode_mon = sim.add_mode_monitor(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(-L, 0, 0),
+            size=mp.Vector3(0, cell_size.y, 0),
+        ),
+    )
+
+    sim.load_minus_flux_data(mode_mon, input_flux_data)
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            50,
+            mp.Ez,
+            mp.Vector3(-L, 0, 0),
+            1e-6,
+        ),
+    )
+
+    res = sim.get_eigenmode_coefficients(
+        mode_mon,
+        bands=[1],
+        eig_parity=eig_parity,
+        kpoint_func=lambda *not_used: k,
+        direction=mp.NO_DIRECTION,
+    )
+
+    # mode coefficient of reflected planewave
+    refl_mode_coeff = res.alpha[0, 0, 1]
+
+    # reflection coefficient
+    refl_coeff = refl_mode_coeff / input_mode_coeff
+
+    # apply phase correction factor
+    refl_coeff /= cmath.exp(1j * k.x * 2 * math.pi * 2 * L)
+
+    return refl_coeff
+
+
+def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex:
+    """Computes the reflection coefficient of a TIR mode using the Fresnel
+       equations.
+
+    Args:
+        pol: polarization of the incident planewave (S or P).
+        theta: angle of the incident planewave (degrees).
+    """
+    if pol.name == "S":
+        refl_coeff = (
+            math.cos(theta) - ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+        ) / (math.cos(theta) + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5)
+    else:
+        refl_coeff = (
+            -((n2 / n1) ** 2) * math.cos(theta)
+            + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+        ) / (
+            (n2 / n1) ** 2 * math.cos(theta)
+            + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+        )
+
+    return refl_coeff
+
+
+if __name__ == "__main__":
+    thetas = [54.3, 48.5]  # angle of incident planewave (degrees)
+    Ls = [0.4, 1.2]  # position of mode monitor relative to flat interface
+    pols = [Polarization.S, Polarization.P]  # polarization of incident planewave
+
+    for pol, theta, L in zip(pols, thetas, Ls):
+        theta_rad = np.radians(theta)
+        rc_m = refl_coeff_meep(pol, theta_rad, L)
+        rc_f = refl_coeff_Fresnel(pol, theta_rad)
+
+        rc_m_str = f"{rc_m.real:.5f}{rc_m.imag:+.5f}j"
+        rc_f_str = f"{rc_f.real:.5f}{rc_f.imag:+.5f}j"
+        print(
+            f"refl-coeff:, {pol.name}, {theta}, {rc_m_str} (Meep), "
+            f"{rc_f_str} (Fresnel)"
+        )
+
+        mag_m = abs(rc_m)
+        mag_f = abs(rc_f)
+        err_mag = abs(mag_m - mag_f) / mag_f
+        print(
+            f"magnitude:, {mag_m:.5f} (Meep), {mag_f:.5f} (Fresnel), "
+            f"{err_mag:.5f} (error)"
+        )
+
+        phase_m = cmath.phase(rc_m)
+        phase_f = cmath.phase(rc_f)
+        err_phase = abs(phase_m - phase_f) / abs(phase_f)
+        print(
+            f"phase:, {phase_m:.5f} (Meep), {phase_f:.5f} (Fresnel), "
+            f"{err_phase:.5f} (error)"
+        )