OpenFUSIONToolkit · matt-pharr · Aug 17, 2026 · Aug 18, 2026 · Aug 18, 2026
diff --git a/docs/development/architecture.md b/docs/development/architecture.md
@@ -57,11 +57,14 @@ Splines are provided by the external `FastInterpolations` package rather than by
    - Solves ideal MHD eigenvalue problem with force-free boundary conditions
    - Identifies singular surfaces where ξ·∇ψ = 0
    - Key files:
-     - `ForceFreeStatesStructs.jl` - Core data structures
+     - `CoreTypes.jl` - Module-wide types (`ForceFreeStatesControl`, `ForceFreeStatesInternal`)
      - `Result.jl` - `ForceFreeStatesResult`, the published solve product every downstream stage reads
-     - `Ode.jl` - ODE solver for Euler-Lagrange equations
-     - `Sing.jl` - Singular point handling and layer analysis
-     - `Fourfit.jl` - Fourier fitting routines
+     - `EulerLagrange.jl` - ODE integration of the Euler-Lagrange equations (`OdeState`, derivative kernel)
+     - `Surfaces/` - Singular-surface finding, Frobenius asymptotics, and GGJ coefficients
+     - `Riccati/` - Chunked fundamental-matrix (STRIDE) driver and Δ' boundary-value problem
+     - `Galerkin/` - RDCON outer-region singular Galerkin Δ' solver
+     - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `resonant_match_rpec`)
+     - `Fourfit.jl` - Fourier fitting routines (`FourFitVars`)
      - `FixedBoundaryStability.jl` - Fixed boundary analysis
      - `Free.jl` - Free boundary stability
    - Status: Stable, core DCON functionality implemented

diff --git a/docs/src/citations.md b/docs/src/citations.md
@@ -36,7 +36,7 @@ The primary reference for the `ForceFreeStates` module. Derives the Euler-Lagran
 > *Physics of Plasmas* **25**, 032507 (2018).
 > DOI: [10.1063/1.5007042](https://doi.org/10.1063/1.5007042)
 
-Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati.jl` and enabled via `integrator = "riccati"` in `[ForceFreeStates]`.
+Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati/` and enabled via `integrator = "riccati"` in `[ForceFreeStates]`.
 
 ---
 

diff --git a/docs/src/conventions.md b/docs/src/conventions.md
@@ -155,7 +155,7 @@ modes on the high side as `delta_mlow` does on the low side.
 
 ### Why Positive ``m`` Is Always Resonant
 
-Resonant surfaces are found by `sing_find!` (`src/ForceFreeStates/Sing.jl`), which locates flux
+Resonant surfaces are found by `sing_find!` (`src/ForceFreeStates/Surfaces/Finding.jl`), which locates flux
 surfaces where ``m = n\,q(\psi)``. Since ``n > 0`` by convention and ``q > 0`` for a standard
 tokamak, the resonant ``m`` is always positive:
 

diff --git a/docs/src/stability.md b/docs/src/stability.md
@@ -292,7 +292,7 @@ The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separa
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
-Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Result.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"]
+Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/Resist.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
 ```
 
 ## Example usage

diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
@@ -22,19 +22,35 @@ using Printf
 using DoubleFloats
 import StaticArrays: @MMatrix
 
-# Include all necessary files
-include("ForceFreeStatesStructs.jl")
-include("Resist.jl")
-include("EulerLagrange.jl")
-include("Sing.jl")
-include("ResistEval.jl")
+# Types with cross-subsystem consumers, loaded before the code that dispatches on them
+include("Surfaces/Types.jl")
+include("CoreTypes.jl")
+include("Riccati/Types.jl")
 include("Fourfit.jl")
+include("Matching/DeltaPrime.jl")
+
+include("EulerLagrange.jl")
+
+# Singular-surface machinery: finding/filtering, Frobenius asymptotics, GGJ coefficients
+include("Surfaces/Finding.jl")
+include("Surfaces/Asymptotics.jl")
+include("Surfaces/Resist.jl")
+include("Surfaces/ResistEval.jl")
+
+# Outer<->inner resistive matching
+include("Matching/ResonantMatch.jl")
+
 include("FixedKineticMatrices.jl")
 include("Kinetic.jl")
 include("FixedBoundaryStability.jl")
 include("Utils.jl")
 include("Free.jl")
-include("Riccati.jl")
+
+# Chunked fundamental-matrix (Riccati/STRIDE) driver
+include("Riccati/Propagators.jl")
+include("Riccati/Crossings.jl")
+include("Riccati/DeltaPrimeBVP.jl")
+include("Riccati/Driver.jl")
 
 # RDCON outer-region singular Galerkin Δ′ solver (gal_solve port)
 include("Galerkin/GalerkinStructs.jl")