From 0ad362d002229a143f8242ef5b5c25381d451c71 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 12 Aug 2026 14:41:26 -0400
Subject: [PATCH 1/5] ForceFreeStates - REFACTOR - Reduce whole-struct passing
 and in-place struct mutation

First installment against #139: make it clear from a call site what a function
changes, by shrinking what the shared mutable structs carry and narrowing who is
allowed to write to them. Further struct cleanups will land on this branch.

Freeze ForceFreeStatesControl:
- ForceFreeStatesControl is now immutable, so ctrl can never appear as the first
  argument of a ! function. ctrl.nn_low/nn_high keep what the user asked for; the
  resolved toroidal range lives on intr.nlow/nhigh, and the Fortran DCON
  delta_mhigh doubling is applied where it is consumed rather than written back
  onto ctrl 130 lines earlier.
- The Riccati dense-xi pass no longer clobbers and restores three ctrl flags
  around a nested integration. The serial branch of eulerlagrange_integration is
  extracted as serial_eulerlagrange_integration and called directly, which also
  removes a mutation of shared state from a function reachable on the threaded
  path.
- sing_lim! normalizes dmlim into a local instead of writing it back to ctrl.

Return values instead of parking them on intr:
- Local stability is returned by the new compute_local_stability rather than
  stored on intr.locstab. The field is removed, along with the zeros-filled
  spline that was built on every run even with local_stability_flag = false;
  write_outputs_to_HDF5 takes locstab as a keyword.

Docstrings:
- eulerlagrange_integration and parallel_eulerlagrange_integration both
  documented "-> OdeState" while returning a 4-tuple. Corrected, and dropped
  stale references to euler.h5 and to setting ctrl fields programmatically.

No numerical change is intended: every edit is a data-flow refactor.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/ForceFreeStates/Ballooning.jl             | 13 ++++
 src/ForceFreeStates/EulerLagrange.jl          | 49 ++++++------
 src/ForceFreeStates/ForceFreeStatesStructs.jl | 11 ++-
 src/ForceFreeStates/Riccati.jl                | 75 +++++++------------
 src/ForceFreeStates/Sing.jl                   | 14 ++--
 src/GeneralizedPerturbedEquilibrium.jl        | 75 ++++++++-----------
 test/runtests_eulerlagrange.jl                | 13 +---
 test/runtests_sing.jl                         |  9 +--
 8 files changed, 122 insertions(+), 137 deletions(-)

diff --git a/src/ForceFreeStates/Ballooning.jl b/src/ForceFreeStates/Ballooning.jl
index d24f9ff81..73492a673 100644
--- a/src/ForceFreeStates/Ballooning.jl
+++ b/src/ForceFreeStates/Ballooning.jl
@@ -78,6 +78,19 @@ function compute_ballooning_stability!(
 
 end
 
+"""
+    compute_local_stability(ctrl, plasma_eq) -> CubicSeriesInterpolant
+
+Local stability profile spline over `plasma_eq.profiles.xs`, with the columns filled by
+[`compute_ballooning_stability!`](@ref): 1 = `D_I·ψ`, 2 = `D_R·ψ`, 4 = ballooning `Δ'`.
+"""
+function compute_local_stability(ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium)
+    xs = plasma_eq.profiles.xs
+    locstab_fs = zeros(Float64, length(xs), 5)
+    compute_ballooning_stability!(ctrl, locstab_fs, plasma_eq)
+    return cubic_interp(xs, Series(locstab_fs); extrap=ExtendExtrap())
+end
+
 """
     resistive_interchange_h(flux_surface_index, plasma_eq)
 
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index 1f8fe3bca..631fea55a 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -137,25 +137,16 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc
 end
 
 """
-    eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
-
-Main driver for integrating the Euler-Lagrange equations across the plasma and detecting singular surfaces.
-Formerly `ode_run`. Has the same functionality as `ode_run` in the Fortran code, with the addition of
-a single dump to the `euler.h5` file at the end of integration instead of multiple dumps
-to `euler.bin` throughout the integration. We have made the control logic more clear
-by pre-computing all integration chunks upfront and using a for loop to iterate through them,
-eliminating the while-loop logic and making integration bounds explicit at each step.
-We now perform significant post-processing after integration including finding the peak dW
-in the edge region and evaluating the stability criterion over the entire integration,
-which were previously done during integration in the Fortran code.
+    eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left)
 
-### TODOs
-
-restype functionality if we decide to do this
+Integrate the Euler-Lagrange equations from the axis to `intr.psilim`, crossing each singular
+surface on the way (Fortran `ode_run`). Dispatches on `ctrl` to the parallel propagator BVP
+(`use_parallel`), the dual Riccati formulation (`use_riccati`), or
+[`serial_eulerlagrange_integration`](@ref).
 
-### Returns
-
-An OdeState struct containing the final state of the ODE solver after integration is complete.
+Only the parallel branch populates `propagators` / `chunks` / `S_left`, which
+`compute_delta_prime_matrix!` consumes for the Δ' BVP; the other two return `nothing` for all
+three.
 """
 function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
@@ -166,6 +157,20 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
     elseif ctrl.use_riccati
         return (riccati_eulerlagrange_integration(ctrl, equil, ffit, intr), nothing, nothing, nothing)
     end
+    return serial_eulerlagrange_integration(ctrl, equil, ffit, intr)
+end
+
+"""
+    serial_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing)
+
+Serial shooting branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk,
+applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.ucrit` and undoing it
+via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. Call
+directly to force this branch regardless of `ctrl.use_parallel` / `ctrl.use_riccati`; `verbose`
+overrides `ctrl.verbose` for progress logging.
+"""
+function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal;
+    verbose::Bool=ctrl.verbose)
 
     # Initialization
     odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
@@ -182,7 +187,7 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
     chunks = chunk_el_integration_bounds(odet, ctrl, intr)
 
     # Print initial integration condition
-    if ctrl.verbose
+    if verbose
         @info "   ψ = $((@sprintf "%.3f" odet.psifac)),  q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))"
     end
 
@@ -190,7 +195,7 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
     for chunk in chunks
         # Integrate this region and display progress
         integrate_el_region!(odet, ctrl, equil, ffit, intr, chunk)
-        if ctrl.verbose
+        if verbose
             @info "   ψ = $((@sprintf "%.3f" odet.psifac)),  q = $((@sprintf "%.3f" odet.q)),  steps = $(odet.total_steps)"
         end
 
@@ -229,20 +234,20 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
             intr.psilim = odet.psi_store[end]
             intr.qlim = odet.q_store[end]
             odet.u .= odet.u_store[:, :, :, end]
-            if ctrl.verbose
+            if verbose
                 @info "Truncating integration at peak edge dW (LEGACY — Δ'/δW unreliable): ψ = $((@sprintf "%.3f" odet.psi_store[odet.step])),  q = $((@sprintf "%.3f" odet.q_store[odet.step]))"
             end
         else
             odet.psifac = saved_psifac
             odet.u .= saved_u
-            if ctrl.verbose
+            if verbose
                 @info "Edge-dW peak (diagnostic): ψ = $((@sprintf "%.3f" odet.psi_store[peak_step])),  q = $((@sprintf "%.3f" odet.q_store[peak_step])); integration domain unchanged"
             end
         end
     end
 
     # Evaluate stability criterion (critical determinant) of saved solutions
-    if ctrl.verbose
+    if verbose
         @info "Evaluating fixed-boundary stability criterion"
     end
     odet.nzero = evaluate_stability_criterion!(odet, equil.profiles)
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index d7569a437..0388b43f1 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -140,8 +140,8 @@ A mutable struct holding internal state variables for stability calculations.
   - `mlow::Int` - Lowest poloidal mode number
   - `mhigh::Int` - Highest poloidal mode number
   - `mpert::Int` - Number of poloidal modes (mhigh - mlow + 1)
-  - `nlow::Int` - Lowest toroidal mode number
-  - `nhigh::Int` - Highest toroidal mode number
+  - `nlow::Int` - Lowest toroidal mode number, resolved from `ctrl.nn_low`/`nn_high`
+  - `nhigh::Int` - Highest toroidal mode number, resolved from `ctrl.nn_low`/`nn_high`
   - `npert::Int` - Number of toroidal modes (nhigh - nlow + 1)
   - `numpert_total::Int` - Total number of modes (mpert × npert)
   - `keq_out::Bool` - Flag to output equilibrium quantities (not yet implemented)
@@ -158,7 +158,6 @@ A mutable struct holding internal state variables for stability calculations.
   - `psilim::Float64` - Flux limit for integration
   - `qlim::Float64` - Safety factor at psilim
   - `q1lim::Float64` - Safety factor derivative at psilim
-  - `locstab::CubicSeriesInterpolant` - Spline for local stability analysis
   - `wall_settings::Vacuum.WallShapeSettings` - Wall shape settings for vacuum calculations
 """
 @kwdef mutable struct ForceFreeStatesInternal
@@ -185,7 +184,6 @@ A mutable struct holding internal state variables for stability calculations.
     psilow::Float64 = 0.0   # lower integration bound; raised above the axis by sing_min! when qlow > qmin (RDCON gal)
     qlim::Float64 = 0.0
     q1lim::Float64 = 0.0
-    locstab::FastInterpolations.CubicSeriesInterpolant = cubic_interp(collect(0.0:0.25:1.0), Series(zeros(5, 5)); bc=ZeroCurvBC())
     debug_settings::DebugSettings = DebugSettings()
     wall_settings::Vacuum.WallShapeSettings = Vacuum.WallShapeSettings()
     """
@@ -217,7 +215,8 @@ end
 """
     ForceFreeStatesControl
 
-A mutable struct containing control parameters for stability analysis, set by the user in gpec.toml.
+An immutable struct containing 'ForceFreeStates' parameters set by the user in
+gpec.toml.
 
 ## Fields
 
@@ -260,7 +259,7 @@ A mutable struct containing control parameters for stability analysis, set by th
   - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline.  This populates `u_store` / `ud_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly.  Without it the parallel path stores only chunk-endpoint Riccati S matrices and zeros for `ud_store` (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable.  Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`.  Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run!` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases).  **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage).  Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case.  Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N).
   - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64.
 """
-@kwdef mutable struct ForceFreeStatesControl
+@kwdef struct ForceFreeStatesControl
     verbose::Bool = true
     local_stability_flag::Bool = false
     vac_flag::Bool = false
diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl
index 1be78f7f2..cb7cfeee8 100644
--- a/src/ForceFreeStates/Riccati.jl
+++ b/src/ForceFreeStates/Riccati.jl
@@ -1279,19 +1279,19 @@ end
 """
     riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) -> OdeState
 
-Main driver for integrating the dual Riccati ODE across the plasma.
-Functionally identical to `eulerlagrange_integration` except:
-
-1. Uses `riccati_integrate_chunk!`: drives `sing_der!` with `riccati_integrator_callback!`
-   which applies `renormalize_riccati_inplace!` (instead of Gaussian reduction) when
-   column norms exceed ucrit
-2. Uses `riccati_cross_ideal_singular_surf!` instead of `cross_ideal_singular_surf!`:
-   skips Gaussian reduction (avoids near-zero pivot issues when S is small near axis)
-   and renormalizes to (S_new, I) in one step
-3. Skips `transform_u!` — S is already the true solution, no Gaussian-reduction undo needed
-
-Enable via `use_riccati = true` in `[ForceFreeStates]` section of gpec.toml, or by
-setting `ctrl.use_riccati = true` programmatically.
+Integrate the dual Riccati ODE S = U₁·U₂⁻¹ across the plasma (Glasser 2018 Phys. Plasmas 25,
+032507). Reduces stiffness relative to [`serial_eulerlagrange_integration`](@ref), which it
+otherwise mirrors, differing in three places:
+
+1. `riccati_integrate_chunk!` drives `sing_der!` with `riccati_integrator_callback!`, which
+   applies `renormalize_riccati_inplace!` rather than Gaussian reduction when column norms
+   exceed `ctrl.ucrit`
+2. `riccati_cross_ideal_singular_surf!` replaces `cross_ideal_singular_surf!`: it skips
+   Gaussian reduction (avoiding near-zero pivots where S is small near the axis) and
+   renormalizes to (S_new, I) in one step
+3. `transform_u!` is skipped — S is already the true solution, so there is no reduction to undo
+
+Enable via `use_riccati = true` in the `[ForceFreeStates]` section of gpec.toml.
 """
 function riccati_eulerlagrange_integration(
     ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium,
@@ -1581,12 +1581,14 @@ function apply_propagator_inverse!(odet::OdeState, prop::ChunkPropagator)
 end
 
 """
-    parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) -> OdeState
+    parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left)
 
-Parallel fundamental matrix (propagator) driver for the EL integration.
+Parallel fundamental matrix (propagator) driver for the EL integration. The trailing three
+return values feed the Δ' BVP in `compute_delta_prime_matrix!`; this is the only branch that
+produces them.
 
-Functionally equivalent to `eulerlagrange_integration`, integrating all bulk chunks
-concurrently using `Threads.@threads`, then re-integrating the outer plasma serially:
+Equivalent to [`serial_eulerlagrange_integration`](@ref), but integrates all bulk chunks
+concurrently using `Threads.@threads`, then re-integrates the outer plasma serially:
 
 1. **Chunk generation**: calls `chunk_el_integration_bounds`, then `balance_integration_chunks`
    to sub-divide chunks for load-balanced parallel execution.
@@ -1602,22 +1604,18 @@ concurrently using `Threads.@threads`, then re-integrating the outer plasma seri
    without renormalization); Riccati integration keeps matrices bounded and provides dense
    checkpoints for `findmax_dW_edge!`.
 
-Enable via `use_parallel = true` in `[ForceFreeStates]` of gpec.toml, or by setting
-`ctrl.use_parallel = true` programmatically. Requires `singfac_min != 0`.
+Enable via `use_parallel = true` in `[ForceFreeStates]` of gpec.toml. Requires `singfac_min != 0`.
 
-**Key differences from standard integration:**
+**Key differences from serial integration:**
 - No Gaussian reduction in the propagator BVP phase (crossings use the
   Riccati-style algorithm, parallel `odet.ifix` stays 0)
 - `transform_u!` is called on the parallel odet but is a no-op (ifix=0)
 - Outer plasma uses serial Riccati integration for numerical stability
-- A serial Euler-Lagrange **dense pass** is appended at the end and
-  replaces the parallel `odet` so that `u_store` / `ud_store` are dense and
-  in axis basis — the only convention the PerturbedEquilibrium downstream
-  code consumes correctly.  Δ' (`singular/delta_prime_matrix`) is computed
-  from the parallel BVP and is bit-identical with vs. without this pass.
-  Toggle off with `ctrl.populate_dense_xi = false` if only Δ' / vacuum /
-  energies are needed and the extra serial-EL cost is unwanted (HDF5
-  `integration/xi_*` will then be sparse / zero).
+- When `ctrl.populate_dense_xi` is set, a serial EL dense pass is appended and replaces the
+  parallel `odet`, so `u_store` / `ud_store` come back in the axis basis that
+  PerturbedEquilibrium requires. Δ' is computed from the parallel BVP either way and is
+  bit-identical between the two. See the `populate_dense_xi` entry in the
+  [`ForceFreeStatesControl`](@ref) docstring for the cost trade-off.
 
 **Bidirectional integration for large-N accuracy:**
 The crossing chunk (nearest to each rational surface singL[j]) is integrated *backward*
@@ -1931,27 +1929,12 @@ function _populate_dense_xi_via_serial_el!(
         ) for s in 1:msing],
     )
 
-    # Temporarily switch dispatch flags so `eulerlagrange_integration`
-    # follows the serial EL branch (axis-basis u_store) for this call.
-    saved_use_parallel = ctrl.use_parallel
-    saved_use_riccati  = ctrl.use_riccati
-    saved_verbose      = ctrl.verbose
-    ctrl.use_parallel = false
-    ctrl.use_riccati  = false
-    ctrl.verbose      = false  # suppress duplicate per-chunk logging
-
-    if saved_verbose
+    if ctrl.verbose
         @info "   S → ξ: serial EL dense pass for HDF5 integration/xi_*"
     end
 
-    local fresh_odet::OdeState
-    try
-        fresh_odet, _, _, _ = eulerlagrange_integration(ctrl, equil, ffit, intr)
-    finally
-        ctrl.use_parallel = saved_use_parallel
-        ctrl.use_riccati  = saved_use_riccati
-        ctrl.verbose      = saved_verbose
-    end
+    # Run the serial branch but suppress logging
+    fresh_odet, _, _, _ = serial_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=false)
 
     # Restore BVP-result fields on `intr`.
     intr.psilim = saved.psilim
diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl
index 211e6f30d..1255c493a 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -122,19 +122,19 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     # cutoff depends on which n is used, so it isn't well-defined. Single-n with nn_low <= 0
     # (e.g. uninitialized default) is also skipped because the formula divides by nn_low.
     # Both cases fall back to qhigh / psihigh truncation with a warning.
-    if ctrl.set_psilim_via_dmlim && ctrl.nn_low != ctrl.nn_high
-        @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(ctrl.nn_low), nn_high=$(ctrl.nn_high)); falling back to qhigh / psihigh truncation."
-    elseif ctrl.set_psilim_via_dmlim && ctrl.nn_low <= 0
-        @warn "set_psilim_via_dmlim = true requires nn_low > 0; got nn_low=$(ctrl.nn_low). Falling back to qhigh / psihigh truncation."
+    if ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh
+        @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation."
+    elseif ctrl.set_psilim_via_dmlim && intr.nlow <= 0
+        @warn "set_psilim_via_dmlim = true requires nn_low > 0; got nn_low=$(intr.nlow). Falling back to qhigh / psihigh truncation."
     elseif ctrl.set_psilim_via_dmlim
         @info "Setting psilim via dmlim: initial qlim = $(@sprintf("%.3f", intr.qlim)), dmlim = $(@sprintf("%.3f", ctrl.dmlim))"
         # Normalize dmlim ∈ [0,1)
-        ctrl.dmlim = mod(ctrl.dmlim, 1.0)
-        intr.qlim = (trunc(Int, ctrl.nn_low * intr.qlim) + ctrl.dmlim) / ctrl.nn_low
+        dmlim = mod(ctrl.dmlim, 1.0)
+        intr.qlim = (trunc(Int, intr.nlow * intr.qlim) + dmlim) / intr.nlow
 
         # Reduce qlim if above qmax
         while intr.qlim > equil.params.qmax
-            intr.qlim -= 1.0 / ctrl.nn_low
+            intr.qlim -= 1.0 / intr.nlow
         end
     end
 
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index f6b1c7b41..8c45231dd 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -64,7 +64,7 @@ include("Rerun.jl")
 # Import ForceFreeStates types and functions needed for main
 using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, VacuumData, OdeState, FourFitVars
 using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry
-using .ForceFreeStates: compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries
+using .ForceFreeStates: compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries
 using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix
 using .ForceFreeStates: find_kinetic_singular_surfaces!
 using .ForceFreeStates: eulerlagrange_integration, free_run!
@@ -185,32 +185,26 @@ function main_from_inputs(
     _drop_deprecated_keys!(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates")
     ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in ffs_table)...)
 
-    # Determine toroidal mode numbers (n >= 1 required; 0 means "not specified").
-    # Validated before equilibrium formation: the two-pass grid refinement needs the
-    # n range to pin rational-surface knots.
-    if ctrl.nn_low == 0 && ctrl.nn_high == 0
+    # Determine toroidal mode numbers (n >= 1 required; 0 means "not specified")
+    intr.nlow, intr.nhigh = ctrl.nn_low, ctrl.nn_high
+    if intr.nlow == 0 && intr.nhigh == 0
         error("Either nn_low or nn_high must be set in [ForceFreeStates] (both are 0)")
-    elseif ctrl.nn_low == 0
-        ctrl.nn_low = ctrl.nn_high
-    elseif ctrl.nn_high == 0
-        ctrl.nn_high = ctrl.nn_low
+    elseif intr.nlow == 0
+        intr.nlow = intr.nhigh
+    elseif intr.nhigh == 0
+        intr.nhigh = intr.nlow
     end
-    if ctrl.nn_low > ctrl.nn_high
-        error("nn_low=$(ctrl.nn_low) cannot be greater than nn_high=$(ctrl.nn_high)")
+    if intr.nlow > intr.nhigh
+        error("nn_low=$(intr.nlow) cannot be greater than nn_high=$(intr.nhigh)")
     end
-    # checks for negative n
-    # note that negative n in fortran had code adding the identitiy matrix to grad Green for n=0
-    # and some n, nu sign switching in vacuum but was not actually supported by DCON sing_find, etc.
-    if ctrl.nn_high < 1
-        error("All requested toroidal modes (n=$(ctrl.nn_low):$(ctrl.nn_high)) are below 1; " *
+    if intr.nhigh < 1
+        error("All requested toroidal modes (n=$(intr.nlow):$(intr.nhigh)) are below 1; " *
               "n < 1 modes are not supported")
     end
-    if ctrl.nn_low < 1
-        @warn "Clamping nn_low from $(ctrl.nn_low) to 1; n < 1 modes are not supported"
-        ctrl.nn_low = 1
+    if intr.nlow < 1
+        @warn "Clamping nn_low from $(intr.nlow) to 1; n < 1 modes are not supported"
+        intr.nlow = 1
     end
-    intr.nlow = ctrl.nn_low
-    intr.nhigh = ctrl.nn_high
     intr.npert = intr.nhigh - intr.nlow + 1
     nstring = intr.npert == 1 ? "$(intr.nlow)" : "$(intr.nlow):$(intr.nhigh)"
 
@@ -246,7 +240,7 @@ function main_from_inputs(
     # kinetic profiles), pin knots on rational surfaces, and re-form on the refined grid
     # from the in-memory input — no file re-read.
     if Equilibrium.wants_two_pass(eq_config)
-        mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high)
+        mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh)
         psi_nodes = Equilibrium.refined_psi_grid(equil;
             tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory)
         rerun_input = if additional_input !== nothing
@@ -321,10 +315,6 @@ function main_from_inputs(
     @info "\n  Force-Free States\n$_SECTION"
     ffs_start = time()
 
-    # Set up variables
-    # TODO: parallel threads logic
-    ctrl.delta_mhigh *= 2 # for consistency with Fortran DCON TODO: why is this present in the Fortran?
-
     # Determine psilim and qlim (where we will integrate to)
     sing_lim!(intr, ctrl, equil)
 
@@ -337,18 +327,15 @@ function main_from_inputs(
         # equil = set_up_equilibrium(equil.config)
     end
 
-    # Compute local stability (if desired). This holds `D_I` from the
-    # ballooning coefficient system and the local ballooning result.
-    profiles_xs = equil.profiles.xs
-    locstab_fs = zeros(Float64, length(profiles_xs), 5)
+    # Compute local stability (if desired). `locstab` holds `D_I` from the ballooning
+    # coefficient system and the local ballooning result; `nothing` when not computed.
+    locstab = nothing
     ballooning_boundary = (psi=Float64[], alpha=Float64[], alpha_critical=Float64[])
     if ctrl.local_stability_flag
-        compute_ballooning_stability!(ctrl, locstab_fs, equil)
+        locstab = compute_local_stability(ctrl, equil)
         # First ballooning stability boundary (α vs ψ_N) for BALOO-style diagnostics.
         ballooning_boundary = ballooning_alpha_boundary(ctrl, equil)
     end
-    # Fit data to splines
-    intr.locstab = cubic_interp(profiles_xs, Series(locstab_fs); extrap=ExtendExtrap())
 
     # Find all singular surfaces in the equilibrium
     sing_find!(intr, equil)
@@ -386,19 +373,21 @@ function main_from_inputs(
     end
 
     # Determine poloidal mode numbers
+    # TODO: delta_mhigh is doubled for consistency with Fortran - why is this present in the Fortran?
+    delta_mhigh = 2 * ctrl.delta_mhigh
     if ctrl.delta_mlow < 0 || ctrl.delta_mhigh < 0
         error("Negative delta_mlow or delta_mhigh not allowed")
     end
     if ctrl.sing_start == 0
         intr.mlow = trunc(Int, min(intr.nlow * equil.params.qmin, 0)) - 4 - ctrl.delta_mlow
-        intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+        intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + delta_mhigh
     else
         intr.mmin = Inf # HUGE in Fortran
         for ising in Int(ctrl.sing_start):intr.msing
             intr.mmin = min(intr.mmin, sing[ising].m)
         end
         intr.mlow = intr.mmin - ctrl.delta_mlow
-        intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+        intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + delta_mhigh
     end
     intr.mpert = intr.mhigh - intr.mlow + 1
     intr.numpert_total = intr.mpert * intr.npert
@@ -504,6 +493,7 @@ function main_from_inputs(
             inputs,
             forcing_modes_snapshot,
             gal_data;
+            locstab=locstab,
             ballooning_boundary=ballooning_boundary
         )
         @info "Results written to $(ctrl.HDF5_filename)"
@@ -709,6 +699,7 @@ function write_outputs_to_HDF5(
     inputs::Union{Nothing,Dict{String,Any}}=nothing,
     forcing_modes::Union{Nothing,Vector{ForcingTerms.ForcingMode}}=nothing,
     gal_data::Union{GalerkinResult,Nothing}=nothing;
+    locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing,
     ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[])
 )
 
@@ -783,17 +774,17 @@ function write_outputs_to_HDF5(
         # locstab/di = Mercier D_I (det(d0bar)); locstab/dr = resistive interchange D_R;
         # locstab/ballooning_Delta_prime = high-n ballooning Δ' (distinct from the Riccati
         # tearing Δ' under perturbed_equilibrium/singular_coupling/delta_prime).
-        if ctrl.local_stability_flag
-            locstab_xs = intr.locstab.cache.x
-            out_h5["locstab/di"] = intr.locstab.y[:, 1] ./ locstab_xs
-            out_h5["locstab/dr"] = intr.locstab.y[:, 2] ./ locstab_xs
+        if locstab !== nothing
+            locstab_xs = locstab.cache.x
+            out_h5["locstab/di"] = locstab.y[:, 1] ./ locstab_xs
+            out_h5["locstab/dr"] = locstab.y[:, 2] ./ locstab_xs
         else
             out_h5["locstab/di"] = Float64[]
             out_h5["locstab/dr"] = Float64[]
         end
-        out_h5["singular/di0"] = (ctrl.local_stability_flag && !isempty(intr.sing)) ?
-                                 [intr.locstab(sing.psifac)[1] / sing.psifac for sing in intr.sing] : Float64[]
-        out_h5["locstab/ballooning_Delta_prime"] = ctrl.local_stability_flag ? intr.locstab.y[:, 4] : Float64[]
+        out_h5["singular/di0"] = (locstab !== nothing && !isempty(intr.sing)) ?
+                                 [locstab(sing.psifac)[1] / sing.psifac for sing in intr.sing] : Float64[]
+        out_h5["locstab/ballooning_Delta_prime"] = locstab !== nothing ? locstab.y[:, 4] : Float64[]
 
         # First ballooning stability boundary: experimental α vs critical α (BALOO-style).
         out_h5["locstab/psi"] = ballooning_boundary.psi
diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl
index 6f3ebc548..358d01bae 100644
--- a/test/runtests_eulerlagrange.jl
+++ b/test/runtests_eulerlagrange.jl
@@ -214,8 +214,7 @@ end
         mpert = 2
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(mpert, 10, 10, 10)
         intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; mpert=mpert)
-        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl()
-        ctrl.ucrit = 10.0
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; ucrit=10.0)
 
         # Case 1: Basic norm computation
         odet.u = zeros(ComplexF64, 2, 2, 2)
@@ -291,9 +290,7 @@ end
 
     @testset "chunk_el_integration_bounds tests" begin
         # Helper to build a minimal control and internal structs
-        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl()
-        ctrl.numsteps_init = 10
-        ctrl.numunorms_init = 5
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; numsteps_init=10, numunorms_init=5, singfac_min=1e-4)
 
         # Case 1: No singular surfaces -> single chunk to edge
         intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; mpert=1, numpert_total=1)
@@ -303,7 +300,6 @@ end
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
         odet.psifac = 0.0
 
-        ctrl.singfac_min = 1e-4
         chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr)
         @test length(chunks) == 1
         @test chunks[1].needs_crossing == false
@@ -323,7 +319,6 @@ end
 
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
         odet.psifac = 0.0
-        ctrl.singfac_min = 1e-4
 
         chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr)
         @test length(chunks) == 2
@@ -343,7 +338,7 @@ end
         intr.mhigh = 1
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
         odet.psifac = 0.0
-        ctrl.singfac_min = 1e-6
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; numsteps_init=10, numunorms_init=5, singfac_min=1e-6)
         chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr)
         @test length(chunks) == 3
         @test all(c.needs_crossing == true for c in chunks[1:2])
@@ -356,7 +351,7 @@ end
         intr.psilim = 1.0
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
         odet.psifac = 0.0
-        ctrl.singfac_min = 0.0
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; numsteps_init=10, numunorms_init=5, singfac_min=0.0)
         chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr)
         @test length(chunks) == 1
         @test chunks[1].needs_crossing == false
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index 26c4a9ee9..ff66bc035 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -90,7 +90,8 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
     @testset "sing_der" begin
 
         equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))
-        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl()
+        # default single toroidal mode used in test data; matches intr.nlow set below
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; nn_low=1)
         intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal()
         intr.numpert_total = 32 # replacing mpert (we set equal to 32). This is the same as msol
         # set mode ranges so sing_der can form singfac_vec consistently
@@ -101,7 +102,6 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         intr.nlow = 1
         intr.nhigh = 1
         intr.npert = intr.nhigh - intr.nlow + 1
-        ctrl.nn_low = intr.nlow
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(; numpert_total=intr.numpert_total,
             numsteps_init=ctrl.numsteps_init, numunorms_init=ctrl.numunorms_init, msing=intr.msing)
 
@@ -171,14 +171,13 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
     # ---------------------------------
     @testset "sing_lim" begin
         equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))
-        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl()
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=equil.params.qmax)
         intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal()
-        ctrl.qhigh = equil.params.qmax
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         @test isapprox(intr.qlim, equil.params.qmax; atol=1e-12)
         @test isapprox(intr.psilim, equil.config.psihigh; atol=1e-12)
 
-        ctrl.qhigh = max(equil.params.qmin + 0.1, equil.params.qmax - 0.5)
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=max(equil.params.qmin + 0.1, equil.params.qmax - 0.5))
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         @test intr.qlim < equil.params.qmax + 1e-12
         @test intr.psilim <= equil.config.psihigh

From 8cf0b98dcbbfcbbf91859b35d0365a74010a59ef Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 12 Aug 2026 14:58:15 -0400
Subject: [PATCH 2/5] Removing dead FFSInternal fields

---
 src/ForceFreeStates/ForceFreeStatesStructs.jl | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index 0388b43f1..5c854c5d4 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -144,10 +144,6 @@ A mutable struct holding internal state variables for stability calculations.
   - `nhigh::Int` - Highest toroidal mode number, resolved from `ctrl.nn_low`/`nn_high`
   - `npert::Int` - Number of toroidal modes (nhigh - nlow + 1)
   - `numpert_total::Int` - Total number of modes (mpert × npert)
-  - `keq_out::Bool` - Flag to output equilibrium quantities (not yet implemented)
-  - `theta_out::Bool` - Flag to output theta coordinate data (not yet implemented)
-  - `xlmda_out::Bool` - Flag to output eigenvalue data (not yet implemented)
-  - `sol_base::Int` - Base index for solution vectors (not yet implemented)
   - `msing::Int` - Number of ideal singular surfaces
   - `kmsing::Int` - Number of kinetic singular surfaces (det(F̄) near-zeros)
   - `sing::Vector{SingType}` - Vector of ideal singular surface data
@@ -169,10 +165,6 @@ A mutable struct holding internal state variables for stability calculations.
     nhigh::Int = 0
     npert::Int = 0
     numpert_total::Int = 0
-    keq_out::Bool = false
-    theta_out::Bool = false
-    xlmda_out::Bool = false
-    sol_base::Int = 50
     msing::Int = 0
     kmsing::Int = 0
     sing::Vector{SingType} = SingType[]

From 90f2bb342306123097202b4e6a4886a815b1b8d1 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 12 Aug 2026 17:49:55 -0400
Subject: [PATCH 3/5] ForceFreeStates - BUGFIX - Pass sing_order explicitly
 instead of mutating a ctrl copy

Freezing ForceFreeStatesControl left one consumer behind: galerkin_solve built
ctrl_gal = deepcopy(ctrl) and assigned ctrl_gal.sing_order per surface, which
throws on an immutable struct. Any run with gal_flag = true and at least one
singular surface failed with

    setfield!: immutable struct of type ForceFreeStatesControl cannot be changed

No unit suite covers the Galerkin path, so this only surfaced under the
regression harness.

compute_sing_asymptotics and compute_sing_mmat! now take a sing_order keyword
defaulting to ctrl.sing_order, and galerkin_solve passes the per-surface value
directly, so the deepcopy is gone. As a side effect SingAsymptotics now records
the order actually used rather than whatever was last written onto the copy.

No numerical change: diiid_n1 48/48, gal_resistive_diiid 10/10 and
gal_resistive_pe 4/4 unchanged against a0cad260.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 11 ++-
 src/ForceFreeStates/Sing.jl                   | 68 ++++++++++---------
 2 files changed, 41 insertions(+), 38 deletions(-)

diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index 0b7dd658e..b551a7a31 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -79,19 +79,18 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
     # right = sig=+1, left = sig=-1, no √det normalization. The Mercier exponent α is a
     # property of the surface, so the left series reuses the right's α (alpha_override). The order is
     # raised to gal_sing_order + ceil(2·Re(α)) for high-Mercier-index surfaces (Fortran sing1_vmat).
-    ctrl_gal = deepcopy(ctrl)
     asymps = GalSingAsymp[]
     for s in sings
-        ctrl_gal.sing_order = ctrl.gal_sing_order
-        ar = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=1.0)
+        sing_order = ctrl.gal_sing_order
+        ar = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=1.0, sing_order=sing_order)
         if ctrl.gal_sing_order_ceiling
             order = ctrl.gal_sing_order + ceil(Int, 2 * real(ar.alpha[1]))
             if order > ctrl.gal_sing_order
-                ctrl_gal.sing_order = order
-                ar = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=1.0)
+                sing_order = order
+                ar = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=1.0, sing_order=sing_order)
             end
         end
-        al = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=-1.0, alpha_override=ar.alpha)
+        al = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=-1.0, alpha_override=ar.alpha, sing_order=sing_order)
         push!(asymps, GalSingAsymp(ar, al))
     end
 
diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl
index 1255c493a..6ddcacb40 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -203,6 +203,8 @@ See equations 41-48 in the Glasser Phys. Plasmas 2016 112506 for the mathematica
 ### Arguments
 
   - `singp::SingType`: Singular surface parameters
+  - `sing_order`: Expansion order, defaulting to `ctrl.sing_order`. The Galerkin path overrides it
+    per surface (`gal_sing_order`, raised for high-Mercier-index surfaces).
 
 ### Returns
 
@@ -215,12 +217,13 @@ function compute_sing_asymptotics(
     ffit::FourFitVars,
     intr::ForceFreeStatesInternal;
     sig::Float64=1.0,
-    alpha_override::Union{Nothing,Vector{ComplexF64}}=nothing
+    alpha_override::Union{Nothing,Vector{ComplexF64}}=nothing,
+    sing_order::Int=ctrl.sing_order
 )
 
     # Allocations
-    vmat = zeros(ComplexF64, intr.numpert_total, 2 * intr.numpert_total, 2, 2 * ctrl.sing_order + 1)
-    mmat = zeros(ComplexF64, intr.numpert_total, 2 * intr.numpert_total, 2, 2 * ctrl.sing_order + 3)
+    vmat = zeros(ComplexF64, intr.numpert_total, 2 * intr.numpert_total, 2, 2 * sing_order + 1)
+    mmat = zeros(ComplexF64, intr.numpert_total, 2 * intr.numpert_total, 2, 2 * sing_order + 3)
     power = zeros(ComplexF64, 2 * intr.numpert_total)
 
     # Compute the resonant (r) and nonresonant (n) indices of the shearing transformation matrix R
@@ -236,7 +239,7 @@ function compute_sing_asymptotics(
     # Compute mmat Taylor coefficients with direction parameter sig.
     # Fortran computes separate mmatl (sig=-1) and mmatr (sig=+1) — the sig flips
     # odd derivatives of all input quantities (q, F, G, K splines).
-    compute_sing_mmat!(mmat, singp, ctrl, equil.profiles, ffit, intr; sig=sig)
+    compute_sing_mmat!(mmat, singp, ctrl, equil.profiles, ffit, intr; sig=sig, sing_order=sing_order)
 
     # Extract direction-specific m0mat from zeroth-order mmat
     m0mat = if length(r1) == 1
@@ -283,7 +286,7 @@ function compute_sing_asymptotics(
     end
 
     # Higher order solutions — sig propagates through the recursion (Fortran STRIDE sing_solve).
-    for k in 1:(2*ctrl.sing_order)
+    for k in 1:(2*sing_order)
         solve_higher_order_vmat!(vmat, mmat, m0mat, alpha, r1, r2, n1, n2, power, intr, k; sig=sig)
     end
 
@@ -303,7 +306,7 @@ function compute_sing_asymptotics(
         msg *= @sprintf("  psifac= %+.12e, r1=%d, ipert0=%d\n", singp.psifac, r1[1], ipert0)
         msg *= @sprintf("  vmat(ip,ip,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0, 2, 1]), imag(vmat[ipert0, ipert0, 2, 1]))
         msg *= @sprintf("  vmat(ip,ip+N,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0+N, 2, 1]), imag(vmat[ipert0, ipert0+N, 2, 1]))
-        for k in 0:(2*ctrl.sing_order)
+        for k in 0:(2*sing_order)
             msg *= @sprintf("  k=%2d vmat(ip,ip,1)=%+.8e %+.8ei vmat(ip,ip,2)=%+.8e %+.8ei\n",
                 k, real(vmat[ipert0, ipert0, 1, k+1]), imag(vmat[ipert0, ipert0, 1, k+1]),
                 real(vmat[ipert0, ipert0, 2, k+1]), imag(vmat[ipert0, ipert0, 2, k+1]))
@@ -314,7 +317,7 @@ function compute_sing_asymptotics(
         msg
     end
 
-    return SingAsymptotics(ctrl.sing_order, alpha, r1, r2, n1, n2, power, vmat, mmat, m0mat)
+    return SingAsymptotics(sing_order, alpha, r1, r2, n1, n2, power, vmat, mmat, m0mat)
 end
 
 """
@@ -355,7 +358,8 @@ Add a spline for F directly instead of the lower triangular factorization to avo
     profiles::Equilibrium.ProfileSplines,
     ffit::FourFitVars,
     intr::ForceFreeStatesInternal;
-    sig::Float64=1.0
+    sig::Float64=1.0,
+    sing_order::Int=ctrl.sing_order
 )
 
     q_spline = profiles.q_spline
@@ -370,13 +374,13 @@ Add a spline for F directly instead of the lower triangular factorization to avo
     f_lower_interp = zeros!(pool, ComplexF64, Npert, Npert, 4)
     g_interp = zeros!(pool, ComplexF64, Npert, Npert, 4)
     k_interp = zeros!(pool, ComplexF64, Npert, Npert, 4)
-    f_lower = zeros!(pool, ComplexF64, Npert, Npert, ctrl.sing_order + 1)
+    f_lower = zeros!(pool, ComplexF64, Npert, Npert, sing_order + 1)
     f0_lower = zeros!(pool, ComplexF64, Npert, Npert)
-    ff_lower = zeros!(pool, ComplexF64, Npert, Npert, ctrl.sing_order + 1)
-    g_lower = zeros!(pool, ComplexF64, Npert, Npert, ctrl.sing_order + 1)
-    k = zeros!(pool, ComplexF64, Npert, Npert, ctrl.sing_order + 1)
+    ff_lower = zeros!(pool, ComplexF64, Npert, Npert, sing_order + 1)
+    g_lower = zeros!(pool, ComplexF64, Npert, Npert, sing_order + 1)
+    k = zeros!(pool, ComplexF64, Npert, Npert, sing_order + 1)
     v = zeros!(pool, ComplexF64, Npert, 2 * Npert, 2)
-    x = zeros!(pool, ComplexF64, Npert, 2 * Npert, 2, ctrl.sing_order + 1)
+    x = zeros!(pool, ComplexF64, Npert, 2 * Npert, 2, sing_order + 1)
     tmp_vec = acquire!(pool, ComplexF64, Npert)
 
     # Evaluate q spline and its derivatives, applying sig to odd derivatives.
@@ -436,34 +440,34 @@ Add a spline for F directly instead of the lower triangular factorization to avo
                 ipert = ipert_m + (ipert_n - 1) * intr.mpert
                 jpert = jpert_m + (ipert_n - 1) * intr.mpert
                 f_lower[ipert, jpert, 1] = singfac[ipert, 1] * f_lower_interp[ipert, jpert, 1]
-                if ctrl.sing_order ≥ 1
+                if sing_order ≥ 1
                     f_lower[ipert, jpert, 2] = singfac[ipert, 1] * f_lower_interp[ipert, jpert, 2] +
                                                singfac[ipert, 2] * f_lower_interp[ipert, jpert, 1]
                 end
-                if ctrl.sing_order ≥ 2
+                if sing_order ≥ 2
                     f_lower[ipert, jpert, 3] =
                         singfac[ipert, 1] * f_lower_interp[ipert, jpert, 3] +
                         2 * singfac[ipert, 2] * f_lower_interp[ipert, jpert, 2] +
                         singfac[ipert, 3] * f_lower_interp[ipert, jpert, 1]
                 end
-                if ctrl.sing_order ≥ 3
+                if sing_order ≥ 3
                     f_lower[ipert, jpert, 4] =
                         singfac[ipert, 1] * f_lower_interp[ipert, jpert, 4] +
                         3 * singfac[ipert, 2] * f_lower_interp[ipert, jpert, 3] +
                         3 * singfac[ipert, 3] * f_lower_interp[ipert, jpert, 2] +
                         singfac[ipert, 4] * f_lower_interp[ipert, jpert, 1]
                 end
-                if ctrl.sing_order ≥ 4
+                if sing_order ≥ 4
                     f_lower[ipert, jpert, 5] =
                         4 * singfac[ipert, 2] * f_lower_interp[ipert, jpert, 4] +
                         6 * singfac[ipert, 3] * f_lower_interp[ipert, jpert, 3] +
                         4 * singfac[ipert, 4] * f_lower_interp[ipert, jpert, 2]
                 end
-                if ctrl.sing_order ≥ 5
+                if sing_order ≥ 5
                     f_lower[ipert, jpert, 6] = 10 * singfac[ipert, 3] * f_lower_interp[ipert, jpert, 4] +
                                                10 * singfac[ipert, 4] * f_lower_interp[ipert, jpert, 3]
                 end
-                if ctrl.sing_order ≥ 6
+                if sing_order ≥ 6
                     f_lower[ipert, jpert, 7] = 20 * singfac[ipert, 4] * f_lower_interp[ipert, jpert, 4]
                 end
             end
@@ -478,7 +482,7 @@ Add a spline for F directly instead of the lower triangular factorization to avo
     # Julia will handle filling the upper half via the Hermitian property
     # internally, just like LAPACK does in Fortran
     fac0 = 1
-    for n in 0:ctrl.sing_order
+    for n in 0:sing_order
         fac1 = 1
         for j in 0:n
             for ipert_n in 1:intr.npert
@@ -507,34 +511,34 @@ Add a spline for F directly instead of the lower triangular factorization to avo
                 ipert = ipert_m + (ipert_n - 1) * intr.mpert
                 jpert = jpert_m + (ipert_n - 1) * intr.mpert
                 k[ipert, jpert, 1] = singfac[ipert, 1] * k_interp[ipert, jpert, 1]
-                if ctrl.sing_order ≥ 1
+                if sing_order ≥ 1
                     k[ipert, jpert, 2] = singfac[ipert, 1] * k_interp[ipert, jpert, 2] +
                                          singfac[ipert, 2] * k_interp[ipert, jpert, 1]
                 end
-                if ctrl.sing_order ≥ 2
+                if sing_order ≥ 2
                     k[ipert, jpert, 3] =
                         singfac[ipert, 1] * k_interp[ipert, jpert, 3] / 2 +
                         singfac[ipert, 2] * k_interp[ipert, jpert, 2] +
                         singfac[ipert, 3] * k_interp[ipert, jpert, 1] / 2
                 end
-                if ctrl.sing_order ≥ 3
+                if sing_order ≥ 3
                     k[ipert, jpert, 4] =
                         singfac[ipert, 1] * k_interp[ipert, jpert, 4] / 6 +
                         singfac[ipert, 2] * k_interp[ipert, jpert, 3] / 2 +
                         singfac[ipert, 3] * k_interp[ipert, jpert, 2] / 2 +
                         singfac[ipert, 4] * k_interp[ipert, jpert, 1] / 6
                 end
-                if ctrl.sing_order ≥ 4
+                if sing_order ≥ 4
                     k[ipert, jpert, 5] =
                         singfac[ipert, 2] * k_interp[ipert, jpert, 4] / 6 +
                         singfac[ipert, 3] * k_interp[ipert, jpert, 3] / 4 +
                         singfac[ipert, 4] * k_interp[ipert, jpert, 2] / 6
                 end
-                if ctrl.sing_order ≥ 5
+                if sing_order ≥ 5
                     k[ipert, jpert, 6] = singfac[ipert, 3] * k_interp[ipert, jpert, 4] / 12 +
                                          singfac[ipert, 4] * k_interp[ipert, jpert, 3] / 12
                 end
-                if ctrl.sing_order ≥ 6
+                if sing_order ≥ 6
                     k[ipert, jpert, 7] = singfac[ipert, 4] * k_interp[ipert, jpert, 4] / 36
                 end
             end
@@ -549,13 +553,13 @@ Add a spline for F directly instead of the lower triangular factorization to avo
                 ipert = ipert_m + (ipert_n - 1) * intr.mpert
                 jpert = jpert_m + (ipert_n - 1) * intr.mpert
                 g_lower[ipert, jpert, 1] = g_interp[ipert, jpert, 1]
-                if ctrl.sing_order ≥ 1
+                if sing_order ≥ 1
                     g_lower[ipert, jpert, 2] = g_interp[ipert, jpert, 2]
                 end
-                if ctrl.sing_order ≥ 2
+                if sing_order ≥ 2
                     g_lower[ipert, jpert, 3] = g_interp[ipert, jpert, 3] / 2
                 end
-                if ctrl.sing_order ≥ 3
+                if sing_order ≥ 3
                     g_lower[ipert, jpert, 4] = g_interp[ipert, jpert, 4] / 6
                 end
             end
@@ -578,7 +582,7 @@ Add a spline for F directly instead of the lower triangular factorization to avo
     @views x[:, :, 1, 1] = UpperTriangular(f0_lower') \ (LowerTriangular(f0_lower) \ x[:, :, 1, 1])
 
     # Higher-order: ∑Fⱼx¹ₙ₋ⱼ = -Kₙv¹ → x¹ₙ = F₀⁻¹(-∑Fⱼxₙ₋ⱼ - Kₙv¹)
-    for i in 1:ctrl.sing_order
+    for i in 1:sing_order
         for isol in 1:(2*intr.numpert_total)
             for j in 1:i
                 @views mul!(tmp_vec, Hermitian(ff_lower[:, :, j+1], :L), x[:, isol, 1, i-j+1])
@@ -591,7 +595,7 @@ Add a spline for F directly instead of the lower triangular factorization to avo
     end
 
     # Solve x²ₙ = (G - K^†F⁻¹K)v¹ + K^†F⁻¹v² = Gₙv¹ + ∑Kⱼ^† x¹ₙ₋ⱼ at each order
-    for i in 0:ctrl.sing_order
+    for i in 0:sing_order
         for isol in 1:(2*intr.numpert_total)
             for j in 0:i
                 @views mul!(tmp_vec, adjoint(k[:, :, j+1]), x[:, isol, 1, i-j+1])
@@ -615,7 +619,7 @@ Add a spline for F directly instead of the lower triangular factorization to avo
     # Start with the S⁻¹LS components
     # Glasser PoP 2023 eq. 39: at each other of L, we get contributions to z^k from RLR,
     # z^k+0.5 from RLA and ALR, and z^k+1 from ALA (where A is the nonresonant part)
-    for i in 0:ctrl.sing_order
+    for i in 0:sing_order
         mmat[r1, r2, :, j+1] .= x[r1, r2, :, i+1]
         mmat[r1, n2, :, j+2] .= x[r1, n2, :, i+1]
         mmat[n1, r2, :, j+2] .= x[n1, r2, :, i+1]

From 2cfa325fbcbb54cadb3ce68901ad94e3b7cb9b3a Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Wed, 12 Aug 2026 19:49:45 -0400
Subject: [PATCH 4/5] Riccati - MINOR - fix docstring accuracy

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

diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl
index 8bb4fc713..f67214c64 100644
--- a/src/ForceFreeStates/Riccati.jl
+++ b/src/ForceFreeStates/Riccati.jl
@@ -1925,7 +1925,7 @@ and does NOT populate `delta_prime`; we keep the parallel pass's values
 which `compute_delta_prime_matrix!` uses).
 
 Called from `parallel_eulerlagrange_integration` when
-`ctrl.populate_dense_xi = true` (default).  Approximate cost: one serial
+`ctrl.populate_dense_xi = true`.  Approximate cost: one serial
 EL integration on top of the parallel BVP phase.  Required to make
 `use_parallel = true` produce DCON eigenfunctions usable by the
 PerturbedEquilibrium downstream pipeline.

From 54d1ecd1bc1828c2fdb7a6daef69b0bd22e6e182 Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Fri, 14 Aug 2026 10:11:20 -0400
Subject: [PATCH 5/5] ForceFreeStates - MINOR - Fix ordering of sing_lim! with
 nlow/nhigh setting

---
 docs/src/stability.md                 |  2 +-
 src/ForceFreeStates/Sing.jl           | 20 +++++++++++---------
 test/runtests_parallel_integration.jl | 14 +++++++-------
 test/runtests_riccati.jl              |  2 +-
 test/runtests_sing.jl                 |  7 +++++--
 5 files changed, 25 insertions(+), 20 deletions(-)

diff --git a/docs/src/stability.md b/docs/src/stability.md
index 9e2d7f8d4..f74a41649 100644
--- a/docs/src/stability.md
+++ b/docs/src/stability.md
@@ -295,8 +295,8 @@ equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(
 intr  = FFS.ForceFreeStatesInternal(; dir_path=ex)
 intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
     (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-FFS.sing_lim!(intr, ctrl, equil)
 intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+FFS.sing_lim!(intr, ctrl, equil)
 FFS.sing_find!(intr, equil)
 intr.mlow  = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
 intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl
index 60bf920b7..194cf58a9 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -91,7 +91,7 @@ function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEqui
 end
 
 """
-    sing_lim!(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
+    sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
 
 Compute and set integration ψ, q, and q' limits by handling cases where user truncates
 before the last singular surface. Performs a similar function to `sing_lim`
@@ -117,15 +117,17 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     intr.q1lim = profiles.q_deriv(profiles.xs[end]; hint=Ref(profiles.npts_minus_1))
     intr.psilim = equil.config.psihigh
 
-    # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent).
-    # Multi-n runs (nn_low != nn_high) are not supported — the "outermost rational + dmlim/n"
-    # cutoff depends on which n is used, so it isn't well-defined. Single-n with nn_low <= 0
-    # (e.g. uninitialized default) is also skipped because the formula divides by nn_low.
-    # Both cases fall back to qhigh / psihigh truncation with a warning.
-    if ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh
+    # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent). The cutoff reads
+    # the *resolved* toroidal range on `intr`, so callers must assign intr.nlow / intr.nhigh
+    # before calling; an unresolved range is an error rather than a silent change of truncation
+    # strategy. Multi-n runs are not supported — the "outermost rational + dmlim/n" cutoff depends
+    # on which n is used — and fall back to qhigh / psihigh truncation with a warning.
+    if ctrl.set_psilim_via_dmlim && intr.nlow <= 0
+        error("sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " *
+              "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " *
+              "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead.")
+    elseif ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh
         @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation."
-    elseif ctrl.set_psilim_via_dmlim && intr.nlow <= 0
-        @warn "set_psilim_via_dmlim = true requires nn_low > 0; got nn_low=$(intr.nlow). Falling back to qhigh / psihigh truncation."
     elseif ctrl.set_psilim_via_dmlim
         @info "Setting psilim via dmlim: initial qlim = $(@sprintf("%.3f", intr.qlim)), dmlim = $(@sprintf("%.3f", ctrl.dmlim))"
         # Normalize dmlim ∈ [0,1)
diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl
index 307615f2d..f9d7ba512 100644
--- a/test/runtests_parallel_integration.jl
+++ b/test/runtests_parallel_integration.jl
@@ -115,10 +115,10 @@ using TOML
             (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
         eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
         equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
-        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
+        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
         intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
         intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
@@ -179,10 +179,10 @@ using TOML
             (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
         eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
         equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
-        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
+        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
         intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
         intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
@@ -242,10 +242,10 @@ using TOML
             equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
             intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
                 (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-            GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
             intr.nlow = ctrl.nn_low;
             intr.nhigh = ctrl.nn_high;
             intr.npert = 1
+            GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
             GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
             intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
             intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
@@ -301,10 +301,10 @@ using TOML
             end
             intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
                 (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-            GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
             intr.nlow = ctrl.nn_low;
             intr.nhigh = ctrl.nn_high;
             intr.npert = 1
+            GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
             GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
             intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
             intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
@@ -350,10 +350,10 @@ using TOML
             (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
         eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
         equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
-        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
+        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
         intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
         intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
@@ -416,10 +416,10 @@ using TOML
             equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
             intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
                 (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-            GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
             intr.nlow = ctrl.nn_low;
             intr.nhigh = ctrl.nn_high;
             intr.npert = 1
+            GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
             GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
             intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
             intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
@@ -509,10 +509,10 @@ using TOML
         end
         intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
             (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
+        GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
         intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
         intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl
index e156a5be0..84194d948 100644
--- a/test/runtests_riccati.jl
+++ b/test/runtests_riccati.jl
@@ -9,10 +9,10 @@ function make_solovev_intr(inputs, ctrl, equil, ex)
     intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-    FFS.sing_lim!(intr, ctrl, equil)
     intr.nlow = ctrl.nn_low;
     intr.nhigh = ctrl.nn_high;
     intr.npert = 1
+    FFS.sing_lim!(intr, ctrl, equil)
     FFS.sing_find!(intr, equil)
     intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
     intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index ff66bc035..0536bf175 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -171,18 +171,21 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
     # ---------------------------------
     @testset "sing_lim" begin
         equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))
-        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=equil.params.qmax)
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=equil.params.qmax, set_psilim_via_dmlim=false)
         intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal()
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         @test isapprox(intr.qlim, equil.params.qmax; atol=1e-12)
         @test isapprox(intr.psilim, equil.config.psihigh; atol=1e-12)
 
-        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=max(equil.params.qmin + 0.1, equil.params.qmax - 0.5))
+        ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=max(equil.params.qmin + 0.1, equil.params.qmax - 0.5), set_psilim_via_dmlim=false)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
         @test intr.qlim < equil.params.qmax + 1e-12
         @test intr.psilim <= equil.config.psihigh
         q_at_psilim = equil.profiles.q_spline(intr.psilim)
         @test isapprox(q_at_psilim, intr.qlim; atol=1e-6)
+        ctrl_dmlim = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=equil.params.qmax, set_psilim_via_dmlim=true)
+        intr_unresolved = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal()
+        @test_throws ErrorException GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr_unresolved, ctrl_dmlim, equil)
     end
 
 end