diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index db60b19e..3f9bbaeb 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -138,7 +138,7 @@ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, pr
         odet.ising_start = searchsortedfirst(getfield.(intr.sing, :psifac), odet.psifac) - 1
     end
 
-    # Initialize solutions with the identity matrix for U_22 as described in [Glasser PoP 2016] Section VI
+    # Initialize solutions with the identity matrix for U_22 [Glasser Phys. Plasmas 2016 112506 Section VI]
     for ipert in 1:intr.numpert_total
         odet.u[ipert, ipert, 2] = 1
     end
diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
index 65a7a24c..84016331 100644
--- a/src/ForceFreeStates/Fourfit.jl
+++ b/src/ForceFreeStates/Fourfit.jl
@@ -130,7 +130,7 @@ end
     make_matrix(metric::MetricData, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) -> FourFitVars
 
 Constructs main ForceFreeStates matrices for a given toroidal mode number and returns
-them as a new `FourFitVars` object. See the appendix of the 2016 Glasser
+them as a new `FourFitVars` object. See the appendix of the Glasser Phys. Plasmas 2016 112506
 DCON paper for details on the matrix definitions. Performs the same function
 as `fourfit_make_matrix` in the Fortran code, except F, G, and K are now
 stored as dense matrices. The matrix F is stored in factorized form with
@@ -281,9 +281,9 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
             end
         end
 
-        # Factorize and build composites
-        # Note: we store the nonsingular forms F̄ and K̄ with F = QF̄Qᴴ, K = QK̄ (eq. 29 in Glasser 2016)
-        # We multiply by Q (singfac) later when performing computations later
+        # Factorize and build composite matrices [Glasser Phys. Plasmas 2016 112506 Appendix A]
+        # Nonsingular forms F̄ and K̄ with F = QF̄Qᴴ, K = QK̄ [Glasser Phys. Plasmas 2016 112506 eq. 29]
+        # We multiply by Q (singfac = m - n*q) later when performing computations
         amat = reshape(amats_flatview, intr.numpert_total, intr.numpert_total)
         cmat = reshape(cmats_flatview, intr.numpert_total, intr.numpert_total)
         dmat = reshape(dmats_flatview, intr.numpert_total, intr.numpert_total)
@@ -294,12 +294,14 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
         gmat = reshape(gmats_flatview, intr.numpert_total, intr.numpert_total)
         # TODO: Fortran threw an error if factorization fails for A/F due to small matrix bandwidth,
         # Add this check back in if we implement banded matrices
+
+        # Schur complement reduction [Glasser Phys. Plasmas 2016 112506 eq. A5-A7]
         amat_fact = cholesky(Hermitian(amat, :L))
-        ldiv!(a_inv_dmat_temp, amat_fact, dmat)
-        ldiv!(a_inv_cmat_temp, amat_fact, cmat)
-        fmat .-= adjoint(dmat) * a_inv_dmat_temp
-        kmat .= emat .- (adjoint(kmat) * a_inv_cmat_temp)
-        gmat .= hmat .- (adjoint(cmat) * a_inv_cmat_temp)
+        ldiv!(a_inv_dmat_temp, amat_fact, dmat)         # A⁻¹D
+        ldiv!(a_inv_cmat_temp, amat_fact, cmat)         # A⁻¹C
+        fmat .-= adjoint(dmat) * a_inv_dmat_temp        # F̃ = F - D†A⁻¹D
+        kmat .= emat .- (adjoint(kmat) * a_inv_cmat_temp)  # K̃ = E - K†A⁻¹C
+        gmat .= hmat .- (adjoint(cmat) * a_inv_cmat_temp)  # G̃ = H - C†A⁻¹C
 
         # Store factorized F matrix (lower triangular only) since we always will need F⁻¹ later
         # and this make computation more efficient via combined forward and back substitution
diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl
index e0cfa5e1..6c4ba56f 100644
--- a/src/ForceFreeStates/Free.jl
+++ b/src/ForceFreeStates/Free.jl
@@ -36,7 +36,7 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr
         # Compute vacuum energy matrix and both Green's functions
         wv_block, vac.grri, vac.grre, vac.xzpts = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
 
-        # Equation 126 in Chance 1997 - scale by (m - n*q)(m' - n*q)
+        # Scale by (m - n*q)(m' - n*q) [Chance Phys. Plasmas 1997 2161 eq. 126]
         singfac = collect(intr.mlow:intr.mhigh) .- (n * intr.qlim)
         @inbounds for ipert in 1:intr.mpert
             @views wv_block[ipert, :] .*= singfac[ipert]
diff --git a/src/ForceFreeStates/Mercier.jl b/src/ForceFreeStates/Mercier.jl
index 79f6621e..f18f235a 100644
--- a/src/ForceFreeStates/Mercier.jl
+++ b/src/ForceFreeStates/Mercier.jl
@@ -1,9 +1,9 @@
 """
     mercier_scan!(locstab_fs::Array{Float64,5}, plasma_eq::PlasmaEquilibrium)
 
-Evaluates Mercier criterion for local stability and modifies results in place
-within the local stability array. Performs the same function as `mercier_scan`
-in the Fortran code.
+Evaluates Mercier criterion for local stability [Glasser Phys. Plasmas 2016 112506]
+and modifies results in place within the local stability array.
+Performs the same function as `mercier_scan` in the Fortran code.
 """
 function mercier_scan!(locstab_fs::Matrix{Float64}, plasma_eq::Equilibrium.PlasmaEquilibrium)
 
diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl
index 0bc8241c..d52665d6 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -138,7 +138,7 @@ Formerly `sing_vmat!`. Returns a `SingAsymptotics` struct with the computed data
 mutating the `SingType` struct. This makes it clear that asymptotics are computed on-demand
 for ideal ForceFreeStates and are not inherent properties of the singular surface.
 
-See equations 41-48 in the 2016 Glasser DCON paper for the mathematical details.
+See equations 41-48 in the Glasser Phys. Plasmas 2016 112506 for the mathematical details.
 
 ### Arguments
 
@@ -157,7 +157,7 @@ function compute_sing_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl,
 
     # Compute the resonant (r) and nonresonant (n) indices of the shearing transformation matrix R
     # 1 indexes along the N*M dimension, and 2 along the 2*N*M dimension
-    # In 2D, see eq. 41 of 2016 Glasser DCON paper
+    # In 2D, see eq. 41 of Glasser Phys. Plasmas 2016 112506
     # TODO: if we remove the 3rd dimension, no need for both r1 and r2
     ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert
     r1 = ipert_res
@@ -689,7 +689,9 @@ end
         psieval::Float64
     )
 
-Evaluate the derivative of the Euler-Lagrange equations, i.e. u' in equation 24 of Glasser 2016.
+Evaluate the derivative of the Euler-Lagrange equations [Glasser Phys. Plasmas 2016 112506 eq. 24].
+This implements du/dψ for the ideal MHD eigenvalue problem.
+
 This function performs the same role as `sing_der` in the Fortran code, with main differences
 coming from hiding LAPACK operations under the hood via Julia's LinearAlgebra package,
 so the code is much more straightforward.
diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl
index fd053e5b..d2c39e15 100644
--- a/src/PerturbedEquilibrium/FieldReconstruction.jl
+++ b/src/PerturbedEquilibrium/FieldReconstruction.jl
@@ -3,6 +3,7 @@ Field reconstruction from eigenmode response.
 
 Converts eigenmode response coefficients to physical displacement and magnetic
 field perturbations in mode space, following the GPEC gpeq module approach.
+[Park Phys. Plasmas 2007 052110]
 
 This module mimics the GPEC Fortran subroutines:
 - gpeq_sol: Get equilibrium solution at each radial point
@@ -18,7 +19,9 @@ Fields are computed using ideal MHD relations in flux coordinates:
 
 where χ₁ = 2π * Ψ₀ (total poloidal flux normalization).
 
-Reference: GPEC gpeq.f lines 100-102
+References:
+- [Park Phys. Plasmas 2007 052110] - 3D equilibrium perturbations
+- GPEC gpeq.f lines 100-102
 """
 
 """
@@ -33,8 +36,8 @@ Reference: GPEC gpeq.f lines 100-102
 Reconstruct displacement and magnetic field from eigenmode response in mode space.
 
 This function mimics GPEC's gpeq module, computing fields using ideal MHD
-algebraic relations in flux coordinates. All fields are returned in mode space
-[npsi, mpert] rather than real space.
+algebraic relations in flux coordinates [Park Phys. Plasmas 2007 052110].
+All fields are returned in mode space [npsi, mpert] rather than real space.
 
 # Process (following GPEC)
 
@@ -86,6 +89,7 @@ function reconstruct_physical_fields(
     )
 
     # Step 2: Compute perturbed field in mode space using ideal MHD relations
+    # [Park Phys. Plasmas 2007 052110 eq. 8-10]
     # This mimics GPEC's gpeq_sol and gpeq_contra
     b_psi_modes, b_theta_modes, b_zeta_modes = compute_perturbed_field_modes(
         xi_psi_modes,
@@ -196,9 +200,9 @@ Compute perturbed magnetic field from displacement using ideal MHD relations in
 
 This function mimics GPEC's gpeq_sol subroutine (gpeq.f lines 100-102), computing
 contravariant field components from covariant displacement using the algebraic
-ideal MHD relations in flux coordinates.
+ideal MHD relations in flux coordinates [Park Phys. Plasmas 2007 052110 eq. 8-10].
 
-# Ideal MHD Relations (from GPEC gpeq.f)
+# Ideal MHD Relations [Park Phys. Plasmas 2007 052110 eq. 8-10]
 
 ```fortran
 bwp_mn = (chi1*singfac*twopi*ifac*xsp_mn)                              ! b^ψ
@@ -298,14 +302,14 @@ function compute_perturbed_field_modes(
             xsp1 = xi_psi1_modes[ipsi, ipert]   # ∂ξ_ψ/∂ψ
             xss = 0.0 + 0.0im                    # ξ_ζ = 0 (not computed from ForceFreeStates)
 
-            # GPEC gpeq.f line 100-102: Compute contravariant field
-            # b^ψ = i * χ₁ * (m - n*q) * ξ_ψ
+            # Contravariant field from ideal MHD [Park Phys. Plasmas 2007 052110 eq. 8-10]
+            # b^ψ = i * χ₁ * (m - n*q) * ξ_ψ  [eq. 8]
             b_psi_modes[ipsi, ipert] = chi1 * singfac * twopi * ifac * xsp
 
-            # b^θ = -i * (χ₁ * ∂ξ_ψ/∂ψ + n * ξ_ζ)
+            # b^θ = -i * (χ₁ * ∂ξ_ψ/∂ψ + n * ξ_ζ)  [eq. 9]
             b_theta_modes[ipsi, ipert] = -(chi1 * xsp1 + twopi * ifac * nn * xss)
 
-            # b^ζ = -i * (χ₁ * (q'*ξ_ψ + q*∂ξ_ψ/∂ψ) + m * ξ_ζ)
+            # b^ζ = -i * (χ₁ * (q'*ξ_ψ + q*∂ξ_ψ/∂ψ) + m * ξ_ζ)  [eq. 10]
             b_zeta_modes[ipsi, ipert] = -(chi1 * (q1 * xsp + q * xsp1) + twopi * ifac * m * xss)
         end
     end
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index abc5d569..46f62445 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -3,6 +3,8 @@ Response matrix construction for perturbed equilibrium calculations.
 
 Based on gpresp.f from GPEC, implementing resp_index=0 (energy-based inductance).
 Uses ForceFreeStates eigenmode solutions and vacuum response data.
+
+Reference: [Park Phys. Plasmas 2009 056115]
 """
 
 # Use FourierTransform utility instead of FFTW for theta ↔ mode transforms
@@ -91,15 +93,17 @@ end
 Compute normal magnetic field at plasma boundary from eigenmode displacements.
 
 This is the key step that converts eigenmode displacements to magnetic field perturbations
-at the plasma surface. Formula from GPEC:
+at the plasma surface. From the ideal MHD constraint [Park Phys. Plasmas 2009 056115 eq. 4]:
+
+    B_n = i * (dΨ/dρ) * (m - n*q) * ξ_ψ
 
-    bwp_mn[i,j] = i * (dΨ/dρ) * (m[i] - n*q_boundary) * ξ_ψ[i,j]
+where ξ_ψ is the radial displacement eigenfunction.
 
-## Physical Interpretation:
+## Physical Interpretation [Park Phys. Plasmas 2007 052110 Section II]:
 - ξ_ψ[i,j]: Displacement of mode i due to eigenmode j
-- singfac[i] = m[i] - n*q: Measures distance from rational surface
-- dΨ/dρ: Converts displacement to flux perturbation
-- Factor of i: Phase relationship for oscillating fields
+- singfac[i] = m[i] - n*q: Singular factor measuring distance from rational surface
+- dΨ/dρ: Converts displacement to flux perturbation (poloidal flux gradient)
+- Factor of i: Phase relationship for oscillating fields in complex representation
 
 ## Arguments
 - `boundary_data`: Output from extract_boundary_displacements()
@@ -126,7 +130,8 @@ function compute_normal_magnetic_field(
     dPsi_drho = boundary_data.dPsi_drho
     q_boundary = boundary_data.q_boundary
 
-    # Compute singular factor for each Fourier mode: singfac[i] = m[i] - n*q_boundary
+    # Compute singular factor for each Fourier mode [Park Phys. Plasmas 2009 056115 eq. 4]
+    # singfac[i] = m[i] - n*q_boundary measures distance from rational surface
     # Mode indexing: modes are ordered as (m, n) pairs
     # Linear index i corresponds to: m = (i-1) % mpert + mlow, n = (i-1) ÷ mpert + nlow
     singfac = zeros(Float64, numpert_total)
@@ -136,7 +141,8 @@ function compute_normal_magnetic_field(
         singfac[i] = m_mode - n_mode * q_boundary
     end
 
-    # Compute normal magnetic field: bwp_mn[i,j] = i * (dΨ/dρ) * singfac[i] * ξ_ψ[i,j]
+    # Compute normal magnetic field [Park Phys. Plasmas 2009 056115 eq. 4]
+    # bwp_mn[i,j] = i * (dΨ/dρ) * singfac[i] * ξ_ψ[i,j]
     for i in 1:numpert_total
         for j in 1:numpert_total
             bwp_mn[i, j] = 1im * dPsi_drho * singfac[i] * ξ_psi[i, j]
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index 924b8ba2..626c8638 100644
--- a/src/PerturbedEquilibrium/SingularCoupling.jl
+++ b/src/PerturbedEquilibrium/SingularCoupling.jl
@@ -8,7 +8,10 @@ Implements GPEC's singular surface analysis (gpout.f lines 585-665, 1700-1810):
 - Chirikov parameters (island overlap)
 - Coupling matrices
 
-Reference: ~/Code/gpec/gpec/gpout.f
+References:
+- [Park Phys. Plasmas 2009 056115] - Plasma response calculation
+- [Park Phys. Rev. Lett. 2007 195003] - RMP control
+- [Glasser Phys. Plasmas 2016 072505] - Resistive Δ' calculation
 """
 
 """
@@ -26,6 +29,7 @@ Compute singular layer coupling metrics and island diagnostics.
 
 Calculates the coupling between forcing modes and resonant surfaces, which determines
 the effectiveness of external perturbations at rational surfaces where q = m/n.
+[Park Phys. Plasmas 2009 056115]
 
 Implements GPEC algorithm from gpout.f:
 
diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl
index 1da6bc1d..07a7816c 100644
--- a/src/Vacuum/Kernel2D.jl
+++ b/src/Vacuum/Kernel2D.jl
@@ -70,7 +70,7 @@ but grad_greenfunction is not since it fills a different block of the
 
   - Uses Simpson's rule for integration away from singular points
   - Uses Gaussian quadrature near singular points for improved accuracy
-  - Implements analytical singularity removal following Chance 1997
+  - Implements analytical singularity removal [Chance Phys. Plasmas 1997 2161]
 """
 function kernel!(
     grad_greenfunction::Matrix{Float64},
@@ -98,7 +98,7 @@ function kernel!(
     # 𝒢ⁿ only needed for plasma as source term (RHS of eqs. 26/27 in Chance 1997)
     populate_greenfunction = source isa PlasmaGeometry
 
-    # S₁ᵢ in Chance 1997, eq.(78)
+    # S₁ᵢ logarithmic correction factors [Chance Phys. Plasmas 1997 2161 eq. 78]
     log_correction_0=16.0*dtheta*(log(2*dtheta)-68.0/15.0)/15.0
     log_correction_1=128.0*dtheta*(log(2*dtheta)-8.0/15.0)/45.0
     log_correction_2=4.0*dtheta*(7.0*log(2*dtheta)-11.0/15.0)/45.0
@@ -132,7 +132,7 @@ function kernel!(
                 greenfunction[j, isrc] += G_n * wsimpson
             end
             grad_greenfunction_block[j, isrc] += gradG_n * wsimpson
-            # Subtract regular integral component of δⱼᵢK⁰ in eq. 83
+            # Subtract regular integral component of δⱼᵢK⁰ [Chance Phys. Plasmas 1997 2161 eq. 83]
             grad_greenfunction_block[j, j] -= gradG_0 * wsimpson
         end
 
@@ -156,19 +156,20 @@ function kernel!(
 
                 # Redefine hardcoded Gaussian weights on the interval [-1, 1] to physical interval with length 2 * dtheta
                 wgauss = GAUSSIANWEIGHTS[ig] * dtheta
-                # Calculate p = θ/Δ = (θⱼ - θ')/Δ
+                # Normalized coordinate p = (θⱼ - θ')/Δθ
                 pgauss=(theta_gauss[ig]-theta_obs)/dtheta
-                # Compute 5-point Lagrange basis polynomials at the Gauss point and multiply by quadrature weight
-                A0 = (pgauss^2-1)*(pgauss^2-4)/4.0 * wgauss
-                A1_plus = -(pgauss+1)*pgauss*(pgauss^2-4)/6.0 * wgauss
-                A1_minus = -(pgauss-1)*pgauss*(pgauss^2-4)/6.0 * wgauss
-                A2_plus = (pgauss^2-1)*pgauss*(pgauss+2)/24.0 * wgauss
-                A2_minus = (pgauss^2-1)*pgauss*(pgauss-2)/24.0 * wgauss
+                # 5-point Lagrange interpolation polynomials [Chance Phys. Plasmas 1997 2161 eq. 76]
+                # Weighted by Gaussian quadrature for accurate singular integral evaluation
+                A0 = (pgauss^2-1)*(pgauss^2-4)/4.0 * wgauss                    # L₀(p) centered at j
+                A1_plus = -(pgauss+1)*pgauss*(pgauss^2-4)/6.0 * wgauss        # L₁(p) at j+1
+                A1_minus = -(pgauss-1)*pgauss*(pgauss^2-4)/6.0 * wgauss       # L₋₁(p) at j-1
+                A2_plus = (pgauss^2-1)*pgauss*(pgauss+2)/24.0 * wgauss        # L₂(p) at j+2
+                A2_minus = (pgauss^2-1)*pgauss*(pgauss-2)/24.0 * wgauss       # L₋₂(p) at j-2
 
                 # First type of singularity: 𝒢ⁿ, occurs plasma as source only (see RHS of Chance eqs. 26/27)
                 if populate_greenfunction
                     if observer isa PlasmaGeometry
-                        # Remove singular behavior by adding on leading-order term, Chance eq.(75)
+                        # Remove singular behavior by adding on leading-order term [Chance Phys. Plasmas 1997 2161 eq. 75]
                         G_n += log((theta_obs - theta_gauss[ig])^2) / x_obs
                     end
                     greenfunction[j, sing_idx[1]] += G_n * A2_minus
@@ -178,7 +179,7 @@ function kernel!(
                     greenfunction[j, sing_idx[5]] += G_n * A2_plus
                 end
 
-                # Second type of singularity: 𝒦ⁿ (Eq. 86: 𝒦ⁿαᵢ - δⱼᵢK⁰)
+                # Second type of singularity: 𝒦ⁿ [Chance Phys. Plasmas 1997 2161 eq. 83, 86]
                 grad_greenfunction_block[j, sing_idx[1]] += gradG_n * A2_minus
                 grad_greenfunction_block[j, sing_idx[2]] += gradG_n * A1_minus
                 grad_greenfunction_block[j, sing_idx[3]] += gradG_n * A0
@@ -189,7 +190,7 @@ function kernel!(
             end
         end
 
-        # Subtract off analytic singular integral from Chance eq.(75) if plasma-plasma block
+        # Subtract off analytic singular integral [Chance Phys. Plasmas 1997 2161 eq. 75] if plasma-plasma block
         if populate_greenfunction && observer isa PlasmaGeometry
             greenfunction[j, sing_idx[1]] -= log_correction_2 / x_obs
             greenfunction[j, sing_idx[2]] -= log_correction_1 / x_obs
@@ -205,16 +206,8 @@ function kernel!(
         grad_greenfunction_block .*= -1
     end
 
-    # Add analytic singular integral (second type) from Table I of Chance 1997 + existing δⱼᵢ in eq. 69
-    # Would need to pass in wall geometry to generalize this to open walls
-    is_closed_toroidal = true
-    if is_closed_toroidal # Chance eq. 89
-        residue = (observer isa WallGeometry) ? 0.0 : (source isa PlasmaGeometry ? 2.0 : -2.0)
-    else # Chance eq. 90
-        # TODO: this line can be gotten rid of if we are never doing open walls
-        residue = (typeof(observer) == typeof(source)) ? 2.0 : 0.0
-    end
-    # Add residue value from eq. 89/90 to block diagonal
+    # Add analytic singular integral (second type) to block diagonal [Chance Phys. Plasmas 1997 2161 Table I, eq. 69, 89]
+    residue = (observer isa WallGeometry) ? 0.0 : (source isa PlasmaGeometry ? 2.0 : -2.0)
     @inbounds for i in 1:mtheta
         grad_greenfunction_block[i, i] += residue
     end
@@ -639,13 +632,13 @@ function green(x_obs::Float64, z_obs::Float64, x_source::Float64, z_source::Floa
 
     ρ2 = x_minus2 + ζ2
 
-    # Chance 1997 eq.(41) ℛ = R
+    # Distance parameter ℛ [Chance Phys. Plasmas 1997 2161 eq. 41]
     R4 = ρ2 * (ρ2 + 4 * x_multiple)
     R2 = sqrt(R4)
     R = sqrt(R2)
     R5 = R4 * R
 
-    # Chance 1997 eq.(42) 𝘴 = s
+    # Argument of Legendre function 𝘴 [Chance Phys. Plasmas 1997 2161 eq. 42]
     s = (x_obs2 + x_source2 + ζ2) / R2
 
     # Legendre functions for
@@ -661,29 +654,30 @@ function green(x_obs::Float64, z_obs::Float64, x_source::Float64, z_source::Floa
     pnp1 = legendre[end]
     pn = legendre[end-1]
 
-    # Chance 1997 eq.(40) 2π𝒢ⁿ = G_n
+    # Green's function 2π𝒢ⁿ = G_n [Chance Phys. Plasmas 1997 2161 eq. 40]
     gg = 2 * sqrt(π) * gamma(0.5 - n) / R
     G_n = gg * pn
 
-    # Chance 1997 eq.(44) (Note this equation in the paper has an erroneous extra factor of 2π)
+    # Gradient factor [Chance Phys. Plasmas 1997 2161 eq. 44]
+    # NOTE: Paper has erroneous extra factor of 2π
     grad_gg = gg / R4 / 2π
 
-    # ∂Gⁿ/∂X' = dG_dX
+    # Derivatives of Green's function [Chance Phys. Plasmas 1997 2161 eq. 36-38]
+    # ∂Gⁿ/∂X' using chain rule: ∂Gⁿ/∂X' = (∂Gⁿ/∂R)(∂R/∂X') + (∂Gⁿ/∂s)(∂s/∂X')
     xterm1 = (n * (x_obs2 + x_source2 + ζ2) * (x_obs2 - x_source2 + ζ2) - x_source2*(x_source2-x_obs2+ζ2)) * pn
     xterm2 = (2.0 * x_source * x_obs * (x_obs2-x_source2+ζ2)) * pnp1
     dG_dX = grad_gg * (xterm1 + xterm2) / x_source
 
-    # ∂Gⁿ/∂Z' = dG_dZ
+    # ∂Gⁿ/∂Z' using chain rule
     zterm1 = (2.0 * n + 1.0) * (x_obs2 + x_source2 + ζ2) * pn
     zterm2 = 4.0 * x_multiple * pnp1
     dG_dZ = grad_gg * (zterm1 + zterm2) * ζ
 
-    # Chance 1997 eq.(51)
-    # 𝒥 ∇'𝒢ⁿ∇'ℒ = aval
-    # ∂X'/∂θ = xtp, ∂Z'/∂θ = ztp
+    # Coupling term 𝒥 ∇'𝒢ⁿ∇'ℒ [Chance Phys. Plasmas 1997 2161 eq. 51]
+    # Jacobian factor from coordinate transformation
     coupling_n = -x_source * (dz_dtheta * dG_dX - dx_dtheta * dG_dZ)
 
-    # for 𝓃⩵0,  aval0 = 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ
+    # Special case for n=0: coupling_0 = 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ
     dG_dX0_R5 = ((2.0 * x_obs * (x_obs2-x_source2+ζ2)) * p1 - x_source * (x_source2-x_obs2+ζ2) * p0)
     dG_dZ0_R5 = ζ * ((x_obs2 + x_source2 + ζ2) * p0 + 4.0 * x_multiple * p1)
     coupling_0 = -x_source * (dz_dtheta * dG_dX0_R5 - dx_dtheta * dG_dZ0_R5) / R5
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 1b848b53..5b57e880 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -221,7 +221,7 @@ For kernelsign < 0 (interior potential), multiply by -1 and add 2 to diagonal.
 function apply_kernelsign!(grad_greenfunction_mat::Matrix{Float64}, kernelsign::Float64, mtheta::Int)
     if kernelsign < 0
         grad_greenfunction_mat .*= kernelsign
-        # Account for factor of 2 in diagonal terms in eq. 90 of Chance
+        # Account for factor of 2 in diagonal terms [Chance Phys. Plasmas 1997 2161 eq. 90]
         for i in 1:(2*mtheta)
             grad_greenfunction_mat[i, i] += 2.0
         end
@@ -307,7 +307,7 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe
 
     # Add cn0 to make grdgre nonsingular for n=0 modes
     cn0 = 1.0 # expose to user if anyone ever actually tries to use this
-    (n == 0 && !wall.nowall && wall.is_closed_toroidal) && begin
+    (n == 0 && !wall.nowall) && begin
         @warn "Adding $cn0 to diagonal of grdgre to regularize n=0 mode; this may affect accuracy of results."
         mth12 = wall.nowall ? mtheta : 2 * mtheta
         for i in 1:mth12, j in 1:mth12
@@ -339,18 +339,18 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe
 
     # There's some logic that computes xpass/zpass and chiwc/chiws here, might eventually be needed?
 
-    wv = zeros(mpert, mpert)
-    xzpts = zeros(Float64, inputs.mtheta, 4)
+    wv = zeros(ComplexF64, mpert, mpert)
+    xzpts = zeros(inputs.mtheta, 4)
     if !green_only
-        # Perform inverse Fourier transforms to get response matrix components (eq. 115-118 of Chance 2007)
+        # Perform inverse Fourier transforms to get response matrix components [Chance Phys. Plasmas 2007 052506 eq. 115-118]
         arr, aii, ari, air = ntuple(_ -> zeros(mpert, mpert), 4)
         fourier_inverse_transform!(arr, grre, cos_mn_basis, PLASMA_ROW_OFFSET, COS_COL_OFFSET)
         fourier_inverse_transform!(aii, grre, sin_mn_basis, PLASMA_ROW_OFFSET, SIN_COL_OFFSET)
         fourier_inverse_transform!(ari, grre, sin_mn_basis, PLASMA_ROW_OFFSET, COS_COL_OFFSET)
         fourier_inverse_transform!(air, grre, cos_mn_basis, PLASMA_ROW_OFFSET, SIN_COL_OFFSET)
 
-        # Final form of vacuum response matrix (eq. 114 of Chance 2007)
-        wv = complex.(arr .+ aii, air .- ari)
+        # Final form of vacuum response matrix [Chance Phys. Plasmas 2007 052506 eq. 114]
+        wv .= complex.(arr .+ aii, air .- ari)
 
         # Force symmetry of response matrix if desired
         force_wv_symmetry && hermitianpart!(wv)
diff --git a/src/Vacuum/VacuumFromEquilibrium.jl b/src/Vacuum/VacuumFromEquilibrium.jl
index 0b0481dd..d9f9fa4b 100644
--- a/src/Vacuum/VacuumFromEquilibrium.jl
+++ b/src/Vacuum/VacuumFromEquilibrium.jl
@@ -35,8 +35,9 @@ From these we compute:
   - R = R₀ + r_minor * cos(θ_cyl)
   - Z = Z₀ + r_minor * sin(θ_cyl)
 
-## Reference    # Allocate output arrays
+## Reference
 
+[Chance Phys. Plasmas 1997 2161]
 Matches GPEC's ahg_write and gpvacuum_flxsurf approach (gpvacuum.f line 291-296)
 """
 function extract_plasma_surface_at_psi(equil::Equilibrium.PlasmaEquilibrium, ψ::Float64)