diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl
index 8419953015..15e93888ea 100644
--- a/benchmarks/benchmark_against_fortran_run.jl
+++ b/benchmarks/benchmark_against_fortran_run.jl
@@ -514,7 +514,7 @@ function setup_equil(p::FortranRunParams, fortran_dir::String)
         "flux_surface_rtol" => 1e-7
     )
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, fortran_dir)
-    return Equilibrium.setup_equilibrium(eq_config)
+    return Equilibrium.PlasmaEquilibrium(eq_config)
 end
 
 # ─── Phi_x comparison (coil mode) ────────────────────────────────────────────
diff --git a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl
index e292c24c0b..8da5616830 100644
--- a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl
+++ b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl
@@ -207,7 +207,7 @@ t_equil = @elapsed begin
         "flux_surface_rtol" => 1e-7,
     )
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, run_dir)
-    equil = Equilibrium.setup_equilibrium(eq_config)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config)
 end
 @printf "    Done in %.1f s  (bt_sign=%d)\n" t_equil equil.params.bt_sign
 
diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl
index 9390c1eb16..94c7dc0eee 100644
--- a/benchmarks/benchmark_delta_prime_methods.jl
+++ b/benchmarks/benchmark_delta_prime_methods.jl
@@ -28,7 +28,7 @@ function setup_and_run_solovev()
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
     eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     FFS.sing_lim!(intr, ctrl, equil)
diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl
index 4fad6fd286..ff1957c16f 100644
--- a/benchmarks/benchmark_q_vs_iota_edge.jl
+++ b/benchmarks/benchmark_q_vs_iota_edge.jl
@@ -23,7 +23,7 @@ function reference_q()
     inputs, _, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
     equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "core_edge_packed", "mpsi" => 1024))
     eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR)
-    equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
+    equil = GPE.Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
     # The dict rides along so the mpsi scan can rebuild configs; EquilibriumConfig is immutable.
     return equil, equil_dict
 end
diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl
index e9d88ef82f..eaaa64f895 100644
--- a/benchmarks/benchmark_riccati_der.jl
+++ b/benchmarks/benchmark_riccati_der.jl
@@ -29,7 +29,7 @@ function setup_solovev()
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
     eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     FFS.sing_lim!(intr, ctrl, equil)
diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl
index 1b2912a70a..8a9a05dc8f 100644
--- a/benchmarks/benchmark_threads.jl
+++ b/benchmarks/benchmark_threads.jl
@@ -17,7 +17,7 @@ function run_ffs(ex; integrator)
     ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
     eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
diff --git a/benchmarks/equil_method_comparison.jl b/benchmarks/equil_method_comparison.jl
index 450423d6a1..27fe93f95e 100644
--- a/benchmarks/equil_method_comparison.jl
+++ b/benchmarks/equil_method_comparison.jl
@@ -50,9 +50,9 @@ for method in methods
     local t2 = 0.0
     local success = true
     try
-        setup_equilibrium(cfg)  # JIT warmup
-        t1 = @elapsed pe = setup_equilibrium(cfg)
-        t2 = @elapsed setup_equilibrium(cfg)
+        PlasmaEquilibrium(cfg)  # JIT warmup
+        t1 = @elapsed pe = PlasmaEquilibrium(cfg)
+        t2 = @elapsed PlasmaEquilibrium(cfg)
     catch e
         @warn "  FAILED: $e"
         success = false
diff --git a/benchmarks/equil_numerical_params.jl b/benchmarks/equil_numerical_params.jl
index f80f0f64ca..e6d8e86168 100644
--- a/benchmarks/equil_numerical_params.jl
+++ b/benchmarks/equil_numerical_params.jl
@@ -56,7 +56,7 @@ end
 # JIT warmup
 println("JIT warmup...")
 warmup_cfg = make_config(config_path, "efit", 64, 128)
-setup_equilibrium(warmup_cfg)
+PlasmaEquilibrium(warmup_cfg)
 
 rows = []
 
@@ -72,7 +72,7 @@ for method in ["efit", "efit_arclength", "efit_by_inversion"]
         local runtime = NaN
         local rt_err = NaN
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             raw = Equilibrium.read_efit(cfg)
             rt_err = roundtrip_error(pe, raw)
@@ -98,7 +98,7 @@ for method in ["efit", "efit_arclength", "efit_by_inversion"]
         local runtime = NaN
         local rt_err = NaN
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             raw = Equilibrium.read_efit(cfg)
             rt_err = roundtrip_error(pe, raw)
@@ -123,7 +123,7 @@ println("=" ^ 65)
 
 # ── Contour grid refinement scan (efit_by_inversion only) ─────────────────────
 # The refine factor in equilibrium_solver_by_inversion is a keyword arg.
-# We call it directly here, bypassing the setup_equilibrium wrapper.
+# We call it directly here, bypassing the PlasmaEquilibrium constructor.
 println("\n" * "=" ^ 65)
 println("Contour grid refinement scan (efit_by_inversion only, mpsi=128, mtheta=256)")
 println("=" ^ 65)
diff --git a/benchmarks/equil_psihigh_scan.jl b/benchmarks/equil_psihigh_scan.jl
index 56f6e9ec6b..9bb89f6f24 100644
--- a/benchmarks/equil_psihigh_scan.jl
+++ b/benchmarks/equil_psihigh_scan.jl
@@ -97,7 +97,7 @@ for method in methods
         local success = false
         local err_msg = ""
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             success = true
         catch e
diff --git a/benchmarks/equil_psilow_scan.jl b/benchmarks/equil_psilow_scan.jl
index ff5902b7b6..882f87ed32 100644
--- a/benchmarks/equil_psilow_scan.jl
+++ b/benchmarks/equil_psilow_scan.jl
@@ -50,7 +50,7 @@ end
 # JIT warmup
 println("JIT warmup...")
 warmup_cfg = make_config(config_path, "efit", 1e-2)
-setup_equilibrium(warmup_cfg)
+PlasmaEquilibrium(warmup_cfg)
 
 println("=" ^ 70)
 println("psilow near-axis accuracy scan")
@@ -67,7 +67,7 @@ for method in methods
         local success = false
         local err_msg = ""
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             success = true
         catch e
diff --git a/benchmarks/equil_spline_comparison.jl b/benchmarks/equil_spline_comparison.jl
index 4cf5337e26..22ffacbc0b 100644
--- a/benchmarks/equil_spline_comparison.jl
+++ b/benchmarks/equil_spline_comparison.jl
@@ -82,7 +82,7 @@ for m in all_methods
     println("--- Running: $m ---")
     cfg = make_config(config_path, m, psihigh_arg)
     try
-        pes[m] = setup_equilibrium(cfg)
+        pes[m] = PlasmaEquilibrium(cfg)
         println("  Done.\n")
     catch e
         println("  FAILED: $e\n")
diff --git a/benchmarks/plot_grid_knot_placement.jl b/benchmarks/plot_grid_knot_placement.jl
index b8a6632340..50dc3ef968 100644
--- a/benchmarks/plot_grid_knot_placement.jl
+++ b/benchmarks/plot_grid_knot_placement.jl
@@ -19,7 +19,7 @@ function build_grids()
     tau = eq_config.psi_accuracy
 
     # Pass 1 (coarse layout) and the two-pass refined grid, as the driver builds them
-    equil1 = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
+    equil1 = GPE.Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
     pass1 = copy(equil1.profiles.xs)
     rationals = GPE.ForceFreeStates.rational_psi_nodes(equil1; nlow=1, nhigh=1)
     auto = GPE.Equilibrium.refined_psi_grid(equil1; tau, mandatory=rationals)
diff --git a/docs/development/architecture.md b/docs/development/architecture.md
index 059539a8b5..5b976b5e7a 100644
--- a/docs/development/architecture.md
+++ b/docs/development/architecture.md
@@ -26,7 +26,7 @@ Splines are provided by the external `FastInterpolations` package rather than by
 ### Core Physics Modules
 
 2. **Equilibrium** (`src/Equilibrium/`) - MHD equilibrium solvers
-   - Main entry point: `setup_equilibrium(path)` or `setup_equilibrium(config)`
+   - Main entry point: `PlasmaEquilibrium(file; eq_type)`, `PlasmaEquilibrium(analytic_config)`, `PlasmaEquilibrium(dd)`, or `PlasmaEquilibrium(config, source)`
    - Supports multiple equilibrium types:
      - `efit` - EFIT g-file format
      - `chease`, `chease2` - CHEASE equilibrium code formats
@@ -149,14 +149,14 @@ Example configuration files are provided in:
 - `examples/Solovev_ideal_example/gpec.toml`
 - `examples/DIIID-like_ideal_example/gpec.toml`
 
-**Note**: Legacy configuration files (`equil.toml`, `vac.in`) are deprecated.
+**Note**: Legacy configuration files (`equil.toml`, `vac.in`) are not read; `gpec.toml` is the only input deck.
 
 ## Data Flow
 
 The complete GPEC analysis pipeline:
 
 1. **Equilibrium Setup**:
-   - `setup_equilibrium(config)` reads configuration from `gpec.toml`
+   - `PlasmaEquilibrium(config, source)` forms the equilibrium from the `[Equilibrium]` section of `gpec.toml`
    - Parses equilibrium data (EFIT, CHEASE, or analytical)
    - Runs Grad-Shafranov solver (direct or inverse)
    - Computes global parameters: q-profile, pressure, current density, β
diff --git a/docs/src/api.md b/docs/src/api.md
index bf734f7621..864cd5d40e 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -28,9 +28,16 @@ ffs  = solve(prob, Riccati())
 The four objects map one-to-one onto the pipeline stages:
 
   - [`PlasmaEquilibrium`](@ref) reads and processes the equilibrium — the `[Equilibrium]`
-    section. Analytic equilibria (`sol`, `lar`, `tj_analytic`) take their parameters from a
-    separate config object and are built with
-    `setup_equilibrium(eq_config, analytic_config)` instead.
+    section. Analytic equilibria take their parameters from their config object, with the
+    `eq_type` implied, and IMAS equilibria from a data dictionary:
+
+    ```julia
+    eq = PlasmaEquilibrium(Equilibrium.SolovevConfig(; elongation=1.5); mpsi=0)
+    eq = PlasmaEquilibrium(dd; imas_cocos=11)
+    ```
+
+    Every equilibrium remembers its source, so with `mpsi=0` `solve` re-forms it on the same
+    refined two-pass grid the deck uses; `Equilibrium.reform(eq, psi_nodes)` re-forms it on any grid.
   - `solve` runs the force-free-states stage — the `[ForceFreeStates]` section — and returns
     a `ForceFreeStatesResult`, exactly the object the TOML driver publishes.
   - [`RMPField`](@ref) describes the external field — the `[ForcingTerms]` section — without
diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md
index 03c2fe8983..99f9584258 100644
--- a/docs/src/equilibrium.md
+++ b/docs/src/equilibrium.md
@@ -98,9 +98,9 @@ built by a **two-pass measured-curvature refinement** driven by the single accur
 Every region's knot count scales as τ^(-1/3), so tightening `psi_accuracy` refines
 the core, pedestal, and edge proportionally. The fixed `grid_type = "core_edge_packed"`
 (sin²-spaced), `"pow1"`, `"uniform"`, and explicit `mpsi > 0` (single-pass, fixed
-layout) are still supported. Library users calling `setup_equilibrium` directly with
-`mpsi = 0` receive the coarse pass-1 grid; use `refined_psi_grid` and the
-`override_psi_nodes` keyword to apply the refinement manually.
+layout) are still supported. A `PlasmaEquilibrium` built directly with `mpsi = 0` keeps
+the coarse pass-1 grid until `solve` re-forms it; to refine it yourself, pass the nodes from
+`refined_psi_grid` to `Equilibrium.reform`.
 
 The packing on the DIII-D-like example (n=1) compared to fixed `core_edge_packed` grids — note the
 coarse spacing across the smooth mid-radius, the spacing dips at each rational surface,
@@ -126,11 +126,8 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 
 ## Important types
 
-- `EquilibriumConfig` — top-level configuration container parsed from a
-	TOML file (outer constructor `EquilibriumConfig(path::String)` is
-	provided).
-- `EquilibriumControl` — low-level control parameters (grid, jacobian
-	type, tolerances, etc.).
+- `EquilibriumConfig` — configuration container, built from the `[Equilibrium]`
+	section of `gpec.toml` or from keywords.
 - `PlasmaEquilibrium` — the runtime structure containing spline fields,
 	geometry, profiles, and computed diagnostics (q-profile, separatrix,
 	etc.).
@@ -140,9 +137,11 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 
 ## Key functions
 
-- `setup_equilibrium(path::String = "equil.toml")`
-	— main entry point that reads configuration, builds the equilibrium,
-	runs the solver, and returns a `PlasmaEquilibrium` instance.
+- `PlasmaEquilibrium(file; eq_type="efit", kwargs...)`, `PlasmaEquilibrium(analytic_config; kwargs...)`
+	and `PlasmaEquilibrium(dd; kwargs...)` — read or solve the equilibrium and return a
+	`PlasmaEquilibrium`; `kwargs` are `EquilibriumConfig` fields.
+- `PlasmaEquilibrium(config::EquilibriumConfig, source=nothing)` — the form every other
+	constructor ends in, with `source` an analytic config, an IMAS `dd`, or `nothing` to read `eq_filename`.
 - `equilibrium_separatrix_find!(pe::PlasmaEquilibrium)` — locate
 	separatrix and related boundary geometry in-place.
 - `equilibrium_global_parameters!(pe::PlasmaEquilibrium)` — populate
@@ -154,15 +153,14 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 
 ## Example usage
 
-Basic example: read a TOML config and build an equilibrium
+Basic example: read a g-file and build an equilibrium
 
 ```julia
 using GeneralizedPerturbedEquilibrium
 
-# Build from a TOML file (searches relative paths if needed)
-pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium("docs/examples/ForceFreeStates.toml")
+pe = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada")
 
-println("Magnetic axis: ", pe.params.r0, ", ", pe.params.z0)
+println("Magnetic axis: ", pe.ro, ", ", pe.zo)
 println("q(0) = ", pe.params.q0)
 
 # Find separatrix (in-place) and inspect results
@@ -175,9 +173,8 @@ Analytic / testing example: construct a large-aspect-ratio model
 ```julia
 using GeneralizedPerturbedEquilibrium
 
-# Create a LAR config from a small TOML fragment or file
 larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(r0=10.0, a=1.0, beta0=1e-3)
-pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(control=Dict("eq_filename"=>"unused","eq_type"=>"lar")), larcfg)
+pe = PlasmaEquilibrium(larcfg)  # eq_type = "lar" is implied by the config type
 
 println("Built LAR equilibrium with a = ", larcfg.a)
 ```
diff --git a/docs/src/set_up.md b/docs/src/set_up.md
index fcebbc0a39..371dbac8e8 100644
--- a/docs/src/set_up.md
+++ b/docs/src/set_up.md
@@ -341,7 +341,7 @@ using GeneralizedPerturbedEquilibrium.Vacuum
 using GeneralizedPerturbedEquilibrium.ForceFreeStates
 
 # Set up equilibrium only
-equil = Equilibrium.setup_equilibrium("path/to/gpec.toml")
+equil = Equilibrium.PlasmaEquilibrium("path/to/g-file"; eq_type="efit")
 
 # Access equilibrium data
 println("q at axis: ", equil.params.q0)
diff --git a/docs/src/stability.md b/docs/src/stability.md
index b056753abf..4c2c9c2234 100644
--- a/docs/src/stability.md
+++ b/docs/src/stability.md
@@ -414,7 +414,7 @@ inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
 
 ctrl  = FFS.ForceFreeStatesControl(;
             (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(
+equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(
             GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex))
 
 intr  = FFS.ForceFreeStatesInternal(; dir_path=ex)
diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl
index 715e3ec76c..244465ab74 100644
--- a/examples/DIIID-like_ideal_example/analyze_example.jl
+++ b/examples/DIIID-like_ideal_example/analyze_example.jl
@@ -58,7 +58,7 @@ for geqdsk_file in geqdsk_cases
     eq_dict = deepcopy(base_toml["Equilibrium"])
     eq_dict["eq_filename"] = geqdsk_file
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, @__DIR__)
-    equil = Equilibrium.setup_equilibrium(eq_config)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config)
     psi_norm = Vector(equil.profiles.xs)
 
     # s and alpha profiles vs normalized flux
diff --git a/examples/LAR_beta_scan/run_scan.jl b/examples/LAR_beta_scan/run_scan.jl
index 79f251059b..4e15fa2139 100644
--- a/examples/LAR_beta_scan/run_scan.jl
+++ b/examples/LAR_beta_scan/run_scan.jl
@@ -16,7 +16,7 @@ using Pkg
 Pkg.activate(joinpath(@__DIR__, "../.."))
 
 using GeneralizedPerturbedEquilibrium
-using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig, setup_equilibrium
+using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig
 using HDF5
 using TOML
 using Printf
diff --git a/examples/LAR_epsilon_scan/run_scan.jl b/examples/LAR_epsilon_scan/run_scan.jl
index 2e6a156ea4..8d364ca199 100644
--- a/examples/LAR_epsilon_scan/run_scan.jl
+++ b/examples/LAR_epsilon_scan/run_scan.jl
@@ -16,7 +16,7 @@ using Pkg
 Pkg.activate(joinpath(@__DIR__, "../.."))
 
 using GeneralizedPerturbedEquilibrium
-using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig, setup_equilibrium
+using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig
 using HDF5
 using TOML
 using Printf
diff --git a/examples/efit_fixedbdy_separatrix_example/README.md b/examples/efit_fixedbdy_separatrix_example/README.md
index 1282056325..de6d8e6a7d 100644
--- a/examples/efit_fixedbdy_separatrix_example/README.md
+++ b/examples/efit_fixedbdy_separatrix_example/README.md
@@ -13,7 +13,7 @@ loads it regardless of box size.
 
 This g-file is the fixture for the `efit_fixedbdy_separatrix` regression case
 (`regression-harness/cases/efit_fixedbdy_separatrix.toml`, a `kind = "computed"` case that calls
-`setup_equilibrium` directly). It is **not** a runnable example deck — there is no `gpec.toml`.
+`PlasmaEquilibrium` directly). It is **not** a runnable example deck — there is no `gpec.toml`.
 
 Free-boundary equilibria (real coil currents) and standard diverted EFITs are unaffected: their
 outboard midplane flux keeps dropping past the LCFS, so the bracket was never ambiguous.
diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl
index 2e7ec54157..2c421654a4 100644
--- a/regression-harness/src/runner.jl
+++ b/regression-harness/src/runner.jl
@@ -151,7 +151,7 @@ end
 # computational box hugs the LCFS (eps=0.05 TokaMaker aspect-scan g-file): outside the prescribed
 # LCFS the coil-vacuum flux turns back above the boundary value before the grid edge, so the old
 # bracketed-Brent separatrix finder could not bracket psi=sibry and equilibrium setup threw. Calls
-# setup_equilibrium directly and writes the leading equilibrium scalars: errors on the buggy code,
+# PlasmaEquilibrium directly and writes the leading equilibrium scalars: errors on the buggy code,
 # passes with the Newton finder. Runtime and environment are written to ARGS[2] by the run-info epilogue.
 const COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE = """
 using Pkg
@@ -161,9 +161,8 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
 using HDF5
 root = dirname(Base.active_project())
 gfile = joinpath(root, "examples", "efit_fixedbdy_separatrix_example", "eq_eps0.0500000_k1.000_d0.000.geqdsk")
-cfg = Equilibrium.EquilibriumConfig(; eq_type="efit", eq_filename=gfile)
 t_start = time()
-pe = Equilibrium.setup_equilibrium(cfg)
+pe = Equilibrium.PlasmaEquilibrium(gfile; eq_type="efit")
 elapsed = time() - t_start
 h5open(ARGS[1], "w") do fid
     fid["Equilibrium/psi_total"]  = pe.psio
diff --git a/src/Equilibrium/DirectEquilibriumByInversion.jl b/src/Equilibrium/DirectEquilibriumByInversion.jl
index 41fe717b53..589c0a8928 100644
--- a/src/Equilibrium/DirectEquilibriumByInversion.jl
+++ b/src/Equilibrium/DirectEquilibriumByInversion.jl
@@ -703,7 +703,7 @@ function equilibrium_solver_by_inversion(
     rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_table; itp_opts2d...)
 
     # Intermediate inverse input; the captured DirectIngest rides on the original eq_input,
-    # which setup_equilibrium forwards onto the equilibrium, so this one carries ingest=nothing.
+    # which PlasmaEquilibrium forwards onto the equilibrium, so this one carries ingest=nothing.
     inv_input = InverseRunInput(raw_profile.config, raw_profile.sq_in,
         rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, nothing, raw_profile.psihigh_resolved)
 
diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl
index 51c92eaf93..79076d0fd6 100644
--- a/src/Equilibrium/Equilibrium.jl
+++ b/src/Equilibrium/Equilibrium.jl
@@ -10,6 +10,7 @@ using FastInterpolations
 using AdaptiveArrayPools
 import StaticArrays: @MMatrix, SVector
 import ..Utilities
+import IMASdd
 
 # --- Internal Module Structure ---
 include("EquilibriumTypes.jl")
@@ -28,7 +29,7 @@ include("GeometryProfiles.jl")
 include("KineticProfiles.jl")
 
 # --- Expose types and functions to the user ---
-export setup_equilibrium, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParameters,
+export reform, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParameters, AnalyticEquilibriumConfig,
     ProfileSplines, GeometryProfileSplines, compute_geometry_profiles,
     KineticProfileSplines, load_kinetic_profiles, shift_exb_rotation,
     KineticProfileData, read_kinetic_file, write_kinetic_h5
@@ -63,8 +64,8 @@ struct AnalyticEqSpec
     run_fn::Function
 end
 
-# Registry of analytic equilibrium kinds. Both `setup_equilibrium` (fresh runs) and the
-# rerun input builder dispatch off this table, so adding a new analytic kind is one new row.
+# Registry of analytic equilibrium kinds. Both `run_input` (fresh runs) and the rerun input
+# builder dispatch off this table, so adding a new analytic kind is one new row.
 const ANALYTIC_EQ = Dict(
     "sol" => AnalyticEqSpec("SOL_INPUT", SolovevConfig, sol_run),
     "lar" => AnalyticEqSpec("LAR_INPUT", LargeAspectRatioConfig, lar_run),
@@ -73,75 +74,79 @@ const ANALYTIC_EQ = Dict(
 )
 
 """
-    setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothing; override_psi_nodes=nothing)
+    run_input(eq_config::EquilibriumConfig, source) -> DirectRunInput or InverseRunInput
 
-Read an equilibrium file, run the appropriate solver, and return the processed
-`PlasmaEquilibrium` with global parameters, q-profile, and GSE diagnostics.
+Build the solver input for `eq_config` from `source`: `nothing` reads `eq_config.eq_filename`, an
+[`AnalyticEquilibriumConfig`](@ref) runs its analytic solver, an `IMASdd.dd` is read by
+[`read_imas`](@ref), and a run input rebuilt from a stored ingest is rebound onto `eq_config`.
+"""
+function run_input(eq_config::EquilibriumConfig, ::Nothing)
+    eq_type = eq_config.eq_type
+    eq_type in EFIT_KINDS && return resolve_psihigh!(read_efit(eq_config))
+    eq_type in ("chease2", "chease_ascii") && return read_chease_ascii(eq_config)
+    eq_type in ("chease", "chease_binary") && return read_chease_binary(eq_config)
+    eq_type == "ifile" && return resolve_psihigh!(read_ifile(eq_config))
+    if haskey(ANALYTIC_EQ, eq_type)
+        spec = ANALYTIC_EQ[eq_type]
+        error("analytic eq_type=\"$eq_type\" needs its $(spec.config_type) as the source (built from the embedded [$(spec.section)] section)")
+    end
+    eq_type == "imas" && error("eq_type=\"imas\" needs an IMASdd.dd as the source")
+    return error("Equilibrium type $(eq_type) is not implemented")
+end
+
+function run_input(eq_config::EquilibriumConfig, cfg::AnalyticEquilibriumConfig)
+    spec = get(ANALYTIC_EQ, eq_config.eq_type, nothing)
+    spec !== nothing && cfg isa spec.config_type ||
+        error("eq_type=\"$(eq_config.eq_type)\" is not an analytic kind solved from a $(typeof(cfg))")
+    return spec.run_fn(eq_config, cfg)
+end
+
+function run_input(eq_config::EquilibriumConfig, dd::IMASdd.dd)
+    eq_config.eq_type == "imas" || error("an IMASdd.dd source needs eq_type=\"imas\", not \"$(eq_config.eq_type)\"")
+    return read_imas(eq_config, dd)
+end
+
+# A run input rebuilt from an ingest is rebound onto the current config, resetting psihigh so a
+# pass-1 clamped value does not shadow a psihigh overridden in the rerun TOML.
+function run_input(eq_config::EquilibriumConfig, input::DirectRunInput)
+    input.config = eq_config
+    input.psihigh_resolved = eq_config.psihigh
+    return eq_config.eq_type in EFIT_KINDS ? resolve_psihigh!(input) : input
+end
+
+function run_input(eq_config::EquilibriumConfig, input::InverseRunInput)
+    input.config = eq_config
+    input.psihigh_resolved = eq_config.psihigh
+    return resolve_psihigh!(input)
+end
+
+"""
+    PlasmaEquilibrium(eq_config::EquilibriumConfig, source=nothing; override_psi_nodes=nothing) -> PlasmaEquilibrium
+
+Form the equilibrium `eq_config` describes from `source` (see [`run_input`](@ref)), run the
+solver, and return it with global parameters, q-profile, and GSE diagnostics. Every other
+constructor ends here.
 
 `override_psi_nodes` bypasses the config-driven radial grid (`grid_type`/`mpsi`) and forms
 the equilibrium on the given strictly increasing ψ_N node vector, whose endpoints must match
-`psilow`/`psihigh`. Used by the two-pass auto-grid refinement in the main driver; standalone
+`psilow`/`psihigh`. Used by the two-pass auto-grid refinement through [`reform`](@ref); standalone
 callers with `mpsi=0` get the coarse first-pass grid unless they refine and re-form.
 """
-function setup_equilibrium(path::String="equil.toml")
-    return setup_equilibrium(EquilibriumConfig(path))
-end
-function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothing; override_psi_nodes::Union{Nothing,Vector{Float64}}=nothing)
-
-    eq_type = eq_config.eq_type
-
-    # Rerun path: caller (e.g. build_inputs_from_h5) already rebuilt the solver input from
-    # the stored raw arrays. Bind it onto the current config and skip the reader.
-    if additional_input isa DirectRunInput
-        eq_input = additional_input
-        eq_input.config = eq_config
-        # Reset before re-resolving: a pass-1 clamped value must not shadow a psihigh
-        # overridden in the rerun TOML.
-        eq_input.psihigh_resolved = eq_config.psihigh
-        eq_type in EFIT_KINDS && resolve_psihigh!(eq_input)
-    elseif additional_input isa InverseRunInput
-        eq_input = additional_input
-        eq_input.config = eq_config
-        eq_input.psihigh_resolved = eq_config.psihigh
-        resolve_psihigh!(eq_input)
-    elseif eq_type in EFIT_KINDS
-        eq_input = resolve_psihigh!(read_efit(eq_config))
-    elseif eq_type in ["chease2", "chease_ascii"]
-        eq_input = read_chease_ascii(eq_config)
-    elseif eq_type in ["chease", "chease_binary"]
-        eq_input = read_chease_binary(eq_config)
-    elseif eq_type == "ifile"
-        eq_input = resolve_psihigh!(read_ifile(eq_config))
-    elseif haskey(ANALYTIC_EQ, eq_type)
-        # Analytic kinds (sol/lar/tj_analytic[_direct]) dispatch off the ANALYTIC_EQ registry.
-        # Their parameters live in the embedded `[*_INPUT]` section and are passed in as the
-        # `*Config` additional_input (built by build_analytic_config on the TOML/rerun paths).
-        spec = ANALYTIC_EQ[eq_type]
-        additional_input isa spec.config_type ||
-            error(
-                "setup_equilibrium: analytic eq_type=\"$eq_type\" requires its $(spec.config_type) " *
-                "passed as additional_input (built from the embedded [$(spec.section)] section)."
-            )
-        eq_input = spec.run_fn(eq_config, additional_input)
-    elseif eq_type == "imas"
-        if additional_input === nothing
-            error("setup_equilibrium: eq_type=\"imas\" requires an IMASdd.dd passed as additional_input")
-        end
-        eq_input = read_imas(eq_config, additional_input)
-    else
-        error("Equilibrium type $(eq_type) is not implemented")
-    end
+function PlasmaEquilibrium(eq_config::EquilibriumConfig, source=nothing; override_psi_nodes::Union{Nothing,Vector{Float64}}=nothing)
+    eq_input = run_input(eq_config, source)
 
-    if eq_type == "efit_by_inversion"
-        plasma_equilibrium = equilibrium_solver_by_inversion(eq_input; override_psi_nodes)
-    elseif eq_type == "efit_arclength"
-        plasma_equilibrium = equilibrium_solver(eq_input, arclength_fieldline_int; override_psi_nodes)
+    plasma_equilibrium = if eq_config.eq_type == "efit_by_inversion"
+        equilibrium_solver_by_inversion(eq_input; override_psi_nodes)
+    elseif eq_config.eq_type == "efit_arclength"
+        equilibrium_solver(eq_input, arclength_fieldline_int; override_psi_nodes)
     else
-        plasma_equilibrium = equilibrium_solver(eq_input; override_psi_nodes)
+        equilibrium_solver(eq_input; override_psi_nodes)
     end
 
     # Forward the captured ingest so the gpec.h5 writer can snapshot it (nothing for analytic).
     plasma_equilibrium.ingest = eq_input.ingest
+    # Analytic configs are still mutable; copy so later edits by the caller cannot change a re-form.
+    source isa AnalyticEquilibriumConfig && (plasma_equilibrium.analytic = deepcopy(source))
 
     equilibrium_global_parameters!(plasma_equilibrium)
     equilibrium_qfind!(plasma_equilibrium)
@@ -150,6 +155,27 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin
     return plasma_equilibrium
 end
 
+"""
+    reform(equil::PlasmaEquilibrium, psi_nodes::Vector{Float64}) -> PlasmaEquilibrium
+
+Re-form `equil` on the ψ_N nodes `psi_nodes` from the source it remembers: the stored
+`analytic` config, else its raw `ingest`, so no file is re-read. Only an equilibrium carrying
+neither (one built with the positional constructor, or by a reader that captures no ingest)
+falls back to re-reading `config.eq_filename`. The node endpoints must match `psilow`/`psihigh`.
+"""
+function reform(equil::PlasmaEquilibrium, psi_nodes::Vector{Float64})
+    source = if equil.analytic !== nothing
+        equil.analytic
+    elseif equil.ingest isa DirectIngest
+        build_direct_from_ingest(equil.config, equil.ingest)
+    elseif equil.ingest isa InverseIngest
+        build_inverse_from_ingest(equil.config, equil.ingest)
+    else
+        nothing
+    end
+    return PlasmaEquilibrium(equil.config, source; override_psi_nodes=psi_nodes)
+end
+
 """
     PlasmaEquilibrium(path::AbstractString; eq_type="efit", kwargs...) -> PlasmaEquilibrium
 
@@ -157,9 +183,8 @@ Read the equilibrium file at `path` and return the processed equilibrium. Conven
 point of the scripting API: `kwargs` are [`EquilibriumConfig`](@ref) fields, so
 `PlasmaEquilibrium("g000001.00001"; jac_type="hamada", mpsi=128)` is the whole setup.
 
-Only file-based equilibria go through this constructor. Analytic kinds (`sol`, `lar`,
-`tj_analytic`) take their parameters from a separate config object and are built with
-`setup_equilibrium(config, analytic_config)` instead.
+Only file-based equilibria go through this constructor. Analytic and IMAS equilibria have
+their own `PlasmaEquilibrium(::AnalyticEquilibriumConfig)` and `PlasmaEquilibrium(::IMASdd.dd)` forms.
 
 ```julia
 eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
@@ -167,8 +192,40 @@ eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
 """
 function PlasmaEquilibrium(path::AbstractString; eq_type::String="efit", kwargs...)
     haskey(ANALYTIC_EQ, eq_type) &&
-        error("$eq_type is an analytic equilibrium: build it with setup_equilibrium(config, $(ANALYTIC_EQ[eq_type].config_type)(...)) instead")
-    return setup_equilibrium(EquilibriumConfig(; eq_type, eq_filename=abspath(path), kwargs...))
+        error("$eq_type is an analytic equilibrium: build it with PlasmaEquilibrium($(ANALYTIC_EQ[eq_type].config_type)(...)) instead")
+    return PlasmaEquilibrium(EquilibriumConfig(; eq_type, eq_filename=abspath(path), kwargs...))
+end
+
+# Default analytic eq_type for a config: the first registry kind, in sorted order, that takes it.
+_analytic_eq_type(cfg::AnalyticEquilibriumConfig) = first(k for k in sort!(collect(keys(ANALYTIC_EQ))) if ANALYTIC_EQ[k].config_type === typeof(cfg))
+
+"""
+    PlasmaEquilibrium(cfg::AnalyticEquilibriumConfig; eq_type=<from ANALYTIC_EQ>, kwargs...) -> PlasmaEquilibrium
+
+Solve the analytic equilibrium described by `cfg` and return the processed equilibrium. The
+`eq_type` is implied by the config type; `TJAnalyticConfig` defaults to `tj_analytic` and takes
+`eq_type="tj_analytic_direct"` for the direct solve. `kwargs` are [`EquilibriumConfig`](@ref) fields.
+
+```julia
+eq = PlasmaEquilibrium(Equilibrium.SolovevConfig(; elongation=1.5); mpsi=0, psi_accuracy=1e-3)
+```
+"""
+function PlasmaEquilibrium(cfg::AnalyticEquilibriumConfig; eq_type::String=_analytic_eq_type(cfg), kwargs...)
+    return PlasmaEquilibrium(EquilibriumConfig(; eq_type, kwargs...), cfg)
+end
+
+"""
+    PlasmaEquilibrium(dd::IMASdd.dd; kwargs...) -> PlasmaEquilibrium
+
+Read the active equilibrium time slice of the IMAS data dictionary `dd` and return the processed
+equilibrium. `kwargs` are [`EquilibriumConfig`](@ref) fields; `imas_cocos` sets the ψ convention.
+
+```julia
+eq = PlasmaEquilibrium(dd; imas_cocos=11, jac_type="hamada")
+```
+"""
+function PlasmaEquilibrium(dd::IMASdd.dd; kwargs...)
+    return PlasmaEquilibrium(EquilibriumConfig(; eq_type="imas", kwargs...), dd)
 end
 
 """
diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index 31012ab0af..6f64faf894 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -32,7 +32,7 @@ specified in the input.
     equilibrium is actually formed on is `EquilibriumParameters.psihigh_resolved`.
   - `mpsi::Int` - Number of radial grid intervals; 0 with grid_type="auto" selects the
     two-pass auto grid: the main driver forms a coarse pass-1 equilibrium, measures its curvature,
-    pins knots on rational surfaces, and re-forms on the refined grid. Standalone `setup_equilibrium`
+    pins knots on rational surfaces, and re-forms on the refined grid. Standalone `PlasmaEquilibrium`
     callers get the coarse pass-1 grid unless they refine via `refined_psi_grid` + `override_psi_nodes`.
   - `psi_accuracy::Float64` - Target relative accuracy τ of splined profile derivatives for the
     two-pass auto grid (knot count scales as τ^(-1/3))
@@ -190,28 +190,13 @@ function EquilibriumConfig(equil_dict::Dict{String,Any}, base_path::String="./")
 end
 
 """
-Outer constructor for EquilibriumConfig that enables a toml file
-interface for specifying the configuration settings
+    AnalyticEquilibriumConfig
 
-DEPRECATED: Use [Equilibrium] section in gpec.toml instead
+Supertype of the parameter structs of the analytic equilibrium kinds. Each subtype is the
+`config_type` of one or more `ANALYTIC_EQ` rows, which pair it with an `eq_type`, its TOML
+section, and its solver.
 """
-function EquilibriumConfig(path::String)
-    raw = TOML.parsefile(path)
-
-    # Extract EQUIL_CONTROL with default fallback
-    config_data = get(raw, "EQUIL_CONTROL", Dict())
-
-    # Check for required fields
-    required_keys = ("eq_filename", "eq_type")
-    missingkeys = filter(k -> !haskey(config_data, k), required_keys)
-
-    if !isempty(missingkeys)
-        error("Missing required key(s) in [EQUIL_CONTROL]: $(join(missingkeys, ", "))")
-    end
-
-    # Construct validated struct
-    return EquilibriumConfig(Dict{String,Any}(config_data), dirname(path))
-end
+abstract type AnalyticEquilibriumConfig end
 
 """
     LargeAspectRatioConfig(...)
@@ -231,7 +216,7 @@ A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equi
   - `ma`: The number of grid points in the radial direction.
   - `zeroth`: If set to true, it neglects the Shafranov shift
 """
-@kwdef mutable struct LargeAspectRatioConfig
+@kwdef mutable struct LargeAspectRatioConfig <: AnalyticEquilibriumConfig
     r0::Float64 = 10.0
     a::Float64 = 1.0
     beta0::Float64 = 1e-3
@@ -275,7 +260,7 @@ profile is p₂(r) = p0·(1-r²)^μ with μ = `pressure_peaking`.
 
 Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ
 """
-@kwdef mutable struct TJAnalyticConfig
+@kwdef mutable struct TJAnalyticConfig <: AnalyticEquilibriumConfig
     r0::Float64 = 10.0         # Major radius R₀ [m]
     a::Float64 = 1.0           # Minor radius a [m] (ε = a/R₀)
     q0::Float64 = 1.5          # On-axis safety factor
@@ -313,7 +298,7 @@ A mutable struct holding parameters for the Solev'ev (SOL) plasma equilibrium mo
   - `b0fac`: scale toroidal field at constant beta (s*Phi,s*f,s^2*P. bt changes. Shape,beta constant)
   - `f0fac`: scale toroidal field at constant pressure (s*f. beta,q changes. Phi,p,bp constant)
 """
-@kwdef mutable struct SolovevConfig
+@kwdef mutable struct SolovevConfig <: AnalyticEquilibriumConfig
     mr::Int = 128      # number of radial grid zones
     mz::Int = 128      # number of axial grid zones
     ma::Int = 128      # number of flux grid zones
@@ -536,7 +521,7 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags
   - `zo::Union{Nothing,Float64}` - Z-coordinate of the magnetic axis [m]
   - `psio::Union{Nothing,Float64}` - Total flux difference |ψ_axis - ψ_boundary| [Wb/rad]
   - `psihigh_resolved::Union{Nothing,Float64}` - Outer flux limit the equilibrium was formed on
-    (the outermost ψ node); the plasma edge downstream of `setup_equilibrium`.
+    (the outermost ψ node); the plasma edge downstream of `PlasmaEquilibrium`.
   - `rsep::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma boundary [m]
   - `zsep::Union{Nothing,Vector{Float64}}` - Z-coordinates of the plasma boundary [m]
   - `rext::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma edge [m]
@@ -894,6 +879,10 @@ This object provides a complete representation of the processed plasma equilibri
   - `ingest::EquilibriumIngest`: raw arrays forwarded from the equilibrium input for the
     `gpec.h5` rerun snapshot — a [`DirectIngest`](@ref)/[`InverseIngest`](@ref) for file-based
     equilibria, or `nothing` for analytic ones (regenerated from their TOML section on replay)
+  - `analytic::Union{Nothing,AnalyticEquilibriumConfig}`: the parameters an analytic kind
+    (`sol`, `lar`, `tj_analytic`, `tj_analytic_direct`) was solved from, so [`reform`](@ref) can
+    re-solve it on a new grid; `nothing` for file-based and IMAS equilibria. A defensive copy, so
+    later edits to the caller's config do not reach it.
 """
 mutable struct PlasmaEquilibrium{P<:ProfileSplines,G<:GeometryProfileSplines,I2D<:FastInterpolations.CubicInterpolantND}
     config::EquilibriumConfig
@@ -921,13 +910,14 @@ mutable struct PlasmaEquilibrium{P<:ProfileSplines,G<:GeometryProfileSplines,I2D
     psio::Float64
 
     ingest::EquilibriumIngest
+    analytic::Union{Nothing,AnalyticEquilibriumConfig}
 end
 
-# Solvers build the equilibrium before setup_equilibrium forwards eq_input.ingest, so allow
-# construction without it; ingest defaults to nothing and is assigned post-construction.
+# Solvers build the equilibrium before PlasmaEquilibrium records its source, so allow
+# construction without it; ingest and analytic default to nothing and are assigned post-construction.
 PlasmaEquilibrium(config, params, profiles, geometry, rzphi_xs, rzphi_ys,
     rzphi_rsquared, rzphi_offset, rzphi_nu, rzphi_jac,
     eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio) =
     PlasmaEquilibrium(config, params, profiles, geometry, rzphi_xs, rzphi_ys,
         rzphi_rsquared, rzphi_offset, rzphi_nu, rzphi_jac,
-        eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio, nothing)
+        eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio, nothing, nothing)
diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl
index cc80c9cebe..9999dcd5f1 100644
--- a/src/Equilibrium/GridRefinement.jl
+++ b/src/Equilibrium/GridRefinement.jl
@@ -5,7 +5,7 @@ Pass 1 forms the equilibrium on a coarse grid; `refined_psi_grid` then derives a
 density from the measured nodal data (1D profiles, 2D geometry channels, and optional
 kinetic profiles) and equidistributes it, pinning mandatory knots (e.g. rational surfaces)
 via `merge_mandatory_nodes`. Pass 2 re-forms the equilibrium on the refined grid through
-the `override_psi_nodes` keyword of `setup_equilibrium`.
+the `override_psi_nodes` keyword of `PlasmaEquilibrium`.
 
 There is ONE sizing rule everywhere — the cubic *derivative* interpolation error model
 err(f′) ≈ h³|f''''|/24 ≤ τ·f_scale — applied to three sources of curvature |f''''|:
diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl
index f8436f26e2..31aa6e34a9 100644
--- a/src/Equilibrium/ReadEquilibrium.jl
+++ b/src/Equilibrium/ReadEquilibrium.jl
@@ -121,7 +121,7 @@ function read_efit(config::EquilibriumConfig)
     psi_in = cubic_interp((psi_in_xs, psi_in_ys), psi_proc; extrap=ExtendExtrap())
 
     # Capture the raw arrays that reconstruct sq_in and psi_in, so the rerun path
-    # (gpec.h5 → setup_equilibrium) can skip the g-file parse entirely.
+    # (gpec.h5 → PlasmaEquilibrium) can skip the g-file parse entirely.
     ingest = DirectIngest(sq_xs, sq_fs_nodes, psi_in_xs, psi_in_ys, psi_proc, rmin, rmax, zmin, zmax, psio, fpol_sign, ip_sign)
 
     return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, fpol_sign, ip_sign, ingest)
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 251e6da84c..ee815f5047 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -138,9 +138,9 @@ end
 
 Build the pipeline inputs from a working directory containing `gpec.toml`. Returns the
 parsed `inputs` dict, the `EquilibriumConfig`, and the `additional_input` consumed by
-`setup_equilibrium` — an analytic `*Config` for
+`PlasmaEquilibrium` — an analytic `*Config` for
 sol/lar/tj equilibria (parameters from the embedded TOML section), the `dd` data dictionary
-for IMAS, or `nothing` for file-based equilibria (efit, chease) that `setup_equilibrium`
+for IMAS, or `nothing` for file-based equilibria (efit, chease) that `PlasmaEquilibrium`
 reads from disk.
 """
 function build_inputs_from_toml(path::String; dd::Union{IMASdd.dd,Nothing}=nothing)
@@ -171,7 +171,7 @@ Shared pipeline body that every input source funnels into. The caller (one of th
 `build_inputs_from_*` builders) is responsible for producing a fully merged
 `inputs::Dict`, an `EquilibriumConfig`, and a ready `additional_input` — this body
 does no source dispatch of its own. `additional_input` is forwarded directly to
-`setup_equilibrium`: a prebuilt `DirectRunInput`/`InverseRunInput` (rerun path), an
+`PlasmaEquilibrium`: a prebuilt `DirectRunInput`/`InverseRunInput` (rerun path), an
 analytic `*Config` or IMAS `dd` (TOML path), or `nothing` for file-based equilibria.
 
 `preloaded_forcing_modes` lets the rerun path inject a `Vector{ForcingMode}`
@@ -220,11 +220,10 @@ function main_from_inputs(
 
     resolve_mode_space!(intr, ctrl)
 
-    equil = Equilibrium.setup_equilibrium(eq_config, additional_input)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
 
     kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil)
-    equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles;
-        kf_ctrl=kf_ctrl, kf_species=kf_species)
+    equil = maybe_reform_equilibrium(equil, intr, ctrl, kinetic_profiles; kf_ctrl=kf_ctrl, kf_species=kf_species)
 
     equil_dt = time() - equil_start
     push!(runtimes, "equilibrium" => equil_dt)
@@ -457,23 +456,22 @@ function load_kinetic_context(
 end
 
 """
-    maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles;
-                             kf_ctrl=nothing, kf_species=nothing) -> equil
+    maybe_reform_equilibrium(equil, intr, ctrl, kinetic_profiles; kf_ctrl=nothing, kf_species=nothing) -> equil
 
 Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, kinetic
-profiles), pin knots on rational surfaces, and re-form on the refined grid from the in-memory
-input — no file re-read. Returns `equil` untouched when the configuration wants a single pass.
+profiles), pin knots on rational surfaces, and re-form on the refined grid with
+`Equilibrium.reform` from the source the equilibrium remembers. Returns `equil` untouched when
+its configuration wants a single pass.
 """
 function maybe_reform_equilibrium(
     equil::Equilibrium.PlasmaEquilibrium,
-    eq_config::Equilibrium.EquilibriumConfig,
-    additional_input,
     intr::ForceFreeStatesInternal,
     ctrl::ForceFreeStatesControl,
     kinetic_profiles;
     kf_ctrl::Union{Nothing,KineticForces.KineticForcesControl}=nothing,
     kf_species=nothing
 )
+    eq_config = equil.config
     Equilibrium.wants_two_pass(eq_config) || return equil
 
     mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh)
@@ -489,19 +487,7 @@ function maybe_reform_equilibrium(
     psi_nodes = Equilibrium.refined_psi_grid(equil;
         tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, pinned=pinned,
         singfac_min=ctrl.singfac_min, n_min=n_min)
-    rerun_input = if additional_input !== nothing
-        # Analytic *Config, IMAS dd, or prebuilt RunInput — all re-formable. The IMAS
-        # path re-runs read_imas, which must resolve the same psihigh both passes;
-        # _validate_psi_nodes errors loudly if it does not.
-        additional_input
-    elseif equil.ingest isa Equilibrium.DirectIngest
-        Equilibrium.build_direct_from_ingest(eq_config, equil.ingest)
-    elseif equil.ingest isa Equilibrium.InverseIngest
-        Equilibrium.build_inverse_from_ingest(eq_config, equil.ingest)
-    else
-        nothing  # fall back to re-reading the input file
-    end
-    equil = Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes)
+    equil = Equilibrium.reform(equil, psi_nodes)
     implied = Equilibrium.implied_knot_count(equil; tau=eq_config.psi_accuracy, kin=kinetic_profiles)
     if implied > 1.5 * (length(psi_nodes) - 1)
         @warn "Two-pass psi grid: refined equilibrium implies $implied knots vs $(length(psi_nodes) - 1) used — " *
@@ -957,12 +943,7 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat
     # placeholder `prepare_force_free_states!` threads into its (unused) callback.
     kf_ctrl = KineticForces.KineticForcesControl()
 
-    if Equilibrium.wants_two_pass(equil.config) && equil.ingest === nothing
-        @warn "Two-pass auto grid needs the equilibrium's raw ingest, which analytic and IMAS equilibria do not carry; " *
-              "solving on the single-pass grid. Set mpsi explicitly to choose the grid."
-    else
-        equil = maybe_reform_equilibrium(equil, equil.config, nothing, intr, ctrl, nothing)
-    end
+    equil = maybe_reform_equilibrium(equil, intr, ctrl, nothing)
 
     locstab, ballooning_boundary = run_local_stability(ctrl, equil)
     metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing)
diff --git a/src/Rerun.jl b/src/Rerun.jl
index 687f729155..a2a5d7012a 100644
--- a/src/Rerun.jl
+++ b/src/Rerun.jl
@@ -8,7 +8,7 @@
 # Analytic equilibria carry their parameters in an embedded TOML section rather than an ingest
 # array; replay regenerates them from that section. The eq_type -> section mapping (and which
 # kinds are analytic at all) lives in the `Equilibrium.ANALYTIC_EQ` registry, the single source
-# of truth shared with `setup_equilibrium`.
+# of truth shared with `Equilibrium.run_input`.
 
 """
     build_analytic_config(eq_type, inputs)
@@ -277,7 +277,7 @@ end
 """
     rebuild_equilibrium_inputs(inputs, ingest, output_dir) -> (eq_config, additional_input)
 
-The equilibrium configuration and the input `setup_equilibrium` consumes, rebuilt from a stored
+The equilibrium configuration and the input `PlasmaEquilibrium` consumes, rebuilt from a stored
 TOML `inputs` dict and the equilibrium `ingest` read from a `gpec.h5`. Analytic kinds regenerate
 from their TOML section; file-based kinds rebuild splines from the stored ingest. A file-based
 run with no ingest can only come from a pre-ingest `gpec.h5` and is an error. Mutates `inputs`
@@ -310,7 +310,7 @@ end
     equilibrium_from_h5(h5path) -> (; equil, inputs, psilim)
 
 Rebuild the `PlasmaEquilibrium` of a finished run from its `gpec.h5` alone — the stored TOML
-and equilibrium ingest, through `setup_equilibrium` on the ψ_N grid the run ended with
+and equilibrium ingest, through `PlasmaEquilibrium` on the ψ_N grid the run ended with
 (`Equilibrium/Geometry/psi`) — without running any stage or writing anything. `inputs` is the
 parsed TOML the run used and `psilim` the control surface the solve integrated to
 (`Info/psilim`), so post-hoc analyses can evaluate new coil geometry on the surface the stored
@@ -329,7 +329,7 @@ function equilibrium_from_h5(h5path::AbstractString)
     end
     inputs = TOML.parse(toml_raw)
     eq_config, additional_input = rebuild_equilibrium_inputs(inputs, ingest, dirname(abspath(h5path)))
-    equil = Equilibrium.setup_equilibrium(eq_config, additional_input; override_psi_nodes=psi_nodes)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config, additional_input; override_psi_nodes=psi_nodes)
     return (; equil, inputs, psilim)
 end
 
diff --git a/test/runtests_coils.jl b/test/runtests_coils.jl
index b6f0c75c24..da2bf9bc55 100644
--- a/test/runtests_coils.jl
+++ b/test/runtests_coils.jl
@@ -166,7 +166,7 @@ end
 
     inputs = TOML.parsefile(joinpath(equil_dir, "gpec.toml"))
     eq_config = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], equil_dir)
-    equil = Equilibrium.setup_equilibrium(eq_config, Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
+    equil = Equilibrium.PlasmaEquilibrium(eq_config, Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
 
     mtheta = 48
     nzeta = 16
@@ -305,7 +305,7 @@ end
 
     inputs2 = TOML.parsefile(joinpath(equil_dir, "gpec.toml"))
     eq_config2 = Equilibrium.EquilibriumConfig(inputs2["Equilibrium"], equil_dir)
-    equil = Equilibrium.setup_equilibrium(eq_config2)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config2)
 
     # Build coil sets manually (standard I-coil currents for RMP studies)
     il_cfg = ForcingTerms.CoilSetConfig(;
@@ -458,7 +458,7 @@ end
 @testset "Analytic coils: window_pane standoff" begin
     equil_dir = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example")
     inputs = TOML.parsefile(joinpath(equil_dir, "gpec.toml"))
-    equil = Equilibrium.setup_equilibrium(Equilibrium.EquilibriumConfig(inputs["Equilibrium"], equil_dir),
+    equil = Equilibrium.PlasmaEquilibrium(Equilibrium.EquilibriumConfig(inputs["Equilibrium"], equil_dir),
         Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
 
     # surface_point_and_normal: outboard midplane (poloidal_angle = 0)
@@ -649,8 +649,8 @@ end
         end
     end
     cfg(path) = Equilibrium.EquilibriumConfig(; eq_filename=path, eq_type="efit", jac_type="hamada", grid_type="core_edge_packed", psilow=0.01, psihigh=0.99, mpsi=32, mtheta=64)
-    eq_pos = Equilibrium.setup_equilibrium(cfg(gfile))
-    eq_neg = Equilibrium.setup_equilibrium(cfg(flipped))
+    eq_pos = Equilibrium.PlasmaEquilibrium(cfg(gfile))
+    eq_neg = Equilibrium.PlasmaEquilibrium(cfg(flipped))
     @test eq_pos.params.ip_sign == 1 && eq_neg.params.ip_sign == -1
     @test eq_pos.params.bt_sign == eq_neg.params.bt_sign == -1
     @test eq_pos.params.crnt > 0 && eq_neg.params.crnt > 0          # the computed current carries no sign
diff --git a/test/runtests_coordinate_invariant.jl b/test/runtests_coordinate_invariant.jl
index 64300b0dde..42787b45e9 100644
--- a/test/runtests_coordinate_invariant.jl
+++ b/test/runtests_coordinate_invariant.jl
@@ -108,7 +108,7 @@ end
     equil_dir = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example")
     inputs = TOML.parsefile(joinpath(equil_dir, "gpec.toml"))
     eq_config = EQ.EquilibriumConfig(inputs["Equilibrium"], equil_dir)
-    equil = EQ.setup_equilibrium(eq_config, EQ.SolovevConfig(inputs["SOL_INPUT"]))
+    equil = EQ.PlasmaEquilibrium(eq_config, EQ.SolovevConfig(inputs["SOL_INPUT"]))
     psi = equil.rzphi_xs[end]
     mtheta = length(equil.rzphi_ys)
     mpert, mlow = 21, -8
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 66e7844450..d21fc477d0 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -15,7 +15,7 @@ using Statistics
             psilow=0.01,
             psihigh=0.994
         )
-        global plasma_eq_efit = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(efit_config)
+        global plasma_eq_efit = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(efit_config)
 
         @test plasma_eq_efit isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium
         @test 6.5 < plasma_eq_efit.ro < 7.5 # Physical sanity check
@@ -30,7 +30,7 @@ using Statistics
             psilow=0.01,
             psihigh=0.994
         )
-        global plasma_eq_arclength = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(arclength_config)
+        global plasma_eq_arclength = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(arclength_config)
 
         @test plasma_eq_arclength isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium
         @test 6.5 < plasma_eq_arclength.ro < 7.5
@@ -54,7 +54,7 @@ using Statistics
             psilow=0.01,
             psihigh=0.994
         )
-        global plasma_eq_inversion = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(inversion_config)
+        global plasma_eq_inversion = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(inversion_config)
 
         @test plasma_eq_inversion isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium
         @test 6.5 < plasma_eq_inversion.ro < 7.5
@@ -95,7 +95,7 @@ using Statistics
             psilow=0.01,
             psihigh=1.0
         )
-        plasma_eq_edge = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(edge_config)
+        plasma_eq_edge = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(edge_config)
         @test edge_config.psihigh == 1.0
         @test plasma_eq_edge.params.psihigh_resolved <= 1.0
         @test plasma_eq_edge.params.psihigh_resolved == plasma_eq_edge.rzphi_xs[end]
@@ -116,19 +116,19 @@ using Statistics
         )
 
         base_config = newq0_config(0)
-        eq_base = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(base_config)
+        eq_base = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(base_config)
 
         # -1 means "flip the sign of the axis extrapolation": f0fac vanishes, so ffac is
         # exactly -1 and the revised q-profile is the negated baseline.
         flip_config = newq0_config(-1)
-        eq_flip = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(flip_config)
+        eq_flip = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(flip_config)
         @test flip_config.newq0 == -1.0
         @test eq_flip.profiles.q_spline.y == -eq_base.profiles.q_spline.y
 
         # A non-integer target used to throw InexactError at construction (newq0 was ::Int)
         target_config = newq0_config(1.05)
         @test target_config.newq0 == 1.05
-        eq_target = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(target_config)
+        eq_target = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(target_config)
         @test all(isfinite, eq_target.profiles.q_spline.y)
         # Extrapolated q(0) meets the target; not exact because ffac uses the first node's F
         # rather than the extrapolated axis value.
@@ -183,13 +183,13 @@ using Statistics
         finite(eq) = all(isfinite, eq.profiles.dVdpsi_spline.y) && all(all(isfinite, s.nodal_derivs.partials) for s in splines(eq))
 
         # original: the file's own ψ nodes between psilow and psihigh, whatever mpsi
-        eq_o = Eq.setup_equilibrium(cfg("original"))
+        eq_o = Eq.PlasmaEquilibrium(cfg("original"))
         xs_in = Eq.read_chease_binary(cfg("original")).rz_in_xs
         @test eq_o.rzphi_xs == [0.01; filter(x -> 0.01 < x < 0.99, xs_in); 0.99]
         @test finite(eq_o)
 
         # rho: starts at psilow (Fortran's rho starts at psihigh/(mpsi+1)²)
-        eq_r = Eq.setup_equilibrium(cfg("rho"))
+        eq_r = Eq.PlasmaEquilibrium(cfg("rho"))
         @test eq_r.rzphi_xs[[1, end]] ≈ [0.01, 0.99]
         @test finite(eq_r)
     end
@@ -200,9 +200,9 @@ using Statistics
         tk_dir = joinpath(@__DIR__, "test_data", "TokaMaker_ifile")
         cfg(type, file; kw...) = Eq.EquilibriumConfig(; eq_type=type, eq_filename=joinpath(tk_dir, file), jac_type="hamada",
             grid_type="core_edge_packed", mpsi=64, psilow=0.01, psihigh=0.99, kw...)
-        eq_i = Eq.setup_equilibrium(cfg("ifile", "i33x65.ifile"))
-        eq_g = Eq.setup_equilibrium(cfg("efit", "g65.geqdsk"))
-        eq_s = Eq.setup_equilibrium(cfg("ifile", "i33x65_single.ifile"))
+        eq_i = Eq.PlasmaEquilibrium(cfg("ifile", "i33x65.ifile"))
+        eq_g = Eq.PlasmaEquilibrium(cfg("efit", "g65.geqdsk"))
+        eq_s = Eq.PlasmaEquilibrium(cfg("ifile", "i33x65_single.ifile"))
 
         # Same equilibrium as the g-file of the same solve
         @test isapprox(eq_i.ro, eq_g.ro; atol=1e-4) && isapprox(eq_i.zo, eq_g.zo; atol=1e-4)
@@ -218,6 +218,16 @@ using Statistics
         @test all(isapprox(eq_s.profiles.q_spline(x), eq_i.profiles.q_spline(x); rtol=1e-6) for x in (0.2, 0.5, 0.8))
         @test maximum(abs.(eq_s.ingest.R_nodes .- eq_i.ingest.R_nodes)) < 1e-6
 
+        # Re-forming from the captured ingest matches re-reading the file on the same nodes
+        @test eq_i.ingest isa Eq.InverseIngest
+        nodes = collect(range(0.01, 0.99; length=33))
+        reformed = Eq.reform(eq_i, nodes)
+        fresh = Eq.PlasmaEquilibrium(cfg("ifile", "i33x65.ifile"); override_psi_nodes=nodes)
+        @test reformed.rzphi_xs == fresh.rzphi_xs == nodes
+        @test map(reformed.profiles.q_spline, nodes) == map(fresh.profiles.q_spline, nodes)
+        @test all(reformed.rzphi_rsquared(pt) == fresh.rzphi_rsquared(pt) && reformed.rzphi_jac(pt) == fresh.rzphi_jac(pt)
+                  for pt in ((0.1, 0.0), (0.5, 0.3), (0.9, 0.7)))
+
         # The trailing FF′ and p′ records are ignored: a file without them reads the same
         table(file) = Eq.read_ifile(cfg("ifile", file)).ingest.sq_fs
         tabulated = table("i33x65.ifile")
@@ -234,7 +244,7 @@ using Statistics
         bumped_file = joinpath(mktempdir(), "bumped_q.ifile")
         write(bumped_file, bumped)
         @test Eq.read_ifile(cfg("ifile", bumped_file)).ingest.sq_fs[:, 3] ≈ 1.1 .* tabulated[:, 3]
-        @test Eq.setup_equilibrium(cfg("ifile", bumped_file)).profiles.q_spline.y == eq_i.profiles.q_spline.y
+        @test Eq.PlasmaEquilibrium(cfg("ifile", bumped_file)).profiles.q_spline.y == eq_i.profiles.q_spline.y
     end
 
     @testset "TokaMaker i-file ψ bounds record" begin
@@ -268,10 +278,10 @@ using Statistics
         @test ing_sb.sq_fs[:, 1:3] == ing_db.sq_fs[:, 1:3] == base_s.sq_fs[:, 1:3]
 
         # psihigh beyond the last surface is clamped to it; a rerun from the ingest keeps the normalisation
-        eq = @test_logs (:warn, r"clamped") match_mode = :any Eq.setup_equilibrium(cfg(db; psihigh=0.995))
+        eq = @test_logs (:warn, r"clamped") match_mode = :any Eq.PlasmaEquilibrium(cfg(db; psihigh=0.995))
         @test eq.params.psihigh_resolved == ing_db.sq_xs[end] ≈ eq.rzphi_xs[end]
         @test all(isfinite, eq.profiles.q_spline.y) && all(isfinite, eq.eqfun_B.nodal_derivs.partials)
-        rerun = Eq.setup_equilibrium(cfg(db; psihigh=0.995), Eq.build_inverse_from_ingest(cfg(db), eq.ingest))
+        rerun = Eq.PlasmaEquilibrium(cfg(db; psihigh=0.995), Eq.build_inverse_from_ingest(cfg(db), eq.ingest))
         @test rerun.rzphi_xs == eq.rzphi_xs && rerun.psio == eq.psio
         @test rerun.profiles.q_spline.y == eq.profiles.q_spline.y
     end
@@ -304,7 +314,7 @@ using Statistics
             r0_norm=6.8,
             b0_norm=7.4
         )
-        global plasma_eq_binary = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(binary_config)
+        global plasma_eq_binary = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(binary_config)
 
         @test plasma_eq_binary isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium
     end
@@ -320,7 +330,7 @@ using Statistics
             r0_norm=6.8,
             b0_norm=7.4
         )
-        global plasma_eq_ascii = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(ascii_config)
+        global plasma_eq_ascii = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(ascii_config)
 
         @test plasma_eq_ascii isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium
     end
@@ -417,7 +427,7 @@ using Statistics
             r0=10.0, a=1.0, q0=1.5, beta0=1e-4,
             mtheta=64, ma=64, zeroth=true
         )
-        plasma_eq_lar = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(lar_config, lar_input)
+        plasma_eq_lar = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(lar_config, lar_input)
 
         @test plasma_eq_lar isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium
 
diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl
index e160032af2..9c6202fc84 100644
--- a/test/runtests_eulerlagrange.jl
+++ b/test/runtests_eulerlagrange.jl
@@ -97,7 +97,7 @@ end
         F = GeneralizedPerturbedEquilibrium.ForceFreeStates
         ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
         inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
-        equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(
             GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex),
             GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
 
@@ -416,7 +416,7 @@ end
             eq_inputs = copy(inputs["Equilibrium"])
             eq_inputs["psilow"] = psilow
             eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(eq_inputs, ex)
-            equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
+            equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
             ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., kwargs...)
             intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
             intr.nlow = ctrl.nn_low
@@ -561,7 +561,7 @@ end
         ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
         eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], example_dir)
         sol_cfg = haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing
-        equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, sol_cfg)
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, sol_cfg)
         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
diff --git a/test/runtests_grid_refinement.jl b/test/runtests_grid_refinement.jl
index 66d34332ad..ed0fd77bdd 100644
--- a/test/runtests_grid_refinement.jl
+++ b/test/runtests_grid_refinement.jl
@@ -153,7 +153,7 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium
             eq_type="sol", grid_type="log_asymptotic", mpsi=0, psi_accuracy=1e-3)
         @test GridRef.wants_two_pass(alias_config)
 
-        equil1 = GridRef.setup_equilibrium(eq_config, sol_config)
+        equil1 = GridRef.PlasmaEquilibrium(eq_config, sol_config)
 
         # Rational nodes land exactly on q = m/n surfaces
         rat = GeneralizedPerturbedEquilibrium.ForceFreeStates.rational_psi_nodes(equil1; nlow=1, nhigh=2)
@@ -179,10 +179,27 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium
         @test minimum(diff(grid)) > 0
 
         # Pass 2 honors the refined grid exactly
-        equil2 = GridRef.setup_equilibrium(eq_config, sol_config; override_psi_nodes=grid)
+        equil2 = GridRef.PlasmaEquilibrium(eq_config, sol_config; override_psi_nodes=grid)
         @test equil2.profiles.xs == grid
         @test equil2.rzphi_xs == grid
 
+        # The equilibrium keeps a private copy of its config, so reform needs no input and
+        # ignores later edits to the caller's config.
+        @test equil1.analytic isa GridRef.SolovevConfig && equil1.analytic !== sol_config
+        sol_config.elongation += 0.1
+        equil2r = GridRef.reform(equil1, grid)
+        sol_config.elongation -= 0.1
+        @test equil2r.rzphi_xs == equil2.rzphi_xs
+        @test map(equil2r.profiles.q_spline, grid) == map(equil2.profiles.q_spline, grid)
+        @test equil2r.rzphi_rsquared((0.5, 0.3)) == equil2.rzphi_rsquared((0.5, 0.3))
+
+        # The one-call constructor implies eq_type from the config type and gives the same pass 1.
+        equil1c = GridRef.PlasmaEquilibrium(sol_config; grid_type="auto", mpsi=0, psi_accuracy=1e-3, psilow=0.01, psihigh=0.995)
+        @test equil1c.config.eq_type == "sol" && equil1c.rzphi_xs == equil1.rzphi_xs
+        @test_throws ErrorException GridRef.PlasmaEquilibrium(sol_config; eq_type="lar")
+        @test_throws ErrorException GridRef.PlasmaEquilibrium(eq_config)  # analytic kind without its config
+        @test GridRef._analytic_eq_type(GridRef.TJAnalyticConfig()) == "tj_analytic"
+
         # Solovev has finite edge q, so the auto grid must classify it as limited: no separatrix
         # log law, on either pass, even with the absolute r² floor disabled (the relative
         # log-vs-linear comparison is what rejects it).
@@ -195,7 +212,7 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium
         using TOML
         dir_path = joinpath(dirname(@__DIR__), "examples", "DIIID-like_ideal_example")
         inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml"))
-        equil = GridRef.setup_equilibrium(GridRef.EquilibriumConfig(inputs["Equilibrium"], dir_path), nothing)
+        equil = GridRef.PlasmaEquilibrium(GridRef.EquilibriumConfig(inputs["Equilibrium"], dir_path), nothing)
         fit = GridRef.edge_q_law(equil)
         @test fit !== nothing
         @test fit.A < 0
diff --git a/test/runtests_imas.jl b/test/runtests_imas.jl
index 7970d6e4ba..de35804fcf 100644
--- a/test/runtests_imas.jl
+++ b/test/runtests_imas.jl
@@ -108,6 +108,26 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         end
     end
 
+    # Re-forming from the captured ingest must match re-reading the dd on the same nodes. The mock's
+    # surfaces close inside its R grid only to psihigh ≈ 0.5, below the auto grid's edge region.
+    @testset "reform from ingest matches a fresh read_imas" begin
+        dd, _ = make_mock_dd()
+        equil = PlasmaEquilibrium(dd; imas_cocos=11, grid_type="core_edge_packed", mpsi=64, psilow=0.01, psihigh=0.5)
+        @test equil.config.eq_type == "imas"
+        @test equil.ingest isa Equilibrium.DirectIngest && equil.analytic === nothing
+        @test_throws ErrorException PlasmaEquilibrium(EquilibriumConfig(; eq_type="efit"), dd)  # dd needs eq_type="imas"
+
+        nodes = Equilibrium.refined_psi_grid(equil; tau=1e-3, mandatory=Float64[], singfac_min=1e-4, n_min=1)
+        reformed = Equilibrium.reform(equil, nodes)
+        fresh = PlasmaEquilibrium(equil.config, dd; override_psi_nodes=nodes)
+        @test reformed.rzphi_xs == fresh.rzphi_xs == nodes
+        @test map(reformed.profiles.q_spline, nodes) == map(fresh.profiles.q_spline, nodes)
+        for pt in ((0.1, 0.0), (0.3, 0.2), (0.45, 0.7))
+            @test reformed.rzphi_rsquared(pt) == fresh.rzphi_rsquared(pt)
+            @test reformed.rzphi_jac(pt) == fresh.rzphi_jac(pt)
+        end
+    end
+
     # Test 4: write_imas — single n, correct metadata and energy_perturbed
     @testset "write_imas: single n_tor" begin
         dd = IMASdd.dd()
diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl
index 0cedfb9e9c..6b1684bbf5 100644
--- a/test/runtests_layer_overlap_fitzpatrick.jl
+++ b/test/runtests_layer_overlap_fitzpatrick.jl
@@ -151,7 +151,7 @@ end
 
     dir_path = joinpath(dirname(@__DIR__), "examples", "Solovev_ideal_example")
     inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml"))
-    equil = Equilibrium.setup_equilibrium(
+    equil = Equilibrium.PlasmaEquilibrium(
         Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path),
         Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
 
@@ -305,7 +305,7 @@ end
 
     trunc_dir = joinpath(dirname(@__DIR__), "examples", "DIIID-like_truncation_example")
     trunc_inputs = TOML.parsefile(joinpath(trunc_dir, "gpec.toml"))
-    diiid(psihigh) = Equilibrium.setup_equilibrium(
+    diiid(psihigh) = Equilibrium.PlasmaEquilibrium(
         Equilibrium.EquilibriumConfig(merge(trunc_inputs["Equilibrium"], Dict("psihigh" => psihigh)), trunc_dir), nothing)
     long, short = diiid(0.9999), diiid(0.995)
 
@@ -345,7 +345,7 @@ end
     # No 2D flux, no X-point: an analytic equilibrium is not extrapolated.
     sol_dir = joinpath(dirname(@__DIR__), "examples", "Solovev_ideal_example")
     sol_inputs = TOML.parsefile(joinpath(sol_dir, "gpec.toml"))
-    solovev = Equilibrium.setup_equilibrium(
+    solovev = Equilibrium.PlasmaEquilibrium(
         Equilibrium.EquilibriumConfig(sol_inputs["Equilibrium"], sol_dir), Equilibrium.SolovevConfig(sol_inputs["SOL_INPUT"]))
     @test separatrix_q_law(solovev) === nothing
 end
diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl
index b427c37958..e868bde7be 100644
--- a/test/runtests_parallel_integration.jl
+++ b/test/runtests_parallel_integration.jl
@@ -114,7 +114,7 @@ using TOML
         ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
             (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)
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
@@ -201,7 +201,7 @@ using TOML
         ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
             (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)
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
@@ -262,7 +262,7 @@ using TOML
             ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
                 (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)
+            equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(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"])...)
             intr.nlow = ctrl.nn_low;
@@ -311,7 +311,7 @@ using TOML
             ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
                 (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)
+            equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
             # Apply the two-pass auto grid exactly as the main driver does (the example ships
             # grid_type="auto", mpsi=0): measured-curvature refinement with rational surfaces
             # pinned as mandatory knots, re-formed from the captured ingest. The pinned values
@@ -320,7 +320,7 @@ using TOML
                 mand = GeneralizedPerturbedEquilibrium.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high)
                 psi_nodes = GeneralizedPerturbedEquilibrium.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand)
                 rerun_input = GeneralizedPerturbedEquilibrium.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest)
-                equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes)
+                equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes)
             end
             intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
                 (Symbol(k) => v for (k, v) in inputs["Wall"])...)
@@ -372,7 +372,7 @@ using TOML
         ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
             (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)
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
         intr.nlow = ctrl.nn_low;
         intr.nhigh = ctrl.nn_high;
         intr.npert = 1
@@ -423,7 +423,7 @@ using TOML
             ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
                 (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
             eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], example_dir)
-            equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
+            equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(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"])...)
             intr.nlow = ctrl.nn_low;
@@ -476,7 +476,7 @@ using TOML
         ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
             (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)
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing)
         # Apply the two-pass auto grid (measured-curvature refinement, rational surfaces pinned
         # as mandatory knots) exactly as the main driver does (see the FM testset above); the
         # pinned values below are for this grid, the production default.
@@ -484,7 +484,7 @@ using TOML
             mand = GeneralizedPerturbedEquilibrium.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high)
             psi_nodes = GeneralizedPerturbedEquilibrium.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand)
             rerun_input = GeneralizedPerturbedEquilibrium.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest)
-            equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes)
+            equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes)
         end
         intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
             (Symbol(k) => v for (k, v) in inputs["Wall"])...)
diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl
index 13172216bb..8d2658b914 100644
--- a/test/runtests_rerun_from_h5.jl
+++ b/test/runtests_rerun_from_h5.jl
@@ -304,7 +304,7 @@ end
     @test bad["ForceFreeStates"]["mode"] == "blorp"
 end
 
-# Drift guard for the analytic-equilibrium registry: setup_equilibrium and the rerun builder both
+# Drift guard for the analytic-equilibrium registry: `run_input` and the rerun builder both
 # dispatch off Equilibrium.ANALYTIC_EQ, so a malformed entry (missing Dict ctor, wrong run_fn,
 # empty section) would silently break one path. Asserting each entry is well-formed and that
 # build_analytic_config resolves it keeps adding a new analytic kind a one-row change.
diff --git a/test/runtests_resist_eval.jl b/test/runtests_resist_eval.jl
index f2181bfea9..2ddb37620b 100644
--- a/test/runtests_resist_eval.jl
+++ b/test/runtests_resist_eval.jl
@@ -15,7 +15,7 @@
     inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml"))
     eq_cfg = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path)
     sol_cfg = Equilibrium.SolovevConfig(inputs["SOL_INPUT"])
-    equil = Equilibrium.setup_equilibrium(eq_cfg, sol_cfg)
+    equil = Equilibrium.PlasmaEquilibrium(eq_cfg, sol_cfg)
 
     # Spitzer-like η at ~1 keV and deuterium ρ at n_e = 5e19 m⁻³.
     eta, rho = 3e-8, 2 * Utilities.PhysicalConstants.M_P * 5e19
diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl
index e5989aed8b..2140b5ac4f 100644
--- a/test/runtests_riccati.jl
+++ b/test/runtests_riccati.jl
@@ -101,7 +101,7 @@ end
 
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(
         GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex),
         GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
 
diff --git a/test/runtests_sfl_sampling.jl b/test/runtests_sfl_sampling.jl
index 332078510d..161a43c12a 100644
--- a/test/runtests_sfl_sampling.jl
+++ b/test/runtests_sfl_sampling.jl
@@ -44,7 +44,7 @@ using GeneralizedPerturbedEquilibrium
 
     @testset "EFIT: efit and efit_arclength tracers agree" begin
         efit_file = joinpath(@__DIR__, "test_data", "CHEASE_test_data", "EQDSK_COCOS_02")
-        efit_eq(kind) = Eq.setup_equilibrium(Eq.EquilibriumConfig(; eq_filename=efit_file, eq_type=kind,
+        efit_eq(kind) = Eq.PlasmaEquilibrium(Eq.EquilibriumConfig(; eq_filename=efit_file, eq_type=kind,
             jac_type="boozer", grid_type="core_edge_packed", psilow=0.01, psihigh=0.994, mpsi=64, mtheta=128))
         direct, arclength = efit_eq("efit"), efit_eq("efit_arclength")
         # Two integrators with different step sequences along the same surfaces. Measured: Jacobian 8.7e-8
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index b77309c9dd..c5f17a0499 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -84,7 +84,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         inputs = TOML.parsefile(gpec_path)
         eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dirname(gpec_path))
         addl = haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing
-        return GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, addl)
+        return GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config, addl)
     end
 
     @testset "sing_der" begin
diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl
index ccb881ecb3..ecc7a06bb4 100644
--- a/test/runtests_slayer_inputs.jl
+++ b/test/runtests_slayer_inputs.jl
@@ -12,7 +12,7 @@
     inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml"))
     eq_cfg = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path)
     sol_cfg = Equilibrium.SolovevConfig(inputs["SOL_INPUT"])
-    equil = Equilibrium.setup_equilibrium(eq_cfg, sol_cfg)
+    equil = Equilibrium.PlasmaEquilibrium(eq_cfg, sol_cfg)
 
     # Synthetic profiles (simple linear-in-ψ temperature decrease)
     psi_pts = collect(0.0:0.1:1.0)
diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl
index 2a238b3a38..184047139b 100644
--- a/test/runtests_solve_api.jl
+++ b/test/runtests_solve_api.jl
@@ -10,7 +10,7 @@ using TOML
     deck = TOML.parsefile(joinpath(template, "gpec.toml"))
 
     # The fixture equilibrium and wall, built exactly as `build_inputs_from_toml` builds them.
-    equil = GPEC.Equilibrium.setup_equilibrium(
+    equil = GPEC.Equilibrium.PlasmaEquilibrium(
         GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template),
         GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"])
     )
diff --git a/test/runtests_tj_analytic.jl b/test/runtests_tj_analytic.jl
index da080440c2..6b32dd90da 100644
--- a/test/runtests_tj_analytic.jl
+++ b/test/runtests_tj_analytic.jl
@@ -2,7 +2,7 @@ using Test
 using Printf
 using GeneralizedPerturbedEquilibrium.Equilibrium
 using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig,
-    setup_equilibrium, tj_analytic_run, tj_analytic_run_direct
+    tj_analytic_run, tj_analytic_run_direct
 
 # Two-path smoke tests for the TJ-analytic equilibrium model
 # (GPEC adaptation of R. Fitzpatrick's TJ code,
@@ -25,7 +25,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium
         eq = EquilibriumConfig(eq_type = "tj_analytic",
                                psilow = 0.01, psihigh = 0.995,
                                mpsi = 64, mtheta = 128, flux_surface_rtol = 1e-7)
-        pe = setup_equilibrium(eq, tj)
+        pe = PlasmaEquilibrium(eq, tj)
 
         # psio is a physical-scale ψ; regressions in the a→a² normalization
         # or the dψ/dr construction would change it by factors of a.
@@ -52,7 +52,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium
         eq = EquilibriumConfig(eq_type = "tj_analytic_direct",
                                psilow = 0.01, psihigh = 0.995,
                                mpsi = 64, mtheta = 128, flux_surface_rtol = 1e-7)
-        pe = setup_equilibrium(eq, tj)
+        pe = PlasmaEquilibrium(eq, tj)
 
         @test pe.psio > 0
         @test isfinite(pe.psio)