From b7c2a1959361de80da622f62d74c866773ac0be1 Mon Sep 17 00:00:00 2001
From: logan-nc <6198372+logan-nc@users.noreply.github.com>
Date: Thu, 12 Feb 2026 21:55:34 -0500
Subject: [PATCH 1/3] DOCUMENTATION - IMPROVEMENT - Standardize citations to
 consistent format across all modules

Add standardized paper citations throughout the codebase using the format:
[Lead Author Last Name, Journal, Year Page]

Changes:
- Vacuum module: Add citations to Chance 1997 and 2007 papers for equations
  and algorithms in Kernel2D.jl, Vacuum.jl, VacuumFromEquilibrium.jl
- ForceFreeStates module: Standardize Glasser 2016 citations in EulerLagrange.jl,
  Free.jl, Mercier.jl, Fourfit.jl, and Sing.jl
- PerturbedEquilibrium module: Add Park et al. citations (2007-2009) to
  ResponseMatrices.jl, SingularCoupling.jl, and FieldReconstruction.jl

This improves traceability between theory papers and implementation, making it
easier for users to understand the mathematical foundation of the algorithms.

Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
---
 src/ForceFreeStates/EulerLagrange.jl          |  2 +-
 src/ForceFreeStates/Fourfit.jl                |  4 +--
 src/ForceFreeStates/Free.jl                   |  2 +-
 src/ForceFreeStates/Mercier.jl                |  6 ++---
 src/ForceFreeStates/Sing.jl                   |  4 +--
 .../FieldReconstruction.jl                    |  9 ++++---
 src/PerturbedEquilibrium/ResponseMatrices.jl  |  4 ++-
 src/PerturbedEquilibrium/SingularCoupling.jl  |  6 ++++-
 src/Vacuum/DataTypes.jl                       |  2 +-
 src/Vacuum/Kernel2D.jl                        | 25 +++++++++----------
 src/Vacuum/Vacuum.jl                          |  6 ++---
 src/Vacuum/VacuumFromEquilibrium.jl           |  2 ++
 12 files changed, 41 insertions(+), 31 deletions(-)

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..335639da 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
@@ -282,7 +282,7 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
         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)
+        # Note: we store the nonsingular forms F̄ and K̄ with F = QF̄Qᴴ, K = QK̄ (eq. 29 in Glasser Phys. Plasmas 2016 112506)
         # We multiply by Q (singfac) later when performing computations later
         amat = reshape(amats_flatview, intr.numpert_total, intr.numpert_total)
         cmat = reshape(cmats_flatview, intr.numpert_total, intr.numpert_total)
diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl
index 9dd14a14..f19e3b9f 100644
--- a/src/ForceFreeStates/Free.jl
+++ b/src/ForceFreeStates/Free.jl
@@ -50,7 +50,7 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr
             @save "vacuum_response_inputs.jld2" benchmark_inputs
         end
 
-        # 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..22f702f1 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
diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl
index fd053e5b..8382dcc7 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)
 
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index abc5d569..13211f3e 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,7 +93,7 @@ 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 [Park Phys. Plasmas 2009 056115]. Formula from GPEC:
 
     bwp_mn[i,j] = i * (dΨ/dρ) * (m[i] - n*q_boundary) * ξ_ψ[i,j]
 
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index 99d25c70..9b0c950c 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:
 1. For each forcing mode, apply unit field and compute plasma response
diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index cfcde7bf..2983fa21 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -304,7 +304,7 @@ function initialize_wall(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_
     end
 
     if any(x_wall .<= 0.0) && !nowall
-        # to add support for x<0 walls, be sure to carefully replicate Chance's fortran code x<0 handling in the kernel function to account for the additional singularities associated with this
+        # to add support for x<0 walls, be sure to carefully replicate [Chance Phys. Plasmas 1997 2161] fortran code x<0 handling in the kernel function to account for the additional singularities associated with this
         error("Wall R-coordinates contain non-physical values (R <= 0). Check wall geometry.")
     end
 
diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl
index bce1b7a7..a3312682 100644
--- a/src/Vacuum/Kernel2D.jl
+++ b/src/Vacuum/Kernel2D.jl
@@ -74,7 +74,7 @@ but grad_greenfunction_mat 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_mat::Matrix{Float64},
@@ -104,12 +104,12 @@ function kernel!(
         error("Length of input arrays (xobs, zobs, xsource, zsce) are different. All length should be the same")
     end
 
-    # 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
 
-    # Used for Z'_θ and X'_θ in eq.(51)
+    # Used for Z'_θ and X'_θ [Chance Phys. Plasmas 1997 2161 eq. 51]
     # Close the loop for periodic BC by appending first point at the end
     theta_closed = vcat(collect(theta_grid), theta_grid[end] + dtheta)
     x_closed = vcat(x_sourcepoints, x_sourcepoints[1])
@@ -146,7 +146,7 @@ function kernel!(
             x_source=x_sourcepoints[isrc]
             z_source=z_sourcepoints[isrc]
 
-            # G_n is 2pi𝒢ⁿ; coupling_n is 𝒥 ∇'𝒢ⁿ∇'ℒ; coupling_0 is 𝒥 ∇'𝒢ⁿ∇'ℒ for n=0
+            # G_n is 2π𝒢ⁿ; coupling_n is 𝒥 ∇'𝒢ⁿ∇'ℒ [Chance Phys. Plasmas 1997 2161 eq. 36-42, 51]
             G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_source, z_source, dx_dtheta[isrc], dz_dtheta[isrc], n)
 
             # Sum contributions to Green's function matrices using Simpson weight
@@ -173,7 +173,7 @@ function kernel!(
                 dz_dtheta_gauss = d1_spline_z(theta_gauss0)
                 G_n, coupling_n, coupling_0 = green(x_obs, z_obs, x_gauss, z_gauss, dx_dtheta_gauss, dz_dtheta_gauss, n)
 
-                # Add logarithm to G_n to analytically isolate the singularity (first type), Chance eq.(75)
+                # Add logarithm to G_n to analytically isolate singularity [Chance Phys. Plasmas 1997 2161 eq. 75]
                 G_n_nonsingular = plasma_plasma_block ? G_n + log((theta_obs-theta_gauss[ig])^2)/x_obs : G_n
 
                 # Redefine hardcoded Gaussian weights on the interval [-1, 1] to physical interval with length 2 * dtheta
@@ -187,7 +187,7 @@ function kernel!(
                 A2_plus = (pgauss^2-1)*pgauss*(pgauss+2)/24.0 * wgauss
                 A2_minus = (pgauss^2-1)*pgauss*(pgauss-2)/24.0 * wgauss
 
-                # First type of singularity: 𝒢ⁿ, occurs plasma as source only (see RHS of Chance eqs. 26/27)
+                # First type of singularity: 𝒢ⁿ [Chance Phys. Plasmas 1997 2161 eq. 26-27]
                 if plasma_is_source
                     greenfunction_mat[j, js[1]] += G_n_nonsingular * A2_minus
                     greenfunction_mat[j, js[2]] += G_n_nonsingular * A1_minus
@@ -196,8 +196,7 @@ function kernel!(
                     greenfunction_mat[j, js[5]] += G_n_nonsingular * A2_plus
                 end
 
-                # Second type of singularity: 𝒦ⁿ
-                # Eq. 86: 𝒦ⁿαᵢ - δⱼᵢK⁰
+                # Second type of singularity: 𝒦ⁿ [Chance Phys. Plasmas 1997 2161 eq. 83, 86]
                 grad_green_work[js[1]] += isgn * coupling_n * A2_minus
                 grad_green_work[js[2]] += isgn * coupling_n * A1_minus
                 grad_green_work[js[3]] += isgn * coupling_n * A0
@@ -208,19 +207,19 @@ function kernel!(
             end
         end
 
-        # Set residue based on logic similar to Table I of Chance 1997 + existing δⱼᵢ in eq. 69
+        # Set residue [Chance Phys. Plasmas 1997 2161 Table I, eq. 69, 89-90]
         # Would need to pass in wall geometry to generalize this to open walls
         is_closed_toroidal = true
         if is_closed_toroidal
-            residue = (j1 == 2.0) ? 0.0 : (j2 == 1 ? 2.0 : -2.0) # Chance eq. 89
+            residue = (j1 == 2.0) ? 0.0 : (j2 == 1 ? 2.0 : -2.0) # eq. 89
         else
             # TODO: this line can be gotten rid of if we are never doing open walls
-            residue = (j1 == j2) ? 2.0 : 0.0 # Chance eq. 90
+            residue = (j1 == j2) ? 2.0 : 0.0 # eq. 90
         end
-        # Subtract regular integral component of δⱼᵢK⁰ in eq. 83 and add residue value in eq. 89/90
+        # Subtract regular integral component of δⱼᵢK⁰ [Chance Phys. Plasmas 1997 2161 eq. 83]
         grad_green_work[j] = grad_green_work[j] - isgn * grad_green_0 + residue
 
-        # 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
             greenfunction_mat[j, js[1]] -= log_correction_2 / x_obs
             greenfunction_mat[j, js[2]] -= log_correction_1 / x_obs
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index a2696d88..adc00fa2 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
@@ -339,7 +339,7 @@ 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?
 
-    # 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 = zeros(mpert, mpert)
     aii = zeros(mpert, mpert)
     ari = zeros(mpert, mpert)
@@ -349,7 +349,7 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe
     fourier_inverse_transform!(ari, grre, sin_ln_basis, 0, 0)
     fourier_inverse_transform!(air, grre, cos_ln_basis, 0, mpert)
 
-    # Final form of vacuum response matrix (eq. 114 of Chance 2007)
+    # 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
diff --git a/src/Vacuum/VacuumFromEquilibrium.jl b/src/Vacuum/VacuumFromEquilibrium.jl
index 296621a6..558dfe37 100644
--- a/src/Vacuum/VacuumFromEquilibrium.jl
+++ b/src/Vacuum/VacuumFromEquilibrium.jl
@@ -49,6 +49,7 @@ From these we compute:
 - delta = rzphi.f[3] / qa
 
 ## 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(
@@ -191,6 +192,7 @@ Much faster than full compute_vacuum_response() since it skips:
 - Mode coupling matrices
 
 ## Reference
+[Chance Phys. Plasmas 1997 2161]
 Mimics GPEC's gpvacuum_flxsurf which only computes Green's functions
 for surface inductance (gpvacuum.f line 279-283)
 """

From 066e7c6a2954de8dcb349dd5e2ed1ba310340314 Mon Sep 17 00:00:00 2001
From: logan-nc <6198372+logan-nc@users.noreply.github.com>
Date: Thu, 12 Feb 2026 22:02:40 -0500
Subject: [PATCH 2/3] DOCUMENTATION - IMPROVEMENT - Add detailed
 equation-specific citations throughout Vacuum, ForceFreeStates, and
 PerturbedEquilibrium modules
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Enhance code documentation with specific equation and section references from key papers:

Vacuum module (Kernel2D.jl):
- Add eq. 36-38, 41-42, 44, 51, 73, 76 citations for Green's function derivatives,
  Legendre parameters, Simpson integration, and Lagrange polynomials
- Improve comments on coordinate transformations and coupling terms
- Clarify physical interpretation of each computation step

ForceFreeStates module:
- Fourfit.jl: Add eq. 29, A5-A7 citations for nonsingular matrix forms and
  Schur complement reduction from [Glasser Phys. Plasmas 2016 112506 Appendix A]
- Sing.jl: Add eq. 24 citation for Euler-Lagrange equation derivative du/dψ

PerturbedEquilibrium module:
- ResponseMatrices.jl: Add eq. 4 citations from [Park Phys. Plasmas 2009 056115]
  for normal magnetic field computation and singular factor definition
- Add Section II citation from [Park Phys. Plasmas 2007 052110] for physical
  interpretation of displacement-to-field conversion
- FieldReconstruction.jl: Add eq. 8-10 citations from [Park Phys. Plasmas 2007 052110]
  for ideal MHD relations computing contravariant field components (b^ψ, b^θ, b^ζ)

These additions make it much easier to trace implementation back to theoretical
foundations and verify correctness of numerical algorithms.

Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
---
 src/ForceFreeStates/Fourfit.jl                | 18 +++++-----
 src/ForceFreeStates/Sing.jl                   |  4 ++-
 .../FieldReconstruction.jl                    | 13 +++----
 src/PerturbedEquilibrium/ResponseMatrices.jl  | 20 ++++++-----
 src/Vacuum/Kernel2D.jl                        | 36 ++++++++++---------
 5 files changed, 51 insertions(+), 40 deletions(-)

diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
index 335639da..84016331 100644
--- a/src/ForceFreeStates/Fourfit.jl
+++ b/src/ForceFreeStates/Fourfit.jl
@@ -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 Phys. Plasmas 2016 112506)
-        # 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/Sing.jl b/src/ForceFreeStates/Sing.jl
index 22f702f1..d52665d6 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -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 8382dcc7..d2c39e15 100644
--- a/src/PerturbedEquilibrium/FieldReconstruction.jl
+++ b/src/PerturbedEquilibrium/FieldReconstruction.jl
@@ -89,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,
@@ -199,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^ψ
@@ -301,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 13211f3e..46f62445 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -93,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 [Park Phys. Plasmas 2009 056115]. Formula from GPEC:
+at the plasma surface. From the ideal MHD constraint [Park Phys. Plasmas 2009 056115 eq. 4]:
 
-    bwp_mn[i,j] = i * (dΨ/dρ) * (m[i] - n*q_boundary) * ξ_ψ[i,j]
+    B_n = i * (dΨ/dρ) * (m - n*q) * ξ_ψ
 
-## Physical Interpretation:
+where ξ_ψ is the radial displacement eigenfunction.
+
+## 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()
@@ -128,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)
@@ -138,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/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl
index a3312682..39dca50a 100644
--- a/src/Vacuum/Kernel2D.jl
+++ b/src/Vacuum/Kernel2D.jl
@@ -178,14 +178,15 @@ 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: 𝒢ⁿ [Chance Phys. Plasmas 1997 2161 eq. 26-27]
                 if plasma_is_source
@@ -648,13 +649,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
@@ -670,29 +671,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

From 9d106935e24ed4c88b466200975802b78e10bb15 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Fri, 13 Feb 2026 16:48:39 -0500
Subject: [PATCH 3/3] VACUUM - BUGFIX - tricky fix that was caught by unit
 tests

---
 src/Vacuum/Vacuum.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 6e034e9e..5b57e880 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -339,7 +339,7 @@ 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)
+    wv = zeros(ComplexF64, mpert, mpert)
     xzpts = zeros(inputs.mtheta, 4)
     if !green_only
         # Perform inverse Fourier transforms to get response matrix components [Chance Phys. Plasmas 2007 052506 eq. 115-118]