diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst
index 5e7007a1ab..6faeb59e19 100644
--- a/docs/source/usersguide/settings.rst
+++ b/docs/source/usersguide/settings.rst
@@ -297,17 +297,17 @@ Tokamak Plasma Sources
 For fusion applications, the :class:`openmc.TokamakSource` class provides a
 native parametric neutron source for tokamak plasmas. Rather than specifying
 spatial, angular, and energy distributions separately, the source is defined by
-the plasma geometry (using a Miller-style flux-surface parameterization) and a
-radial emission profile. Source sites are sampled directly from the plasma
-volume without rejection.
+the plasma geometry (using a `Miller-style flux-surface parameterization
+<https://doi.org/10.1063/1.872666>`_) and a radial emission profile. Source
+sites are sampled directly from the plasma volume without rejection.
 
 The plasma shape is described by the major radius :math:`R_0`, minor radius
 :math:`a`, elongation :math:`\kappa`, triangularity :math:`\delta`, and
 Shafranov shift :math:`\Delta`. The neutron birth profile is given as an
 emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that
 runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of
-the profile matters, since it is normalized internally. The emission density
-is linearly interpolated between the supplied points and refined internally for
+the profile matters, since it is normalized internally. The emission density is
+linearly interpolated between the supplied points and refined internally for
 radial sampling. For example::
 
   import numpy as np
@@ -333,8 +333,8 @@ The ``energy`` argument accepts either a single
 sequence with one distribution per ``r_over_a`` grid point to model a
 radially-varying neutron spectrum (energies are then sampled by stochastic
 interpolation between the two distributions bracketing the sampled radius). A
-time distribution can be given with the ``time`` argument; by default,
-particles are born at :math:`t=0`. The toroidal extent can be restricted with
+time distribution can be given with the ``time`` argument; by default, particles
+are born at :math:`t=0`. The toroidal extent can be restricted with
 ``phi_start`` and ``phi_extent`` to model a sector of the plasma, and
 ``vertical_shift`` translates the plasma center along the z-axis.
 
diff --git a/openmc/source.py b/openmc/source.py
index 4e9f84dfc0..13b1492f67 100644
--- a/openmc/source.py
+++ b/openmc/source.py
@@ -899,19 +899,27 @@ class FileSource(SourceBase):
 
 
 class TokamakSource(SourceBase):
-    """A source representing neutron emission from a tokamak plasma.
+    r"""A source representing neutron emission from a tokamak plasma.
 
     This source samples neutron positions from a tokamak plasma geometry using
     Miller-style flux surface parameterization. The user provides an emission
     profile S(r/a) as a function of normalized minor radius, along with one or
     more energy distributions.
 
-    The flux surface parameterization is:
-        R = R0 + r*cos(α + δ*sin(α)) + Δ*(1 - (r/a)²)
-        Z = κ * r * sin(α)
+    The flux surface parameterization is
 
-    where R0 is major radius, a is minor radius, κ is elongation, δ is
-    triangularity, and Δ is the Shafranov shift.
+    .. math::
+
+        \begin{aligned}
+        R &= R_0 + r \cos\left(\alpha + \delta \sin\alpha\right)
+             + \Delta \left[1 - \left(\frac{r}{a}\right)^2\right] \\
+        Z &= Z_\mathrm{shift} + \kappa r \sin\alpha
+        \end{aligned}
+
+    where :math:`R_0` is major radius, :math:`a` is minor radius,
+    :math:`\kappa` is elongation, :math:`\delta` is triangularity,
+    :math:`\Delta` is the Shafranov shift, and :math:`Z_\mathrm{shift}` is
+    the vertical shift.
 
     .. versionadded:: 0.15.4
 
diff --git a/src/source.cpp b/src/source.cpp
index 8106ac1780..1a537e77b8 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -850,8 +850,8 @@ void TokamakSource::precompute_sampling_cdfs()
   radial_poly_b_ = 0.375 * c1 * eps; // 3/8 * c1 * eps
   radial_poly_c_ = 2.0 * eps * Dt;
 
-  // Build a refined radial grid that retains the user-provided knots. The
-  // emission density is interpreted as linear-linear between those knots.
+  // Build a refined radial grid that retains the user-specified grid points.
+  // The emission density is interpreted as linear-linear between those points.
   constexpr int MIN_SUBINTERVALS = 8;
   constexpr double MAX_GRID_SPACING = 1.0e-3;
   vector<double> radial_grid {r_over_a_.front()};