Citations

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

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

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]

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)
 

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.

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

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]

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: