From 5ae268cd422936b76de1dd70f49da27b980292e0 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 1 Jul 2026 12:22:33 -0400
Subject: [PATCH 1/4] PE - MINOR - cleaning up functions/comments as I work
 through the code

---
 src/PerturbedEquilibrium/Response.jl         |  10 +-
 src/PerturbedEquilibrium/ResponseMatrices.jl | 312 +++----------------
 2 files changed, 42 insertions(+), 280 deletions(-)

diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl
index ba2827246..1ca4d08a8 100644
--- a/src/PerturbedEquilibrium/Response.jl
+++ b/src/PerturbedEquilibrium/Response.jl
@@ -30,10 +30,10 @@ function compute_plasma_response!(
     end
 
     # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode]
-    flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, vac_data, ffs_intr)
+    flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, ffs_intr)
 
-    # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default)
-    plasma_inductance = calc_plasma_inductance(vac_data, ffs_intr, equil.psio)
+    # Compute plasma inductance
+    plasma_inductance = calc_plasma_inductance(ffs_intr, vac_data.wt0, equil.psio)
 
     # Surface inductance L from Green's functions at psilim.
     # Requires a 2D (nzvac=1) vacuum response so rows are theta points only,
@@ -45,7 +45,7 @@ function compute_plasma_response!(
     grre_2d = Matrix{ComplexF64}(grre_2d_raw)
     ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν
     surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac)
-    permeability = calc_permeability(plasma_inductance, surface_inductance)
+    permeability = plasma_inductance / surface_inductance
 
     # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats).
     # Λ (plasma inductance) is not Hermitian — its anti-Hermitian part is the dissipative/torque response — so the adjoint matters.
@@ -74,7 +74,7 @@ function compute_plasma_response!(
     # forcing arrives as Φ_x, the field reconstruction below consumes Φ_tot, and the b̃ spectra are
     # formed via the conform operator R = S·A (Φ = R·b̃).
     forcing_flux = map_forcing_to_eigenmodes(intr.forcing_modes, ffs_intr)
-    response_flux = compute_plasma_response_vector(permeability, forcing_flux)
+    response_flux = permeability * forcing_flux
 
     # Output forcing/response in the three Pharr field representations (all tesla):
     #   b̃ (root-area-weighted) = R⁻¹·Φ,  b (bare) = Σ⁻¹·b̃,  b̄ (area-weighted) = S·b̃.
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index a4b798a6d..f635bd95b 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -11,283 +11,72 @@ Reference: [Park Phys. Plasmas 2009 056115]
 using ..Utilities.FourierTransforms
 
 """
-    extract_boundary_displacements(
-        equil::Equilibrium.PlasmaEquilibrium,
-        ForceFreeStates_results::OdeState,
-        intr::ForceFreeStatesInternal
-    )::NamedTuple
-
-Extract eigenmode displacements and equilibrium quantities at the plasma boundary.
-
-This function extracts the data needed to compute the normal magnetic field at the
-plasma surface from ForceFreeStates eigenmode solutions.
-
-## What's extracted:
-
- 1. **Boundary displacement**: ξ_ψ from `u_store[:, :, 1, end]`
-
-      + This is the radial (normal) component of the eigenmode displacement
-      + At the last radial integration point (plasma edge)
-      + Dimensions: [numpert_total, numpert_total]
-
- 2. **Flux surface spacing**: dΨ/dρ at boundary
-
-      + From equilibrium bicubic spline evaluation
-      + Needed to convert displacement to magnetic field
-
- 3. **Safety factor**: q at boundary
-
-      + Used to compute singular factors (m - n*q)
-      + Identifies resonant surfaces
+    build_flux_matrix(equil, ffs_results, intr) -> Matrix{ComplexF64}
 
-## Arguments
+Calculate the vacuum poloidal flux matrix at the plasma boundary:
+bwp_mn[i,j] = 1im * (dΨ/dρ) * (m - n * q) * ξ_ψ(i,j)
 
-  - `equil`: Equilibrium solution containing flux surfaces and q-profile
-  - `ForceFreeStates_results`: ODE integration results containing u_store with eigenmodes
-  - `intr`: ForceFreeStates internal state with boundary location (psilim)
+Arguments:
 
-## Returns
+  - `equil`: Plasma equilibrium
+  - `ffs_results`: ForceFreeStates eigenmode ODE results
+  - `intr`: ForceFreeStates internal state with mode information
 
-Named tuple with:
+Returns:
 
-  - `ξ_psi_boundary`: Boundary displacement [numpert_total, numpert_total]
-  - `dPsi_drho`: Flux surface spacing at boundary (scalar)
-  - `q_boundary`: Safety factor at boundary (scalar)
-  - `psi_boundary`: Normalized flux at boundary (scalar)
+  - `bwp_mn`: Complex flux matrix [mode_i, eigenmode_j]
 """
-function extract_boundary_displacements(
-    equil::Equilibrium.PlasmaEquilibrium,
-    ForceFreeStates_results::OdeState,
-    intr::ForceFreeStatesInternal
-)
-    # Extract boundary displacement (normal component)
-    # u_store dimensions: [numpert_total, numpert_total, 2, numsteps]
-    # Index 1 in 3rd dimension is ξ_ψ (radial displacement)
-    # Last index in 4th dimension is the boundary
-    ξ_psi_boundary = ForceFreeStates_results.u_store[:, :, 1, ForceFreeStates_results.step]
-
-    # Get boundary location in normalized flux coordinates
-    psi_boundary = ForceFreeStates_results.psi_store[ForceFreeStates_results.step]
-
-    # Evaluate equilibrium quantities at boundary
-    # Safety factor at boundary
-    q_boundary = ForceFreeStates_results.q_store[ForceFreeStates_results.step]
-
-    # FFS ODE integrates in ψ (normalized flux), so bwp_mn = chi1·singfac·2πi·ξ_ψ
-    # where chi1 = 2π·psio  (Fortran idcon.f: chi1 = twopi*psio)
-    # Combined flux factor = chi1·2π = (2π)²·psio  (gpeq.f: bwp_mn = chi1·singfac·twopi·ifac·xsp)
-    dPsi_drho = (2π)^2 * equil.psio
+function build_flux_matrix(equil::Equilibrium.PlasmaEquilibrium, ffs_results::OdeState, intr::ForceFreeStatesInternal)::Matrix{ComplexF64}
 
-    return (
-        ξ_psi_boundary=ξ_psi_boundary,
-        dPsi_drho=dPsi_drho,
-        q_boundary=q_boundary,
-        psi_boundary=psi_boundary
-    )
-end
-
-"""
-    compute_normal_magnetic_field(
-        boundary_data::NamedTuple,
-        intr::ForceFreeStatesInternal
-    )::Matrix{ComplexF64}
-
-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. From the ideal MHD constraint [Park Phys. Plasmas 2009 056115 eq. 4]:
-
-    B_n = i * (dΨ/dρ) * (m - n*q) * ξ_ψ
-
-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: 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
+    # Extract boundary displacements and equilibrium quantities
+    ξ_psi_boundary = ffs_results.u_store[:, :, 1, ffs_results.step]
+    q_boundary = ffs_results.q_store[ffs_results.step]
+    dPsi_drho = (2π)^2 * equil.psio
 
-  - `boundary_data`: Output from extract_boundary_displacements()
-      + ξ_psi_boundary: Boundary displacement [numpert_total, numpert_total]
-      + dPsi_drho: Flux surface spacing at boundary (scalar)
-      + q_boundary: Safety factor at boundary (scalar)
-      + psi_boundary: Normalized flux at boundary (scalar)
-  - `intr`: ForceFreeStates internal state with mode arrays (mlow, mhigh, nlow, etc.)
+    # Compute singular factor for each Fourier mode
+    singfac = vec((intr.mlow:intr.mhigh) .- q_boundary .* (intr.nlow:intr.nhigh)')
 
-## Returns
-
-  - `bwp_mn[numpert_total, numpert_total]`: Normal magnetic field matrix where bwp_mn[i,j]
-    is the normal field of Fourier mode i in response to eigenmode j
-"""
-function compute_normal_magnetic_field(
-    boundary_data::NamedTuple,
-    intr::ForceFreeStatesInternal
-)::Matrix{ComplexF64}
-
-    numpert_total = intr.numpert_total
-    bwp_mn = zeros(ComplexF64, numpert_total, numpert_total)
-
-    # Extract boundary data
-    ξ_psi = boundary_data.ξ_psi_boundary
-    dPsi_drho = boundary_data.dPsi_drho
-    q_boundary = boundary_data.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)
-    for i in 1:numpert_total
-        m_mode = (i - 1) % intr.mpert + intr.mlow
-        n_mode = (i - 1) ÷ intr.mpert + intr.nlow
-        singfac[i] = m_mode - n_mode * q_boundary
-    end
-
-    # 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]
-        end
-    end
+    # Compute normal magnetic field at plasma boundary [Park Phys. Plasmas 2009 056115 eq. 4]
+    bwp_mn = 1im * dPsi_drho .* singfac .* ξ_psi_boundary
 
     return bwp_mn
 end
 
 """
-    build_flux_matrix(
-        equil::Equilibrium.PlasmaEquilibrium,
-        ForceFreeStates_results::OdeState,
-        vac_data::VacuumData,
-        intr::ForceFreeStatesInternal
-    )::Matrix{ComplexF64}
-
-Build vacuum poloidal flux matrix from ForceFreeStates eigenmode solutions.
-
-This extracts the vacuum flux response for each eigenmode at the plasma boundary.
-In GPEC, this comes from `bwp_mn` (boundary normal field) computed from eigenmode
-displacements.
-
-The flux matrix relates eigenmode displacements to vacuum poloidal flux:
-
- 1. Extract eigenmode displacement at plasma boundary from u_store
- 2. Compute normal magnetic field: B_ψ = i×(dΨ/dρ)×(m - n×q)×ξ_ψ
- 3. Result is flux[mode_i, eigenmode_j] = bwp_mn[i,j]
-
-## Arguments
-
-  - `equil`: Equilibrium solution containing flux surfaces and q-profile
-  - `ForceFreeStates_results`: ForceFreeStates ODE integration results containing eigenmodes
-  - `vac_data`: Vacuum response data from free boundary calculation
-  - `intr`: ForceFreeStates internal state with mode information
-
-## Returns
-
-  - `flxmats[numpert_total, numpert_total]`: Complex flux matrix where flxmats[i,j] is the
-    vacuum flux of mode i in response to eigenmode j
-"""
-function build_flux_matrix(
-    equil::Equilibrium.PlasmaEquilibrium,
-    ForceFreeStates_results::OdeState,
-    vac_data::VacuumData,
-    intr::ForceFreeStatesInternal
-)::Matrix{ComplexF64}
-
-    # Step 1: Extract boundary displacements and equilibrium quantities
-    boundary_data = extract_boundary_displacements(equil, ForceFreeStates_results, intr)
-
-    # Step 2: Compute normal magnetic field at plasma boundary
-    # This is the actual implementation of GPEC's bwp_mn calculation
-    # bwp_mn[i,j] = i * (dΨ/dρ) * (m[i] - n*q_boundary) * ξ_ψ[i,j]
-    flxmats = compute_normal_magnetic_field(boundary_data, intr)
-
-    return flxmats
-end
-
-"""
-    calc_plasma_inductance(
-        vac_data::VacuumData,
-        ffs_intr::ForceFreeStatesInternal,
-        psio::Float64
-    )::Matrix{ComplexF64}
-
-Calculate plasma inductance matrix Λ using the wt0-based energy formula
-(matches Fortran `gpresp_induct` with `resp_induct_flag=TRUE`).
-
-    Λ = inv(2·t₁·wt0·t₂)
+    calc_plasma_inductance(ffs_intr, wt0, psio) -> Matrix{ComplexF64}
 
-where t₁ = im/(χ₁·s_i·2π), t₂ = -im/(χ₁·s_j·2π), s_i = m_i - n·q_lim.
+Compute the plasma inductance Λ from the displacement-space energy matrix ``W_0 = W_p + W_v``
+(Fortran `gpresp.f`, `resp_induct_flag=TRUE`):
 
-Note: `vac_data.wt0` already contains singfac² factors (s_i·s_j) baked into the
-vacuum term via the scaling in `free_run!`. The t₁/t₂ factors divide by s_i·s_j,
-correctly recovering the properly-normalized inductance.
+ 1. Normalize `wt0` to SI units (`ψ₀²/(μ₀·2)`, as in `idcon.f`).
+ 2. Convert ξ → Φ via ``Λ⁻¹ = 2·T₁·W₀·T₂`` with ``T_{ii} ∝ 1/(χ₁ s_i 2π)`` and
+    ``s_i = m_i - n·q_lim``; the ``s_i s_j`` factor undoes the `wv` scaling in `free_run!`.
+ 3. Invert to obtain Λ
 
 ## Arguments
 
-  - `vac_data`: Vacuum data containing wt0 (total energy matrix before eigenvector sorting)
-  - `ffs_intr`: ForceFreeStates internal state with mode info (mlow, mpert, nlow, qlim)
-  - `psio`: Total toroidal flux [Wb/rad] from equilibrium (equil.psio)
+  - `ffs_intr`: ForceFreeStates internal state (`mlow`, `mhigh`, `nlow`, `nhigh`, `qlim`)
+  - `wt0`: Total energy matrix `wp + wv` [numpert_total × numpert_total]
+  - `psio`: Toroidal flux [Wb/rad]
 
 ## Returns
 
-  - Plasma inductance matrix Lambda [numpert_total × numpert_total]
+  - Plasma inductance Λ matrix [numpert_total × numpert_total]
 """
-function calc_plasma_inductance(
-    vac_data::VacuumData,
-    ffs_intr::ForceFreeStatesInternal,
-    psio::Float64
-)::Matrix{ComplexF64}
-
-    mpert = ffs_intr.numpert_total
-    chi1 = 2π * psio          # = Fortran's chi1 = twopi*psio
-    n = ffs_intr.nlow
-    qlim = ffs_intr.qlim      # q at psilim
+function calc_plasma_inductance(ffs_intr::ForceFreeStatesInternal, wt0::Matrix{ComplexF64}, psio::Float64)::Matrix{ComplexF64}
 
     # Singular factors s_i = m_i - n*qlim  (same as Fortran: mfac(i) - nn*qlim)
-    s = [((i-1) % ffs_intr.mpert + ffs_intr.mlow) - n * qlim for i in 1:mpert]
+    singfac = vec((ffs_intr.mlow:ffs_intr.mhigh) .- ffs_intr.qlim .* (ffs_intr.nlow:ffs_intr.nhigh)')
 
-    # Fortran idcon_norm: wt0 = wt0/(mu0*2)*psio^2
-    # Julia's vac_data.wt0 is raw wp+wv; Fortran additionally scales by psio^2/(mu0*2)
+    # Convert to metric units
     mu0 = 4π * 1e-7
-    wt0_norm = vac_data.wt0 .* (psio^2 / (mu0 * 2))
-
-    # Build temp2[i,j] = 2·t1_i·wt0[i,j]·t2_j (matches Fortran gpresp_induct)
-    temp2 = Matrix{ComplexF64}(undef, mpert, mpert)
-    for i in 1:mpert, j in 1:mpert
-        t1 = im / (chi1 * s[i] * 2π)
-        t2 = -im / (chi1 * s[j] * 2π)
-        temp2[i, j] = 2 * t1 * wt0_norm[i, j] * t2
-    end
-
-    return inv(temp2)
-end
-
-"""
-    calc_permeability(
-        plasma_inductance::Matrix{ComplexF64},
-        surface_inductance::Matrix{ComplexF64}
-    )::Matrix{ComplexF64}
+    wt0_norm = wt0 .* (psio^2 / (mu0 * 2))
 
-Calculate permeability matrix P = Λ·L⁻¹ (matches Fortran `gpresp_permeab`).
+    # Convert from displacement to flux space using 1 / (singfac * chi1 * 2π) factor
+    chi1 = 2π * psio
+    wt0_norm .*= 2 ./ (2π * chi1)^2 ./ (singfac' .* singfac)
 
-## Arguments
-
-  - `plasma_inductance`: Plasma inductance matrix Lambda
-  - `surface_inductance`: Surface inductance matrix L
-
-## Returns
-
-  - Permeability matrix P = Lambda * L^{-1} [mpert, mpert]
-"""
-function calc_permeability(
-    plasma_inductance::Matrix{ComplexF64},
-    surface_inductance::Matrix{ComplexF64}
-)::Matrix{ComplexF64}
-    # P = Lambda * L^{-1}  (right-division solves for P s.t. P*L = Lambda)
-    return plasma_inductance / surface_inductance
+    return inv(wt0_norm)
 end
 
 """
@@ -424,30 +213,3 @@ function map_forcing_to_eigenmodes(
 
     return forcing_vector
 end
-
-"""
-    compute_plasma_response_vector(
-        permeability::Matrix{ComplexF64},
-        forcing_vector::Vector{ComplexF64}
-    )::Vector{ComplexF64}
-
-Compute plasma response to external forcing.
-
-Response = Permeability * Forcing
-
-## Arguments
-
-  - `permeability`: Permeability matrix
-  - `forcing_vector`: External forcing in eigenmode basis
-
-## Returns
-
-  - Plasma response vector in eigenmode basis
-"""
-function compute_plasma_response_vector(
-    permeability::Matrix{ComplexF64},
-    forcing_vector::Vector{ComplexF64}
-)::Vector{ComplexF64}
-
-    return permeability * forcing_vector
-end

From eb0bc12c80b39511ab46f07558d2c2ae361e98ee Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 1 Jul 2026 14:30:31 -0400
Subject: [PATCH 2/4] PE - MINOR - some more cosmetic cleanups

---
 src/GeneralizedPerturbedEquilibrium.jl        | 31 +++++------------
 .../PerturbedEquilibrium.jl                   | 11 +++----
 src/PerturbedEquilibrium/Response.jl          | 33 +++++++------------
 src/PerturbedEquilibrium/ResponseMatrices.jl  |  3 +-
 src/PerturbedEquilibrium/SingularCoupling.jl  | 25 +++++---------
 5 files changed, 35 insertions(+), 68 deletions(-)

diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 51bd422b4..7c2379efc 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -396,12 +396,6 @@ function main_from_inputs(
                 find_kinetic_singular_surfaces!(ffit, equil, intr)
             end
         end
-
-        # NOTE: Asymptotic calculations for ideal ForceFreeStates are now computed on-demand during
-        # singular surface crossings in cross_ideal_singular_surf!. This makes it clear that
-        # asymptotics are only needed for ideal ForceFreeStates and are not inherent properties of
-        # the singular surface.
-
     end
 
     # Integrate Euler-Lagrange Equation
@@ -422,12 +416,10 @@ function main_from_inputs(
             @info "Computing free boundary energies ($wall_desc)"
         end
         vac_data = free_run!(odet, ctrl, equil, ffit, intr)
-        if real(vac_data.et[1]) < 0
-            if ctrl.verbose
+        if ctrl.verbose
+            if real(vac_data.et[1]) < 0
                 @warn "Free-boundary mode unstable for n = $nstring"
-            end
-        else
-            if ctrl.verbose
+            else
                 @info "All free-boundary modes stable for n = $nstring"
             end
         end
@@ -479,6 +471,12 @@ function main_from_inputs(
         return
     end
 
+    # No perturbed equilibrium calculations if free-boundary mode is unstable
+    if real(vac_data.et[1]) < 0
+        @warn "Since a free-boundary mode is unstable, perturbed equilibrium calculations will not run."
+        return
+    end
+
     # ----------------------------------------------------------------
     # Perturbed Equilibrium
     # ----------------------------------------------------------------
@@ -585,15 +583,8 @@ function main_from_inputs(
         @info "KineticForces completed in $(@sprintf("%.3f", time() - kf_start)) s"
     end
 
-    # ----------------------------------------------------------------
-    # Done
-    # ----------------------------------------------------------------
     @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER"
 
-    # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found
-
-    return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, vac_data=ctrl.vac_flag ? vac_data : nothing)
-
 end
 
 """
@@ -604,10 +595,6 @@ This combines the functionality of several pieces of the Fortran code in `ode_ou
 primarily `ode_output_open` and the various `bin_euler` writes that occur throughout the
 integration. Some parameters are only dumped in their respective flags are true, e.g.
 vacuum data if `vac_flag` is true.
-
-### TODOs
-
-Combine spline unpacking if possible, too many extra lines
 """
 function write_outputs_to_HDF5(
     ctrl::ForceFreeStatesControl,
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
index 5f1fe83b5..669c51835 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
@@ -1,5 +1,7 @@
 module PerturbedEquilibrium
 
+const μ0 = 4π * 1e-7
+
 # Imports
 using HDF5
 using Printf
@@ -78,7 +80,7 @@ function compute_perturbed_equilibrium(
 
     state = PerturbedEquilibriumState()
 
-    # Step 0: Initialize mode arrays for convenient indexing
+    # Initialize mode arrays for convenient indexing
     initialize_mode_arrays!(intr, ffs_intr)
 
     # Load forcing data. On the gpec.h5 replay path the caller preloads
@@ -113,7 +115,7 @@ function compute_perturbed_equilibrium(
         end
     end
 
-    # Step 2: Compute plasma response
+    # Compute plasma response
     if ctrl.compute_response
         if vac_data === nothing
             @warn "Vacuum data not available. Skipping plasma response calculation. Set vac_flag=true in [ForceFreeStates] section."
@@ -122,7 +124,7 @@ function compute_perturbed_equilibrium(
         end
     end
 
-    # Step 3: Compute singular coupling metrics
+    # Compute singular coupling metrics
     if ctrl.compute_singular_coupling
         if vac_data === nothing
             @warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section."
@@ -131,9 +133,6 @@ function compute_perturbed_equilibrium(
         end
     end
 
-    # Step 4: Output eigenmode fields (integrated into HDF5 output)
-    # This is handled by write_outputs_to_HDF5 in main()
-
     return state
 end
 
diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl
index 1ca4d08a8..b4e8938f4 100644
--- a/src/PerturbedEquilibrium/Response.jl
+++ b/src/PerturbedEquilibrium/Response.jl
@@ -25,9 +25,7 @@ function compute_plasma_response!(
     metric::MetricData,
     ffit::FourFitVars
 )
-    if ctrl.verbose
-        @info "Computing plasma response (wt0-based inductance)"
-    end
+    ctrl.verbose && @info "Computing plasma response"
 
     # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode]
     flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, ffs_intr)
@@ -35,30 +33,25 @@ function compute_plasma_response!(
     # Compute plasma inductance
     plasma_inductance = calc_plasma_inductance(ffs_intr, vac_data.wt0, equil.psio)
 
-    # Surface inductance L from Green's functions at psilim.
-    # Requires a 2D (nzvac=1) vacuum response so rows are theta points only,
+    # Surface inductance L from Green's functions at psilim - requires no wall limit
     nn = ffs_intr.nlow
     vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn])
     wall_nowall = Vacuum.WallShapeSettings(; shape="nowall")
-    _, grri_2d_raw, grre_2d_raw, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
-    grri_2d = Matrix{ComplexF64}(grri_2d_raw)
-    grre_2d = Matrix{ComplexF64}(grre_2d_raw)
+    _, grri_2d, grre_2d, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
     ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν
     surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac)
+
+    # Compute permeability P = Λ·L⁻¹ and store in internal state for singular coupling / field reconstruction.
     permeability = plasma_inductance / surface_inductance
+    intr.plasma_response = permeability
 
-    # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats).
-    # Λ (plasma inductance) is not Hermitian — its anti-Hermitian part is the dissipative/torque response — so the adjoint matters.
+    # Compute reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹
+    # Λ is not Hermitian — its anti-Hermitian part is the dissipative/torque response — so the adjoint matters.
     L_inv = inv(surface_inductance)
     reluctance = L_inv * (plasma_inductance' - surface_inductance) * L_inv
 
-    # Store permeability in internal state for singular coupling / field reconstruction.
-    # These consumers operate on the physical control-surface flux Φ_x, so the internal
-    # copy stays in flux space; only the stored/output quantities are conformed to fields below.
-    intr.plasma_response = permeability
-
     # Conform the control-surface matrices to the coordinate-invariant root-area-weighted
-    # field (b̃) space for output (issue #233 / Pharr 2026). Store the b̃→b̄ operator S = Σ/√A
+    # field (b̃) space for output (Pharr 2026). Store the b̃→b̄ operator S = Σ/√A
     # and the scalar surface area A so users can recover the area-weighted field (b̄ = S·b̃) or
     # flux (Φ = A·b̄) — see Utils.jl output docs.
     rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs_intr)
@@ -70,9 +63,7 @@ function compute_plasma_response!(
     state.rootarea_to_area_weight = rootarea_to_area_weight
     state.surface_area = surface_area
 
-    # Forcing and response on the control surface. Flux Φ appears only as a brief internal bridge:
-    # forcing arrives as Φ_x, the field reconstruction below consumes Φ_tot, and the b̃ spectra are
-    # formed via the conform operator R = S·A (Φ = R·b̃).
+    # Compute actual flux from external flux and permiability, Φ = P Φ^x
     forcing_flux = map_forcing_to_eigenmodes(intr.forcing_modes, ffs_intr)
     response_flux = permeability * forcing_flux
 
@@ -119,7 +110,5 @@ function compute_plasma_response!(
     state.b_n_modes = b_n_modes
     state.xi_n_modes = xi_n_modes
 
-    if ctrl.verbose
-        @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))"
-    end
+    ctrl.verbose && @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))"
 end
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index f635bd95b..8319928a7 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -69,8 +69,7 @@ function calc_plasma_inductance(ffs_intr::ForceFreeStatesInternal, wt0::Matrix{C
     singfac = vec((ffs_intr.mlow:ffs_intr.mhigh) .- ffs_intr.qlim .* (ffs_intr.nlow:ffs_intr.nhigh)')
 
     # Convert to metric units
-    mu0 = 4π * 1e-7
-    wt0_norm = wt0 .* (psio^2 / (mu0 * 2))
+    wt0_norm = wt0 .* (psio^2 / (μ0 * 2))
 
     # Convert from displacement to flux space using 1 / (singfac * chi1 * 2π) factor
     chi1 = 2π * psio
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index eb41cd76d..89a003585 100644
--- a/src/PerturbedEquilibrium/SingularCoupling.jl
+++ b/src/PerturbedEquilibrium/SingularCoupling.jl
@@ -361,7 +361,6 @@ function compute_current_density(
     psi::Float64
 )::Float64
     # Physical constants
-    μ₀ = 4π * 1e-7
     chi1 = 2π * equil.psio
     twopi = 2π
 
@@ -412,8 +411,8 @@ function compute_current_density(
     # Trapezoidal rule end correction (subtract half of last point contribution)
     integral -= last_jac * last_delpsi * last_jcfun / mthsurf
 
-    # Final normalization: j_c = (1/integral) * χ₁² * q / μ₀
-    j_c = (1.0 / integral) * chi1^2 * q / μ₀
+    # Final normalization: j_c = (1/integral) * χ₁² * q / μ0
+    j_c = (1.0 / integral) * chi1^2 * q / μ0
 
     return j_c
 end
@@ -456,11 +455,9 @@ Surface inductance matrix [mpert × mpert]
 )::Matrix{ComplexF64}
     mpert = ffs_intr.mpert
     mtheta = length(ν)
-    μ₀ = 4π * 1e-7
 
     ft = FourierTransforms.FourierTransform(mtheta, mpert, ffs_intr.mlow)
 
-    flux_matrix = zeros!(pool, ComplexF64, mpert, mpert)
     current_matrix = zeros!(pool, ComplexF64, mpert, mpert)
 
     kax = zeros!(pool, ComplexF64, mtheta)
@@ -471,17 +468,15 @@ Surface inductance matrix [mpert × mpert]
     phase = cis.(-nn .* ν)
 
     for i in 1:mpert
-        flux_matrix[i, i] = 1.0
-
         # Complex grri/e stores exp(i(mθ-nν)) projection, need conjugate for exp(-i(mθ-nν))
-        kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ₀ * (2π)^2)
+        kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ0 * (2π)^2)
 
-        # Port of Fortran gpvacuum_flxsurf: apply toroidal phase, reverse theta, forward-DFT.
+        # Apply toroidal phase, reverse theta, forward-DFT.
         g_phased = kax .* phase
         current_matrix[:, i] = ft(_reverse_theta(g_phased))
     end
 
-    # Compute surface inductance: L_surf = flux * inv(current) = inv(current)
+    # Compute surface inductance: L_surf = flux * inv(current) = inv(current) when flux is the identity matrix
     L_surf = zeros(ComplexF64, mpert, mpert)
 
     current_mag = maximum(abs.(current_matrix))
@@ -489,21 +484,19 @@ Surface inductance matrix [mpert × mpert]
     if current_mag < 1e-15
         @warn "Current matrix is all zeros! Cannot compute surface inductance." maxlog=1
         for i in 1:mpert
-            L_surf[i, i] = μ₀ * 1e-6
+            L_surf[i, i] = μ0 * 1e-6
         end
     else
         try
             regularization = 1e-12 * current_mag
             current_reg = current_matrix + regularization * I
 
-            L_surf = flux_matrix * inv(current_reg)
-
-            # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1))))
-            L_surf = 0.5 * (L_surf + L_surf')
+            L_surf = inv(current_reg)
+            hermitianpart!(L_surf)
         catch e
             @warn "Surface inductance inversion failed: $e" maxlog=1
             for i in 1:mpert
-                L_surf[i, i] = μ₀ * 1e-6
+                L_surf[i, i] = μ0 * 1e-6
             end
         end
     end

From 629f62d298cbd940ebd8a7370ccb82512ad05b1b Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Fri, 10 Jul 2026 08:53:26 -0400
Subject: [PATCH 3/4] PE - MINOR - more cleanups

---
 src/GeneralizedPerturbedEquilibrium.jl        |  18 +++
 .../PerturbedEquilibrium.jl                   |  12 +-
 src/PerturbedEquilibrium/Response.jl          |   6 +-
 src/PerturbedEquilibrium/ResponseMatrices.jl  |  18 +--
 src/PerturbedEquilibrium/SingularCoupling.jl  | 105 +++---------------
 5 files changed, 40 insertions(+), 119 deletions(-)

diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 7c2379efc..ec8a1462b 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -471,12 +471,30 @@ function main_from_inputs(
         return
     end
 
+    # No perturbed equilibrium calculations if vacuum data is not available
+    if !ctrl.vac_flag
+        @warn "Vacuum data not available. Skipping perturbed equilibrium calculations. Set vac_flag=true in [ForceFreeStates] section."
+        return
+    end
+
     # No perturbed equilibrium calculations if free-boundary mode is unstable
     if real(vac_data.et[1]) < 0
         @warn "Since a free-boundary mode is unstable, perturbed equilibrium calculations will not run."
         return
     end
 
+    # No perturbed equilibrium calculations if wall shape is not nowall
+    if intr.wall_settings.shape !== "nowall"
+        @warn "Perturbed equilibrium calculations are only supported for nowall cases - skipping perturbed equilibrium calculations."
+        return
+    end
+
+    # Perturbed equilibrium calculations still need multi-n support added
+    if intr.npert > 1
+        @warn "Perturbed equilibrium calculations are only supported for single-n cases - skipping perturbed equilibrium calculations."
+        return
+    end
+
     # ----------------------------------------------------------------
     # Perturbed Equilibrium
     # ----------------------------------------------------------------
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
index 669c51835..6fe7465a9 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
@@ -117,20 +117,12 @@ function compute_perturbed_equilibrium(
 
     # Compute plasma response
     if ctrl.compute_response
-        if vac_data === nothing
-            @warn "Vacuum data not available. Skipping plasma response calculation. Set vac_flag=true in [ForceFreeStates] section."
-        else
-            compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit)
-        end
+        compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit)
     end
 
     # Compute singular coupling metrics
     if ctrl.compute_singular_coupling
-        if vac_data === nothing
-            @warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section."
-        else
-            compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl)
-        end
+        compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl)
     end
 
     return state
diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl
index b4e8938f4..c9ac7e40d 100644
--- a/src/PerturbedEquilibrium/Response.jl
+++ b/src/PerturbedEquilibrium/Response.jl
@@ -33,13 +33,11 @@ function compute_plasma_response!(
     # Compute plasma inductance
     plasma_inductance = calc_plasma_inductance(ffs_intr, vac_data.wt0, equil.psio)
 
-    # Surface inductance L from Green's functions at psilim - requires no wall limit
+    # Surface inductance L from Green's functions
     nn = ffs_intr.nlow
     vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn])
-    wall_nowall = Vacuum.WallShapeSettings(; shape="nowall")
-    _, grri_2d, grre_2d, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
     ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν
-    surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac)
+    surface_inductance = compute_surface_inductance_from_greens(vac_data.grri, vac_data.grre, ffs_intr, nn, ν_vac)
 
     # Compute permeability P = Λ·L⁻¹ and store in internal state for singular coupling / field reconstruction.
     permeability = plasma_inductance / surface_inductance
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index 8319928a7..fcca26d9b 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -192,21 +192,11 @@ function map_forcing_to_eigenmodes(
     numpert_total = intr.mpert * intr.npert
     forcing_vector = zeros(ComplexF64, numpert_total)
 
-    # Create mode index map: (m,n) -> linear index
+    # Combine all forcing modes into a single vector
     for forcing_mode in forcing_modes
-        # Find matching mode in eigenmode basis
-        for i in 1:numpert_total
-            # Calculate m and n for this index
-            # Using 0-based indexing converted to 1-based:
-            # m = (i-1) % mpert + mlow
-            # n = (i-1) ÷ mpert + nlow
-            m_mode = (i - 1) % intr.mpert + intr.mlow
-            n_mode = (i - 1) ÷ intr.mpert + intr.nlow
-
-            if m_mode == forcing_mode.m && n_mode == forcing_mode.n
-                forcing_vector[i] = forcing_mode.amplitude
-                break
-            end
+        i = (forcing_mode.m - intr.mlow) + (forcing_mode.n - intr.nlow) * intr.mpert + 1
+        if 1 <= i <= numpert_total
+            forcing_vector[i] = forcing_mode.amplitude
         end
     end
 
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index 89a003585..1c6b4b784 100644
--- a/src/PerturbedEquilibrium/SingularCoupling.jl
+++ b/src/PerturbedEquilibrium/SingularCoupling.jl
@@ -99,7 +99,6 @@ function compute_singular_coupling_metrics!(
     end
 
     chi1 = 2π * equil.psio
-    twopi = 2π
     mtheta = vac_data.mthvac
     wall_settings = Vacuum.WallShapeSettings(; shape="nowall")
 
@@ -162,9 +161,7 @@ function compute_singular_coupling_metrics!(
 
         # Compute Green's functions at this surface for this n (once per pair)
         vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn])
-        _, grri_raw, grre_raw, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings)
-        grri = Matrix{ComplexF64}(grri_raw)
-        grre = Matrix{ComplexF64}(grre_raw)
+        _, grri, grre, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings)
         ffs_intr.sing[s].grri = grri
         ffs_intr.sing[s].grre = grre
 
@@ -177,11 +174,8 @@ function compute_singular_coupling_metrics!(
         L_mm = L_surf[m_idx, m_idx]
 
         j_c = compute_current_density(equil, sing_surf.psifac)
-        area = compute_surface_area(equil, sing_surf.psifac)
-        # Matches Fortran gpout_resp: shear = m*dq/dψ / q² = n*dq/dψ / q (since m=n*q).
-        # Uses abs(nn) because island_half_width = sqrt(abs(island_width_sq)), so the sign
-        # of shear only affects the sign of C_island_width_sq, not the physical island width.
-        shear = abs(nn) * sing_surf.q1 / sing_surf.q
+        area = Equilibrium.flux_surface_area(equil, sing_surf.psifac, length(equil.rzphi_ys) - 1)
+        shear = nn * sing_surf.q1 / sing_surf.q # m*dq/dψ / q² = n*dq/dψ / q (since m=n*q)
 
         # Evaluate bwp1_mn = ∂b^ψ/∂ψ at lpsi and rpsi using permeability-weighted eigenstates.
         # Matches Fortran gpout_resp: evaluate bwp1_mn at lpsi/rpsi via gpeq_sol
@@ -254,14 +248,14 @@ function compute_singular_coupling_metrics!(
         #  - resonant (shielding) current: j_c already integrates jac·|∇ψ| over the surface, so the
         #    Jacobian weighting is carried inside j_c — no separate area factor needed.
         #  - resonant flux → field: Φ^r/A^r [T], invariant [Park 2008; Pharr 2026].
-        state.C_delta_prime[row, :] = jump_vec ./ (twopi * chi1)
-        state.C_resonant_current[row, :] = jump_vec .* (-j_c / (twopi * m_res))
+        state.C_delta_prime[row, :] = jump_vec ./ (2π * chi1)
+        state.C_resonant_current[row, :] = jump_vec .* (-j_c / (2π * m_res))
         # Matches Fortran gpout_resp: singflx = L·fkaxmn, resonant area-weighted field = singflx/area,
         # islandhwids = 4·singflx/(2π·shear·q·chi1)
-        singflx_pre = (L_mm / (twopi * nn)) .* state.C_resonant_current[row, :]
+        singflx_pre = (L_mm / (2π * nn)) .* state.C_resonant_current[row, :]
         state.C_resonant_area_weighted_field[row, :] = singflx_pre ./ area
         if abs(shear) > 1e-10
-            state.C_island_width_sq[row, :] = (4.0 / (twopi * shear * sing_surf.q * chi1)) .* singflx_pre
+            state.C_island_width_sq[row, :] = abs.(4.0 / (2π * shear * sing_surf.q * chi1)) .* singflx_pre
         end
 
         state.rational_psi[row] = sing_surf.psifac
@@ -329,24 +323,6 @@ j_c = χ₁² * q / (μ₀ * integral)
 where the integral is computed via flux surface integration:
 integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ
 
-## GPEC Formula
-
-```fortran
-DO itheta=0,mthsurf
-   CALL bicube_eval(rzphi,respsi,theta(itheta),1)
-   rfac=SQRT(rzphi%f(1))
-   jac=rzphi%f(4)
-   w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac
-   w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac)
-   delpsi(itheta)=SQRT(w(1,1)**2+w(1,2)**2)
-   sqreqb(itheta)=(sq%f(1)**2+chi1**2*delpsi(itheta)**2)/(twopi*r(itheta))**2
-   jcfun(itheta)=sqreqb(itheta)/(delpsi(itheta)**3)
-   j_c(ising)=j_c(ising)+jac*delpsi(itheta)*jcfun(itheta)/mthsurf
-ENDDO
-j_c(ising)=j_c(ising)-jac*delpsi(mthsurf)*jcfun(mthsurf)/mthsurf  ! trapezoidal rule
-j_c(ising)=1.0/j_c(ising)*chi1**2*sq%f(4)/mu0
-```
-
 ## Implementation
 
 Uses trapezoidal rule integration around the flux surface with metric quantities
@@ -362,14 +338,11 @@ function compute_current_density(
 )::Float64
     # Physical constants
     chi1 = 2π * equil.psio
-    twopi = 2π
 
     # Get equilibrium quantities at this surface
     F_tor = equil.profiles.F_spline(psi)  # Toroidal field function (2π·R·B_tor in GPEC convention)
     q = equil.profiles.q_spline(psi)      # Safety factor
 
-    ro = equil.ro
-
     # Number of theta points for integration
     # Match GPEC's mthsurf (typically 101 points from theta=0 to theta=1)
     mthsurf = length(equil.rzphi_ys) - 1
@@ -392,7 +365,7 @@ function compute_current_density(
         delpsi = m.delpsi  # flux gradient magnitude |∇ψ|
 
         # sqreqb = (F² + χ₁²|∇ψ|²) / (2πR)²  where F = R·B_tor (Fortran sq%f(1))
-        sqreqb = (F_tor^2 + chi1^2 * delpsi^2) / (twopi * m.r)^2
+        sqreqb = (F_tor^2 + chi1^2 * delpsi^2) / (2π * m.r)^2
 
         # Integrand function
         jcfun = sqreqb / (delpsi^3)
@@ -469,10 +442,12 @@ Surface inductance matrix [mpert × mpert]
 
     for i in 1:mpert
         # Complex grri/e stores exp(i(mθ-nν)) projection, need conjugate for exp(-i(mθ-nν))
+        # Eq. 10 of Park 2007
         kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ0 * (2π)^2)
 
         # Apply toroidal phase, reverse theta, forward-DFT.
         g_phased = kax .* phase
+        # Eq. 21b of Park 2007
         current_matrix[:, i] = ft(_reverse_theta(g_phased))
     end
 
@@ -487,67 +462,15 @@ Surface inductance matrix [mpert × mpert]
             L_surf[i, i] = μ0 * 1e-6
         end
     else
-        try
-            regularization = 1e-12 * current_mag
-            current_reg = current_matrix + regularization * I
-
-            L_surf = inv(current_reg)
-            hermitianpart!(L_surf)
-        catch e
-            @warn "Surface inductance inversion failed: $e" maxlog=1
-            for i in 1:mpert
-                L_surf[i, i] = μ0 * 1e-6
-            end
-        end
+        # Add a small regularization to the current matrix to avoid division by zero
+        current_matrix += 1e-12 * current_mag * I
+        L_surf .= inv(current_matrix)
+        hermitianpart!(L_surf)
     end
 
     return L_surf
 end
 
-"""
-    compute_surface_area(
-        equil::Equilibrium.PlasmaEquilibrium,
-        psi::Float64
-    )::Float64
-
-Compute flux surface area at given ψ.
-
-Implements GPEC's area calculation (Fortran `gpout_respinfo`):
-area = ∫ jac * |∇ψ| dθ
-
-where the integral is computed around the flux surface.
-
-## GPEC Formula
-
-```fortran
-DO itheta=0,mthsurf
-   CALL bicube_eval(rzphi,respsi,theta(itheta),1)
-   rfac=SQRT(rzphi%f(1))
-   jac=rzphi%f(4)
-   w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac
-   w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac)
-   delpsi(itheta)=SQRT(w(1,1)**2+w(1,2)**2)
-   area(ising)=area(ising)+jac*delpsi(itheta)/mthsurf
-ENDDO
-area(ising)=area(ising)-jac*delpsi(mthsurf)/mthsurf  ! trapezoidal rule
-```
-
-## Implementation
-
-Uses trapezoidal rule integration around the flux surface with:
-
-  - jac: Jacobian of flux coordinates from rzphi
-  - |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor
-"""
-function compute_surface_area(
-    equil::Equilibrium.PlasmaEquilibrium,
-    psi::Float64
-)::Float64
-    # mthsurf matches GPEC's flux-surface theta resolution
-    mthsurf = length(equil.rzphi_ys) - 1
-    return Equilibrium.flux_surface_area(equil, psi, mthsurf)
-end
-
 """
     compute_island_diagnostics!(state::PerturbedEquilibriumState, n_rational::Int)
 

From a9edf4a31d74ce50078291f597d537963934ab93 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 30 Jul 2026 08:38:03 -0400
Subject: [PATCH 4/4] PE - MINOR - improving plasma inductance comment

---
 src/PerturbedEquilibrium/ResponseMatrices.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index fcca26d9b..2298916b3 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -72,6 +72,7 @@ function calc_plasma_inductance(ffs_intr::ForceFreeStatesInternal, wt0::Matrix{C
     wt0_norm = wt0 .* (psio^2 / (μ0 * 2))
 
     # Convert from displacement to flux space using 1 / (singfac * chi1 * 2π) factor
+    # The factor of 2 comes from the 1/2 + inverse in eq. 31 of Park Phys. Plasmas 2007
     chi1 = 2π * psio
     wt0_norm .*= 2 ./ (2π * chi1)^2 ./ (singfac' .* singfac)