diff --git a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md
index d4733ea2d..1a3bc1fc1 100644
--- a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md
+++ b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md
@@ -4,3 +4,4 @@
 - [Known Physics Issues in PerturbedEquilibrium](known_issues_perteq.md) — xss_mn not implemented, xi_modes placeholders
 - [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify
 - [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify
+- [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS
diff --git a/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md b/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md
new file mode 100644
index 000000000..c916e3b00
--- /dev/null
+++ b/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md
@@ -0,0 +1,37 @@
+---
+name: galerkin-assembly-map
+description: RDCON singular Galerkin Δ′ solver — gal.f ↔ GalerkinAssembly.jl element-level weak-form correspondence and audit notes
+metadata:
+  type: project
+---
+
+# RDCON Galerkin Δ′ assembly: gal.f ↔ GalerkinAssembly.jl
+
+File: `src/ForceFreeStates/Galerkin/GalerkinAssembly.jl` ↔ `GPEC_fortran/rdcon/gal.f`.
+On branch `feature/rdcon-galerkin`. Reviewed 2026-06-16.
+
+## Correspondence (verified, not just docstring-trusted)
+- gal_hermite (222-262) ↔ gal_hermite: pb/qb EXACT.
+- gal_get_fkg (1013-1071) ↔ gal_get_fkg: F=Fbar·sf[i]·sf[j], K=K·sf[i], G=G. singfac=mlow-nn*q+(0..mpert-1) = DIRECT (m-nq), not reciprocal. Correct.
+- gal_gauss_quad (1079-1138) ↔ gal_gauss_quad!: F·qq + K·qp + conj(K[jpert,ipert])·pq + G·pp. swap pb(2)/pb(3) (0-based) = Julia pb[3]/pb[4]. Correct.
+- gal_extension (1146-1322) ↔ gal_extension!: emat/ediag/rhs/erhs all match incl surface terms 1277-1300 and w1/w2.
+- gal_resonant! ↔ gal_lsode_int(780-900)+gal_lsode_der(908-954)+gal_make_arrays signs(1366-1369). QuadGK replaces LSODE.
+- gal_assemble_mat (694-772) ↔ gal_assemble_mat!: LU offset=kl+ku+1, chol offset=1, conj(e) Hermitian placement. Correct.
+- gal_assemble_rhs (642-686) ↔ gal_assemble_rhs!: isol++ at ext2-left / res-right. Correct.
+- gal_set_boundary (1554-1620) ↔ gal_set_boundary!: axis idx, edge 3-way (rpec/vac/fixed). DOF index 3(0-based)=Julia 4.
+
+## Resonant sign chain (the trap — verified correct)
+Fortran: LSODE integrates x0→x_lsode giving raw u; then line 882 `IF right: u=-u`; then make_arrays:
+erhs=-u_res(1), ediag=+u_res(2), rhs=-u_hermite1, emat=+u_hermite2.
+Net with s=+1(left)/-1(right): erhs=-s·res1, ediag=s·res2, rhs=-s·hbig, emat=s·hsmall.
+Julia sgn_u = s (=-1 right/+1 left): erhs=-sgn_u·raw1, ediag=sgn_u·raw2, rhs=-sgn_u·hbig, emat=sgn_u·hsmall. MATCHES.
+QuadGK limits x0→x_lsode preserve LSODE orientation so raw == pre-negation u. Correct.
+du_res conj-weight uses the SMALL solution (isol(2)) in both. Correct.
+
+## Only open item (MINOR, out-of-file)
+sing_get_ua_gal/dua_gal two-sided z<0 → sig=-1 series at -z with derivative ×(-1) is a REIMPLEMENTATION
+of sing.f's internal singular-surface handling, not a line port. The ×(-1) chain is internally consistent
+but worth a dedicated check against sing.f sing_get_ua/dua if Δ′ ever disagrees with Fortran. Lives outside
+GalerkinAssembly.jl.
+
+## Verdict: PASS WITH REQUIRED ANNOTATIONS (citation comments only; no physics defects found).
diff --git a/benchmarks/compare_gal_vs_el.jl b/benchmarks/compare_gal_vs_el.jl
new file mode 100644
index 000000000..f3facb532
--- /dev/null
+++ b/benchmarks/compare_gal_vs_el.jl
@@ -0,0 +1,78 @@
+# Overlay the IDEAL gal matched ξ(ψ) against the EL total-energy eigenmode ξ(ψ), for the same edge
+# eigenvector w. The ideal gal solution should reproduce the EL (DCON) ideal solution.
+#
+# EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, integration/xi_psi)
+# gal:  ξ_gal(ψ) = U_gal(ψ) · w                       (identity-at-edge ⇒ coefficient is w itself)
+# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability/XiNorm).
+#
+# Needs ONE gpec.h5 from a run with populate_dense_xi=true (EL dense u_store), gal_match_flag=true,
+# gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator).
+# Usage: julia --project=. benchmarks/compare_gal_vs_el.jl [gpec.h5] [out.png] [mode]
+#   mode = "highest" (default, most stable), "lowest" (most unstable), or an integer eigenmode index.
+
+using HDF5, LinearAlgebra, Plots, Printf
+
+h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
+outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "compare_gal_vs_el.png")
+ksel = length(ARGS) >= 3 ? ARGS[3] : "highest"
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
+    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+        to_c(read(f["integration/xi_psi"])), read(f["integration/psi"]),
+        to_c(read(f["galerkin/match/xi"])), read(f["galerkin/solution/psi"]),
+        Bool.(read(f["galerkin/solution/issing"])), read(f["info/mlow"]),
+        read(f["galerkin/sing_psi"]))
+end
+
+mpert = size(u1, 1)
+k = ksel == "highest" ? argmax(real.(et)) :
+    ksel == "lowest" ? argmin(real.(et)) :
+    parse(Int, ksel)
+w = wt[:, k]
+@printf("comparing eigenmode k=%d   δW = %.4e %+.4ei\n", k, real(et[k]), imag(et[k]))
+
+# EL profile
+cEL = u1[:, :, end] \ w
+xiE = reduce(hcat, (u1[:, :, ip] * cEL for ip in 1:size(u1, 3)))   # (mpert, nE)
+
+# gal ideal profile (identity-at-edge ⇒ coefficient = w); drop on-surface points
+keep = .!issing
+psiGk = psiG[keep]
+xiG = reduce(hcat, (gxi[:, ip, :] * w for ip in findall(keep)))    # (mpert, nGk)
+
+ms = mlow .+ (0:mpert-1)
+peak = [maximum(abs, @view xiE[i, :]) for i in 1:mpert]
+order = sortperm(peak; rev=true)
+ndom = min(5, mpert)
+
+# quantitative bulk agreement on the dominant harmonic (linear-interp EL onto gal grid, off rationals/edge)
+lininterp(xq, x, y) = (j = clamp(searchsortedlast(x, xq), 1, length(x) - 1);
+t = (xq - x[j]) / (x[j+1] - x[j]); y[j] * (1 - t) + y[j+1] * t)
+idom = order[1]
+inbulk(p) = 0.1 <= p <= 0.92 && all(abs(p - ps) > 8e-3 for ps in sing_psi)
+sel = [ip for ip in eachindex(psiGk) if inbulk(psiGk[ip])]
+num = sum(abs2(xiG[idom, ip] - lininterp(psiGk[ip], psiE, xiE[idom, :])) for ip in sel)
+den = sum(abs2(xiG[idom, ip]) for ip in sel)
+@printf("dominant m=%d:  bulk relative L2 ‖gal−EL‖/‖gal‖ = %.3e  (n=%d, off-layer/off-edge)\n",
+    ms[idom], sqrt(num / den), length(sel))
+
+plt = plot(; size=(1000, 650), xlabel="ψ_N", ylabel="|ξ^ψ_m(ψ)|",
+    title=@sprintf("Ideal gal vs EL total-energy eigenmode  (δW=%.2e)  — solid: EL, dashed: gal", real(et[k])),
+    legend=:topleft, left_margin=12Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+palette = theme_palette(:auto)
+for (c, i) in enumerate(order[1:ndom])
+    col = palette[mod1(c, length(palette))]
+    plot!(plt, psiE, abs.(@view xiE[i, :]); color=col, lw=2.2, label="EL m=$(ms[i])")
+    plot!(plt, psiGk, abs.(@view xiG[i, :]); color=col, lw=1.6, ls=:dash, label="gal m=$(ms[i])")
+end
+for (j, ps) in enumerate(sing_psi)
+    vline!(plt, [ps]; color=:black, ls=:dot, lw=1, label=(j == 1 ? "rational q" : ""))
+end
+
+display(plt)
+savefig(plt, outpng)
+println("dominant harmonics (m): ", [ms[i] for i in order[1:ndom]])
+println("saved: ", abspath(outpng))
diff --git a/benchmarks/compare_jbgradpsi_m2.jl b/benchmarks/compare_jbgradpsi_m2.jl
new file mode 100644
index 000000000..b8da06fcc
--- /dev/null
+++ b/benchmarks/compare_jbgradpsi_m2.jl
@@ -0,0 +1,61 @@
+# Compare the area-normalized b^ψ (perturbed_equilibrium/response/psi_area = b^ψ/⟨J·|∇ψ|⟩_θ) for one
+# poloidal harmonic between two GPEC runs that are identical except for which ξ feeds PerturbedEquilibrium:
+#   (1) IDEAL galerkin matched ξ   (gal_match_flag=true, gal_ideal_flag=true)
+#   (2) SHOOTING ξ                  (gal_match_flag=false)
+#
+# PE writes no ψ grid, so it's reconstructed: gal-ideal → galerkin/solution/psi minus issing points;
+# shooting → integration/psi.
+# Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m]
+
+using HDF5, Plots, Printf
+
+gal_h5 = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
+sh_h5 = length(ARGS) >= 2 ? ARGS[2] : "/tmp/shooting_test/gpec.h5"
+outpng = length(ARGS) >= 3 ? ARGS[3] : joinpath(@__DIR__, "cmp_jbgradpsi_m2.png")
+mtarget = length(ARGS) >= 4 ? parse(Int, ARGS[4]) : 2
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+# gal-ideal run: PE grid = gal solution grid with the on-surface (issing) points dropped
+pa_g, psi_g, mlow, sing_psi, sing_m = h5open(gal_h5) do f
+    pa = to_c(read(f["perturbed_equilibrium/response/psi_area"]))   # [npsi, mpert]
+    iss = Bool.(read(f["galerkin/solution/issing"]))
+    (pa, read(f["galerkin/solution/psi"])[.!iss], read(f["info/mlow"]),
+        read(f["galerkin/sing_psi"]), read(f["galerkin/sing_m"]))
+end
+# shooting run: PE grid = integration/psi
+pa_s, psi_s = h5open(sh_h5) do f
+    (to_c(read(f["perturbed_equilibrium/response/psi_area"])), read(f["integration/psi"]))
+end
+
+size(pa_g, 1) == length(psi_g) || error("gal grid mismatch: npsi=$(size(pa_g,1)) vs grid=$(length(psi_g))")
+size(pa_s, 1) == length(psi_s) || error("shoot grid mismatch: npsi=$(size(pa_s,1)) vs grid=$(length(psi_s))")
+
+col = mtarget - mlow + 1
+@printf("m=%d → column %d (mlow=%d, mpert=%d)\n", mtarget, col, mlow, size(pa_g, 2))
+g = @view pa_g[:, col]
+s = @view pa_s[:, col]
+
+# resonant surface for this m (q = m/n)
+psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
+# value at the grid point nearest the resonant surface, each run
+if !isnan(psi_res)
+    ig = argmin(abs.(psi_g .- psi_res))
+    is = argmin(abs.(psi_s .- psi_res))
+    @printf("m=%d resonant surface ψ=%.4f\n", mtarget, psi_res)
+    @printf("  nearest |b^ψ_area|:  gal-ideal=%.4e (ψ=%.4f)   shooting=%.4e (ψ=%.4f)   ratio=%.3f\n",
+        abs(g[ig]), psi_g[ig], abs(s[is]), psi_s[is], abs(g[ig]) / abs(s[is]))
+end
+
+plt = plot(; size=(1000, 620), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩|  (area-normalized)",
+    title="m=$mtarget perturbed normal field — gal-ideal vs shooting PE", legend=:topleft,
+    left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+plot!(plt, psi_g, abs.(g); lw=2.2, color=1, label="gal-ideal ξ → PE")
+plot!(plt, psi_s, abs.(s); lw=1.8, color=2, ls=:dash, label="shooting ξ → PE")
+for (j, ps) in enumerate(sing_psi)
+    ism2 = sing_m[j] == mtarget
+    vline!(plt, [ps]; color=(ism2 ? :red : :black), ls=:dot, lw=(ism2 ? 1.8 : 1),
+        label=(ism2 ? "q=$mtarget surface" : (j == 1 ? "rational q" : "")))
+end
+savefig(plt, outpng)
+println("saved: ", abspath(outpng))
diff --git a/benchmarks/gal_validation/gal_validation/gal_beta_scan.jl b/benchmarks/gal_validation/gal_validation/gal_beta_scan.jl
new file mode 100644
index 000000000..76ce10887
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/gal_beta_scan.jl
@@ -0,0 +1,111 @@
+#!/usr/bin/env julia
+# Gal-enabled TJ-analytic beta (pc) scan: at each pc, run the full GPEC pipeline with gal_flag=true
+# and capture three Δ′ measures per resonant surface (m=2 q=2, m=3 q=3):
+#   gal      = galerkin/pest3_Delta diagonal      (RDCON singular-Galerkin, this port)
+#   stride   = singular/delta_prime_matrix diag   (STRIDE BVP, pr178)
+#   cajump   = (ca_r-ca_l)/((2pi)^2 psio)          (ca-jump, pr178 primary)
+# Reuses examples/LAR_beta_scan/gpec.toml. Writes /tmp/gal_beta_scan_results.h5 + CSV.
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", "Projects", "GPEC_dev", "docs", "GPEC_julia"))
+const REPO = "/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia"
+Pkg.activate(REPO)
+using GeneralizedPerturbedEquilibrium
+using HDF5, TOML, Printf
+
+_warped_grid(x0, x1, N; p=2.0) = [x0 + (x1 - x0) * (1 - (1 - i / (N - 1))^p) for i in 0:N-1]
+const PC_FULL = _warped_grid(0.001, 0.1735, 40; p=2.0)
+const PC_TEST = [0.001, 0.05, 0.10, 0.15, 0.17]
+const BASE = TOML.parsefile(joinpath(REPO, "examples", "LAR_beta_scan", "gpec.toml"))
+
+function extract(h5)
+    h5open(h5, "r") do f
+        out = Dict{String,Any}()
+        out["q0"] = read(f, "equil/q0"); out["qmax"] = read(f, "equil/qmax")
+        out["dW_total"] = real(read(f, "vacuum/et")[1])
+        psio = read(f, "equil/psio"); denom = (2pi)^2 * psio
+        mlow = read(f, "info/mlow"); mpert = read(f, "info/mpert"); nlow = read(f, "info/nlow")
+        # surface m-values (STRIDE side)
+        m_s = vec(read(f, "singular/m")); n_s = vec(read(f, "singular/n")); msing = length(m_s)
+        # STRIDE delta_prime_matrix
+        dpm = haskey(f["singular"], "delta_prime_matrix") ? read(f, "singular/delta_prime_matrix") : nothing
+        cal = read(f, "singular/ca_left"); car = read(f, "singular/ca_right")
+        # gal
+        gm = haskey(f, "galerkin") && haskey(f["galerkin"], "sing_m") ? vec(read(f, "galerkin/sing_m")) : Int[]
+        gp = haskey(f, "galerkin") && haskey(f["galerkin"], "pest3_Delta") ? read(f, "galerkin/pest3_Delta") : nothing
+        for (s, m) in enumerate(m_s)
+            key = "m$(m)"
+            # ca-jump
+            ipr = 1 + m - mlow + (n_s[s] - nlow) * mpert
+            out["cajump_$key"] = (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom
+            # stride
+            out["stride_$key"] = dpm !== nothing && s <= size(dpm, 1) ? dpm[s, s] : NaN + NaN*im
+            # gal (match by m)
+            gi = findfirst(==(m), gm)
+            out["gal_$key"] = (gp !== nothing && gi !== nothing) ? gp[gi, gi] : NaN + NaN*im
+        end
+        return out
+    end
+end
+
+function run_point(pc)
+    dir = mktempdir(; prefix="gal_beta_")
+    try
+        c = deepcopy(BASE)
+        c["TJ_ANALYTIC_INPUT"]["pc"] = pc
+        c["ForceFreeStates"]["gal_flag"] = true
+        c["ForceFreeStates"]["gal_solver"] = "LU"
+        c["ForceFreeStates"]["HDF5_filename"] = joinpath(dir, "gpec.h5")
+        open(joinpath(dir, "gpec.toml"), "w") do io; TOML.print(io, c); end
+        GeneralizedPerturbedEquilibrium.main([dir])
+        return extract(joinpath(dir, "gpec.h5"))
+    catch e
+        @warn "pc=$pc failed" exception=(e, catch_backtrace())
+        return nothing
+    finally
+        rm(dir; force=true, recursive=true)
+    end
+end
+
+pcs = ("--test" in ARGS) ? PC_TEST : PC_FULL
+@info "Gal beta scan: $(length(pcs)) points (ε=0.2, qc=1.5, qa=3.6)"
+results = Tuple{Float64,Dict{String,Any}}[]
+t0 = time()
+for (i, pc) in enumerate(pcs)
+    ti = time()
+    r = run_point(pc)
+    el = time() - ti; tot = time() - t0
+    if r !== nothing
+        push!(results, (pc, r))
+        g2 = get(r, "gal_m2", NaN); s2 = get(r, "stride_m2", NaN)
+        g3 = get(r, "gal_m3", NaN); s3 = get(r, "stride_m3", NaN)
+        @printf("[%2d/%2d] pc=%.5f  q=2: gal=%+9.3f stride=%+9.3f | q=3: gal=%+9.3f stride=%+9.3f  dW=%+.4f  (%.0fs, tot %.0fs)\n",
+            i, length(pcs), pc, real(g2), real(s2), real(g3), real(s3), get(r, "dW_total", NaN), el, tot)
+        flush(stdout)
+    end
+end
+
+# Save
+out_h5 = "/tmp/gal_beta_scan_results.h5"
+isfile(out_h5) && rm(out_h5)
+h5open(out_h5, "w") do f
+    f["pc"] = [r[1] for r in results]
+    for key in ("gal_m2", "stride_m2", "cajump_m2", "gal_m3", "stride_m3", "cajump_m3", "dW_total")
+        vals = [get(r[2], key, NaN+NaN*im) for r in results]
+        if eltype(vals) <: Complex
+            f["$(key)_re"] = real.(vals); f["$(key)_im"] = imag.(vals)
+        else
+            f[key] = Float64.(vals)
+        end
+    end
+end
+open("/tmp/gal_beta_scan_results.csv", "w") do io
+    println(io, "pc,gal_m2,stride_m2,cajump_m2,gal_m3,stride_m3,cajump_m3,dW_total")
+    for (pc, r) in results
+        @printf(io, "%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f\n", pc,
+            real(get(r,"gal_m2",NaN)), real(get(r,"stride_m2",NaN)), real(get(r,"cajump_m2",NaN)),
+            real(get(r,"gal_m3",NaN)), real(get(r,"stride_m3",NaN)), real(get(r,"cajump_m3",NaN)),
+            get(r,"dW_total",NaN))
+    end
+end
+@info "Saved /tmp/gal_beta_scan_results.h5 and .csv ($(length(results)) points)"
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/gal_epsilon_scan.jl b/benchmarks/gal_validation/gal_validation/gal_epsilon_scan.jl
new file mode 100644
index 000000000..e701bee39
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/gal_epsilon_scan.jl
@@ -0,0 +1,112 @@
+#!/usr/bin/env julia
+# Gal-enabled TJ-analytic beta (pc) scan: at each pc, run the full GPEC pipeline with gal_flag=true
+# and capture three Δ′ measures per resonant surface (m=2 q=2, m=3 q=3):
+#   gal      = galerkin/pest3_Delta diagonal      (RDCON singular-Galerkin, this port)
+#   stride   = singular/delta_prime_matrix diag   (STRIDE BVP, pr178)
+#   cajump   = (ca_r-ca_l)/((2pi)^2 psio)          (ca-jump, pr178 primary)
+# Reuses examples/LAR_beta_scan/gpec.toml. Writes /tmp/gal_epsilon_scan_results.h5 + CSV.
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", "Projects", "GPEC_dev", "docs", "GPEC_julia"))
+const REPO = "/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia"
+Pkg.activate(REPO)
+using GeneralizedPerturbedEquilibrium
+using HDF5, TOML, Printf
+
+_warped_grid(x0, x1, N; p=2.0) = [x0 + (x1 - x0) * (1 - (1 - i / (N - 1))^p) for i in 0:N-1]
+const PC_FULL = _warped_grid(0.125, 0.660, 56; p=2.0)
+const PC_TEST = [0.2, 0.4, 0.55]
+const BASE = TOML.parsefile(joinpath(REPO, "examples", "LAR_epsilon_scan", "gpec.toml"))
+
+function extract(h5)
+    h5open(h5, "r") do f
+        out = Dict{String,Any}()
+        out["q0"] = read(f, "equil/q0"); out["qmax"] = read(f, "equil/qmax")
+        out["dW_total"] = real(read(f, "vacuum/et")[1])
+        psio = read(f, "equil/psio"); denom = (2pi)^2 * psio
+        mlow = read(f, "info/mlow"); mpert = read(f, "info/mpert"); nlow = read(f, "info/nlow")
+        # surface m-values (STRIDE side)
+        m_s = vec(read(f, "singular/m")); n_s = vec(read(f, "singular/n")); msing = length(m_s)
+        # STRIDE delta_prime_matrix
+        dpm = haskey(f["singular"], "delta_prime_matrix") ? read(f, "singular/delta_prime_matrix") : nothing
+        cal = read(f, "singular/ca_left"); car = read(f, "singular/ca_right")
+        # gal
+        gm = haskey(f, "galerkin") && haskey(f["galerkin"], "sing_m") ? vec(read(f, "galerkin/sing_m")) : Int[]
+        gp = haskey(f, "galerkin") && haskey(f["galerkin"], "pest3_Delta") ? read(f, "galerkin/pest3_Delta") : nothing
+        for (s, m) in enumerate(m_s)
+            key = "m$(m)"
+            # ca-jump
+            ipr = 1 + m - mlow + (n_s[s] - nlow) * mpert
+            out["cajump_$key"] = (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom
+            # stride
+            out["stride_$key"] = dpm !== nothing && s <= size(dpm, 1) ? dpm[s, s] : NaN + NaN*im
+            # gal (match by m)
+            gi = findfirst(==(m), gm)
+            out["gal_$key"] = (gp !== nothing && gi !== nothing) ? gp[gi, gi] : NaN + NaN*im
+        end
+        return out
+    end
+end
+
+function run_point(pc)
+    dir = mktempdir(; prefix="gal_eps_")
+    try
+        c = deepcopy(BASE)
+        c["TJ_ANALYTIC_INPUT"]["lar_r0"] = c["TJ_ANALYTIC_INPUT"]["lar_a"] / pc
+        c["Equilibrium"]["eq_type"] = "tj_analytic_direct"
+        c["ForceFreeStates"]["gal_flag"] = true
+        c["ForceFreeStates"]["gal_solver"] = "LU"
+        c["ForceFreeStates"]["HDF5_filename"] = joinpath(dir, "gpec.h5")
+        open(joinpath(dir, "gpec.toml"), "w") do io; TOML.print(io, c); end
+        GeneralizedPerturbedEquilibrium.main([dir])
+        return extract(joinpath(dir, "gpec.h5"))
+    catch e
+        @warn "pc=$pc failed" exception=(e, catch_backtrace())
+        return nothing
+    finally
+        rm(dir; force=true, recursive=true)
+    end
+end
+
+pcs = ("--test" in ARGS) ? PC_TEST : PC_FULL
+@info "Gal epsilon scan: $(length(pcs)) points (ε=0.2, qc=1.5, qa=3.6)"
+results = Tuple{Float64,Dict{String,Any}}[]
+t0 = time()
+for (i, pc) in enumerate(pcs)
+    ti = time()
+    r = run_point(pc)
+    el = time() - ti; tot = time() - t0
+    if r !== nothing
+        push!(results, (pc, r))
+        g2 = get(r, "gal_m2", NaN); s2 = get(r, "stride_m2", NaN)
+        g3 = get(r, "gal_m3", NaN); s3 = get(r, "stride_m3", NaN)
+        @printf("[%2d/%2d] pc=%.5f  q=2: gal=%+9.3f stride=%+9.3f | q=3: gal=%+9.3f stride=%+9.3f  dW=%+.4f  (%.0fs, tot %.0fs)\n",
+            i, length(pcs), pc, real(g2), real(s2), real(g3), real(s3), get(r, "dW_total", NaN), el, tot)
+        flush(stdout)
+    end
+end
+
+# Save
+out_h5 = "/tmp/gal_epsilon_scan_results.h5"
+isfile(out_h5) && rm(out_h5)
+h5open(out_h5, "w") do f
+    f["pc"] = [r[1] for r in results]
+    for key in ("gal_m2", "stride_m2", "cajump_m2", "gal_m3", "stride_m3", "cajump_m3", "dW_total")
+        vals = [get(r[2], key, NaN+NaN*im) for r in results]
+        if eltype(vals) <: Complex
+            f["$(key)_re"] = real.(vals); f["$(key)_im"] = imag.(vals)
+        else
+            f[key] = Float64.(vals)
+        end
+    end
+end
+open("/tmp/gal_epsilon_scan_results.csv", "w") do io
+    println(io, "pc,gal_m2,stride_m2,cajump_m2,gal_m3,stride_m3,cajump_m3,dW_total")
+    for (pc, r) in results
+        @printf(io, "%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f\n", pc,
+            real(get(r,"gal_m2",NaN)), real(get(r,"stride_m2",NaN)), real(get(r,"cajump_m2",NaN)),
+            real(get(r,"gal_m3",NaN)), real(get(r,"stride_m3",NaN)), real(get(r,"cajump_m3",NaN)),
+            get(r,"dW_total",NaN))
+    end
+end
+@info "Saved /tmp/gal_epsilon_scan_results.h5 and .csv ($(length(results)) points)"
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/gal_scan_plot.jl b/benchmarks/gal_validation/gal_validation/gal_scan_plot.jl
new file mode 100644
index 000000000..2e7cf3f15
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/gal_scan_plot.jl
@@ -0,0 +1,98 @@
+# Plot gal vs STRIDE Δ′ across the TJ-analytic beta and epsilon scans.
+# Reads /tmp/gal_beta_scan_results.csv and /tmp/gal_epsilon_scan_results.csv.
+using Pkg
+Pkg.activate("/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia")
+using Plots, DelimitedFiles, Printf
+gr()
+
+_median(v) = (w = sort(v); n = length(w); iseven(n) ? (w[n÷2] + w[n÷2+1]) / 2 : w[(n+1)÷2])
+
+function loadcsv(path)
+    isfile(path) || return nothing
+    data, hdr = readdlm(path, ','; header=true)
+    return Dict(strip(h) => Float64.(data[:, i]) for (i, h) in enumerate(vec(hdr)))
+end
+
+# Linear-extrapolate δW_total → 0 from the last k points → marginal-stability estimate.
+function marginal_point(cols, xname; k=6)
+    x = cols[xname]; dw = cols["dW_total"]
+    ord = sortperm(x); x = x[ord]; dw = dw[ord]
+    k = min(k, length(x))
+    xt = x[end-k+1:end]; yt = dw[end-k+1:end]
+    coef = hcat(xt, ones(k)) \ yt
+    return coef[1] == 0 ? NaN : -coef[2] / coef[1]
+end
+
+# Δ′ vs scan parameter; gal (solid+markers) vs STRIDE (dashed). Log y if all-positive.
+function panel(cols, xname, xlabel, title, xmarg)
+    x = cols[xname]
+    # log-y Δ′; non-positive values (q=3 crosses zero at low ε) are masked to NaN → shown as gaps.
+    poslog(v) = [y > 0 ? y : NaN for y in v]
+    p = plot(; xlabel=xlabel, ylabel="Δ′ (PEST-3)", title=title, legend=:topleft,
+        yscale=:log10,
+        left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=4Plots.mm)
+    for (lbl, col, c) in (("q=2", "m2", :dodgerblue), ("q=3", "m3", :crimson))
+        plot!(p, x, poslog(cols["gal_$col"]); lw=2, color=c, marker=:circle, ms=3, label="gal $lbl")
+        plot!(p, x, poslog(cols["stride_$col"]); lw=2, ls=:dash, color=c, label="STRIDE $lbl")
+    end
+    if isfinite(xmarg)
+        # extend x-range so the asymptote line (just past the last point) is visible
+        plot!(p; xlims=(minimum(x), max(maximum(x), xmarg) * 1.002))
+        vline!(p, [xmarg]; color=:black, ls=:dot, lw=2,
+            label=@sprintf("marginal (δW=0): %s≈%.4f", xname == "pc" ? "pₐ" : "ε", xmarg))
+    end
+    return p
+end
+
+# δW_total vs scan parameter, with a linear extrapolation of the last points to δW=0
+# (the ideal marginal-stability / kink-pole estimate).
+function dwpanel(cols, xname, xlabel)
+    x = cols[xname]; dw = cols["dW_total"]
+    ord = sortperm(x); x = x[ord]; dw = dw[ord]
+    p = plot(; xlabel=xlabel, ylabel="δW_total (Re et[1])", legend=:topright,
+        left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=4Plots.mm)
+    plot!(p, x, dw; lw=2, color=:black, marker=:circle, ms=3, label="δW_total")
+    hline!(p, [0.0]; color=:gray, ls=:dot, label="marginal (δW=0)")
+    x0 = marginal_point(cols, xname)
+    if isfinite(x0)
+        k = min(6, length(x)); xt = x[end-k+1:end]
+        coef = hcat(xt, ones(length(xt))) \ dw[end-k+1:end]
+        xext = range(xt[1], x0; length=20)
+        plot!(p, xext, coef[1] .* xext .+ coef[2]; lw=1.5, ls=:dash, color=:red,
+            label=@sprintf("extrap: δW=0 at %s≈%.4f", xname == "pc" ? "pₐ" : "ε", x0))
+        scatter!(p, [x0], [0.0]; color=:red, ms=6, marker=:star5, label="")
+        @printf("  marginal-point estimate (%s): δW=0 at %.5f\n", xname, x0)
+    end
+    return p
+end
+
+# Relative |gal-STRIDE|/|STRIDE| in %.
+function relpanel(cols, xname, xlabel)
+    x = cols[xname]
+    p = plot(; xlabel=xlabel, ylabel="|gal−STRIDE|/|STRIDE| (%)", yscale=:log10, legend=:topright,
+        left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=4Plots.mm)
+    for (lbl, col, c) in (("q=2", "m2", :dodgerblue), ("q=3", "m3", :crimson))
+        rel = 100 .* abs.(cols["gal_$col"] .- cols["stride_$col"]) ./ max.(abs.(cols["stride_$col"]), 1e-12)
+        rel = max.(rel, 1e-4)
+        plot!(p, x, rel; lw=2, color=c, marker=:circle, ms=3, label=lbl)
+    end
+    return p
+end
+
+for (name, path, xn, xl, ttl) in (
+        ("beta", "/tmp/gal_beta_scan_results.csv", "pc", "pₐ (on-axis pressure factor)", "TJ-analytic β scan (ε=0.2, qc=1.5, qa=3.6)"),
+        ("eps", "/tmp/gal_epsilon_scan_results.csv", "pc", "ε = a/R₀", "TJ-analytic ε scan (pc=0.001, qc=1.5, qa=3.6)"))
+    cols = loadcsv(path)
+    cols === nothing && continue
+    xmarg = marginal_point(cols, xn)
+    fig = plot(panel(cols, xn, xl, ttl, xmarg), dwpanel(cols, xn, xl), relpanel(cols, xn, xl); layout=(3, 1), size=(820, 1320))
+    out = "/tmp/gal_$(name)_comparison.png"
+    savefig(fig, out)
+    println("Saved figure: ", out)
+    println("=== $name: gal-vs-STRIDE relative diff (%) ===")
+    for col in ("m2", "m3")
+        rel = 100 .* abs.(cols["gal_$col"] .- cols["stride_$col"]) ./ max.(abs.(cols["stride_$col"]), 1e-12)
+        @printf("  %s: median=%.3f%%  max=%.2f%%  (n=%d)\n", col, _median(rel), maximum(rel), length(rel))
+    end
+end
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d.jl b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d.jl
new file mode 100644
index 000000000..a4a9a39f2
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d.jl
@@ -0,0 +1,53 @@
+# 2D scan: gal-vs-STRIDE q=3 tracking over (pc, pfac). STRIDE is pfac-independent
+# (recorded per run as a consistency check). Output /tmp/gal_pfac_scan2d.csv.
+using GeneralizedPerturbedEquilibrium
+using HDF5, TOML, Printf
+const REPO = "/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia"
+const BASE = TOML.parsefile(joinpath(REPO, "examples", "LAR_beta_scan", "gpec.toml"))
+const PCS = [0.005, 0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.070, 0.085, 0.100]
+const PFACS = [5e-4, 1e-3, 3e-3, 1e-2, 3e-2, 1e-1]
+
+function extract(h5)
+    h5open(h5, "r") do f
+        mlow = read(f, "info/mlow"); mpert = read(f, "info/mpert"); nlow = read(f, "info/nlow")
+        m_s = vec(read(f, "singular/m"))
+        dpm = read(f, "singular/delta_prime_matrix")
+        gm = vec(read(f, "galerkin/sing_m")); gp = read(f, "galerkin/pest3_Delta")
+        g3 = NaN; s3 = NaN; g2 = NaN; s2 = NaN
+        for (s, m) in enumerate(m_s)
+            gi = findfirst(==(m), gm)
+            if m == 2; g2 = real(gp[gi, gi]); s2 = real(dpm[s, s]); end
+            if m == 3; g3 = real(gp[gi, gi]); s3 = real(dpm[s, s]); end
+        end
+        return (g2, s2, g3, s3)
+    end
+end
+
+io = open("/tmp/gal_pfac_scan2d.csv", "w")
+println(io, "pc,pfac,gal_q2,stride_q2,gal_q3,stride_q3"); flush(io)
+t0 = time(); n = 0; ntot = length(PCS) * length(PFACS)
+for pc in PCS, pfac in PFACS
+    global n += 1
+    dir = mktempdir(; prefix="gal2d_")
+    c = deepcopy(BASE)
+    c["TJ_ANALYTIC_INPUT"]["pc"] = pc
+    c["ForceFreeStates"]["gal_flag"] = true
+    c["ForceFreeStates"]["gal_solver"] = "LU"
+    c["ForceFreeStates"]["gal_pfac"] = pfac
+    c["ForceFreeStates"]["HDF5_filename"] = joinpath(dir, "gpec.h5")
+    open(joinpath(dir, "gpec.toml"), "w") do f; TOML.print(f, c); end
+    try
+        GeneralizedPerturbedEquilibrium.main([dir])
+        g2, s2, g3, s3 = extract(joinpath(dir, "gpec.h5"))
+        @printf(io, "%.6f,%.6e,%.6f,%.6f,%.6f,%.6f\n", pc, pfac, g2, s2, g3, s3); flush(io)
+        @printf("[%3d/%3d] pc=%.4f pfac=%.0e  q3: gal=%+8.3f stride=%+8.3f  relerr=%5.1f%%  tot %.0fs\n",
+            n, ntot, pc, pfac, g3, s3, 100abs(g3 - s3) / abs(s3), time() - t0); flush(stdout)
+    catch e
+        @warn "pc=$pc pfac=$pfac failed" exception=(e, catch_backtrace())
+    finally
+        rm(dir; force=true, recursive=true)
+    end
+end
+close(io)
+println("Saved /tmp/gal_pfac_scan2d.csv")
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d_big_plot.jl b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d_big_plot.jl
new file mode 100644
index 000000000..1ec4b7f72
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d_big_plot.jl
@@ -0,0 +1,41 @@
+# Plot the fine (50 pc x 30 pfac) gal-vs-STRIDE q=3 tracking scan: heatmap + per-pfac curves.
+using Pkg; Pkg.activate("/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia")
+using Plots, DelimitedFiles, Printf; gr()
+
+data, hdr = readdlm("/tmp/gal_pfac_scan2d_big.csv", ','; header=true)
+col(name) = data[:, findfirst(==(name), strip.(vec(hdr)))]
+pc = col("pc"); pfac = col("pfac"); gq3 = col("gal_q3"); sq3 = col("stride_q3")
+relerr = 100 .* abs.(gq3 .- sq3) ./ max.(abs.(sq3), 1e-12)
+
+pcs = sort(unique(round.(pc; digits=6))); pfacs = sort(unique(pfac))
+Z = fill(NaN, length(pfacs), length(pcs))
+for k in eachindex(pc)
+    i = findfirst(≈(pfac[k]), pfacs); j = findfirst(≈(round(pc[k]; digits=6)), pcs)
+    (i === nothing || j === nothing) && continue
+    Z[i, j] = relerr[k]
+end
+Zc = clamp.(Z, 0.01, 100)
+
+p1 = heatmap(pcs, pfacs, log10.(Zc); yscale=:log10, xlabel="pₐ", ylabel="packing pfac",
+    title="gal q=3 tracking vs STRIDE: log₁₀ |gal−STRIDE|/STRIDE (%)   [50×30]",
+    colorbar_title="log₁₀(% err)", c=:viridis,
+    left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=10Plots.mm)
+
+# per-pfac error vs pc, subset of pfacs for legibility
+p2 = plot(; xlabel="pₐ", ylabel="|gal−STRIDE|/STRIDE (%)", yscale=:log10, legend=:topright,
+    left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=4Plots.mm,
+    title="q=3 tracking error per packing ratio")
+sel = [findmin(abs.(log10.(pfacs) .- t))[2] for t in (-4, -3, -2, -1.5, -1, 0)]
+for i in unique(sel)
+    plot!(p2, pcs, [isfinite(z) ? max(z, 1e-3) : NaN for z in Z[i, :]]; lw=2, label=@sprintf("pfac=%.1e", pfacs[i]))
+end
+fig = plot(p1, p2; layout=(2, 1), size=(1000, 1050))
+savefig(fig, "/tmp/gal_pfac_scan2d_big.png")
+println("Saved /tmp/gal_pfac_scan2d_big.png")
+
+med(v) = (w = sort(filter(isfinite, v)); isempty(w) ? NaN : (n = length(w); iseven(n) ? (w[n÷2]+w[n÷2+1])/2 : w[(n+1)÷2]))
+println("\npfac        median%   max%")
+for (i, pf) in enumerate(pfacs)
+    @printf("%.2e    %7.3f  %7.1f\n", pf, med(Z[i, :]), maximum(filter(isfinite, Z[i, :]); init=NaN))
+end
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d_plot.jl b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d_plot.jl
new file mode 100644
index 000000000..9f4f2e67b
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_scan2d_plot.jl
@@ -0,0 +1,45 @@
+# Plot the (pc x pfac) gal-vs-STRIDE q=3 tracking scan: heatmap + per-pfac curves.
+using Pkg; Pkg.activate("/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia")
+using Plots, DelimitedFiles, Printf; gr()
+
+data, hdr = readdlm("/tmp/gal_pfac_scan2d.csv", ','; header=true)
+col(name) = data[:, findfirst(==(name), strip.(vec(hdr)))]
+pc = col("pc"); pfac = col("pfac"); gq3 = col("gal_q3"); sq3 = col("stride_q3")
+relerr = 100 .* abs.(gq3 .- sq3) ./ abs.(sq3)
+
+pcs = sort(unique(pc)); pfacs = sort(unique(pfac))
+# Build relerr grid Z[i_pfac, j_pc]
+Z = fill(NaN, length(pfacs), length(pcs))
+for k in eachindex(pc)
+    i = findfirst(≈(pfac[k]), pfacs); j = findfirst(≈(pc[k]), pcs)
+    Z[i, j] = relerr[k]
+end
+Zc = clamp.(Z, 0.01, 100)  # for log color
+
+# Panel 1: heatmap, log color
+p1 = heatmap(pcs, pfacs, log10.(Zc);
+    yscale=:log10, xlabel="pₐ", ylabel="packing pfac",
+    title="gal q=3 tracking error vs STRIDE:  log₁₀ |gal−STRIDE|/STRIDE (%)",
+    colorbar_title="log₁₀(% err)", c=:viridis,
+    left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=8Plots.mm)
+
+# Panel 2: per-pfac tracking curves
+p2 = plot(; xlabel="pₐ", ylabel="|gal−STRIDE|/STRIDE (%)", yscale=:log10, legend=:topright,
+    left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=4Plots.mm,
+    title="q=3 tracking error per packing ratio")
+for (i, pf) in enumerate(pfacs)
+    plot!(p2, pcs, max.(Z[i, :], 1e-3); lw=2, marker=:circle, ms=3, label=@sprintf("pfac=%.0e", pf))
+end
+
+fig = plot(p1, p2; layout=(2, 1), size=(900, 1000))
+savefig(fig, "/tmp/gal_pfac_scan2d.png")
+println("Saved /tmp/gal_pfac_scan2d.png")
+
+# Per-pfac median/max summary
+println("\npfac      median%   max%   (q=3 over the pc grid)")
+med(v)=(w=sort(filter(isfinite,v));isempty(w) ? NaN : (n=length(w); iseven(n) ? (w[n÷2]+w[n÷2+1])/2 : w[(n+1)÷2]))
+for (i, pf) in enumerate(pfacs)
+    r = Z[i, :]
+    @printf("%.0e   %7.2f  %6.1f\n", pf, med(r), maximum(filter(isfinite, r)))
+end
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_worker.jl b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_worker.jl
new file mode 100644
index 000000000..84197210f
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/pfac_study/gal_pfac_worker.jl
@@ -0,0 +1,80 @@
+# Worker: for ONE pc, run the full pipeline once, then loop galerkin_solve over MANY pfacs.
+# Args: pc, out_csv_path. Appends rows "pc,pfac,gal_q2,stride_q2,gal_q3,stride_q3".
+using GeneralizedPerturbedEquilibrium
+const GPE = GeneralizedPerturbedEquilibrium
+const FFS = GPE.ForceFreeStates
+const EQ = GPE.Equilibrium
+using TOML, Printf, LinearAlgebra
+import FastInterpolations: cubic_interp, Series, ExtendExtrap
+LinearAlgebra.BLAS.set_num_threads(1)   # avoid oversubscription across parallel workers
+
+const REPO = "/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia"
+const PFACS = exp10.(range(-4, 0; length=30))   # 30 packing ratios, 1e-4 .. 1e0
+
+const PCLIST = parse.(Float64, split(ARGS[1], ','))   # comma-separated pc values
+const OUTCSV = ARGS[2]
+
+function run_pc(pc)
+inputs = TOML.parsefile(joinpath(REPO, "examples", "LAR_beta_scan", "gpec.toml"))
+inputs["TJ_ANALYTIC_INPUT"]["pc"] = pc
+ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
+ctrl.gal_flag = true; ctrl.gal_solver = "LU"; ctrl.verbose = false
+dir = mktempdir(; prefix="galw_"); intr = FFS.ForceFreeStatesInternal(; dir_path=dir)
+eq_config = EQ.EquilibriumConfig(inputs["Equilibrium"], dir)
+tj = EQ.TJAnalyticConfig(; (Symbol(k) => v for (k, v) in inputs["TJ_ANALYTIC_INPUT"])...)
+equil = EQ.setup_equilibrium(eq_config, tj)
+
+# --- replicate main() pipeline (mirrors GeneralizedPerturbedEquilibrium.jl) ---
+FFS.sing_lim!(intr, ctrl, equil)
+xs = equil.profiles.xs; ls = zeros(length(xs), 5); FFS.mercier_scan!(ls, equil)
+intr.locstab = cubic_interp(xs, Series(ls); extrap=ExtendExtrap())
+intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+FFS.sing_find!(intr, equil)
+intr.mlow = trunc(Int, min(intr.nlow * equil.params.qmin, 0)) - 4 - ctrl.delta_mlow
+intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+intr.mpert = intr.mhigh - intr.mlow + 1
+intr.mband = min(max(intr.mpert - 1 - ctrl.delta_mband, 0), intr.mpert - 1)
+intr.numpert_total = intr.mpert * intr.npert
+metric = FFS.make_metric(equil; mband=intr.mband, fft_flag=ctrl.fft_flag)
+ffit = FFS.make_matrix(equil, intr, metric)
+odet, fm_propagators, fm_chunks, fm_S_left = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr)
+vac_data = FFS.free_run!(odet, ctrl, equil, ffit, intr)
+if intr.msing > 0 && fm_propagators !== nothing
+    FFS.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks;
+        wv=vac_data.wv, psio=equil.psio, debug=false,
+        S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit)
+end
+dpm = intr.delta_prime_matrix
+m_s = [s.m[1] for s in intr.sing if s.psifac < intr.psilim]
+sidx2 = findfirst(==(2), m_s); sidx3 = findfirst(==(3), m_s)
+s2 = sidx2 === nothing ? NaN : real(dpm[sidx2, sidx2])
+s3 = sidx3 === nothing ? NaN : real(dpm[sidx3, sidx3])
+
+open(OUTCSV, "a") do io
+    for pfac in PFACS
+        ctrl.gal_pfac = pfac
+        g2 = NaN; g3 = NaN
+        try
+            res = FFS.galerkin_solve(ctrl, equil, ffit, intr; vac_data=vac_data)
+            gm = res.sing_m
+            gi2 = findfirst(==(2), gm); gi3 = findfirst(==(3), gm)
+            g2 = gi2 === nothing ? NaN : real(res.Deltap[gi2, gi2])
+            g3 = gi3 === nothing ? NaN : real(res.Deltap[gi3, gi3])
+        catch e
+            @warn "gal failed pc=$pc pfac=$pfac" exception=e
+        end
+        @printf(io, "%.6f,%.6e,%.6f,%.6f,%.6f,%.6f\n", pc, pfac, g2, s2, g3, s3); flush(io)
+    end
+end
+rm(dir; force=true, recursive=true)
+return nothing
+end
+
+for pc in PCLIST
+    try
+        run_pc(pc); println("done pc=$pc"); flush(stdout)
+    catch e
+        @warn "pc=$pc failed" exception=(e, catch_backtrace())
+    end
+end
+println("WORKER DONE")
diff --git a/benchmarks/gal_validation/gal_validation/tokamaker_beta/gal_tok_plot.jl b/benchmarks/gal_validation/gal_validation/tokamaker_beta/gal_tok_plot.jl
new file mode 100644
index 000000000..9e2d657c1
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/tokamaker_beta/gal_tok_plot.jl
@@ -0,0 +1,53 @@
+# Plot the TokaMaker beta-scan pfac convergence study. Data: /tmp/gal_tok_scan.csv
+# (beta,pfac,m,q,gal,stride). Produces: per-pfac error summary, a heatmap (beta x pfac)
+# of median-over-surfaces tracking error, and per-q tracking-vs-beta at the chosen pfac.
+using Pkg; Pkg.activate("/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia")
+using Plots, DelimitedFiles, Printf; gr()
+
+d, h = readdlm("/tmp/gal_tok_scan.csv", ','; header=true)
+col(n) = d[:, findfirst(==(n), strip.(vec(h)))]
+beta = col("beta"); pfac = col("pfac"); m = Int.(col("m")); q = col("q"); g = col("gal"); s = col("stride")
+# Exclude q=1: the internal-kink surface where the outer-region Δ′ construction does not apply
+# (gal and STRIDE diverge structurally, pfac-independent). Keep the tearing surfaces q≥2.
+keep = round.(q; digits=0) .>= 2
+beta = beta[keep]; pfac = pfac[keep]; m = m[keep]; q = q[keep]; g = g[keep]; s = s[keep]
+finite = isfinite.(g) .& isfinite.(s) .& (abs.(s) .> 1e-9)
+relerr = fill(NaN, length(g))
+relerr[finite] = 100 .* abs.(g[finite] .- s[finite]) ./ abs.(s[finite])
+
+betas = sort(unique(beta)); pfacs = sort(unique(pfac))
+med(v) = (w = sort(filter(isfinite, v)); isempty(w) ? NaN : (n = length(w); iseven(n) ? (w[n÷2]+w[n÷2+1])/2 : w[(n+1)÷2]))
+qq(v, p) = (w = sort(filter(isfinite, v)); isempty(w) ? NaN : w[clamp(ceil(Int, p*length(w)), 1, length(w))])
+
+# --- per-pfac summary over ALL (beta, surface) ---
+println("pfac        median%   p90%    finite-max%   n")
+for pf in pfacs
+    r = relerr[pfac .≈ pf]
+    @printf("%.2e    %7.3f  %7.2f  %9.1f   %d\n", pf, med(r), qq(r, 0.9), maximum(filter(isfinite, r); init=NaN), count(isfinite, r))
+end
+
+# --- heatmap: median over surfaces, per (beta, pfac) ---
+Z = fill(NaN, length(pfacs), length(betas))
+for (i, pf) in enumerate(pfacs), (j, b) in enumerate(betas)
+    Z[i, j] = med(relerr[(pfac .≈ pf) .& (beta .≈ b)])
+end
+Zc = clamp.(Z, 0.01, 100)
+p1 = heatmap(betas, pfacs, log10.(Zc); yscale=:log10, xlabel="β_N", ylabel="packing pfac",
+    title="TokaMaker β scan: gal tracking vs STRIDE (median over q≥2 surfaces)\nlog₁₀ |gal−STRIDE|/STRIDE (%)",
+    colorbar_title="log₁₀(% err)", c=:viridis,
+    left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=10Plots.mm)
+
+# --- per-pfac error curves vs beta (median over surfaces) ---
+p2 = plot(; xlabel="β_N", ylabel="median |gal−STRIDE|/STRIDE (%)", yscale=:log10, legend=:topright,
+    left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=4Plots.mm,
+    title="median tracking error vs β per packing ratio")
+sel = [findmin(abs.(log10.(pfacs) .- t))[2] for t in (-4, -3, -2, -1.5, -1, 0)]
+for i in unique(sel)
+    yv = [max(Z[i, j], 1e-3) for j in eachindex(betas)]
+    plot!(p2, betas, yv; lw=2, marker=:circle, ms=2, label=@sprintf("pfac=%.1e", pfacs[i]))
+end
+fig = plot(p1, p2; layout=(2, 1), size=(1000, 1050))
+savefig(fig, "/tmp/gal_tok_scan.png")
+println("Saved /tmp/gal_tok_scan.png")
+println("surfaces present (q): ", sort(unique(round.(q; digits=1))))
+println("DONE")
diff --git a/benchmarks/gal_validation/gal_validation/tokamaker_beta/gal_tok_worker.jl b/benchmarks/gal_validation/gal_validation/tokamaker_beta/gal_tok_worker.jl
new file mode 100644
index 000000000..01d75b276
--- /dev/null
+++ b/benchmarks/gal_validation/gal_validation/tokamaker_beta/gal_tok_worker.jl
@@ -0,0 +1,122 @@
+# Worker: for ONE TokaMaker geqdsk, run the full GPEC pipeline once, then loop galerkin_solve
+# over many pfacs. Records EVERY rational surface. Args: geqdsk_path, beta, out_csv.
+# Appends rows: beta,pfac,m,q,gal,stride
+using GeneralizedPerturbedEquilibrium
+const GPE = GeneralizedPerturbedEquilibrium
+const FFS = GPE.ForceFreeStates
+const EQ = GPE.Equilibrium
+using TOML, Printf, LinearAlgebra
+import FastInterpolations: cubic_interp, Series, ExtendExtrap
+LinearAlgebra.BLAS.set_num_threads(1)
+
+const PFACS = exp10.(range(-4, 0; length=30))
+const GEQDSK = ARGS[1]
+const BETA = parse(Float64, ARGS[2])
+const OUTCSV = ARGS[3]
+
+# DIIID-like stability settings (efit, free-boundary nowall), no PE/forcing.
+function base_inputs(dir)
+    Dict(
+        "Equilibrium" => Dict("eq_filename" => basename(GEQDSK), "eq_type" => "efit",
+            "jac_type" => "hamada", "power_bp" => 0, "power_b" => 0, "power_r" => 0,
+            "grid_type" => "log_asymptotic", "psilow" => 1e-4, "psihigh" => 0.999,
+            "mpsi" => 0, "psi_accuracy" => 0.001, "mtheta" => 256, "newq0" => 0,
+            "etol" => 1e-10, "force_termination" => false),
+        "Wall" => Dict("shape" => "nowall", "a" => 0.2415),
+        "ForceFreeStates" => Dict(
+            "bal_flag" => false, "mat_flag" => true, "ode_flag" => true,
+            "vac_flag" => true, "mer_flag" => true,
+            "qlow" => 1.02, "qhigh" => 1e3, "sing_start" => 0,
+            "nn_low" => 1, "nn_high" => 1, "delta_mlow" => 8, "delta_mhigh" => 8, "delta_mband" => 0,
+            "mthvac" => 512, "thmax0" => 1, "kinetic_factor" => 0.0,
+            "eulerlagrange_tolerance" => 1e-10, "singfac_min" => 1e-4, "ucrit" => 1e4,
+            "use_parallel" => true, "parallel_threads" => 1, "populate_dense_xi" => false,
+            "truncate_at_dW_peak" => false, "set_psilim_via_dmlim" => true, "dmlim" => 0.2,
+            "gal_flag" => true, "gal_solver" => "LU"))
+end
+
+# Run setup + the full stability/STRIDE pipeline at a given psihigh. High-β limited TokaMaker
+# equilibria need a smaller psihigh (the LCFS sits at the grid edge; tracing too far out fails),
+# so the caller retries down a ladder.
+function run_pipeline(psihigh)
+    dir = mktempdir(; prefix="galtok_")
+    cp(GEQDSK, joinpath(dir, basename(GEQDSK)))
+    inputs = base_inputs(dir)
+    inputs["Equilibrium"]["psihigh"] = psihigh
+    ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
+    ctrl.verbose = false
+    intr = FFS.ForceFreeStatesInternal(; dir_path=dir)
+    eq_config = EQ.EquilibriumConfig(inputs["Equilibrium"], dir)
+    equil = EQ.setup_equilibrium(eq_config, nothing)
+    FFS.sing_lim!(intr, ctrl, equil)
+    xs = equil.profiles.xs; ls = zeros(length(xs), 5); FFS.mercier_scan!(ls, equil)
+    intr.locstab = cubic_interp(xs, Series(ls); extrap=ExtendExtrap())
+    intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+    FFS.sing_find!(intr, equil)
+    # Replicate main()'s post-find surface filter: drop surfaces outside [qlow, qlim].
+    # Without this the worker keeps the q=1 internal-kink surface (and mis-aligns the STRIDE
+    # matrix indexing), which main() excludes — the source of the spurious β≈1.6 error cliff.
+    if intr.msing > 0
+        qmin_int = max(ctrl.qlow, equil.params.qmin)
+        keep = [j for j in 1:intr.msing if intr.sing[j].q >= qmin_int && intr.sing[j].psifac <= intr.psilim]
+        intr.sing = intr.sing[keep]
+        intr.msing = length(keep)
+    end
+    intr.mlow = trunc(Int, min(intr.nlow * equil.params.qmin, 0)) - 4 - ctrl.delta_mlow
+    intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+    intr.mpert = intr.mhigh - intr.mlow + 1
+    intr.mband = min(max(intr.mpert - 1 - ctrl.delta_mband, 0), intr.mpert - 1)
+    intr.numpert_total = intr.mpert * intr.npert
+    metric = FFS.make_metric(equil; mband=intr.mband, fft_flag=ctrl.fft_flag)
+    ffit = FFS.make_matrix(equil, intr, metric)
+    odet, fm_propagators, fm_chunks, fm_S_left = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr)
+    vac_data = FFS.free_run!(odet, ctrl, equil, ffit, intr)
+    if intr.msing > 0 && fm_propagators !== nothing
+        FFS.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks;
+            wv=vac_data.wv, psio=equil.psio, debug=false,
+            S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit)
+    end
+    return (ctrl=ctrl, equil=equil, ffit=ffit, intr=intr, vac_data=vac_data, dir=dir, psihigh=psihigh)
+end
+
+P = nothing
+for ph in (0.999, 0.99, 0.98, 0.97, 0.95)
+    try
+        global P = run_pipeline(ph); break
+    catch e
+        @warn "pipeline failed at psihigh=$ph for $(basename(GEQDSK)); retrying lower" exception=e
+    end
+end
+P === nothing && error("all psihigh retries failed for $(basename(GEQDSK))")
+ctrl = P.ctrl; equil = P.equil; ffit = P.ffit; intr = P.intr; vac_data = P.vac_data; dir = P.dir
+@printf("  using psihigh=%.3f for %s\n", P.psihigh, basename(GEQDSK))
+dpm = intr.delta_prime_matrix
+m_s = [s.m[1] for s in intr.sing if s.psifac < intr.psilim]
+q_s = [s.q for s in intr.sing if s.psifac < intr.psilim]
+# STRIDE delta_prime_matrix may have fewer rows than m_s (it drops surfaces beyond its own
+# active-domain cut). Match by position but guard the access so extra surfaces just get NaN.
+ndpm = (dpm === nothing || isempty(dpm)) ? 0 : size(dpm, 1)
+stride_diag = [i <= ndpm ? real(dpm[i, i]) : NaN for i in 1:length(m_s)]
+
+open(OUTCSV, "a") do io
+    for pfac in PFACS
+        ctrl.gal_pfac = pfac
+        gdiag = fill(NaN, length(m_s)); gm = Int[]
+        try
+            res = FFS.galerkin_solve(ctrl, equil, ffit, intr; vac_data=vac_data)
+            gm = res.sing_m
+            for (i, m) in enumerate(m_s)
+                gi = findfirst(==(m), gm)
+                gdiag[i] = gi === nothing ? NaN : real(res.Deltap[gi, gi])
+            end
+        catch e
+            @warn "gal failed beta=$BETA pfac=$pfac" exception=e
+        end
+        for i in 1:length(m_s)
+            @printf(io, "%.4f,%.6e,%d,%.4f,%.6f,%.6f\n", BETA, pfac, m_s[i], q_s[i], gdiag[i], stride_diag[i])
+        end
+        flush(io)
+    end
+end
+rm(dir; force=true, recursive=true)
+@printf("WORKER DONE beta=%.4f  msing=%d  surfaces q=%s\n", BETA, length(m_s), string(round.(q_s; digits=2)))
diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl
new file mode 100644
index 000000000..842e0cc5c
--- /dev/null
+++ b/benchmarks/plot_xi_eigenmode.jl
@@ -0,0 +1,61 @@
+# Plot the ξ(ψ) profile of the eigenmode with the highest total-energy eigenvalue.
+#
+# The total energy operator W = W_plasma + W_vacuum (free_run!, Free.jl); its eigenvectors are the
+# free-boundary edge displacement patterns and the eigenvalues are δW. This picks the eigenmode with the
+# largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix
+# (integration/xi_psi) onto that edge eigenvector:  c = U_edge \ w,  ξ(ψ) = U(ψ)·c.
+#
+# Requires a run with populate_dense_xi=true so integration/xi_psi is the dense axis-basis fundamental
+# matrix (not the Riccati S-matrices).
+# Usage: julia --project=. benchmarks/plot_xi_eigenmode.jl [path/to/gpec.h5] [out.png]
+
+using HDF5, LinearAlgebra, Plots, Printf
+
+h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_ideal_example/gpec.h5"
+outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "xi_eigenmode.png")
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f
+    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+        to_c(read(f["integration/xi_psi"])),
+        read(f["integration/psi"]),
+        read(f["info/mlow"]),
+        haskey(f, "galerkin/sing_psi") ? read(f["galerkin/sing_psi"]) : Float64[])
+end
+
+mpert, _, nstep = size(u1)
+k = argmax(real.(et))                  # eigenmode with the highest total-energy eigenvalue
+@printf("highest-eigenvalue mode: k=%d   δW = %.4e %+.4ei   (of %d modes)\n",
+    k, real(et[k]), imag(et[k]), length(et))
+
+# Reconstruct ξ(ψ) for the eigenmode whose edge displacement is the eigenvector w = wt[:,k].
+w = wt[:, k]
+U_edge = u1[:, :, nstep]               # edge fundamental matrix (psi inner→edge, last index = edge)
+c = U_edge \ w                         # U_edge·c = w
+xi = Array{ComplexF64}(undef, mpert, nstep)
+for ip in 1:nstep
+    @views mul!(xi[:, ip], u1[:, :, ip], c)
+end
+
+ms = mlow .+ (0:mpert-1)
+peak = [maximum(abs, @view xi[i, :]) for i in 1:mpert]
+order = sortperm(peak; rev=true)       # dominant harmonics first
+
+plt = plot(; size=(950, 620), xlabel="ψ_N", ylabel="|ξ^ψ_m(ψ)|",
+    title=@sprintf("Highest total-energy eigenmode   δW = %.3e", real(et[k])),
+    legend=:topleft, left_margin=12Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+for i in 1:mpert                        # all harmonics, faint
+    plot!(plt, psi, abs.(@view xi[i, :]); color=:gray85, lw=0.6, label="")
+end
+for i in order[1:min(8, mpert)]         # dominant harmonics, bold + labeled
+    plot!(plt, psi, abs.(@view xi[i, :]); lw=2, label="m=$(ms[i])")
+end
+for (j, ps) in enumerate(sing_psi)      # rational surfaces
+    vline!(plt, [ps]; color=:black, ls=:dash, lw=1, label=(j == 1 ? "rational q" : ""))
+end
+
+savefig(plt, outpng)
+println("dominant harmonics (m): ", [ms[i] for i in order[1:min(8, mpert)]])
+println("saved: ", abspath(outpng))
diff --git a/benchmarks/scan_resistivity_m2.jl b/benchmarks/scan_resistivity_m2.jl
new file mode 100644
index 000000000..5ed0931dc
--- /dev/null
+++ b/benchmarks/scan_resistivity_m2.jl
@@ -0,0 +1,57 @@
+# Plot the m=2 area-normalized b^ψ (perturbed_equilibrium/response/psi_area) across a resistivity scan
+# of the RESISTIVE gal matched PE runs (gal_match_flag=true, gal_ideal_flag=false), one curve per η.
+# Overlays the shooting (ideal, η→0) reference. The η-scan dirs are produced by the bash loop over
+# /tmp/etascan_<factor> (each a copy of the 0.993 config with gal_eta scaled).
+# Usage: julia --project=. benchmarks/scan_resistivity_m2.jl [out.png] [m]
+
+using HDF5, Plots, Printf, TOML
+
+outpng = length(ARGS) >= 1 ? ARGS[1] : joinpath(@__DIR__, "scan_resistivity_m2.png")
+mtarget = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 2
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+# read m=target area-normalized b^ψ on the run's PE grid
+function read_m2(h5; gal::Bool)
+    h5open(h5) do f
+        pa = to_c(read(f["perturbed_equilibrium/response/psi_area"]))   # [npsi, mpert]
+        col = mtarget - read(f["info/mlow"]) + 1
+        psi = gal ? read(f["galerkin/solution/psi"])[.!Bool.(read(f["galerkin/solution/issing"]))] :
+              read(f["integration/psi"])
+        (psi, pa[:, col])
+    end
+end
+
+# collect scan dirs that finished
+scandirs = filter(d -> isfile(joinpath(d, "gpec.h5")),
+    ["/tmp/etascan_0p01", "/tmp/etascan_0p1", "/tmp/etascan_1", "/tmp/etascan_10", "/tmp/etascan_100"])
+etas = [TOML.parsefile(joinpath(d, "gpec.toml"))["ForceFreeStates"]["gal_eta"][1] for d in scandirs]
+ord = sortperm(etas)
+scandirs, etas = scandirs[ord], etas[ord]
+eta_ref = 8e-8
+@printf("%d scan runs: η = %s  (×η_ref where η_ref=%.0e)\n", length(etas),
+    join((@sprintf("%.0e", e) for e in etas), ", "), eta_ref)
+
+# rational surface for m=target
+sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f
+    (read(f["galerkin/sing_psi"]), read(f["galerkin/sing_m"]))
+end
+psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
+
+cols = cgrad(:viridis, max(length(etas), 2); categorical=true)
+plt = plot(; size=(1000, 640), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩|  (area-normalized)",
+    title="m=$mtarget perturbed normal field — resistivity scan (resistive gal-matched PE)",
+    legend=:topleft, left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+for (i, (d, e)) in enumerate(zip(scandirs, etas))
+    psi, v = read_m2(joinpath(d, "gpec.h5"); gal=true)
+    plot!(plt, psi, abs.(v); color=cols[i], lw=2, label=@sprintf("η = %g×η_ref", round(e / eta_ref; sigdigits=2)))
+end
+# shooting (ideal, η→0) reference
+if isfile("/tmp/shooting_test/gpec.h5")
+    psis, vs = read_m2("/tmp/shooting_test/gpec.h5"; gal=false)
+    plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="shooting (ideal, η→0)")
+end
+isnan(psi_res) || vline!(plt, [psi_res]; color=:red, ls=:dot, lw=1.6, label="q=$mtarget surface")
+
+savefig(plt, outpng)
+println("saved: ", abspath(outpng))
diff --git a/benchmarks/scan_rotation_m2.jl b/benchmarks/scan_rotation_m2.jl
new file mode 100644
index 000000000..63cf2f04d
--- /dev/null
+++ b/benchmarks/scan_rotation_m2.jl
@@ -0,0 +1,51 @@
+# Plot the m=2 area-normalized b^ψ (perturbed_equilibrium/response/psi_area) across a ROTATION scan
+# of the resistive gal matched PE runs (gal_match_flag=true, gal_ideal_flag=false), fixed η=8e-8,
+# rotation f = 1,2,4,8,16 Hz (forced eigenvalue γ_s = 2πi·n·f). One curve per rotation; overlays the
+# shooting (ideal) reference. Scan dirs produced by the bash loop over /tmp/rotscan_<f>.
+# Usage: julia --project=. benchmarks/scan_rotation_m2.jl [out.png] [m]
+
+using HDF5, Plots, Printf, TOML
+
+outpng = length(ARGS) >= 1 ? ARGS[1] : joinpath(@__DIR__, "scan_rotation_m2.png")
+mtarget = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 2
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+function read_m2(h5; gal::Bool)
+    h5open(h5) do f
+        pa = to_c(read(f["perturbed_equilibrium/response/psi_area"]))
+        col = mtarget - read(f["info/mlow"]) + 1
+        psi = gal ? read(f["galerkin/solution/psi"])[.!Bool.(read(f["galerkin/solution/issing"]))] :
+              read(f["integration/psi"])
+        (psi, pa[:, col])
+    end
+end
+
+scandirs = filter(d -> isfile(joinpath(d, "gpec.h5")),
+    ["/tmp/rotscan_1", "/tmp/rotscan_2", "/tmp/rotscan_4", "/tmp/rotscan_8", "/tmp/rotscan_16"])
+rots = [TOML.parsefile(joinpath(d, "gpec.toml"))["ForceFreeStates"]["gal_rotation"][1] for d in scandirs]
+ord = sortperm(rots)
+scandirs, rots = scandirs[ord], rots[ord]
+@printf("%d scan runs: rotation f = %s Hz  (η fixed = 8e-8)\n", length(rots), join((@sprintf("%g", r) for r in rots), ", "))
+
+sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f
+    (read(f["galerkin/sing_psi"]), read(f["galerkin/sing_m"]))
+end
+psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
+
+cols = cgrad(:plasma, max(length(rots), 2); categorical=true)
+plt = plot(; size=(1000, 640), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩|  (area-normalized)",
+    title="m=$mtarget perturbed normal field — rotation scan (η=8e-8, resistive gal-matched PE)",
+    legend=:topleft, left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+for (i, (d, r)) in enumerate(zip(scandirs, rots))
+    psi, v = read_m2(joinpath(d, "gpec.h5"); gal=true)
+    plot!(plt, psi, abs.(v); color=cols[i], lw=2, label=@sprintf("f = %g Hz", r))
+end
+if isfile("/tmp/shooting_test/gpec.h5")
+    psis, vs = read_m2("/tmp/shooting_test/gpec.h5"; gal=false)
+    plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="shooting (ideal)")
+end
+isnan(psi_res) || vline!(plt, [psi_res]; color=:red, ls=:dot, lw=1.6, label="q=$mtarget surface")
+
+savefig(plt, outpng)
+println("saved: ", abspath(outpng))
diff --git a/benchmarks/verify_gal_ideal.jl b/benchmarks/verify_gal_ideal.jl
new file mode 100644
index 000000000..1c0d0ad41
--- /dev/null
+++ b/benchmarks/verify_gal_ideal.jl
@@ -0,0 +1,30 @@
+# Verify the gal_ideal_flag path: the matched solution is the bare ideal coil column (cout=0).
+#   1. cout / cin / deltar all zero (no resistive plasma combination, no inner layer)
+#   2. matched ξ_j == the gal coil column sols(:,:,2·msing+j)   (and ξ′ likewise)
+# Usage: julia --project=. benchmarks/verify_gal_ideal.jl [path/to/gpec.h5]
+using HDF5, Printf, LinearAlgebra
+
+h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+cout, deltar, mxi, mdxi, sols, sols_d, sing_psi = h5open(h5path) do f
+    (to_c(read(f["galerkin/match/cout"])), to_c(read(f["galerkin/match/deltar"])),
+        to_c(read(f["galerkin/match/xi"])), to_c(read(f["galerkin/match/xi_deriv"])),
+        to_c(read(f["galerkin/solution/xi"])), to_c(read(f["galerkin/solution/xi_deriv"])),
+        read(f["galerkin/sing_psi"]))
+end
+msing = length(sing_psi)
+mpert, ngrid, mcoil = size(mxi)
+
+@printf("[1] ‖cout‖ = %.2e, ‖deltar‖ = %.2e   %s\n", norm(cout), norm(deltar),
+    (norm(cout) == 0 && norm(deltar) == 0) ? "✓ no resistive combination (ideal)" : "✗")
+
+# matched ξ_j should equal the coil column sols(:,:,2msing+j)
+err = maximum(1:mcoil) do j
+    csol = 2msing + j
+    max(norm(mxi[:, :, j] - sols[:, :, csol]), norm(mdxi[:, :, j] - sols_d[:, :, csol]))
+end
+@printf("[2] max‖ξ_matched − coil column‖ = %.2e   %s\n", err,
+    err < 1e-12 ? "✓ matched ξ == bare ideal coil column" : "✗")
+
+println((norm(cout) == 0 && err < 1e-12) ? "\n✓ gal_ideal_flag path verified" : "\n✗ investigate")
diff --git a/benchmarks/verify_gal_match.jl b/benchmarks/verify_gal_match.jl
new file mode 100644
index 000000000..86f8cfffe
--- /dev/null
+++ b/benchmarks/verify_gal_match.jl
@@ -0,0 +1,49 @@
+# Piece 2 verification: RPEC outer↔inner matched solution (galerkin/match/*).
+#   1. linear-solve residual ‖mat·cof − rmat‖/‖rmat‖
+#   2. matched ξ / ξ′ finiteness
+#   3. edge column == identity basis: each coil drive j must give ξ_edge = e_j (the j-th harmonic),
+#      since plasma columns vanish at the edge and coil column j carries the unit edge source
+#   4. inner-layer Δ(Q) per surface (finite, nonzero); forced eigenvalues γ_s = 2πi·n·f
+# Usage: julia --project=. benchmarks/verify_gal_match.jl [path/to/gpec.h5]
+using HDF5, Printf, LinearAlgebra
+
+h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_example/gpec.h5"
+@info "Reading $h5path"
+
+xi, dxi, cout, cin, deltar, eig, resid, sing_psi = h5open(h5path) do f
+    (read(f["galerkin/match/xi"]), read(f["galerkin/match/xi_deriv"]),
+        read(f["galerkin/match/cout"]), read(f["galerkin/match/cin"]),
+        read(f["galerkin/match/deltar"]), read(f["galerkin/match/rpec_eig"]),
+        read(f["galerkin/match/residual"]), read(f["galerkin/sing_psi"]))
+end
+# HDF5 stores ComplexF64 as a compound (re,im); convert if needed
+to_c(a) = eltype(a) <: Complex ? a : map(x -> ComplexF64(x.re, x.im), a)
+xi = to_c(xi); dxi = to_c(dxi); cout = to_c(cout); cin = to_c(cin); deltar = to_c(deltar); eig = to_c(eig)
+
+mpert, ngrid, mcoil = size(xi)
+msing = size(deltar, 1)
+@printf("matched solution: mpert=%d  ngrid=%d  mcoil=%d   msing=%d\n", mpert, ngrid, mcoil, msing)
+
+@printf("\n[1] linear-solve residual  ‖mat·cof − rmat‖/‖rmat‖ = %.3e   %s\n",
+    resid, resid < 1e-8 ? "✓" : "✗ (large)")
+
+nbad = count(!isfinite, xi) + count(!isfinite, dxi)
+@printf("[2] finiteness: non-finite entries = %d   |ξ|max=%.3e  |ξ′|max=%.3e   %s\n",
+    nbad, maximum(abs, xi), maximum(abs, dxi), nbad == 0 ? "✓" : "✗")
+
+# [3] edge column == identity (last grid point = psihigh edge)
+edge = xi[:, ngrid, :]                       # (mpert mode, mcoil drive)
+id_err = norm(edge - Matrix{ComplexF64}(I, mpert, mcoil)) / sqrt(mpert)
+@printf("[3] edge basis: ‖ξ(edge) − I‖/√mpert = %.3e   %s\n",
+    id_err, id_err < 1e-8 ? "✓ identity-at-edge" : "✗")
+
+# [4] inner-layer Δ and forced eigenvalues
+@printf("\n[4] per-surface inner-layer Δ(Q) and forced eigenvalue γ = 2πi·n·f:\n")
+@printf("    %-4s %-10s  %-26s %-26s\n", "s", "ψ_s", "Δ₁", "Δ₂")
+for s in 1:msing
+    @printf("    %-4d %-10.5f  %+.4e%+.4ei  %+.4e%+.4ei   γ=%+.3f%+.3fi\n",
+        s, sing_psi[s], real(deltar[s, 1]), imag(deltar[s, 1]),
+        real(deltar[s, 2]), imag(deltar[s, 2]), real(eig[s]), imag(eig[s]))
+end
+allgood = resid < 1e-8 && nbad == 0 && id_err < 1e-8 && all(isfinite, deltar) && all(!iszero, deltar)
+println(allgood ? "\n✓ Piece 2 matching checks PASS" : "\n✗ some checks failed — investigate")
diff --git a/benchmarks/verify_gal_solution.jl b/benchmarks/verify_gal_solution.jl
new file mode 100644
index 000000000..e6c02b653
--- /dev/null
+++ b/benchmarks/verify_gal_solution.jl
@@ -0,0 +1,77 @@
+# Piece 1 verification: reconstructed gal ξ(ψ) and analytic ξ′(ψ).
+#   1. shapes / finiteness sanity of galerkin/solution arrays
+#   2. analytic ξ′ vs centered finite-difference of ξ — agree in the smooth interior, diverge at the
+#      packed edge (the spline-endpoint-derivative artifact we deliberately avoid)
+# Usage: julia --project=. verify_gal_solution.jl [path/to/gpec.h5]
+using HDF5, Printf, Statistics
+
+h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_example/gpec.h5"
+@info "Reading $h5path"
+
+psi, q, issing, xi, dxi, sing_psi = h5open(h5path) do f
+    (read(f["galerkin/solution/psi"]), read(f["galerkin/solution/q"]),
+        read(f["galerkin/solution/issing"]), read(f["galerkin/solution/xi"]),
+        read(f["galerkin/solution/xi_deriv"]), read(f["galerkin/sing_psi"]))
+end
+issing = Bool.(issing)
+mpert, ngrid, nsol = size(xi)
+@printf("grid: mpert=%d  ngrid=%d  nsol=%d   psi∈[%.4f, %.4f]\n", mpert, ngrid, nsol, psi[1], psi[end])
+@printf("singular surfaces at psi = %s\n", join((@sprintf("%.4f", p) for p in sing_psi), ", "))
+
+# --- finiteness (skip on-surface points, which are intentionally left zero) ---
+good = .!issing
+nbad = count(!isfinite, xi[:, good, :]) + count(!isfinite, dxi[:, good, :])
+@printf("non-finite entries (off-surface): %d   |xi|max=%.3e  |dxi|max=%.3e\n",
+    nbad, maximum(abs, xi[:, good, :]), maximum(abs, dxi[:, good, :]))
+
+# --- analytic ξ′ vs centered finite difference of ξ ---
+# For each column, centered diff at interior grid points (using off-surface neighbours), compared to the
+# analytic xi_deriv. Report median relative error binned by distance to the nearest singular surface and
+# to the packed edge.
+dist_sing(p) = isempty(sing_psi) ? Inf : minimum(abs.(p .- sing_psi))
+relerr = Float64[]
+dsurf = Float64[]
+dedge = Float64[]
+for isol in 1:nsol
+    for ip in 2:(ngrid-1)
+        (issing[ip] || issing[ip-1] || issing[ip+1]) && continue
+        h1 = psi[ip] - psi[ip-1]
+        h2 = psi[ip+1] - psi[ip]
+        (h1 <= 0 || h2 <= 0) && continue
+        for m in 1:mpert
+            # nonuniform centered difference
+            fd = (xi[m, ip+1, isol] - xi[m, ip-1, isol]) / (h2 + h1)
+            an = dxi[m, ip, isol]
+            scale = max(abs(an), abs(fd), 1e-30)
+            push!(relerr, abs(fd - an) / scale)
+            push!(dsurf, dist_sing(psi[ip]))
+            push!(dedge, psi[end] - psi[ip])
+        end
+    end
+end
+
+@printf("\nFinite-diff(ξ) vs analytic ξ′   (%d samples)\n", length(relerr))
+println("  binned by distance to nearest singular surface:")
+for (lo, hi) in [(0.0, 1e-3), (1e-3, 1e-2), (1e-2, 1e-1), (1e-1, Inf)]
+    sel = (dsurf .>= lo) .& (dsurf .< hi)
+    n = count(sel)
+    n == 0 && continue
+    @printf("    |ψ-ψ_s|∈[%-7.0e,%-7.0e)  n=%-7d  median relerr=%.2e  p90=%.2e\n",
+        lo, hi, n, median(relerr[sel]), quantile(relerr[sel], 0.9))
+end
+println("  binned by distance to packed edge ψ=psihigh:")
+for (lo, hi) in [(0.0, 1e-4), (1e-4, 1e-3), (1e-3, 1e-2), (1e-2, Inf)]
+    sel = (dedge .>= lo) .& (dedge .< hi)
+    n = count(sel)
+    n == 0 && continue
+    @printf("    edge dist∈[%-7.0e,%-7.0e)  n=%-7d  median relerr=%.2e  p90=%.2e\n",
+        lo, hi, n, median(relerr[sel]), quantile(relerr[sel], 0.9))
+end
+
+# Interior (away from both layer and edge) should agree well; near edge/layer the finite diff degrades.
+interior = (dsurf .> 1e-2) .& (dedge .> 1e-2)
+@printf("\nINTERIOR (>1e-2 from layer AND edge): median relerr = %.2e  (n=%d)\n",
+    median(relerr[interior]), count(interior))
+println(median(relerr[interior]) < 1e-3 ?
+        "  ✓ analytic ξ′ matches finite-diff in the smooth interior" :
+        "  ✗ interior mismatch — investigate")
diff --git a/docs/make.jl b/docs/make.jl
index fbbbb3725..74fb7e57d 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -28,6 +28,7 @@ makedocs(;
             "Vacuum" => "vacuum.md",
             "Equilibrium" => "equilibrium.md",
             "Stability Analysis" => "stability.md",
+            "Galerkin Solver" => "galerkin.md",
             "KineticForces" => "kinetic_forces.md",
             "Forcing Terms" => "forcing_terms.md",
             "Perturbed Equilibrium" => "perturbed_equilibrium.md",
diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
new file mode 100644
index 000000000..540bde2c3
--- /dev/null
+++ b/docs/src/galerkin.md
@@ -0,0 +1,34 @@
+# Galerkin Δ′ Solver (RDCON outer region)
+
+The Galerkin solver is the Julia port of RDCON's `gal.f` — Dewar's singular Galerkin method for
+the outer-region Δ′ computation.  Instead of integrating the Euler-Lagrange ODE radially, it
+discretizes the displacement on a packed radial grid of Hermite-cubic finite elements and solves
+a single global banded system.  Cells adjacent to each rational surface ("resonant" and
+"extension" cells) carry the large/small power-law asymptotic solutions analytically, so the
+singular behavior is built into the basis rather than resolved numerically.
+
+The solve produces the inter-surface Δ′ matrix and the PEST-3 matching blocks
+(``A'``, ``B'``, ``\Gamma'``, ``\Delta'``), written to the HDF5 output under the `galerkin/`
+group.  These are the outer-region inputs to resistive matched-asymptotic stability analysis
+[Glasser 2016, Phys. Plasmas **23**, 072505].
+
+The implementation lives in `src/ForceFreeStates/Galerkin/`:
+
+| File | Content |
+|------|---------|
+| `GalerkinStructs.jl` | Cell, interval, and workspace types mirroring the Fortran derived types |
+| `GalerkinGrid.jl` | Packed grid construction and local→global DOF mapping |
+| `GalerkinAssembly.jl` | Element-level assembly: Hermite basis, Gauss-Lobatto stiffness, resonant and extension cells, boundary conditions |
+| `GalerkinSolve.jl` | Top-level driver `galerkin_solve`, banded solve, Δ′ extraction, PEST-3 blocks, HDF5 output |
+
+## API Reference
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
+Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolve.jl"]
+```
+
+## See also
+
+- `docs/src/stability.md` — the ForceFreeStates module this solver belongs to (EL system, metric matrices, singular surfaces)
+- `docs/src/inner_layer.md` — inner-layer solutions to be matched against the outer-region Δ′
diff --git a/docs/src/stability.md b/docs/src/stability.md
index b294125a3..1d0d3671c 100644
--- a/docs/src/stability.md
+++ b/docs/src/stability.md
@@ -222,8 +222,12 @@ environment variable; it is not a runtime parameter.
 
 ## API Reference
 
+The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separately in
+`docs/src/galerkin.md`.
+
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
+Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Mercier.jl", "Resist.jl", "Bal.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"]
 ```
 
 ## Example usage
@@ -306,6 +310,7 @@ end
 
 ## See also
 
+- `docs/src/galerkin.md` — RDCON outer-region Galerkin Δ′ solver (part of this module)
 - `docs/src/equilibrium.md` — build the `PlasmaEquilibrium` object required by this module
 - `docs/src/vacuum.md` — vacuum response computed from the EL solution in `free_run!`
 - `docs/src/perturbed_equilibrium.md` — downstream singular coupling analysis using Δ'
diff --git a/examples/DIIID-like_gal_resistive_example/TkMkr_D3Dlike_Hmode.geqdsk b/examples/DIIID-like_gal_resistive_example/TkMkr_D3Dlike_Hmode.geqdsk
new file mode 120000
index 000000000..4200141b3
--- /dev/null
+++ b/examples/DIIID-like_gal_resistive_example/TkMkr_D3Dlike_Hmode.geqdsk
@@ -0,0 +1 @@
+../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk
\ No newline at end of file
diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml
new file mode 100644
index 000000000..b52110c22
--- /dev/null
+++ b/examples/DIIID-like_gal_resistive_example/gpec.toml
@@ -0,0 +1,88 @@
+[Equilibrium]
+eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
+eq_type = "efit"               # Type of the input 2D equilibrium file
+jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_arc)
+power_bp = 0                   # Poloidal field power exponent for Jacobian
+power_b = 0                    # Toroidal field power exponent for Jacobian
+power_r = 0                    # Major radius power exponent for Jacobian
+grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
+psilow = 1e-4                  # Lower limit of normalized flux coordinate
+psihigh = 0.993                # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here)
+mpsi = 128                     # Number of radial grid points (0 = auto-compute from psi_accuracy)
+psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+mtheta = 256                   # Number of poloidal grid points
+newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+etol = 1e-10                   # Error tolerance for equilibrium solver
+force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
+
+[Wall]
+shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a = 0.2415                     # Distance from plasma (conformal) or shape parameter
+aw = 0.05                      # Half-thickness parameter for Dee-shaped walls
+bw = 1.5                       # Elongation parameter for wall shapes
+cw = 0                         # Offset of wall center from major radius
+dw = 0.5                       # Triangularity parameter for wall shapes
+tw = 0.05                      # Sharpness of wall corners (try 0.05 as initial value)
+equal_arc_wall = true          # Equal arc length distribution of nodes on wall
+
+[ForceFreeStates]
+bal_flag = false               # Ideal MHD ballooning criterion for short wavelengths
+mat_flag = true                # Construct coefficient matrices for diagnostic purposes (required by gal_flag)
+ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
+vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
+mer_flag = true                # Evaluate the Mercier criterion
+
+psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
+qlow = 1.02                    # Integration initiated at q determined by min(q0, qlow)...
+qhigh = 1e3                    # Integration terminated at q limit determined by min(qa, qhigh)...
+sing_start = 0                 # Start integration at the sing_start'th rational from the axis (psilow)
+
+nn_low = 1                     # Smallest toroidal mode number to include
+nn_high = 1                    # Largest toroidal mode number to include
+delta_mlow = 8                 # Expands lower bound of Fourier harmonics
+delta_mhigh = 4                # Expands upper bound of Fourier harmonics
+delta_mband = 0                # Integration keeps only this wide a band of solutions along the m,m' diagonal (0 = full)
+mthvac = 512                   # Number of points used in splines over poloidal angle at plasma-vacuum interface.
+thmax0 = 1                     # Linear multiplier on the automatic choice of theta integration bounds
+
+kinetic_source = "fixed"       # Kinetic matrix source ("fixed" test matrices, or "calculated" from PENTRC)
+kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
+eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations
+save_interval = 3              # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
+singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+ucrit = 1e4                    # Maximum fraction of solutions allowed before re-normalized
+
+# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details)
+use_parallel        = true     # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix)
+parallel_threads    = 2        # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates)
+populate_dense_xi   = false    # Dense axis-basis ξ for the FFS HDF5 output. Not needed here: no [PerturbedEquilibrium] section, and the gal-matched path builds its own dense ξ. Set true only for a shooting-fed PE run.
+truncate_at_dW_peak = false    # Edge-dW scan stays diagnostic; integration domain set by qhigh / psihigh / dmlim
+set_psilim_via_dmlim = false   # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim)
+dmlim               = 0.2      # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true)
+
+# RDCON outer-region singular Galerkin Δ′ solver with rpec coil columns (matched to Fortran &GAL_INPUT).
+# This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
+gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve (requires mat_flag = true)
+gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
+gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
+gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
+gal_dx0                = 5e-4  # Resonant-element integration truncation distance from the rational (× 1/|n q'|)
+gal_dx1                = 1e-3  # Resonant-element size (× 1/|n q'|)
+gal_dx2                = 1e-3  # Extension-element size (× 1/|n q'|)
+gal_cutoff             = 10    # Number of elements carrying the large (driving) solution as a source term
+gal_tol                = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tolerance
+gal_gnstep             = 5000  # Maximum resonant-integration steps
+gal_dx1dx2_flag        = true  # Enable the special dx1/dx2 element sizing for resonant/extension cells
+gal_sing_order         = 6     # Base power-series order for the Galerkin singular asymptotics
+gal_sing_order_ceiling = true  # Auto-raise the series order by ceil(2·Re(α)) on high-Mercier-index surfaces
+gal_rpec_flag          = true  # Append mpert coil-response columns (unit edge sources) to the solve → delta_coil block
+
+# DRIVEN (RPEC) outer↔inner matching. Per-surface resistive inputs (core→edge), mirroring the Fortran
+# rmatch DIII-D example: eta=8e-8, massden=3.3e-7 (ni≈1e20 m^-3), rotation=1 Hz, Γ=5/3. 4 surfaces.
+gal_match_flag         = true  # Solve the coil-driven matched ξ(ψ) from the gal Δ′ + the inner-layer Δ(Q). Requires gal_rpec_flag = true.
+gal_ideal_flag         = false # true → build the IDEAL solution (skip the inner layer, bare coil columns; ≡ DCON). false → the resistive matching below.
+gal_eta                = [8e-8, 8e-8, 8e-8, 8e-8]          # Per-surface resistivity η (core→edge); Fortran rmatch `eta`
+gal_rho                = [3.3e-7, 3.3e-7, 3.3e-7, 3.3e-7]  # Per-surface mass density ρ [kg/m³] (core→edge); Fortran rmatch `massden`
+gal_rotation           = [1.0, 1.0, 1.0, 1.0]             # Per-surface rotation frequency f [Hz] (core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation`
+gal_gamma              = 1.6666666666666667               # Ratio of specific heats Γ for the resistive-layer coefficients (5/3)
diff --git a/examples/DIIID-like_gal_resistive_pe_example/TkMkr_D3Dlike_Hmode.geqdsk b/examples/DIIID-like_gal_resistive_pe_example/TkMkr_D3Dlike_Hmode.geqdsk
new file mode 120000
index 000000000..4200141b3
--- /dev/null
+++ b/examples/DIIID-like_gal_resistive_pe_example/TkMkr_D3Dlike_Hmode.geqdsk
@@ -0,0 +1 @@
+../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk
\ No newline at end of file
diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
new file mode 100644
index 000000000..b52110c22
--- /dev/null
+++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -0,0 +1,88 @@
+[Equilibrium]
+eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
+eq_type = "efit"               # Type of the input 2D equilibrium file
+jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_arc)
+power_bp = 0                   # Poloidal field power exponent for Jacobian
+power_b = 0                    # Toroidal field power exponent for Jacobian
+power_r = 0                    # Major radius power exponent for Jacobian
+grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
+psilow = 1e-4                  # Lower limit of normalized flux coordinate
+psihigh = 0.993                # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here)
+mpsi = 128                     # Number of radial grid points (0 = auto-compute from psi_accuracy)
+psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+mtheta = 256                   # Number of poloidal grid points
+newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+etol = 1e-10                   # Error tolerance for equilibrium solver
+force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
+
+[Wall]
+shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a = 0.2415                     # Distance from plasma (conformal) or shape parameter
+aw = 0.05                      # Half-thickness parameter for Dee-shaped walls
+bw = 1.5                       # Elongation parameter for wall shapes
+cw = 0                         # Offset of wall center from major radius
+dw = 0.5                       # Triangularity parameter for wall shapes
+tw = 0.05                      # Sharpness of wall corners (try 0.05 as initial value)
+equal_arc_wall = true          # Equal arc length distribution of nodes on wall
+
+[ForceFreeStates]
+bal_flag = false               # Ideal MHD ballooning criterion for short wavelengths
+mat_flag = true                # Construct coefficient matrices for diagnostic purposes (required by gal_flag)
+ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
+vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
+mer_flag = true                # Evaluate the Mercier criterion
+
+psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
+qlow = 1.02                    # Integration initiated at q determined by min(q0, qlow)...
+qhigh = 1e3                    # Integration terminated at q limit determined by min(qa, qhigh)...
+sing_start = 0                 # Start integration at the sing_start'th rational from the axis (psilow)
+
+nn_low = 1                     # Smallest toroidal mode number to include
+nn_high = 1                    # Largest toroidal mode number to include
+delta_mlow = 8                 # Expands lower bound of Fourier harmonics
+delta_mhigh = 4                # Expands upper bound of Fourier harmonics
+delta_mband = 0                # Integration keeps only this wide a band of solutions along the m,m' diagonal (0 = full)
+mthvac = 512                   # Number of points used in splines over poloidal angle at plasma-vacuum interface.
+thmax0 = 1                     # Linear multiplier on the automatic choice of theta integration bounds
+
+kinetic_source = "fixed"       # Kinetic matrix source ("fixed" test matrices, or "calculated" from PENTRC)
+kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
+eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations
+save_interval = 3              # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
+singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+ucrit = 1e4                    # Maximum fraction of solutions allowed before re-normalized
+
+# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details)
+use_parallel        = true     # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix)
+parallel_threads    = 2        # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates)
+populate_dense_xi   = false    # Dense axis-basis ξ for the FFS HDF5 output. Not needed here: no [PerturbedEquilibrium] section, and the gal-matched path builds its own dense ξ. Set true only for a shooting-fed PE run.
+truncate_at_dW_peak = false    # Edge-dW scan stays diagnostic; integration domain set by qhigh / psihigh / dmlim
+set_psilim_via_dmlim = false   # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim)
+dmlim               = 0.2      # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true)
+
+# RDCON outer-region singular Galerkin Δ′ solver with rpec coil columns (matched to Fortran &GAL_INPUT).
+# This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
+gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve (requires mat_flag = true)
+gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
+gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
+gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
+gal_dx0                = 5e-4  # Resonant-element integration truncation distance from the rational (× 1/|n q'|)
+gal_dx1                = 1e-3  # Resonant-element size (× 1/|n q'|)
+gal_dx2                = 1e-3  # Extension-element size (× 1/|n q'|)
+gal_cutoff             = 10    # Number of elements carrying the large (driving) solution as a source term
+gal_tol                = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tolerance
+gal_gnstep             = 5000  # Maximum resonant-integration steps
+gal_dx1dx2_flag        = true  # Enable the special dx1/dx2 element sizing for resonant/extension cells
+gal_sing_order         = 6     # Base power-series order for the Galerkin singular asymptotics
+gal_sing_order_ceiling = true  # Auto-raise the series order by ceil(2·Re(α)) on high-Mercier-index surfaces
+gal_rpec_flag          = true  # Append mpert coil-response columns (unit edge sources) to the solve → delta_coil block
+
+# DRIVEN (RPEC) outer↔inner matching. Per-surface resistive inputs (core→edge), mirroring the Fortran
+# rmatch DIII-D example: eta=8e-8, massden=3.3e-7 (ni≈1e20 m^-3), rotation=1 Hz, Γ=5/3. 4 surfaces.
+gal_match_flag         = true  # Solve the coil-driven matched ξ(ψ) from the gal Δ′ + the inner-layer Δ(Q). Requires gal_rpec_flag = true.
+gal_ideal_flag         = false # true → build the IDEAL solution (skip the inner layer, bare coil columns; ≡ DCON). false → the resistive matching below.
+gal_eta                = [8e-8, 8e-8, 8e-8, 8e-8]          # Per-surface resistivity η (core→edge); Fortran rmatch `eta`
+gal_rho                = [3.3e-7, 3.3e-7, 3.3e-7, 3.3e-7]  # Per-surface mass density ρ [kg/m³] (core→edge); Fortran rmatch `massden`
+gal_rotation           = [1.0, 1.0, 1.0, 1.0]             # Per-surface rotation frequency f [Hz] (core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation`
+gal_gamma              = 1.6666666666666667               # Ratio of specific heats Γ for the resistive-layer coefficients (5/3)
diff --git a/regression-harness/cases/gal_resistive_diiid.toml b/regression-harness/cases/gal_resistive_diiid.toml
new file mode 100644
index 000000000..9ec61109e
--- /dev/null
+++ b/regression-harness/cases/gal_resistive_diiid.toml
@@ -0,0 +1,107 @@
+[case]
+name = "gal_resistive_diiid"
+description = "DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with rpec coil columns (delta_coil)"
+example_dir = "examples/DIIID-like_gal_resistive_example"
+
+# Number of resonant surfaces in the gal domain
+[quantities.gal_msing]
+h5path = "galerkin/msing"
+type = "int_scalar"
+extract = "value"
+label = "gal # singular surfaces"
+noise_threshold = 0
+order = 10
+
+[quantities.gal_sing_q]
+h5path = "galerkin/sing_q"
+type = "real_vector"
+extract = "all_real"
+label = "gal singular q values"
+noise_threshold = 1e-8
+order = 11
+
+# PEST-3 Δ matching matrix — per-surface diagonal (the physics-meaningful tearing Δ′)
+[quantities.gal_pest3_delta_diag]
+h5path = "galerkin/pest3_Delta"
+type = "complex_matrix"
+extract = "diagonal_complex"
+label = "gal PEST3 Δ diagonal"
+noise_threshold = 1e-6
+order = 20
+
+# Full outer Δ′ matrix (nsol × 2·msing) — Frobenius norm catches any element drift
+[quantities.gal_delta_norm]
+h5path = "galerkin/delta"
+type = "complex_matrix"
+extract = "norm"
+label = "||gal Δ′ matrix||"
+noise_threshold = 1e-6
+order = 21
+
+# rpec coil-response block (mpert × 2·msing) — Frobenius norm
+[quantities.gal_delta_coil_norm]
+h5path = "galerkin/delta_coil"
+type = "complex_matrix"
+extract = "norm"
+label = "||gal Δ_coil block||"
+noise_threshold = 1e-6
+order = 22
+
+# Mercier index per surface
+[quantities.gal_di]
+h5path = "galerkin/di"
+type = "real_vector"
+extract = "all_real"
+label = "gal D_I per surface"
+noise_threshold = 1e-8
+order = 30
+
+# Resonant exponents α per surface
+[quantities.gal_alpha]
+h5path = "galerkin/alpha"
+type = "complex_vector"
+extract = "all_complex"
+label = "gal α per surface"
+noise_threshold = 1e-8
+order = 31
+
+# --- DRIVEN (RPEC) inner-layer matching (gal_match_flag) ---
+# Only finite, well-defined matching quantities are tracked. The reconstructed ξ/ξ′ are deliberately
+# NOT regressed via norms: ξ′ diverges at the rational surfaces, so any whole-array norm is dominated by
+# the (grid- and ULP-fragile) near-singular points — an ill-posed benchmark. The producer is validated
+# separately by benchmarks/verify_gal_solution.jl (finite-difference self-consistency).
+
+# Inner-layer matching data Δ(Q) per surface (resist_eval geometry + GGJ inner solver). (msing × 2)
+[quantities.gal_match_deltar_norm]
+h5path = "galerkin/match/deltar"
+type = "complex_matrix"
+extract = "norm"
+label = "||gal inner-layer Δ||"
+noise_threshold = 1e-6
+order = 41
+
+# Outer-region matched coefficients cout (2·msing × mcoil) — the 4·msing matching assembly + solve.
+[quantities.gal_match_cout_norm]
+h5path = "galerkin/match/cout"
+type = "complex_matrix"
+extract = "norm"
+label = "||gal match cout||"
+noise_threshold = 1e-6
+order = 42
+
+# Matching linear-solve residual ‖mat·cof − rmat‖/‖rmat‖ — health check (should stay ~machine eps).
+[quantities.gal_match_residual]
+h5path = "galerkin/match/residual"
+type = "real_scalar"
+extract = "value"
+label = "gal match residual"
+noise_threshold = 1e-10
+order = 44
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
diff --git a/regression-harness/cases/gal_resistive_pe.toml b/regression-harness/cases/gal_resistive_pe.toml
new file mode 100644
index 000000000..7a26d7563
--- /dev/null
+++ b/regression-harness/cases/gal_resistive_pe.toml
@@ -0,0 +1,87 @@
+[case]
+name = "gal_resistive_pe"
+description = "DIII-D-like, n=1, DRIVEN/RPEC: gal-matched resistive ξ → PerturbedEquilibrium (coil-driven singular coupling)"
+example_dir = "examples/DIIID-like_gal_resistive_pe_example"
+
+# Full-pipeline coverage for gal_match_flag → gal_matched_odestate → PerturbedEquilibrium. Tracks the
+# FINITE, per-rational-surface singular-coupling scalars (the physics-meaningful driven quantities) — not
+# any whole-profile ξ/ξ′ norm (those diverge at the rationals and are ill-posed as benchmarks).
+
+# Surface locations / count (sanity)
+[quantities.pe_rational_q]
+h5path = "perturbed_equilibrium/singular_coupling/rational_q"
+type = "real_vector"
+extract = "all_real"
+label = "rational q values"
+noise_threshold = 1e-8
+order = 10
+
+# Island half-widths per rational surface — the headline driven quantity
+[quantities.pe_island_half_width]
+h5path = "perturbed_equilibrium/singular_coupling/island_half_width"
+type = "real_vector"
+extract = "all_real"
+label = "island half-widths"
+noise_threshold = 1e-8
+order = 20
+
+# Chirikov overlap parameter per surface
+[quantities.pe_chirikov]
+h5path = "perturbed_equilibrium/singular_coupling/chirikov_parameter"
+type = "real_vector"
+extract = "all_real"
+label = "Chirikov parameter"
+noise_threshold = 1e-8
+order = 21
+
+# Resonant flux Φ_res per surface (applied) — norm + per-surface complex values
+[quantities.pe_resonant_flux_norm]
+h5path = "perturbed_equilibrium/singular_coupling/resonant_flux"
+type = "complex_vector"
+extract = "norm"
+label = "||resonant flux||"
+noise_threshold = 1e-6
+order = 30
+
+[quantities.pe_resonant_flux_all]
+h5path = "perturbed_equilibrium/singular_coupling/resonant_flux"
+type = "complex_vector"
+extract = "all_complex"
+label = "resonant flux Phi_res"
+noise_threshold = 1e-8
+order = 31
+
+# Penetrated (reconnected) resonant field per surface — the resistive shielding signature
+[quantities.pe_penetrated_field]
+h5path = "perturbed_equilibrium/singular_coupling/penetrated_field"
+type = "complex_vector"
+extract = "all_complex"
+label = "penetrated field"
+noise_threshold = 1e-8
+order = 32
+
+# Per-surface Δ' from the PE singular-coupling analysis
+[quantities.pe_delta_prime]
+h5path = "perturbed_equilibrium/singular_coupling/delta_prime"
+type = "complex_vector"
+extract = "all_complex"
+label = "PE Δ' per surface"
+noise_threshold = 1e-8
+order = 33
+
+# Resonant-flux coupling matrix (n_rational × numpert_total) — Frobenius norm, catches broader drift
+[quantities.pe_C_resonant_flux_norm]
+h5path = "perturbed_equilibrium/singular_coupling/C_resonant_flux"
+type = "complex_matrix"
+extract = "norm"
+label = "||C resonant flux||"
+noise_threshold = 1e-6
+order = 40
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 5e7821fe0..9201db5ff 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -200,11 +200,30 @@ function direct_position!(raw_profile::DirectRunInput)
     # ψ = 0 at the separatrix (after renormalization), and ψ changes sign between the
     # magnetic axis (ψ > 0) and the region outside the plasma (ψ < 0), so Brent is
     # globally convergent within the bracket (start_r, end_r) and needs no restarts.
+    # For limited equilibria whose grid box hugs the LCFS (e.g. some TokaMaker gEQDSKs),
+    # ψ can stay positive out to the grid edge with no sign change to bracket. The separatrix
+    # (ψ = 0) then sits just beyond the grid; extend the bracket outward via the extrapolating
+    # ψ spline until ψ changes sign (capped), so rs is a faithful ψ=0 crossing rather than the
+    # arbitrary grid edge (which lands in the degenerate ψ→0, Bp→0 region and breaks tracing).
+    psi_at(r) = (direct_get_bfield!(bfield, r, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio; derivs=0); bfield.psi)
     function find_separatrix_crossing(start_r, end_r, label)
-        r_sol = find_zero(
-            r -> (direct_get_bfield!(bfield, r, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio; derivs=0); bfield.psi),
-            (start_r, end_r), Roots.Brent()
-        )
+        hi = end_r
+        if psi_at(start_r) * psi_at(hi) > 0
+            span = abs(end_r - start_r)
+            step = 0.02 * span
+            hi_max = end_r + sign(end_r - start_r) * 0.5 * span   # cap extrapolation at 50% of the half-width
+            while psi_at(start_r) * psi_at(hi) > 0 && abs(hi - end_r) < abs(hi_max - end_r)
+                hi += sign(end_r - start_r) * step
+            end
+            if psi_at(start_r) * psi_at(hi) > 0
+                @warn "$label: ψ does not change sign even extrapolating to R = $(@sprintf("%.3f", hi)); " *
+                      "using R = $(@sprintf("%.3f", hi)) as the boundary."
+                return hi
+            end
+            @warn "$label: ψ positive to grid edge R = $(@sprintf("%.3f", end_r)) (limited equilibrium); " *
+                  "located ψ=0 by extrapolation in R ∈ (R_edge, $(@sprintf("%.3f", hi)))."
+        end
+        r_sol = find_zero(psi_at, (start_r, hi), Roots.Brent())
         @info "$label separatrix found at R = $(@sprintf("%.3f", r_sol))"
         return r_sol
     end
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
index 4af40c716..53677e2ad 100644
--- a/src/ForceFreeStates/ForceFreeStates.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -11,10 +11,13 @@ using JLD2
 using FastInterpolations
 using AdaptiveArrayPools
 using Roots
+using FastGaussQuadrature: gausslobatto
+using QuadGK: quadgk
 
 import ..Equilibrium
 import ..Utilities
 import ..Vacuum
+import ..InnerLayer
 using Printf
 using DoubleFloats
 import StaticArrays: @MMatrix
@@ -22,6 +25,7 @@ import StaticArrays: @MMatrix
 # Include all necessary files
 include("ForceFreeStatesStructs.jl")
 include("Mercier.jl")
+include("Resist.jl")
 include("Bal.jl")
 include("EulerLagrange.jl")
 include("Sing.jl")
@@ -34,6 +38,14 @@ include("PowerNorm.jl")
 include("Free.jl")
 include("Riccati.jl")
 
+# RDCON outer-region singular Galerkin Δ′ solver (gal_solve port)
+include("Galerkin/GalerkinStructs.jl")
+include("Galerkin/GalerkinGrid.jl")
+include("Galerkin/GalerkinAssembly.jl")
+include("Galerkin/GalerkinSolution.jl")
+include("Galerkin/GalerkinMatch.jl")
+include("Galerkin/GalerkinSolve.jl")
+
 # These are used for various small tolerances and root finders throughout ForceFreeStates
 global eps = 1e-10
 global itmax = 50
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index 852d8893e..69ce552df 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -183,6 +183,7 @@ A mutable struct holding internal state variables for stability calculations.
     kinsing_scan_cond::Vector{Float64} = Float64[]
     kinsing_scan_threshold::Float64 = 0.0
     psilim::Float64 = 0.0
+    psilow::Float64 = 0.0   # lower integration bound; raised above the axis by sing_min! when qlow > qmin (RDCON gal)
     qlim::Float64 = 0.0
     q1lim::Float64 = 0.0
     locstab::FastInterpolations.CubicSeriesInterpolant = cubic_interp(collect(0.0:0.25:1.0), Series(zeros(5, 5)); bc=ZeroCurvBC())
@@ -195,6 +196,9 @@ A mutable struct holding internal state variables for stability calculations.
     raw 2msing×2msing BVP solution to produce the PEST3-compatible tearing parameter.
     """
     delta_prime_matrix::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0)
+    delta_raw_matrix::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0)
+    delta_coil_matrix::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0)
+      # raw BVP Δ' matrix before PEST3 combination
 end
 
 """
@@ -296,6 +300,31 @@ A mutable struct containing control parameters for stability analysis, set by th
     use_parallel::Bool = true    # Default on: unlocks singular/delta_prime_matrix (STRIDE BVP Δ' matrix) used by SLAYER/GGJ downstream.
     populate_dense_xi::Bool = false  # When use_parallel=true, set to true ONLY if a PerturbedEquilibrium pipeline will consume dense ξ. Default false avoids the ~1× parallel-BVP serial-EL re-run for non-PE runs (Δ'/vacuum/ideal-stability only). See ForceFreeStatesControl docstring for the full trade-off (et[1] convention differs by ~0.12% on DIIID between populate=true vs false).
     extended_precision_bvp::Bool = true   # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases).
+
+    # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port) ---
+    gal_flag::Bool = false          # enable outer-region Galerkin Δ′ solve
+    gal_solver::String = "LU"       # "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs)
+    gal_nx::Int = 256               # elements per interval between singular surfaces
+    gal_nq::Int = 6                 # Gauss-Lobatto quadrature order per element
+    gal_pfac::Float64 = 0.001       # grid packing ratio near singular surfaces
+    gal_dx0::Float64 = 5e-4         # resonant-element integration truncation distance (×1/|n q'|)
+    gal_dx1::Float64 = 1e-3         # resonant-element size (×1/|n q'|)
+    gal_dx2::Float64 = 1e-3         # extension-element size (×1/|n q'|)
+    gal_cutoff::Int = 10            # # of elements carrying the large solution as driving term
+    gal_tol::Float64 = 1e-10        # resonant-quadrature (QuadGK) tolerance
+    gal_gnstep::Int = 20000         # max resonant-quadrature evaluations (QuadGK maxevals in gal_resonant!)
+    gal_dx1dx2_flag::Bool = true    # enable special dx1/dx2 treatment for resonant/extension elements
+    gal_sing_order::Int = 6         # base power-series order for the Galerkin singular asymptotics
+    gal_sing_order_ceiling::Bool = true  # auto-raise order by ceil(2·Re(α)) per surface (high Mercier index)
+    gal_rpec_flag::Bool = false     # append mpert coil-response columns to the Δ′ solve (RDCON rpec_flag): unit boundary sources whose plasma response is recorded; needed for the driven (resistive perturbed-equilibrium) Δ_gw
+    gal_edge_onesided::Bool = false # pack the two end intervals (`[psilow,first rational]`, `[last rational,psihigh]`) one-sided toward their single rational end instead of the Fortran-faithful symmetric "both" pack. Leaves the regular domain boundary at coarse spacing, removing the gratuitous fine edge cell that inflates cond(A) (entry spread ~1/h²) and seeds the downstream edge-derivative artifact. Default false = faithful to gal.f.
+    # --- DRIVEN (RPEC) outer↔inner asymptotic matching (rmatch match_rpec port) ---
+    gal_match_flag::Bool = false    # enable the RPEC inner-layer matching: solve the coil-driven matched ξ(ψ) from the gal Δ′ + the inner-layer Δ(Q). Requires gal_rpec_flag=true.
+    gal_ideal_flag::Bool = false    # within the match (gal_match_flag), build the IDEAL solution instead of the resistive one: skip the inner-layer Δ and use the bare coil columns (cout=0). Mirrors Fortran rmatch coil%ideal_flag; "physically equivalent to regular DCON" — the EL reference. eta/rho/rotation are ignored.
+    gal_eta::Vector{Float64} = Float64[]      # per-surface resistivity η (length msing, core→edge); Fortran rmatch `eta`
+    gal_rho::Vector{Float64} = Float64[]      # per-surface mass density ρ [kg/m³] (length msing, core→edge); Fortran rmatch `massden`
+    gal_rotation::Vector{Float64} = Float64[] # per-surface rotation frequency f [Hz] (length msing, core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation`
+    gal_gamma::Float64 = 5 / 3       # ratio of specific heats Γ for the resistive-layer coefficients (resist_eval G term)
     fixed_axis::Bool = false
 end
 
@@ -323,6 +352,7 @@ end
     hmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
     fmats_lower::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
     fmats_prim::S = _empty_series_interp_complex(numpert_total^2, itp_opts)  # primitive F before Schur complement (for kinetic)
+    fmats_gal::S = _empty_series_interp_complex(numpert_total^2, itp_opts)   # reduced Hermitian F̄ (un-factored) for the Galerkin solver
     kmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
     gmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
 
diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
index 45b07369a..1efb40690 100644
--- a/src/ForceFreeStates/Fourfit.jl
+++ b/src/ForceFreeStates/Fourfit.jl
@@ -365,6 +365,7 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
     hmats_flat = zeros(ComplexF64, mpsi, intr.numpert_total^2)
     fmats_lower_flat = zeros(ComplexF64, mpsi, intr.numpert_total^2)
     fmats_prim_flat = zeros(ComplexF64, mpsi, intr.numpert_total^2)  # primitive F before Schur complement (for kinetic)
+    fmats_gal_flat = zeros(ComplexF64, mpsi, intr.numpert_total^2)   # reduced Hermitian F̄ (before Cholesky) for the Galerkin solver
     gmats_flat = zeros(ComplexF64, mpsi, intr.numpert_total^2)
     kmats_flat = zeros(ComplexF64, mpsi, intr.numpert_total^2)
     g11 = Vector{ComplexF64}(undef, 2 * intr.mband + 1)
@@ -497,6 +498,11 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
         kmat .= emat .- (adjoint(kmat) * a_inv_cmat_temp)  # K̃ = E - K†A⁻¹C
         gmat .= hmat .- (adjoint(cmat) * a_inv_cmat_temp)  # G̃ = H - C†A⁻¹C
 
+        # Capture the reduced Hermitian F̄ (full matrix) before it is overwritten by its Cholesky
+        # factor below. The Galerkin solver (gal_get_fkg) needs F̄ directly: F = Q F̄ Qᴴ with
+        # Q = diag(singfac). Mirror of Fortran fmats_gal (fourfit.f).
+        @views fmats_gal_flat[ipsi, :] .= vec(fmat)
+
         # Store factorized F matrix (lower triangular only) since we always will need F⁻¹ later
         # and this make computation more efficient via combined forward and back substitution
         # TODO: does F stay Hermitian in the 3D case, allowing us to use the lower representation?
@@ -518,6 +524,7 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
     ffit.hmats = cubic_interp(metric.xs, Series(hmats_flat); ffit.itp_opts...)
     ffit.fmats_lower = cubic_interp(metric.xs, Series(fmats_lower_flat); ffit.itp_opts...)
     ffit.fmats_prim = cubic_interp(metric.xs, Series(fmats_prim_flat); ffit.itp_opts...)
+    ffit.fmats_gal = cubic_interp(metric.xs, Series(fmats_gal_flat); ffit.itp_opts...)
     ffit.gmats = cubic_interp(metric.xs, Series(gmats_flat); ffit.itp_opts...)
     ffit.kmats = cubic_interp(metric.xs, Series(kmats_flat); ffit.itp_opts...)
 
diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
new file mode 100644
index 000000000..94b6ddb35
--- /dev/null
+++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
@@ -0,0 +1,458 @@
+# GalerkinAssembly.jl
+#
+# Element-level assembly for the outer-region Galerkin solver. Ports:
+#   gal_hermite       (gal.f)   Hermite cubic basis
+#   gal_get_fkg       (gal.f)  F, K, G at a point from F̄/K̄/Ḡ + singfac
+#   gal_gauss_quad    (gal.f) nonresonant stiffness via Gauss-Lobatto quadrature
+#   gal_extension     (gal.f) extension-cell matrix/rhs (big-solution driving + surface terms)
+#   gal_resonant      (gal.f)   resonant-cell integral (QuadGK in place of Fortran LSODE)
+#   gal_assemble_mat  (gal.f)   scatter cell blocks into the banded global matrix (LU/Cholesky)
+#   gal_assemble_rhs  (gal.f)   scatter cell rhs into the global RHS
+#   gal_set_boundary  (gal.f) axis identity + edge (vacuum or fixed) BCs
+#
+# Hermite DOF index ip (Fortran 0:3) ↔ Julia array index ip+1 throughout.
+
+"""Resonant flat mode index `1 + (m-mlow) + (n-nlow)*mpert` for a singular surface (single resonance)."""
+@inline function _ipert_res(sing::SingType, intr::ForceFreeStatesInternal)
+    return 1 + (sing.m[1] - intr.mlow) + (sing.n[1] - intr.nlow) * intr.mpert
+end
+
+# Two-sided asymptotic evaluation (z = ψ - ψ_s). Right (z≥0) uses the sig=+1 series at z; left (z<0)
+# uses the sig=-1 series at the positive distance |z| (real √, no renorm), matching Fortran sing.f.
+# The derivative on the left picks up d|z|/dψ = -1.
+@inline sing_get_ua_gal(g::GalSingAsymp, z::Float64) =
+    z >= 0 ? sing_get_ua(g.right, z) : sing_get_ua(g.left, -z)
+@inline sing_get_dua_gal(g::GalSingAsymp, z::Float64) =
+    z >= 0 ? sing_get_dua(g.right, z) : (-1) .* sing_get_dua(g.left, -z)
+
+"""
+    gal_hermite(x, x0, x1) -> (pb, qb)
+
+Hermite cubic basis values `pb` and derivatives `qb` on `[x0, x1]` (length 4, Fortran `0:3`).
+Port of `gal_hermite` (gal.f).
+"""
+function gal_hermite(x::Real, x0::Real, x1::Real)
+    dx = x1 - x0
+    t0 = (x - x0) / dx
+    t1 = 1 - t0
+    t02 = t0 * t0
+    t12 = t1 * t1
+    pb = (t12 * (1 + 2t0), t12 * t0 * dx, t02 * (1 + 2t1), -t02 * t1 * dx)
+    qb = (-6t0 * t1 / dx, t0 * (3t0 - 4) + 1, 6t1 * t0 / dx, t0 * (3t0 - 2))
+    return pb, qb
+end
+
+"""
+    gal_get_fkg(ffit, intr, x, q) -> (F, K, G)
+
+Evaluate the `mpert×mpert` matrices `F = Q F̄ Qᴴ`, `K = Q K̄`, `G = Ḡ` at flux `x` with safety factor
+`q`, where `Q = diag(singfac)` and `singfac = m - n q` (direct). Port of `gal_get_fkg` (gal.f).
+Uses the un-factored reduced `ffit.fmats_gal` (F̄), `ffit.kmats` (K̄), `ffit.gmats` (Ḡ). F/K/G are the
+ideal-MHD Euler–Lagrange coefficient matrices of the outer-region weak form (Glasser 2016, PoP 23, 112506).
+"""
+function gal_get_fkg(ffit::FourFitVars, intr::ForceFreeStatesInternal, x::Float64, q::Float64)
+    N = intr.numpert_total
+    sf = vec((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)')
+
+    Fbar = Matrix{ComplexF64}(undef, N, N)
+    K = Matrix{ComplexF64}(undef, N, N)
+    G = Matrix{ComplexF64}(undef, N, N)
+    ffit.fmats_gal(vec(Fbar), x; hint=ffit._hint)
+    ffit.kmats(vec(K), x; hint=ffit._hint)
+    ffit.gmats(vec(G), x; hint=ffit._hint)
+
+    F = Matrix{ComplexF64}(undef, N, N)
+    @inbounds for j in 1:N, i in 1:N
+        F[i, j] = Fbar[i, j] * sf[i] * sf[j]
+        K[i, j] = K[i, j] * sf[i]
+    end
+    return F, K, G
+end
+
+"""
+    gal_gauss_quad!(cell, ffit, profiles, intr, nodes, weights, swap_edge)
+
+Accumulate the nonresonant Hermite stiffness block `cell.mat` by Gauss-Lobatto quadrature. Port of
+`gal_gauss_quad` (gal.f). When `swap_edge` (the final cell of the last interval), the
+right-node value/slope DOFs (Fortran pb(2)↔pb(3)) are swapped so the edge-value DOF lands at index 3,
+matching the free-boundary BC (gal.f).
+"""
+function gal_gauss_quad!(cell::GalCell, ffit::FourFitVars, profiles, intr::ForceFreeStatesInternal,
+    nodes::Vector{Float64}, weights::Vector{Float64}, swap_edge::Bool)
+
+    N = intr.numpert_total
+    np = GAL_NP
+    fill!(cell.mat, 0)
+    x1, x2 = cell.x
+    x0c = (x2 + x1) / 2
+    dxc = (x2 - x1) / 2
+    qhint = Ref(1)
+
+    @inbounds for iq in eachindex(nodes)
+        x = x0c + dxc * nodes[iq]
+        w = dxc * weights[iq]
+        q = profiles.q_spline(x; hint=qhint)
+        F, K, G = gal_get_fkg(ffit, intr, x, q)
+        pbt, qbt = gal_hermite(x, x1, x2)
+        pb = collect(pbt)
+        qb = collect(qbt)
+        if swap_edge
+            pb[3], pb[4] = pb[4], pb[3]
+            qb[3], qb[4] = qb[4], qb[3]
+        end
+        for ip in 0:np, ipert in 1:N, jp in 0:np, jpert in 1:N
+            cell.mat[ipert, jpert, ip+1, jp+1] += w * (
+                F[ipert, jpert] * qb[ip+1] * qb[jp+1] +
+                K[ipert, jpert] * qb[ip+1] * pb[jp+1] +
+                conj(K[jpert, ipert]) * pb[ip+1] * qb[jp+1] +
+                G[ipert, jpert] * pb[ip+1] * pb[jp+1])
+        end
+    end
+    return cell
+end
+
+# Resolve (jsing, ψ_s, sign) for an extension/resonant cell from its side.
+@inline function _cell_jsing(cell::GalCell, ising::Int)
+    return cell.extra == GAL_SIDE_LEFT ? ising + 1 : ising
+end
+
+"""
+    gal_extension!(cell, ising, ffit, profiles, intr, asymps, sings, nn, nodes, weights)
+
+Build `cell.emat`, `cell.ediag`, `cell.rhs`, `cell.erhs` for an extension cell (`ext`/`ext1`/`ext2`).
+Port of `gal_extension` (gal.f). The big solution (`isol` column) drives the RHS; for `ext`
+cells the small solution (`isol+N` column) also forms the resonant coupling `emat`/`ediag`.
+"""
+function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles,
+    intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType},
+    nn::Int, nodes::Vector{Float64}, weights::Vector{Float64})
+
+    N = intr.numpert_total
+    np = GAL_NP
+    x1, x2 = cell.x
+    qhint = Ref(1)
+
+    jsing = _cell_jsing(cell, ising)
+    asymp = asymps[jsing]
+    psi_s = sings[jsing].psifac
+    isol = _ipert_res(sings[jsing], intr)
+    ua_at(x) = sing_get_ua_gal(asymp, x - psi_s)
+    dua_at(x) = sing_get_dua_gal(asymp, x - psi_s)
+
+    # Boundary node and weights (gal.f)
+    if cell.extra == GAL_SIDE_LEFT
+        jp = 2          # Fortran jp (0-based) -> Julia jp+1
+        xb = x2
+        ws = -1.0
+    else
+        jp = 0
+        xb = x1
+        ws = 1.0
+    end
+    if cell.etype == GCT_EXT || cell.etype == GCT_EXT1
+        w1 = 1.0
+        w2 = 1.0
+    else # GCT_EXT2
+        w1 = (cell.extra == GAL_SIDE_LEFT) ? 0.0 : 1.0
+        w2 = (cell.extra == GAL_SIDE_LEFT) ? 1.0 : 0.0
+    end
+
+    fill!(cell.emat, 0)
+    cell.ediag = zero(ComplexF64)
+
+    # --- emat and ediag (ext cells only) (gal.f) ---
+    if cell.etype == GCT_EXT
+        uab = ua_at(xb)
+        duab = dua_at(xb)
+        u = uab[:, isol+N, 1]      # small solution
+        du = duab[:, isol+N, 1]
+        for ipert in 1:N
+            @views cell.emat .+= cell.mat[:, ipert, :, jp+1] .* u[ipert] .+
+                                 cell.mat[:, ipert, :, jp+2] .* du[ipert]
+        end
+        for ipert in 1:N
+            cell.ediag += conj(u[ipert]) * cell.emat[ipert, jp+1] +
+                          conj(du[ipert]) * cell.emat[ipert, jp+2]
+        end
+        # surface term at xb
+        qb = profiles.q_spline(xb; hint=qhint)
+        Fb, Kb, _ = gal_get_fkg(ffit, intr, xb, qb)
+        pbt, _ = gal_hermite(xb, x1, x2)
+        Fdu_Ku = Fb * du .+ Kb * u
+        for ip in 0:np
+            @views cell.emat[:, ip+1] .+= Fdu_Ku .* (pbt[ip+1] * ws)
+        end
+        cell.ediag += sum(conj.(u) .* Fdu_Ku) * ws
+    end
+
+    # --- rhs (gal.f) ---
+    ua1 = ua_at(x1)[:, isol, 1]    # big solution at left node
+    du1 = dua_at(x1)[:, isol, 1]
+    ua2 = ua_at(x2)[:, isol, 1]    # big solution at right node
+    du2 = dua_at(x2)[:, isol, 1]
+    fill!(cell.rhs, 0)
+
+    if cell.etype == GCT_EXT || cell.etype == GCT_EXT1
+        cell.erhs = zero(ComplexF64)
+        x0c = (x2 + x1) / 2
+        dxc = (x2 - x1) / 2
+        for iq in eachindex(nodes)
+            x = x0c + dxc * nodes[iq]
+            w = dxc * weights[iq]
+            q = profiles.q_spline(x; hint=qhint)
+            F, K, G = gal_get_fkg(ffit, intr, x, q)
+            pbt, qbt = gal_hermite(x, x1, x2)
+            uax = ua_at(x)
+            ub = uax[:, isol, 1]
+            dub = dua_at(x)[:, isol, 1]
+            term1 = F * dub .+ K * ub
+            term2 = adjoint(K) * dub .+ G * ub
+            for ip in 0:np
+                @views cell.rhs[:, ip+1] .-= (term1 .* qbt[ip+1] .+ term2 .* pbt[ip+1]) .* w
+            end
+        end
+        if cell.etype == GCT_EXT
+            for ipert in 1:N
+                cell.erhs += conj(ua_at(xb)[ipert, isol+N, 1]) * cell.rhs[ipert, jp+1] +
+                             conj(dua_at(xb)[ipert, isol+N, 1]) * cell.rhs[ipert, jp+2]
+            end
+        end
+    elseif cell.etype == GCT_EXT2
+        for ipert in 1:N
+            @views cell.rhs .-= w1 .* (cell.mat[:, ipert, :, 0+1] .* ua1[ipert] .+
+                                       cell.mat[:, ipert, :, 1+1] .* du1[ipert])
+            @views cell.rhs .-= w2 .* (cell.mat[:, ipert, :, 2+1] .* ua2[ipert] .+
+                                       cell.mat[:, ipert, :, 3+1] .* du2[ipert])
+        end
+    end
+
+    # --- surface terms (always; gal.f) ---
+    q_l = profiles.q_spline(x1; hint=qhint)
+    Fl, Kl, _ = gal_get_fkg(ffit, intr, x1, q_l)
+    pbt, _ = gal_hermite(x1, x1, x2)
+    surf_l = Fl * du1 .+ Kl * ua1
+    for ip in 0:np
+        @views cell.rhs[:, ip+1] .-= surf_l .* (pbt[ip+1] * w1)
+    end
+    if cell.etype == GCT_EXT && cell.extra == GAL_SIDE_RIGHT
+        cell.erhs -= sum(conj.(ua_at(x1)[:, isol+N, 1]) .* surf_l)
+    end
+
+    q_r = profiles.q_spline(x2; hint=qhint)
+    Fr, Kr, _ = gal_get_fkg(ffit, intr, x2, q_r)
+    pbt, _ = gal_hermite(x2, x1, x2)
+    surf_r = Fr * du2 .+ Kr * ua2
+    for ip in 0:np
+        @views cell.rhs[:, ip+1] .+= surf_r .* (pbt[ip+1] * w2)
+    end
+    if cell.etype == GCT_EXT && cell.extra == GAL_SIDE_LEFT
+        cell.erhs += sum(conj.(ua_at(x2)[:, isol+N, 1]) .* surf_r)
+    end
+    return cell
+end
+
+"""
+    gal_resonant!(cell, ising, ffit, profiles, intr, asymps, sings, nn, gal_tol, gal_gnstep, verbose)
+
+Compute the resonant-cell contributions `cell.erhs`, `cell.ediag`, `cell.rhs`, `cell.emat` by adaptive
+quadrature (QuadGK) of the residual-operator integrand. Port of `gal_lsode_int`/`gal_lsode_der`
+(gal.f): the Fortran LSODE accumulation is replaced by `quadgk` over `[x_bdy, x_lsode]`.
+"""
+function gal_resonant!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles,
+    intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType},
+    nn::Int, gal_tol::Float64, gal_gnstep::Int, verbose::Bool)
+
+    N = intr.numpert_total
+    np = GAL_NP
+    x1, x2 = cell.x
+    qhint = Ref(1)
+
+    jsing = _cell_jsing(cell, ising)
+    asymp = asymps[jsing]
+    psi_s = sings[jsing].psifac
+    isol_big = _ipert_res(sings[jsing], intr)
+    isol_small = isol_big + N
+    cols = [isol_big, isol_small]
+
+    # Integration limits (gal.f): right → [x2, x_lsode]; left → [x1, x_lsode]
+    x0 = (cell.extra == GAL_SIDE_LEFT) ? x1 : x2
+    x1l = cell.x_lsode
+    sgn_u = (cell.extra == GAL_SIDE_RIGHT) ? -1.0 : 1.0
+
+    L = 2 + 2 * N * (np + 1)
+    function integrand(x)
+        z = x - psi_s
+        ua = sing_get_ua_gal(asymp, z)
+        dua = sing_get_dua_gal(asymp, z)
+        q = profiles.q_spline(x; hint=qhint)
+        ua2 = ua[:, cols, :]
+        dua2 = dua[:, cols, :]
+        mv = sing_matvec(ffit, intr, x, q, ua2, dua2)   # N×2 (col1=big, col2=small)
+        pbt, _ = gal_hermite(x, x1, x2)
+        w = conj.(@view ua2[:, 2, 1])                   # conj small solution, qty1
+        out = Vector{ComplexF64}(undef, L)
+        out[1] = sum(w .* @view mv[:, 1])               # res1
+        out[2] = sum(w .* @view mv[:, 2])               # res2
+        idx = 3
+        for ip in 0:np, ipert in 1:N
+            out[idx] = pbt[ip+1] * mv[ipert, 1]         # hbig (ipert fastest, then ip)
+            idx += 1
+        end
+        for ip in 0:np, ipert in 1:N
+            out[idx] = pbt[ip+1] * mv[ipert, 2]         # hsmall
+            idx += 1
+        end
+        return out
+    end
+
+    # Numerical-method deviation from gal_lsode_int: adaptive QuadGK (Gauss-Kronrod) replaces Fortran's
+    # LSODE accumulation over the same integrand/limits, so this path is NOT bit-exact with gal.f. Match
+    # Fortran's LSODE tolerance (gal_tol) and cap the evaluations at gal_gnstep (Fortran's gnstep limit)
+    # so a near-singular integrand cannot hang.
+    raw, qerr = quadgk(integrand, x0, x1l; rtol=gal_tol, atol=1e-30, maxevals=gal_gnstep)
+    if verbose
+        @info "  resonant jsing=$jsing side=$(cell.extra) qerr=$(qerr) res1=$(raw[1]) res2=$(raw[2])"
+    end
+    hbig = reshape(@view(raw[3:2+N*(np+1)]), N, np + 1)
+    hsmall = reshape(@view(raw[3+N*(np+1):end]), N, np + 1)
+
+    # Apply the right-side negation and the gal_make_arrays signs (gal.f gal_make_arrays)
+    cell.erhs = -sgn_u * raw[1]
+    cell.ediag = sgn_u * raw[2]
+    cell.rhs .= (-sgn_u) .* hbig
+    cell.emat .= sgn_u .* hsmall
+    return cell
+end
+
+"""
+    gal_assemble_mat!(ws, mpert)
+
+Scatter every cell's `mat`/`emat`/`ediag` into the banded global matrix `ws.mat`. Port of
+`gal_assemble_mat` (gal.f). Supports both `"LU"` (full band, `offset = kl+ku+1`) and
+`"cholesky"` (lower band only, `offset = 1`).
+"""
+function gal_assemble_mat!(ws::GalWorkspace, mpert::Int)
+    np = ws.np
+    msing = length(ws.intvl) - 1
+    lu = ws.solver == "LU"
+    offset = lu ? ws.kl + ws.ku + 1 : 1
+
+    for ising in 0:msing
+        for cell in ws.intvl[ising+1].cells
+            res = cell.etype == GCT_RES || cell.etype == GCT_EXT
+            for ip in 0:np, ipert in 1:mpert
+                i = cell.map[ipert, ip+1]
+                # nonresonant stiffness
+                for jp in 0:np, jpert in 1:mpert
+                    j = cell.map[jpert, jp+1]
+                    if lu
+                        ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1]
+                    elseif j <= i
+                        ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1]
+                    end
+                end
+                # resonant coupling (res/ext only)
+                if res
+                    j = cell.emap
+                    e = cell.emat[ipert, ip+1]
+                    if lu
+                        ws.mat[offset+i-j, j] += e
+                        ws.mat[offset+j-i, i] += conj(e)
+                    elseif j < i
+                        ws.mat[offset+i-j, j] += e
+                    elseif j > i
+                        ws.mat[offset+j-i, i] += conj(e)
+                    end
+                end
+            end
+            if res
+                ws.mat[offset, cell.emap] += cell.ediag
+            end
+        end
+    end
+    return ws
+end
+
+"""
+    gal_assemble_rhs!(ws, mpert)
+
+Scatter every cell's `rhs`/`erhs` into the global RHS `ws.rhs`. Port of `gal_assemble_rhs`
+(gal.f). A new RHS column (`isol`) begins at each big-solution driver: an `ext2`-left cell or a
+`res`-right cell.
+"""
+function gal_assemble_rhs!(ws::GalWorkspace, mpert::Int)
+    np = ws.np
+    msing = length(ws.intvl) - 1
+    isol = 0
+    for ising in 0:msing
+        for cell in ws.intvl[ising+1].cells
+            cell.etype == GCT_NONE && continue
+            if cell.etype == GCT_EXT2 && cell.extra == GAL_SIDE_LEFT
+                isol += 1
+            elseif cell.etype == GCT_RES && cell.extra == GAL_SIDE_RIGHT
+                isol += 1
+            end
+            for ip in 0:np, ipert in 1:mpert
+                imap = cell.map[ipert, ip+1]
+                ws.rhs[imap, isol] += cell.rhs[ipert, ip+1]
+            end
+            if cell.etype == GCT_EXT || cell.etype == GCT_RES
+                ws.rhs[cell.emap, isol] += cell.erhs
+            end
+        end
+    end
+    return ws
+end
+
+"""
+    gal_set_boundary!(ws, mpert, wv_edge)
+
+Apply the axis (`u(0)=0` identity) and edge boundary conditions to the boundary cells' `mat`, plus the
+rpec coil-source RHS injection, before global assembly. Port of `gal_set_boundary` (gal.f),
+including its full three-way edge branch (rpec count `ncoil = ws.nsol - 2*msing`):
+  - rpec (`ncoil > 0`): identity edge **and** a unit source at the edge value DOF of coil column
+    `2*msing+ipert` for each poloidal mode `ipert` (gal.f). `ws.rhs` is already allocated and
+    `gal_assemble_rhs!` only fills columns `1:2*msing`, so the injection is not clobbered.
+  - free boundary (`wv_edge !== nothing`): add the vacuum block `wvac·psio²` at the edge value DOF
+    (gal.f).
+  - fixed boundary (else): identity edge (gal.f).
+`wv_edge` is `nothing` for the rpec and fixed cases.
+"""
+function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}})
+    chol = ws.solver == "cholesky"
+    msing = length(ws.intvl) - 1
+    ncoil = ws.nsol - 2 * msing
+
+    # axis u(0): first cell of interval 0, Hermite value DOF ip=0 (Julia index 1)
+    cell = ws.intvl[1].cells[1]
+    cell.mat[:, :, 1, :] .= 0
+    chol && (cell.mat[:, :, :, 1] .= 0)
+    for idx in 1:mpert
+        cell.mat[idx, idx, 1, 1] = 1
+    end
+
+    # edge u(1): last cell of last interval, edge-value DOF ip=3 (Julia index 4, see gauss_quad swap)
+    cell = ws.intvl[msing+1].cells[ws.nx]
+    if ncoil > 0
+        # rpec: identity edge + unit coil sources (gal.f)
+        cell.mat[:, :, 4, :] .= 0
+        chol && (cell.mat[:, :, :, 4] .= 0)
+        for idx in 1:mpert
+            cell.mat[idx, idx, 4, 4] = 1
+        end
+        for ipert in 1:mpert
+            ws.rhs[cell.map[ipert, 4], 2 * msing + ipert] = 1
+        end
+    elseif wv_edge !== nothing
+        # free boundary: vacuum block wvac·psio² (gal.f)
+        @views cell.mat[:, :, 4, 4] .+= wv_edge
+    else
+        # fixed boundary: identity edge (gal.f)
+        cell.mat[:, :, 4, :] .= 0
+        chol && (cell.mat[:, :, :, 4] .= 0)
+        for idx in 1:mpert
+            cell.mat[idx, idx, 4, 4] = 1
+        end
+    end
+    return ws
+end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinGrid.jl b/src/ForceFreeStates/Galerkin/GalerkinGrid.jl
new file mode 100644
index 000000000..88bb26782
--- /dev/null
+++ b/src/ForceFreeStates/Galerkin/GalerkinGrid.jl
@@ -0,0 +1,232 @@
+# GalerkinGrid.jl
+#
+# Grid construction and local→global DOF mapping for the outer-region Galerkin solver.
+# Ports gal_pack (gal.f), gal_make_grid (gal.f), and gal_make_map (gal.f).
+#
+# Node arrays use Fortran 0-based `x(0:nx)` mapped to Julia `xn[i+1]` (so Fortran `x(i)` is `xn[i+1]`).
+# Cells are 1-based in both. Intervals `0:msing` (Fortran) are `ws.intvl[ising+1]` in Julia.
+# The optional `gal_xmin_flag`/`sing1_xmin` adaptive-width branch is intentionally not ported.
+
+"""
+    gal_pack(nx::Int, pfac::Float64, side::String) -> Vector{Float64}
+
+Packed logical grid on `[-1, 1]`. Port of Fortran `gal_pack` (gal.f). For `side="both"`
+returns `2*nx+1` nodes (Fortran `-nx:nx`) clustered toward both ends when `pfac < 1` (strong packing
+near the bracketing singular surfaces). `pfac > 1` packs toward the centre; `pfac == 1` is uniform.
+"""
+function gal_pack(nx::Int, pfac::Float64, side::String)
+    if side == "left"
+        xi = collect(-2nx:0) ./ (2nx)
+    elseif side == "right"
+        xi = collect(0:2nx) ./ (2nx)
+    elseif side == "both"
+        xi = collect(-nx:nx) ./ nx
+    else
+        error("gal_pack: unrecognized side = $side")
+    end
+
+    if pfac > 1
+        λ = sqrt(1 - 1 / pfac)
+        x = @. log((1 + λ * xi) / (1 - λ * xi)) / log((1 + λ) / (1 - λ))
+    elseif pfac == 1
+        x = copy(xi)
+    elseif pfac > 0
+        λ = sqrt(1 - pfac)
+        a = log((1 + λ) / (1 - λ))
+        x = @. (exp(a * xi) - 1) / (exp(a * xi) + 1) / λ
+    else
+        x = @. xi^(-pfac)
+    end
+
+    if side == "left"
+        @. x = 2x + 1
+    elseif side == "right"
+        @. x = 2x - 1
+    end
+    return x
+end
+
+"""
+    gal_make_grid!(intvl::GalInterval, ising::Int, msing::Int, sings::Vector{SingType},
+                   nn::Int, psilow::Float64, psihigh::Float64, ctrl::ForceFreeStatesControl)
+
+Set the cell types (`res`/`ext`/`ext1`/`ext2`/`none`), sides, resonant integration limits, and packed
+node positions for interval `ising` (0-based, `0:msing`). Port of Fortran `gal_make_grid`
+(gal.f), `dx1dx2_flag=true` path. `intvl.x`/`intvl.dx` have length `nx+1`, `intvl.cells` length
+`nx` (preallocated by the caller). `sings` is the filtered list of resonant surfaces (1-based,
+`sings[ising]` = Fortran `sing(ising)`).
+
+Note (faithful to Fortran): the lower-surface width scale uses `|nn·q1|` (gal.f) but the
+upper-surface scale uses `|q1|` without `nn` (gal.f). Identical for `nn = 1`.
+"""
+function gal_make_grid!(intvl::GalInterval, ising::Int, msing::Int, sings::Vector{SingType},
+    nn::Int, psilow::Float64, psihigh::Float64, ctrl::ForceFreeStatesControl)
+
+    nx = ctrl.gal_nx
+    dx0 = ctrl.gal_dx0
+    dx1 = ctrl.gal_dx1
+    dx2 = ctrl.gal_dx2
+    pfac = ctrl.gal_pfac
+    cutoff = ctrl.gal_cutoff
+    xn = intvl.x
+    cells = intvl.cells
+
+    # Reset cell flags
+    for c in cells
+        c.etype = GCT_NONE
+        c.extra = GAL_SIDE_NONE
+    end
+
+    # --- lower bound (gal.f) ---
+    local ixmin::Int
+    if ising == 0
+        ixmin = 0
+        xn[0+1] = psilow
+    else
+        x0 = sings[ising].psifac
+        nq1 = abs(nn * sings[ising].q1)
+        xn[0+1] = x0
+        cells[1].x_lsode = x0 + dx0 / nq1
+        if ctrl.gal_dx1dx2_flag
+            xn[1+1] = x0 + dx1 / nq1
+            xn[2+1] = x0 + (dx1 + dx2) / nq1
+            ixmin = 2
+        else
+            ixmin = 0
+        end
+        cells[1].extra = GAL_SIDE_RIGHT
+        cells[2].extra = GAL_SIDE_RIGHT
+        cells[1].etype = GCT_RES
+        cells[2].etype = GCT_EXT
+        for ix in 3:(cutoff+1)
+            cells[ix].etype = GCT_EXT1
+            cells[ix].extra = GAL_SIDE_RIGHT
+        end
+        ix_ext2 = max(cutoff + 2, 3)   # Fortran loop-exit value of ix (gal.f)
+        (ix_ext2 > nx || cells[ix_ext2].etype != GCT_NONE) &&
+            error("gal_make_grid: too many elements include the big solution (right side)")
+        cells[ix_ext2].etype = GCT_EXT2
+        cells[ix_ext2].extra = GAL_SIDE_RIGHT
+    end
+
+    # --- upper bound (gal.f) ---
+    local ixmax::Int
+    if ising == msing
+        ixmax = nx
+        xn[nx+1] = psihigh
+    else
+        x1 = sings[ising+1].psifac
+        nq1 = abs(sings[ising+1].q1)   # NOTE: no nn factor, faithful to gal.f
+        xn[nx+1] = x1
+        cells[nx].x_lsode = x1 - dx0 / nq1
+        if ctrl.gal_dx1dx2_flag
+            xn[(nx-1)+1] = x1 - dx1 / nq1
+            xn[(nx-2)+1] = x1 - (dx1 + dx2) / nq1
+            ixmax = nx - 2
+        else
+            ixmax = nx
+        end
+        cells[nx-1].extra = GAL_SIDE_LEFT
+        cells[nx].extra = GAL_SIDE_LEFT
+        cells[nx-1].etype = GCT_EXT
+        cells[nx].etype = GCT_RES
+        for ix in (nx-2):-1:(nx-cutoff)
+            cells[ix].etype = GCT_EXT1
+            cells[ix].extra = GAL_SIDE_LEFT
+        end
+        ix_ext2 = min(nx - cutoff - 1, nx - 2)   # Fortran loop-exit value of ix (gal.f)
+        (ix_ext2 < 1 || cells[ix_ext2].etype != GCT_NONE) &&
+            error("gal_make_grid: too many elements include the big solution (left side)")
+        cells[ix_ext2].etype = GCT_EXT2
+        cells[ix_ext2].extra = GAL_SIDE_LEFT
+    end
+
+    # --- interior packed grid (gal.f) ---
+    x0 = xn[ixmin+1]
+    x1 = xn[ixmax+1]
+    xm = (xn[ixmax+1] + xn[ixmin+1]) / 2
+    dxh = (xn[ixmax+1] - xn[ixmin+1]) / 2
+    mx = (ixmax - ixmin) ÷ 2
+    # Default: symmetric "both" pack (faithful to gal.f). With gal_edge_onesided, the two end
+    # intervals pack only toward their single rational end, leaving the regular boundary (psilow at
+    # interval 0, psihigh at interval msing) at coarse spacing — removes the gratuitous fine edge cell.
+    side = "both"
+    if ctrl.gal_edge_onesided && msing > 0
+        if ising == 0
+            side = "right"   # rational sits at the upper (ixmax) end
+        elseif ising == msing
+            side = "left"    # rational sits at the lower (ixmin) end
+        end
+    end
+    packed = gal_pack(mx, pfac, side)          # length 2*mx+1 = (ixmax-ixmin)+1
+    for k in 0:(ixmax-ixmin)
+        xn[ixmin+1+k] = xm + dxh * packed[k+1]
+    end
+    xn[ixmin+1] = x0
+    xn[ixmax+1] = x1
+
+    intvl.dx[0+1] = 0.0
+    for i in 1:nx
+        intvl.dx[i+1] = xn[i+1] - xn[i]
+    end
+    for ix in 1:nx
+        cells[ix].x = (xn[ix], xn[ix+1])       # Fortran (x(ix-1), x(ix))
+    end
+    return intvl
+end
+
+"""
+    gal_make_map!(ws::GalWorkspace, mpert::Int)
+
+Build the local→global DOF numbering with the `2*mpert` C¹ overlap between adjacent Hermite cells and
+the interleaved extra (small resonant) DOFs (`emap`). Port of Fortran `gal_make_map` (gal.f).
+Sets `ws.ndim`. The `skip`/`emap` bookkeeping reserves a global slot for each resonant coefficient and
+shares it between the `res` cell and its neighbouring `ext` cell.
+"""
+function gal_make_map!(ws::GalWorkspace, mpert::Int)
+    np = ws.np
+    msing = length(ws.intvl) - 1
+    imap = 1
+    skip = false
+    emap = 0
+    for ising in 0:msing
+        intvl = ws.intvl[ising+1]
+        for ix in 1:ws.nx
+            cell = intvl.cells[ix]
+            # Nonresonant Hermite DOFs (with 2*mpert overlap onto the previous cell)
+            if imap > 1
+                imap -= 2 * mpert
+            end
+            for ip in 0:np
+                for ipert in 1:mpert
+                    cell.map[ipert, ip+1] = imap
+                    imap += 1
+                end
+                if ip == 1 && skip
+                    imap += 1
+                    skip = false
+                end
+            end
+            # Extra (small resonant) DOF, shared between res and adjacent ext cell
+            if cell.extra == GAL_SIDE_LEFT
+                if cell.etype == GCT_EXT
+                    emap = imap
+                    cell.emap = emap
+                    skip = true
+                elseif cell.etype == GCT_RES
+                    cell.emap = emap
+                end
+            elseif cell.extra == GAL_SIDE_RIGHT
+                if cell.etype == GCT_RES
+                    emap = imap
+                    cell.emap = emap
+                    skip = true
+                elseif cell.etype == GCT_EXT
+                    cell.emap = emap
+                end
+            end
+        end
+    end
+    ws.ndim = imap - 1
+    return ws
+end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl
new file mode 100644
index 000000000..2bdb1d2d2
--- /dev/null
+++ b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl
@@ -0,0 +1,182 @@
+# GalerkinMatch.jl
+#
+# DRIVEN (RPEC) outer↔inner asymptotic matching. Port of rmatch `match_rpec` (match.f) and the
+# outer total-solution construction of `match_output_solution` (match.f).
+#
+# Per rational surface: the forced eigenvalue γ_s = 2πi·n·f_s (rpec, match.f), the inner-layer
+# matching data Δ(Q) via the Galerkin GGJ solver (InnerLayer.solve_inner), then the 4·msing matching
+# system mat·[cout;cin] = rmat is assembled (rmat = −transpose of the gal Δ′ coil block) and solved for
+# the per-coil outer/inner coefficients. The matched outer solution for coil drive j is
+#   ξ_j = Σ_{isol=1}^{2 msing} cout[isol,j]·sols[:,:,isol] + sols[:,:,2 msing+j]
+# and identically for ξ′ with sols_deriv. Stacking over j gives the matched fundamental matrix in the
+# identity-at-edge basis (coil column j has ξ_edge = e_j).
+
+"""
+    gal_match_rpec(ctrl, equil, intr, gal_result) -> GalMatchResult
+
+Solve the coil-driven RPEC matching from the outer Δ′ (`gal_result`) and the per-surface inner-layer Δ.
+Requires `gal_result.solution` (reconstructed outer ξ/ξ′) and the rpec coil block `gal_result.delta_coil`.
+
+The resistive path uses per-surface inputs `ctrl.gal_eta`/`gal_rho`/`gal_rotation` (length `msing`) +
+`ctrl.gal_gamma`. When `ctrl.gal_ideal_flag`, the inner layer is skipped and `cout=0`, so the matched
+solution is the bare ideal coil column (Fortran rmatch `coil%ideal_flag`) — the DCON/EL reference.
+"""
+function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStatesInternal,
+    gal_result::GalerkinResult)
+
+    msing = gal_result.msing
+    mpert = intr.numpert_total
+    nn = intr.nlow
+    mcoil = mpert   # rpec: one coil column per poloidal harmonic
+
+    gal_result.solution !== nothing ||
+        error("gal_match_rpec: gal_result.solution is missing (need the gal-reconstructed ξ/ξ′)")
+    isempty(gal_result.delta_coil) &&
+        error("gal_match_rpec: delta_coil is empty — gal_rpec_flag must be true for RPEC matching")
+
+    if ctrl.gal_ideal_flag
+        # Ideal limit (Fortran rmatch coil%ideal_flag, match.f): skip the inner layer entirely
+        # and set the resistive plasma combination to zero. The gal outer solve is already fully ideal, so
+        # the shared construction below collapses to the bare ideal coil column sols(:,:,csol).
+        cout = zeros(ComplexF64, 2msing, mcoil)
+        cin = zeros(ComplexF64, 2msing, mcoil)
+        deltar = zeros(ComplexF64, msing, 2)
+        rpec_eig = zeros(ComplexF64, msing)
+        residual = 0.0
+    else
+        for (name, v) in (("gal_eta", ctrl.gal_eta), ("gal_rho", ctrl.gal_rho), ("gal_rotation", ctrl.gal_rotation))
+            length(v) == msing || error("gal_match_rpec: $name has length $(length(v)), expected msing=$msing (one value per surface, core→edge)")
+        end
+
+        # Re-derive the gal singular-surface set (must match galerkin_solve's filter) to get the SingType
+        # objects resist_eval needs; verify alignment with the stored gal Δ′ surfaces.
+        psilow = intr.psilow > 0 ? intr.psilow : equil.profiles.xs[1]
+        psihigh = intr.psilim
+        sings = [s for s in intr.sing if psilow < s.psifac < psihigh && intr.mlow <= s.m[1] <= intr.mhigh]
+        length(sings) == msing && isapprox([s.psifac for s in sings], gal_result.sing_psi; rtol=1e-8) ||
+            error("gal_match_rpec: re-derived surface set does not match gal Δ′ surfaces (filter drift?)")
+
+        # --- inner-layer matching data Δ(Q) per surface (deltac_run; match.f) ---
+        inner = InnerLayer.GGJModel(solver=:galerkin)
+        deltar = zeros(ComplexF64, msing, 2)
+        rpec_eig = zeros(ComplexF64, msing)
+        for i in 1:msing
+            params = resist_eval(sings[i], equil, intr; eta=ctrl.gal_eta[i], rho=ctrl.gal_rho[i],
+                gamma=ctrl.gal_gamma, ising=i)
+            γ = 2π * im * nn * ctrl.gal_rotation[i]    # rpec forced eigenvalue 2πi·n·f
+            rpec_eig[i] = γ
+            Δ = InnerLayer.solve_inner(inner, params, γ)   # (Δ₁, Δ₂) in deltac.f convention
+            deltar[i, 1] = Δ[1]
+            deltar[i, 2] = Δ[2]
+        end
+
+        # --- assemble the 4·msing matching system (match.f) ---
+        delta_out = gal_result.delta[1:2msing, 1:2msing]   # outer Δ′ plasma block
+        mat = zeros(ComplexF64, 4msing, 4msing)
+        rmat = zeros(ComplexF64, 4msing, mcoil)
+        @views mat[2msing+1:4msing, 1:2msing] .= transpose(delta_out)          # Δ_out
+        @views rmat[2msing+1:4msing, :] .= .-transpose(gal_result.delta_coil)  # −Δ_coil source
+        for ising in 1:msing
+            idx1 = 2ising - 1
+            idx2 = 2ising
+            idx3 = idx1 + 2msing
+            idx4 = idx2 + 2msing
+            delta1 = deltar[ising, 1]
+            delta2 = deltar[ising, 2]
+            mat[idx1, idx1] = 1
+            mat[idx2, idx2] = 1
+            mat[idx1, idx3] = -1
+            mat[idx1, idx4] = 1
+            mat[idx2, idx3] = -1
+            mat[idx2, idx4] = -1
+            mat[idx3, idx3] = -delta1   # Δ_in
+            mat[idx3, idx4] = delta2
+            mat[idx4, idx3] = -delta1
+            mat[idx4, idx4] = -delta2
+        end
+
+        # --- solve mat·cof = rmat for the outer/inner coefficients (match.f) ---
+        cof = mat \ rmat
+        residual = norm(mat * cof - rmat) / max(norm(rmat), 1e-300)
+        cout = cof[1:2msing, :]
+        cin = cof[2msing+1:4msing, :]
+    end
+
+    # --- matched outer ξ/ξ′ per coil drive (match.f); ideal: cout=0 ⇒ bare coil column ---
+    sols = gal_result.solution.xi          # (mpert, ngrid, nsol)
+    sols_deriv = gal_result.solution.xi_deriv
+    ngrid = size(sols, 2)
+    xi = zeros(ComplexF64, mpert, ngrid, mcoil)
+    xi_deriv = zeros(ComplexF64, mpert, ngrid, mcoil)
+    for j in 1:mcoil
+        csol = 2msing + j                  # this coil's particular-solution column
+        @views xi[:, :, j] .= sols[:, :, csol]
+        @views xi_deriv[:, :, j] .= sols_deriv[:, :, csol]
+        for isol in 1:2msing
+            @views xi[:, :, j] .+= cout[isol, j] .* sols[:, :, isol]
+            @views xi_deriv[:, :, j] .+= cout[isol, j] .* sols_deriv[:, :, isol]
+        end
+    end
+
+    return GalMatchResult(cout, cin, xi, xi_deriv, deltar, rpec_eig, residual)
+end
+
+"""
+    gal_matched_odestate(gal_result, ffit, intr) -> OdeState
+
+Pack the RPEC-matched outer solution into an `OdeState` shaped exactly like the shooting integrator's,
+so `PerturbedEquilibrium` consumes it unchanged. Mirrors Fortran `idcon_build`'s gal branch
+(idcon.f) and `globalsol.bin` (the on-surface `issing` points are dropped, match.f):
+
+  - `u_store[:,:,1,ip] = ξ_ψ`  (matched fundamental matrix, mode×coil-drive, identity-at-edge basis)
+  - `ud_store[:,:,1,ip] = dξ_ψ/dψ`  (analytic)
+  - `ud_store[:,:,2,ip] = ξ_s = −A⁻¹(B·ξ′ + C·ξ)`  via `ffit` (same outer ideal-MHD relation as `sing_der!`)
+  - `u_store[:,:,2] = 0`  (PE never reads it; matches Fortran's unused u2 in the gal path)
+
+The grid is the gal-native grid (inner→edge); `step` indexes the edge so `build_flux_matrix` derives the
+edge BC from `u_store[:,:,1,step]`.
+"""
+function gal_matched_odestate(gal_result::GalerkinResult, ffit::FourFitVars, intr::ForceFreeStatesInternal)
+    gal_result.match !== nothing || error("gal_matched_odestate: no match result (run with gal_match_flag=true)")
+    sol = gal_result.solution
+    m = gal_result.match
+    mpert = intr.numpert_total
+
+    # Drop the on-surface (issing) grid points — zero placeholders where the resonant series diverges.
+    keep = .!sol.issing
+    psi_f = sol.psi[keep]
+    q_f = sol.q[keep]
+    xi_f = m.xi[:, keep, :]          # (mpert, ngrid_f, mcoil)
+    dxi_f = m.xi_deriv[:, keep, :]
+    ngrid_f = length(psi_f)
+
+    u_store = zeros(ComplexF64, mpert, mpert, 2, ngrid_f)
+    ud_store = zeros(ComplexF64, mpert, mpert, 2, ngrid_f)
+
+    amat = Matrix{ComplexF64}(undef, mpert, mpert)
+    bmat = Matrix{ComplexF64}(undef, mpert, mpert)
+    cmat = Matrix{ComplexF64}(undef, mpert, mpert)
+    tmp = Matrix{ComplexF64}(undef, mpert, mpert)
+    hint = Ref(1)
+    for ip in 1:ngrid_f
+        ξ = @view xi_f[:, ip, :]
+        ξ′ = @view dxi_f[:, ip, :]
+        @views u_store[:, :, 1, ip] .= ξ
+        @views ud_store[:, :, 1, ip] .= ξ′
+        # ξ_s = −A⁻¹(B·ξ′ + C·ξ), exactly as the ideal path of sing_der! (Sing.jl:1015-1049)
+        ffit.amats(vec(amat), psi_f[ip]; hint=hint)
+        ffit.bmats(vec(bmat), psi_f[ip]; hint=hint)
+        ffit.cmats(vec(cmat), psi_f[ip]; hint=hint)
+        LinearAlgebra.LAPACK.potrf!('U', amat)
+        LinearAlgebra.LAPACK.potrs!('U', amat, bmat)   # bmat ← A⁻¹ B
+        LinearAlgebra.LAPACK.potrs!('U', amat, cmat)   # cmat ← A⁻¹ C
+        xs = @view ud_store[:, :, 2, ip]
+        mul!(tmp, bmat, ξ′)
+        xs .= .-tmp
+        mul!(tmp, cmat, ξ)
+        xs .-= tmp
+    end
+
+    return OdeState(; numpert_total=mpert, numunorms_init=1, msing=gal_result.msing, numsteps_init=ngrid_f,
+        step=ngrid_f, total_steps=ngrid_f, psi_store=psi_f, q_store=q_f, u_store=u_store, ud_store=ud_store)
+end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolution.jl b/src/ForceFreeStates/Galerkin/GalerkinSolution.jl
new file mode 100644
index 000000000..dfed6dfe7
--- /dev/null
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolution.jl
@@ -0,0 +1,179 @@
+# GalerkinSolution.jl
+#
+# Reconstruct the outer-region radial displacement ξ(ψ) AND its analytic derivative ξ′(ψ) from the
+# solved Galerkin coefficients. Ports:
+#   gal_get_solution    (gal.f)  ξ (and optional ξ′) at one ψ for one solution column
+#   gal_output_solution (gal.f)  evaluate (ξ, ξ′) on a gal-native packed grid, all columns
+#
+# The whole point of carrying ξ′ here is to evaluate it analytically from the Hermite derivative basis
+# (gal_hermite qb) and the asymptotic-series derivative (sing_get_dua_gal) rather than ever
+# spline-differentiating ξ at the packed edge — the spline-endpoint derivative was the root cause of
+# GPEC's fake ~100% shielding (fixed in Fortran by RDCON_singcoup commit c3e8c5bc).
+#
+# Fortran module flags restore_uh/us/ul default .TRUE. (gal.f); the full reconstruction is always
+# performed. The cut path (gal.f, used only for the matching diagnostic) is NOT ported here.
+
+# Dense / coarse sampling per cell, matching Fortran interp_np_res / interp_np (gal.f).
+const GAL_INTERP_NP = 3
+const GAL_INTERP_NP_RES = 60
+
+# Advance the monotone (iintvl, icell) cursor to the cell containing x. Port of gal.f.
+@inline function gal_locate_cell(ws::GalWorkspace, x::Float64, iintvl::Int, icell::Int)
+    msing = length(ws.intvl) - 1
+    for _ in 1:((msing+1)*ws.nx+1)
+        cell = ws.intvl[iintvl+1].cells[icell]
+        if cell.x[1] <= x <= cell.x[2]
+            return iintvl, icell, cell
+        end
+        icell += 1
+        if icell > ws.nx
+            icell = 1
+            iintvl += 1
+        end
+        if iintvl > msing
+            iintvl = 0
+        end
+    end
+    error("gal_locate_cell: x=$x not found in any cell")
+end
+
+"""
+    gal_get_solution(ws, asymps, sings, intr, x, iintvl, icell, isol; want_d=true)
+        -> (sol, dsol, iintvl, icell)
+
+Reconstruct the displacement `sol` (length `mpert`) — and, when `want_d`, its analytic radial
+derivative `dsol = d(sol)/dψ` — at flux `x` for solution column `isol`, from the solved Galerkin
+coefficients `ws.sol`. Port of `gal_get_solution` (gal.f, non-cut path). `iintvl`/`icell`
+are a monotone cell cursor advanced and returned for the next call. `dsol` is `nothing` if `!want_d`.
+"""
+function gal_get_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings::Vector{SingType},
+    intr::ForceFreeStatesInternal, x::Float64, iintvl::Int, icell::Int, isol::Int; want_d::Bool=true)
+
+    msing = length(ws.intvl) - 1
+    mpert = intr.numpert_total
+    np = GAL_NP
+    iintvl, icell, cell = gal_locate_cell(ws, x, iintvl, icell)
+
+    sol = zeros(ComplexF64, mpert)
+    dsol = want_d ? zeros(ComplexF64, mpert) : nothing
+
+    # --- non-resonant (Hermite) part: same coefficients, value basis pb → ξ, derivative basis qb → ξ′ ---
+    pbt, qbt = gal_hermite(x, cell.x[1], cell.x[2])
+    pb = collect(pbt)
+    qb = collect(qbt)
+    if iintvl == msing && icell == ws.nx          # last edge cell: swap node-2/3 DOFs (gal.f)
+        pb[3], pb[4] = pb[4], pb[3]
+        qb[3], qb[4] = qb[4], qb[3]
+    end
+    for ip in 0:np
+        u = @view ws.sol[cell.map[:, ip+1], isol]
+        @views sol .+= u .* pb[ip+1]
+        want_d && (@views dsol .+= u .* qb[ip+1])
+    end
+
+    cell.etype == GCT_NONE && return sol, dsol, iintvl, icell
+
+    # --- small & large resonant/extension solutions (gal.f) ---
+    jsing = _cell_jsing(cell, iintvl)             # left → iintvl+1, right → iintvl (gal_get_solution)
+    asymp = asymps[jsing]
+    psi_s = sings[jsing].psifac
+    ipert0 = _ipert_res(sings[jsing], intr)       # resonant (big) row; small row = ipert0+mpert
+
+    # Edge Hermite weights (fresh, unswapped) for ext/ext2 cells (gal.f)
+    pbe, qbe = gal_hermite(x, cell.x[1], cell.x[2])
+    if cell.extra == GAL_SIDE_LEFT
+        epb1, epb2 = pbe[3], pbe[4]
+        eqb1, eqb2 = qbe[3], qbe[4]
+        xext = cell.x[2]
+    else
+        epb1, epb2 = pbe[1], pbe[2]
+        eqb1, eqb2 = qbe[1], qbe[2]
+        xext = cell.x[1]
+    end
+
+    need_in = cell.etype == GCT_RES || cell.etype == GCT_EXT || cell.etype == GCT_EXT1
+    need_ext = cell.etype == GCT_EXT || cell.etype == GCT_EXT2
+    ua = need_in ? sing_get_ua_gal(asymp, x - psi_s) : nothing
+    dua = (need_in && want_d) ? sing_get_dua_gal(asymp, x - psi_s) : nothing
+    uaext = need_ext ? sing_get_ua_gal(asymp, xext - psi_s) : nothing
+    duaext = need_ext ? sing_get_dua_gal(asymp, xext - psi_s) : nothing
+
+    # small solution Δ-coefficient restore (ext/res only), coeff = ws.sol[emap,isol] (gal.f)
+    if cell.etype == GCT_RES || cell.etype == GCT_EXT
+        delta = ws.sol[cell.emap, isol]
+        if cell.etype == GCT_RES
+            @views sol .+= delta .* ua[:, ipert0+mpert, 1]
+            want_d && (@views dsol .+= delta .* dua[:, ipert0+mpert, 1])
+        else  # GCT_EXT: series held at the fixed edge xext, x-dependence only via Hermite weights
+            @views sol .+= delta .* (epb1 .* uaext[:, ipert0+mpert, 1] .+ epb2 .* duaext[:, ipert0+mpert, 1])
+            want_d && (@views dsol .+= delta .* (eqb1 .* uaext[:, ipert0+mpert, 1] .+ eqb2 .* duaext[:, ipert0+mpert, 1]))
+        end
+    end
+
+    # large solution (coeff 1) only for the column driving this surface/side (gal.f)
+    if (isol == 2jsing - 1 && cell.extra == GAL_SIDE_LEFT) ||
+       (isol == 2jsing && cell.extra == GAL_SIDE_RIGHT)
+        if cell.etype == GCT_RES || cell.etype == GCT_EXT || cell.etype == GCT_EXT1
+            @views sol .+= ua[:, ipert0, 1]
+            want_d && (@views dsol .+= dua[:, ipert0, 1])
+        elseif cell.etype == GCT_EXT2
+            @views sol .+= epb1 .* uaext[:, ipert0, 1] .+ epb2 .* duaext[:, ipert0, 1]
+            want_d && (@views dsol .+= eqb1 .* uaext[:, ipert0, 1] .+ eqb2 .* duaext[:, ipert0, 1])
+        end
+    end
+
+    return sol, dsol, iintvl, icell
+end
+
+"""
+    gal_output_solution(ws, asymps, sings, intr, profiles, psihigh) -> GalerkinSolution
+
+Build the gal-native packed radial grid (inner→edge; `GAL_INTERP_NP_RES` points per resonant/extension
+cell, `GAL_INTERP_NP` elsewhere, plus the edge point ψ=psihigh) and evaluate ξ AND the analytic ξ′ over
+it for every solution column. Port of `gal_output_solution` (gal.f); the binary/ASCII writes
+and the `b_flag` ξ→b^ψ conversion (off by default) are not ported — we keep ξ itself.
+"""
+function gal_output_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings::Vector{SingType},
+    intr::ForceFreeStatesInternal, profiles, psihigh::Float64)
+
+    mpert = intr.numpert_total
+    msing = length(ws.intvl) - 1
+
+    # --- grid (gal.f): per-cell left-inclusive sampling, then the edge point ---
+    psi = Float64[]
+    issing = Bool[]
+    for iintvl in 0:msing
+        for icell in 1:ws.nx
+            cell = ws.intvl[iintvl+1].cells[icell]
+            npts = (cell.etype == GCT_RES || cell.etype == GCT_EXT) ? GAL_INTERP_NP_RES : GAL_INTERP_NP
+            x1, x2 = cell.x
+            for k in 0:(npts-1)
+                push!(psi, x1 + (x2 - x1) * (k / npts))
+                # res-right cell's first point sits exactly on the surface (gal.f)
+                push!(issing, k == 0 && cell.extra == GAL_SIDE_RIGHT && cell.etype == GCT_RES)
+            end
+        end
+    end
+    push!(psi, psihigh)
+    push!(issing, false)
+
+    qhint = Ref(1)
+    q = [profiles.q_spline(p; hint=qhint) for p in psi]
+
+    # --- evaluate (gal.f) ---
+    ngrid = length(psi)
+    xi = zeros(ComplexF64, mpert, ngrid, ws.nsol)
+    xi_deriv = zeros(ComplexF64, mpert, ngrid, ws.nsol)
+    for isol in 1:ws.nsol
+        iintvl = 0
+        icell = 1
+        for ip in 1:ngrid
+            issing[ip] && continue
+            sol, dsol, iintvl, icell = gal_get_solution(ws, asymps, sings, intr, psi[ip], iintvl, icell, isol)
+            @views xi[:, ip, isol] .= sol
+            @views xi_deriv[:, ip, isol] .= dsol
+        end
+    end
+
+    return GalerkinSolution(psi, q, issing, xi, xi_deriv)
+end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
new file mode 100644
index 000000000..01b4670c0
--- /dev/null
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -0,0 +1,273 @@
+# GalerkinSolve.jl
+#
+# Top-level driver for the RDCON outer-region singular Galerkin Δ′ solve, plus the per-cell assembly
+# orchestration, the banded solve, Δ′ extraction, PEST-3 blocks, and HDF5 output.
+# Ports gal_make_arrays (gal.f), gal_solve (gal.f), and gal_write_pest3_data
+# (gal.f). The DRIVEN/RPEC inner-layer matching is wired in via gal_match_rpec (GalerkinMatch.jl).
+
+"""
+    gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge)
+
+Assemble the global banded matrix and RHS. For each cell: Gauss-Lobatto Hermite stiffness
+(`gal_gauss_quad!`), then the resonant (`gal_resonant!`) or extension (`gal_extension!`) contributions;
+then the boundary conditions and the scatter into `ws.mat`/`ws.rhs`. Port of `gal_make_arrays`.
+"""
+function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
+    intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType},
+    nn::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}})
+
+    mpert = intr.numpert_total
+    msing = length(ws.intvl) - 1
+    profiles = equil.profiles
+    nodes, weights = gausslobatto(ws.nq + 1)
+
+    for ising in 0:msing
+        for (ix, cell) in enumerate(ws.intvl[ising+1].cells)
+            swap_edge = (ising == msing && ix == ws.nx)
+            gal_gauss_quad!(cell, ffit, profiles, intr, nodes, weights, swap_edge)
+            if cell.etype == GCT_RES
+                gal_resonant!(cell, ising, ffit, profiles, intr, asymps, sings, nn, ctrl.gal_tol, ctrl.gal_gnstep, ctrl.verbose)
+            elseif cell.etype == GCT_EXT || cell.etype == GCT_EXT1 || cell.etype == GCT_EXT2
+                gal_extension!(cell, ising, ffit, profiles, intr, asymps, sings, nn, nodes, weights)
+            end
+        end
+    end
+
+    gal_set_boundary!(ws, mpert, wv_edge)
+    gal_assemble_mat!(ws, mpert)
+    gal_assemble_rhs!(ws, mpert)
+    return ws
+end
+
+"""
+    galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
+                   intr::ForceFreeStatesInternal; vac_data=nothing) -> GalerkinResult
+
+Compute the outer-region Δ′ matching matrix by the singular Galerkin method. Port of `gal_solve`
+(gal.f). Single toroidal mode only (`intr.npert == 1`). Returns a `GalerkinResult`; if there
+are no resonant surfaces in the domain it returns an empty result.
+
+`vac_data` (a `VacuumData` from `free_run!`) supplies the free-boundary edge term
+`wv_edge = vac_data.wv · psio²`; pass `nothing` for a fixed-boundary edge.
+"""
+function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
+    intr::ForceFreeStatesInternal; vac_data=nothing)
+
+    intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported")
+    ctrl.gal_solver in ("LU", "cholesky") || error("galerkin_solve: gal_solver must be \"LU\" or \"cholesky\"")
+    nn = intr.nlow
+    mpert = intr.numpert_total
+    np = GAL_NP
+
+    # Resonant surfaces inside the integration domain, in increasing ψ (Fortran `sing(1:msing)`).
+    # psilow is raised above the axis by sing_min! when qlow > qmin (excludes the q<qlow core);
+    # fall back to the equilibrium axis bound if sing_min! was not run (psilow still 0).
+    psilow = intr.psilow > 0 ? intr.psilow : equil.profiles.xs[1]
+    psihigh = intr.psilim
+    sings = [s for s in intr.sing if psilow < s.psifac < psihigh && intr.mlow <= s.m[1] <= intr.mhigh]
+    msing = length(sings)
+    if msing == 0
+        ctrl.verbose && @info "galerkin_solve: no resonant surfaces in domain; skipping Δ′ solve"
+        empty2 = Matrix{ComplexF64}(undef, 0, 0)
+        return GalerkinResult(Matrix{ComplexF64}(undef, 0, 0), empty2, empty2, empty2, empty2, 0,
+            Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], empty2, nothing, nothing)
+    end
+
+    ctrl.verbose && @info "Starting outer-region Galerkin Δ′ solve (msing=$msing, solver=$(ctrl.gal_solver))"
+
+    # Per-surface two-sided asymptotic series (pr178's validated convention, matching Fortran sing.f
+    # vmatr/vmatl): right = sig=+1, left = sig=-1, no √det normalization. The Mercier exponent α is a
+    # property of the surface, so the left series reuses the right's α (alpha_override). The order is
+    # raised to gal_sing_order + ceil(2·Re(α)) for high-Mercier-index surfaces (Fortran sing1_vmat).
+    ctrl_gal = deepcopy(ctrl)
+    asymps = GalSingAsymp[]
+    for s in sings
+        ctrl_gal.sing_order = ctrl.gal_sing_order
+        ar = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=1.0)
+        if ctrl.gal_sing_order_ceiling
+            order = ctrl.gal_sing_order + ceil(Int, 2 * real(ar.alpha[1]))
+            if order > ctrl.gal_sing_order
+                ctrl_gal.sing_order = order
+                ar = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=1.0)
+            end
+        end
+        al = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=-1.0, alpha_override=ar.alpha)
+        push!(asymps, GalSingAsymp(ar, al))
+    end
+
+    # --- allocate workspace ---
+    nx = ctrl.gal_nx
+    kl = mpert * (np + 1)
+    ku = kl
+    ldab = ctrl.gal_solver == "LU" ? 2kl + ku + 1 : kl + 1
+    # rpec_flag (RDCON, gal.f): append mpert coil-response columns. Each is a unit source at
+    # the edge value DOF in one poloidal mode; the recorded plasma response is the coil block of Δ_gw.
+    # Cholesky's lower-only edge zeroing can't represent the rpec identity edge, so require LU.
+    ncoil = ctrl.gal_rpec_flag ? mpert : 0
+    if ncoil > 0 && ctrl.gal_solver != "LU"
+        error("galerkin_solve: gal_rpec_flag=true requires gal_solver=\"LU\" (coil edge BC needs the full-band path)")
+    end
+    nsol = 2 * msing + ncoil
+    intvl = [GalInterval(zeros(Float64, nx + 1), zeros(Float64, nx + 1), [GalCell(mpert) for _ in 1:nx])
+             for _ in 0:msing]
+    ws = GalWorkspace(ctrl.gal_solver, nx, ctrl.gal_nq, np, 0, kl, ku, ldab, nsol,
+        intvl, Matrix{ComplexF64}(undef, 0, 0), Matrix{ComplexF64}(undef, 0, 0), Matrix{ComplexF64}(undef, 0, 0))
+
+    for ising in 0:msing
+        gal_make_grid!(ws.intvl[ising+1], ising, msing, sings, nn, psilow, psihigh, ctrl)
+    end
+    gal_make_map!(ws, mpert)
+
+    ws.mat = zeros(ComplexF64, ldab, ws.ndim)
+    ws.rhs = zeros(ComplexF64, ws.ndim, nsol)
+    ws.sol = zeros(ComplexF64, ws.ndim, nsol)
+
+    # Free-boundary edge term wvac·psio² (vac_data.wv is already singfac-scaled at qlim in free_run!).
+    # rpec passes wv_edge=nothing; gal_set_boundary! then applies the identity edge AND injects the coil
+    # unit sources (its three-way branch). The vacuum block is only built for the non-rpec free case.
+    wv_edge = nothing
+    if ctrl.vac_flag && vac_data !== nothing && ncoil == 0
+        wv_edge = Matrix{ComplexF64}(vac_data.wv .* equil.psio^2)
+    end
+
+    gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge)
+
+    if ctrl.verbose
+        offdbg = ws.solver == "LU" ? ws.kl + ws.ku + 1 : 1
+        @info "Galerkin assembly: ndim=$(ws.ndim) nsol=$nsol kl=$(ws.kl) ldab=$(ws.ldab) |rhs|=$(norm(ws.rhs))"
+        for ising in 0:msing
+            for cell in ws.intvl[ising+1].cells
+                if cell.etype == GCT_RES || cell.etype == GCT_EXT
+                    @info "  cell ising=$ising etype=$(cell.etype) side=$(cell.extra) emap=$(cell.emap) ediag=$(cell.ediag) erhs=$(cell.erhs) |rhs|=$(norm(cell.rhs)) |emat|=$(norm(cell.emat)) matdiag[emap]=$(ws.mat[offdbg, cell.emap])"
+                end
+            end
+        end
+    end
+
+    # --- banded solve (gal.f) ---
+    ws.sol .= ws.rhs
+    if ws.solver == "LU"
+        ctrl.verbose && @info "Galerkin LU banded factorization + solve"
+        ab, ipiv = LinearAlgebra.LAPACK.gbtrf!(ws.kl, ws.ku, ws.ndim, ws.mat)
+        LinearAlgebra.LAPACK.gbtrs!('N', ws.kl, ws.ku, ws.ndim, ab, ipiv, ws.sol)
+    else
+        ctrl.verbose && @info "Galerkin Cholesky banded factorization + solve"
+        LinearAlgebra.LAPACK.pbtrf!('L', ws.kl, ws.mat)
+        LinearAlgebra.LAPACK.pbtrs!('L', ws.kl, ws.mat, ws.sol)
+    end
+
+    # --- extract Δ′ from the small resonant coefficients (gal.f) ---
+    delta = zeros(ComplexF64, nsol, 2 * msing)
+    jsol = 0
+    for ising in 0:msing
+        for cell in ws.intvl[ising+1].cells
+            if cell.etype == GCT_RES
+                jsol += 1
+                @views delta[:, jsol] .= ws.sol[cell.emap, :]
+            end
+        end
+    end
+
+    Ap, Bp, Gammap, Deltap = gal_pest3_blocks(delta, msing)
+
+    # Mercier index D_I = -Re(α²); α is the small/large solution exponent (a surface property, taken
+    # from the right series — the left series reuses the same α via alpha_override above).
+    di = [real(-asymps[i].right.alpha[1]^2) for i in 1:msing]
+    alpha = [asymps[i].right.alpha[1] for i in 1:msing]
+    # Coil-response block (rpec_flag): rows 2*msing+1 : 2*msing+mpert of delta (empty otherwise).
+    delta_coil = ncoil > 0 ? delta[(2*msing+1):(2*msing+ncoil), :] : Matrix{ComplexF64}(undef, 0, 0)
+    if ncoil > 0
+        @info "Galerkin Δ′ solve: coil-response block delta_coil = $(size(delta_coil))"
+        @info "Galerkin Δ′ solve: coil-response 1,1 = $(delta_coil[1, 1])"
+    end
+
+    # Reconstruct ξ(ψ) AND analytic ξ′(ψ) on the gal-native grid (gal_output_solution).
+    ctrl.verbose && @info "Reconstructing outer-region ξ and analytic ξ′ on the gal grid"
+    solution = gal_output_solution(ws, asymps, sings, intr, equil.profiles, psihigh)
+
+    sing_psi = [s.psifac for s in sings]
+    sing_q = [s.q for s in sings]
+    sing_m = [s.m[1] for s in sings]
+    sing_n = [s.n[1] for s in sings]
+    result = GalerkinResult(delta, Ap, Bp, Gammap, Deltap, msing,
+        sing_psi, sing_q, sing_m, sing_n, di, alpha, delta_coil, solution, nothing)
+
+    # DRIVEN (RPEC) outer↔inner matching: build the coil-driven matched ξ/ξ′ (gal_match_rpec).
+    ctrl.gal_match_flag || return result
+    ctrl.gal_rpec_flag || error("galerkin_solve: gal_match_flag=true requires gal_rpec_flag=true")
+    ctrl.verbose && @info(ctrl.gal_ideal_flag ?
+                          "RPEC matching: IDEAL solution (inner layer skipped, bare coil columns)" :
+                          "RPEC matching: inner-layer Δ(Q) + outer↔inner solve for the coil-driven ξ")
+    match = gal_match_rpec(ctrl, equil, intr, result)
+    ctrl.gal_ideal_flag || (ctrl.verbose && @info "RPEC matching: linear-solve residual = $(match.residual)")
+    return GalerkinResult(delta, Ap, Bp, Gammap, Deltap, msing,
+        sing_psi, sing_q, sing_m, sing_n, di, alpha, delta_coil, solution, match)
+end
+
+"""
+    gal_pest3_blocks(delta, msing) -> (Ap, Bp, Gammap, Deltap)
+
+PEST-3 matching blocks (each `msing×msing`) as ± combinations of the left/right Δ′ entries. Port of
+`gal_write_pest3_data`. Row index `2i-1`/`2i` = left/right of surface `i`.
+"""
+function gal_pest3_blocks(delta::Matrix{ComplexF64}, msing::Int)
+    Ap = zeros(ComplexF64, msing, msing)
+    Bp = zeros(ComplexF64, msing, msing)
+    Gammap = zeros(ComplexF64, msing, msing)
+    Deltap = zeros(ComplexF64, msing, msing)
+    for i in 1:msing, j in 1:msing
+        d11 = delta[2i-1, 2j-1]
+        d12 = delta[2i-1, 2j]
+        d21 = delta[2i, 2j-1]
+        d22 = delta[2i, 2j]
+        Ap[i, j] = d22 + d21 + d12 + d11
+        Bp[i, j] = d22 - d21 + d12 - d11
+        Gammap[i, j] = d22 + d21 - d12 - d11
+        Deltap[i, j] = d22 - d21 - d12 + d11
+    end
+    return Ap, Bp, Gammap, Deltap
+end
+
+"""
+    write_galerkin!(out_h5, result::GalerkinResult)
+
+Write the Galerkin Δ′ outputs into the open HDF5 file under the `galerkin/` group. Replaces the Fortran
+`delta_gw`/`pest3_data` ASCII/binary outputs.
+"""
+function write_galerkin!(out_h5, result::GalerkinResult)
+    out_h5["galerkin/msing"] = result.msing
+    result.msing == 0 && return nothing
+    out_h5["galerkin/delta"] = result.delta
+    out_h5["galerkin/pest3_A"] = result.Ap
+    out_h5["galerkin/pest3_B"] = result.Bp
+    out_h5["galerkin/pest3_Gamma"] = result.Gammap
+    out_h5["galerkin/pest3_Delta"] = result.Deltap
+    out_h5["galerkin/sing_psi"] = result.sing_psi
+    out_h5["galerkin/sing_q"] = result.sing_q
+    out_h5["galerkin/sing_m"] = result.sing_m
+    out_h5["galerkin/sing_n"] = result.sing_n
+    out_h5["galerkin/di"] = result.di
+    out_h5["galerkin/alpha"] = result.alpha
+    if !isempty(result.delta_coil)
+        out_h5["galerkin/delta_coil"] = result.delta_coil
+    end
+    if result.solution !== nothing
+        sol = result.solution
+        out_h5["galerkin/solution/psi"] = sol.psi
+        out_h5["galerkin/solution/q"] = sol.q
+        out_h5["galerkin/solution/issing"] = collect(sol.issing)
+        out_h5["galerkin/solution/xi"] = sol.xi
+        out_h5["galerkin/solution/xi_deriv"] = sol.xi_deriv
+    end
+    if result.match !== nothing
+        m = result.match
+        out_h5["galerkin/match/cout"] = m.cout
+        out_h5["galerkin/match/cin"] = m.cin
+        out_h5["galerkin/match/xi"] = m.xi
+        out_h5["galerkin/match/xi_deriv"] = m.xi_deriv
+        out_h5["galerkin/match/deltar"] = m.deltar
+        out_h5["galerkin/match/rpec_eig"] = m.rpec_eig
+        out_h5["galerkin/match/residual"] = m.residual
+    end
+    return nothing
+end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
new file mode 100644
index 000000000..b1b3819d3
--- /dev/null
+++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
@@ -0,0 +1,195 @@
+# GalerkinStructs.jl
+#
+# Data structures for the RDCON outer-region singular Galerkin solver — the Julia port of
+# `gal.f` (Dewar's Galerkin method for singular modes). These mirror the Fortran derived types
+# `cell_type`, `interval_type`, and `gal_type` (gal.f).
+#
+# Fortran uses 0-based dimensions `0:np` and `0:nx`/`0:msing`; here those become Julia 1-based with a
+# `+1` index shift (Hermite DOF `ip` ↔ array index `ip+1`; interval `ising` ↔ `intvl[ising+1]`).
+#
+# Reference: A. H. Glasser, Z. R. Wang & J.-K. Park, Phys. Plasmas 23, 112506 (2016).
+
+# Hermite-cubic order: 4 DOFs per node-pair (value + slope at each end), indexed 0:np in Fortran.
+const GAL_NP = 3
+
+"""
+Cell type flanking a singular surface (Fortran `cell%etype`, gal.f).
+"""
+@enum GalCellType GCT_NONE GCT_RES GCT_EXT GCT_EXT1 GCT_EXT2
+
+"""
+Which side of the adjacent singular surface a cell lies on (Fortran `cell%extra`).
+"""
+@enum GalSide GAL_SIDE_NONE GAL_SIDE_LEFT GAL_SIDE_RIGHT
+
+"""
+    GalCell
+
+One grid cell of the Galerkin discretization. Port of Fortran `cell_type` (gal.f).
+
+## Fields
+
+  - `etype::GalCellType` — none / resonant / extension(1,2)
+  - `extra::GalSide` — none / left / right (side of the neighbouring singular surface)
+  - `emap::Int` — global index of the extra (small resonant) DOF; 0 if none
+  - `map::Matrix{Int}` — `(mpert, np+1)` local→global DOF map (Fortran `cell%map(mpert, 0:np)`)
+  - `x::NTuple{2,Float64}` — `[x_left, x_right]` of the cell
+  - `x_lsode::Float64` — resonant-integration truncation point (res cells only)
+  - `erhs::ComplexF64`, `ediag::ComplexF64` — resonant RHS and diagonal contributions
+  - `rhs::Matrix{ComplexF64}` — `(mpert, np+1)` cell RHS (driving term from the big solution)
+  - `emat::Matrix{ComplexF64}` — `(mpert, np+1)` coupling of the extra DOF to the Hermite DOFs
+  - `mat::Array{ComplexF64,4}` — `(mpert, mpert, np+1, np+1)` cell stiffness block
+"""
+mutable struct GalCell
+    etype::GalCellType
+    extra::GalSide
+    emap::Int
+    map::Matrix{Int}
+    x::NTuple{2,Float64}
+    x_lsode::Float64
+    erhs::ComplexF64
+    ediag::ComplexF64
+    rhs::Matrix{ComplexF64}
+    emat::Matrix{ComplexF64}
+    mat::Array{ComplexF64,4}
+end
+
+"""
+Construct an empty `GalCell` for `mpert` poloidal modes and Hermite order `np`.
+"""
+function GalCell(mpert::Int, np::Int=GAL_NP)
+    return GalCell(GCT_NONE, GAL_SIDE_NONE, 0,
+        zeros(Int, mpert, np + 1), (0.0, 0.0), 0.0,
+        zero(ComplexF64), zero(ComplexF64),
+        zeros(ComplexF64, mpert, np + 1),
+        zeros(ComplexF64, mpert, np + 1),
+        zeros(ComplexF64, mpert, mpert, np + 1, np + 1))
+end
+
+"""
+    GalInterval
+
+One interval between consecutive singular surfaces (or between a surface and a domain edge).
+Port of Fortran `interval_type` (gal.f). `x`/`dx` are length `nx+1` (Fortran `0:nx`).
+"""
+struct GalInterval
+    x::Vector{Float64}
+    dx::Vector{Float64}
+    cells::Vector{GalCell}
+end
+
+"""
+    GalWorkspace
+
+Global workspace for the Galerkin solve. Port of Fortran `gal_type` (gal.f).
+
+`intvl` has length `msing+1` (Fortran `intvl(0:msing)`, accessed as `intvl[ising+1]`).
+`kl = ku = mpert*(np+1)`. The banded `mat` storage layout depends on `solver` (gal.f):
+
+  - `"LU"`       → `ldab = 2*kl + ku + 1`, full band (LAPACK `gbtrf!`/`gbtrs!`)
+  - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `pbtrf!`/`pbtrs!`, `'L'`)
+"""
+mutable struct GalWorkspace
+    solver::String
+    nx::Int
+    nq::Int
+    np::Int
+    ndim::Int
+    kl::Int
+    ku::Int
+    ldab::Int
+    nsol::Int
+    intvl::Vector{GalInterval}
+    mat::Matrix{ComplexF64}
+    rhs::Matrix{ComplexF64}
+    sol::Matrix{ComplexF64}
+end
+
+"""
+    GalSingAsymp
+
+Two-sided singular asymptotics for the Galerkin solver: separate right (sig=+1) and left (sig=−1)
+power series, mirroring Fortran sing.f `vmatr`/`vmatl`. Each is evaluated at a positive distance
+`|ψ−ψ_s|` (real √, real pfac, no z<0 renormalization), so the directly-extracted Δ′ matches gal.f.
+"""
+struct GalSingAsymp
+    right::SingAsymptotics
+    left::SingAsymptotics
+end
+
+"""
+    GalerkinSolution
+
+Reconstructed outer-region radial solution on the gal-native grid (ξ and analytic ξ′ consumed by
+PerturbedEquilibrium). Populated by `gal_output_solution` (GalerkinSolution.jl).
+
+  - `psi::Vector{Float64}`, `q::Vector{Float64}` — radial grid (inner→edge) and its safety factor.
+  - `issing::Vector{Bool}` — grid points sitting on a singular surface (skipped; left zero).
+  - `xi::Array{ComplexF64,3}`, `xi_deriv::Array{ComplexF64,3}` — `(mpert, ngrid, nsol)` ξ and dξ/dψ.
+"""
+struct GalerkinSolution
+    psi::Vector{Float64}
+    q::Vector{Float64}
+    issing::Vector{Bool}
+    xi::Array{ComplexF64,3}
+    xi_deriv::Array{ComplexF64,3}
+end
+
+"""
+    GalMatchResult
+
+Coil-driven RPEC matched solution from the outer↔inner asymptotic matching (Fortran rmatch
+`match_rpec`). Populated by `gal_match_rpec` (GalerkinMatch.jl) when `gal_match_flag`.
+
+  - `cout::Matrix{ComplexF64}` — `(2·msing, mcoil)` outer-region plasma-solution coefficients.
+  - `cin::Matrix{ComplexF64}` — `(2·msing, mcoil)` inner-region coefficients.
+  - `xi::Array{ComplexF64,3}`, `xi_deriv::Array{ComplexF64,3}` — `(mpert, ngrid, mcoil)` matched ξ(ψ) and
+    analytic ξ′(ψ) on the gal grid, one column per coil drive (identity-at-edge basis).
+  - `deltar::Matrix{ComplexF64}` — `(msing, 2)` inner-layer matching data `(Δ₁, Δ₂)` per surface.
+  - `rpec_eig::Vector{ComplexF64}` — forced eigenvalues `γ_s = 2πi·n·f_s` per surface.
+  - `residual::Float64` — relative linear-solve residual `‖mat·cof − rmat‖/‖rmat‖`.
+"""
+struct GalMatchResult
+    cout::Matrix{ComplexF64}
+    cin::Matrix{ComplexF64}
+    xi::Array{ComplexF64,3}
+    xi_deriv::Array{ComplexF64,3}
+    deltar::Matrix{ComplexF64}
+    rpec_eig::Vector{ComplexF64}
+    residual::Float64
+end
+
+"""
+    GalerkinResult
+
+Outputs of the outer-region Galerkin solve.
+
+## Fields
+
+  - `delta::Matrix{ComplexF64}` — Δ′ matrix, `(nsol, 2*msing)` with `nsol = 2*msing`; the small
+    resonant coefficients (Fortran `delta`, gal.f).
+  - `Ap, Bp, Gammap, Deltap::Matrix{ComplexF64}` — PEST-3 matching blocks, each `(msing, msing)`
+    (Fortran `gal_write_pest3_data`, gal.f).
+  - `msing::Int` — number of resonant singular surfaces included.
+  - `sing_psi, sing_q::Vector{Float64}`, `sing_m, sing_n::Vector{Int}` — per-surface identifiers.
+  - `di::Vector{Float64}`, `alpha::Vector{ComplexF64}` — Mercier index and exponent per surface.
+  - `solution::Union{Nothing,GalerkinSolution}` — reconstructed radial ξ(ψ) and analytic ξ′(ψ) on the
+    gal-native grid; `nothing` if no resonant surfaces.
+"""
+struct GalerkinResult
+    delta::Matrix{ComplexF64}
+    Ap::Matrix{ComplexF64}
+    Bp::Matrix{ComplexF64}
+    Gammap::Matrix{ComplexF64}
+    Deltap::Matrix{ComplexF64}
+    msing::Int
+    sing_psi::Vector{Float64}
+    sing_q::Vector{Float64}
+    sing_m::Vector{Int}
+    sing_n::Vector{Int}
+    di::Vector{Float64}
+    alpha::Vector{ComplexF64}
+    delta_coil::Matrix{ComplexF64}   # (mpert, 2*msing) coil-response block (rpec_flag); empty if not computed
+    solution::Union{Nothing,GalerkinSolution}
+    match::Union{Nothing,GalMatchResult}   # RPEC matched solution (gal_match_flag); nothing otherwise
+end
diff --git a/src/ForceFreeStates/Resist.jl b/src/ForceFreeStates/Resist.jl
new file mode 100644
index 000000000..ecfbcd0f8
--- /dev/null
+++ b/src/ForceFreeStates/Resist.jl
@@ -0,0 +1,124 @@
+# Resist.jl
+#
+# Per-rational-surface resistive inner-layer parameters (Glasser–Greene–Johnson).
+# Julia port of the RDCON `resist_eval` (resist.f): evaluates the flux-surface-averaged
+# curvature/mass coefficients E, F, H, M, G, K and the local Alfvén / resistive time scales
+# τ_A, τ_R at a singular surface.
+#
+# Returns an `InnerLayer.GGJParameters` directly (the inner-layer matching input type), so a
+# resistive ForceFreeStates run can hand its per-surface coefficients straight to the GGJ inner
+# solver. InnerLayer is included before ForceFreeStates precisely so this type is available here.
+#
+# Resistivity η and mass density ρ are REQUIRED inputs (per surface) — there is no built-in
+# resistivity/density model. We deliberately do NOT port the Fortran `resist.f` Spitzer block
+# (hardcoded ne/Te → η, ρ); the caller supplies η and ρ from whatever transport/profile model is
+# appropriate. τ_A and τ_R are then formed from physically separable pieces: a purely geometric
+# factor times the explicit η and ρ (no RMATCH-style divide-out / re-multiply).
+#
+# Differences from the Fortran reference (deliberate):
+#   - Rational surfaces lie off the equilibrium grid (q = m/n between nodes), so the geometry
+#     is evaluated with the callable bicubic interpolants at the exact (ψ_s, θ) — using
+#     DerivOp(0,1) for θ-derivatives — rather than the nodal-grid access mercier_scan! uses.
+#   - η and ρ are caller-supplied (no Spitzer model, no hardcoded ne/Te/mi defaults).
+#
+# Reference: Glasser, Greene & Johnson, Phys. Fluids 18, 875 (1975); Glasser PoP 23, 072505 (2016).
+
+"""
+    resist_eval(sing, equil, intr; eta, rho, gamma, ising=0) -> InnerLayer.GGJParameters
+
+Compute the GGJ resistive inner-layer parameters at the rational surface `sing`.
+Port of RDCON `resist_eval` (resist.f), but with η and ρ supplied by the caller (no
+built-in Spitzer/density model). Returns an `InnerLayer.GGJParameters` carrying the curvature/mass
+coefficients E, F, G, H, K, M, the local time scales τ_A (`taua`), τ_R (`taur`), `v1` (V′ normalized
+by total volume), and the surface index `ising`.
+
+# Required keyword arguments
+  - `eta::Real` — plasma resistivity η at this surface [Ω·m]
+  - `rho::Real` — mass density ρ at this surface [kg/m³]
+  - `gamma::Real` — ratio of specific heats (enters G; physical thermodynamic value, e.g. 5/3)
+
+# Physics
+The six flux-surface averages (Fortran `avg(1:6)`) are
+`⟨B²/|∇ψ|²⟩, ⟨1/|∇ψ|²⟩, ⟨1/B²⟩, ⟨1/(B²|∇ψ|²)⟩, ⟨B²⟩, ⟨|∇ψ|²/B²⟩`, each weighted by `J/V′`
+and integrated over θ ∈ [0,1) — exactly the resist.f integrands. The timescales are
+`τ_A = √(ρ·M·μ₀) / |2π n q₁ χ₁ / V′|` and `τ_R = (⟨B²/|∇ψ|²⟩/⟨B²⟩)·μ₀/η`, with η and ρ entering
+once, explicitly.
+"""
+function resist_eval(sing::SingType, equil::Equilibrium.PlasmaEquilibrium,
+    intr::ForceFreeStatesInternal; eta::Real, rho::Real, gamma::Real, ising::Int=0)
+
+    profiles = equil.profiles
+    psifac = sing.psifac
+    nn = sing.n[1]
+
+    # --- 1D surface quantities at the (off-grid) rational ψ ---
+    hint = Ref(1)
+    twopif = profiles.F_spline(psifac; hint=hint)
+    p = profiles.P_spline(psifac; hint=hint)
+    p1 = profiles.P_deriv(psifac; hint=hint)
+    v1 = profiles.dVdpsi_spline(psifac; hint=hint)
+    v2 = profiles.dVdpsi_deriv(psifac; hint=hint)
+    q = sing.q
+    q1 = sing.q1
+    chi1 = 2π * equil.psio
+
+    # --- θ-integrate the six geometric averages (resist.f) ---
+    nth = length(equil.rzphi_ys)
+    ff = zeros(nth, 6)
+    hint2d = (Ref(1), Ref(1))
+    for itheta in 1:nth
+        theta = equil.rzphi_ys[itheta]
+        f1 = equil.rzphi_rsquared((psifac, theta); hint=hint2d)
+        f2 = equil.rzphi_offset((psifac, theta); hint=hint2d)
+        jac = equil.rzphi_jac((psifac, theta); hint=hint2d)
+        fy1 = equil.rzphi_rsquared((psifac, theta); deriv=DerivOp(0, 1), hint=hint2d)
+        fy2 = equil.rzphi_offset((psifac, theta); deriv=DerivOp(0, 1), hint=hint2d)
+        fy3 = equil.rzphi_nu((psifac, theta); deriv=DerivOp(0, 1), hint=hint2d)
+
+        rfac = sqrt(f1)
+        eta_geo = 2π * (theta + f2)   # geometric poloidal angle (local; not the resistivity η)
+        r = equil.ro + rfac * cos(eta_geo)
+
+        v21 = fy1 / (2.0 * rfac * jac)
+        v22 = (1.0 + fy2) * 2π * rfac / jac
+        v23 = fy3 * r / jac
+        v33 = 2π * r / jac
+        bsq = chi1^2 * (v21^2 + v22^2 + (v23 + q * v33)^2)
+        dpsisq = (2π * r)^2 * (v21^2 + v22^2)
+
+        ff[itheta, 1] = bsq / dpsisq
+        ff[itheta, 2] = 1.0 / dpsisq
+        ff[itheta, 3] = 1.0 / bsq
+        ff[itheta, 4] = 1.0 / (bsq * dpsisq)
+        ff[itheta, 5] = bsq
+        ff[itheta, 6] = dpsisq / bsq
+        @views ff[itheta, :] .*= jac / v1
+    end
+    itp = cubic_interp(equil.rzphi_ys, Series(ff); bc=PeriodicBC())
+    avg = FastInterpolations.integrate(itp)
+
+    # --- curvature terms E, F, H (resist.f) ---
+    E = p1 * v1 / (q1 * chi1^2)^2 * avg[1] * (twopif * q1 * chi1 / avg[5] - v2)
+    F = (p1 * v1 / (q1 * chi1^2))^2 * (avg[1] * avg[3] + (twopif / chi1)^2 * (avg[1] * avg[4] - avg[2]^2))
+    H = twopif * p1 * v1 / (q1 * chi1^3) * (avg[2] - avg[1] / avg[5])
+
+    # --- mass factor M and derived G, K (resist.f) ---
+    M = avg[1] * (avg[6] + (twopif / chi1)^2 * (avg[3] - 1 / avg[5]))
+    G = avg[5] / (M * gamma * p)
+    K = (q1 * chi1^2 / (p1 * v1))^2 * avg[5] / (M * avg[1])
+
+    # --- Alfvén and resistive times from caller-supplied η, ρ (resist.f) ---
+    # τ_A = √(ρ·M·μ₀) / |2π n q₁ χ₁ / V′|   — geometric factor × √ρ.
+    # τ_R = (⟨B²/|∇ψ|²⟩ / ⟨B²⟩) · μ₀ / η     — geometric factor × 1/η.
+    # η, ρ enter once, explicitly (no RMATCH divide-out / re-multiply).
+    mu0 = 4π * 1e-7
+    geom_taua = sqrt(M * mu0) / abs(2π * nn * q1 * chi1 / v1)
+    taua = sqrt(rho) * geom_taua
+    geom_taur = avg[1] / avg[5] * mu0
+    taur = geom_taur / eta
+
+    v1norm = equil.params.volume === nothing ? v1 : v1 / equil.params.volume
+
+    return InnerLayer.GGJParameters(; E=E, F=F, G=G, H=H, K=K, M=M,
+        taua=taua, taur=taur, v1=v1norm, ising=ising)
+end
diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl
index f134b686b..bc8d304d5 100644
--- a/src/ForceFreeStates/Riccati.jl
+++ b/src/ForceFreeStates/Riccati.jl
@@ -315,34 +315,46 @@ function compute_delta_prime_matrix!(
     debug && _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua,
                             Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing)
 
-    if use_S_axis
+if use_S_axis
         uShootR, uShootL, uAxis = _build_S_axis_shooting_propagators(
             propagators, chunks, i_crossings, sing, msing, N,
             T_left_mats, T_right_mats, has_ua, ctrl, equil, ffit, intr, debug)
         debug && _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis,
                                                   S_at_surface_left, T_left_mats,
                                                   ipert_all, has_ua, msing, N)
-        M, nMat, col_edge = _assemble_bvp_S_axis(
-            uShootR, uShootL, uAxis, ipert_all, msing, N, wv, psio)
+        M, nMat, col_edge, edge_drive_rows = _assemble_bvp_S_axis(
+            uShootR, uShootL, uAxis, ipert_all, msing, N, wv, psio;
+            rpec = (ctrl !== nothing && ctrl.gal_rpec_flag))
     else
         M, nMat, col_edge = _assemble_bvp_FM_axis(
             Phi_L_mats, Phi_R_mats, ipert_all, msing, N,
             T_left_inv, T_right_inv, has_ua, wv, psio)
+        edge_drive_rows = nothing
     end
-
     if debug
         @info "Δ' BVP: nMat=$nMat, rank(M)=$(rank(M)), cond(M)=$(@sprintf("%.2e", cond(M)))"
     end
 
-    intr.delta_prime_matrix = _solve_bvp_and_combine_pest3(
-        M, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug)
+    # Coil solve only on the S-axis path (only path that returns edge_drive_rows) and
+    # only when the coil gate is on.
+    do_coil = ctrl !== nothing && ctrl.gal_rpec_flag && use_S_axis
+
+    intr.delta_prime_matrix, intr.delta_raw_matrix, intr.delta_coil_matrix =
+        _solve_bvp_and_combine_pest3(
+            M, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug;
+            compute_coil = do_coil, edge_drive_rows = edge_drive_rows)
 end
 
 # Column index helpers for the BVP matrix. j is the 1-based singular-surface index,
 # N is numpert_total. Layout: c_axis(N), c_left[1](2N), c_right[1](2N), ..., c_edge(N).
+
 _col_left(j::Int, N::Int)  = (N + 4N*(j-1) + 1):(N + 4N*(j-1) + 2N)
 _col_right(j::Int, N::Int) = (N + 4N*(j-1) + 2N + 1):(N + 4N*j)
 
+#N = number of poloidal modes = N_coil in RDCON (ViaW)
+#j = singular surface index (1-based) (ViaW)
+#its figuring out where in M is the reg. Delta matrix (ViaW)
+
 # Multi-resonance surfaces (one q value satisfying multiple (m,n) tuples in a multi-n run)
 # are not yet handled by the inter-surface BVP. Returns true if any surface has >1 modes;
 # emits a warning as a side effect. The stub per-surface delta_prime is unaffected.
@@ -570,12 +582,13 @@ end
 # Assemble the BVP matrix M with S-based axis BC. The Riccati S matrix at surface 1's left
 # boundary encodes the axis BC (U₁ = S·U₂) in a well-conditioned form (cond ~ 10⁶), avoiding
 # the catastrophically ill-conditioned axis FM. Fortran-matched structure with
-# nMat = (2 + 4·msing)·N. Returns (M, nMat, col_edge).
+# nMat = (2 + 4·msing)·N. Returns (M, nMat, col_edge). (ViaW edit: THIS section is boundary conditions, through 630)
 function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
                               uShootL::Vector{Matrix{ComplexF64}},
                               uAxis::Matrix{ComplexF64}, ipert_all::Vector{Int},
                               msing::Int, N::Int,
-                              wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64)
+                              wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64;
+                              rpec::Bool=false)
     # STRIDE global BVP block structure [Glasser-Kolemen 2018 PoP 25, 032501 Eq. 37].
     nMat = (2 + 4 * msing) * N
     col_axis = 1:N
@@ -583,10 +596,15 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
     M = zeros(ComplexF64, nMat, nMat)
 
     # Axis matching: uShootL[1] · c_left[1] = uAxis · c_axis  (2N equations)
-    M[1:2N, _col_left(1, N)] .= uShootL[1]
+    M[1:2N, _col_left(1, N)] .= uShootL[1]#(ViaW: axis condition where q(psi_axis)=0)
     M[1:2N, col_axis]        .= -uAxis
     row_offset = 2N
 
+    # ViaW code for testing boundary conditions
+    #M[1:2N, col_edge] .= 0.0
+    #M[N+1:2N, col_axis] .= -I(N) #U_2 = I, U_1 = 0 -> for conditions 0 and 1
+
+    edge_drive_rows = Int[]    
     for j in 1:msing
         ipert_j = ipert_all[j]
         # Crossing: non-resonant modes continuity (asymptotic basis = identity)
@@ -598,15 +616,26 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
             end
         end
 
-        junc_rows = (row_offset + 1):(row_offset + 2N)
+        junc_rows = (row_offset + 1):(row_offset + 2N) #(Via edit: these are your boundary conditions)
+        @info "junc_rows = $junc_rows, j=$j, msing=$msing, N=$N"
         if j < msing
             # Midpoint matching between consecutive surfaces
             M[junc_rows, _col_right(j, N)]   .= -uShootR[j]
             M[junc_rows, _col_left(j+1, N)]  .=  uShootL[j+1]
         else
-            # Edge junction
+            # Edge junction (ViaW: edge condition where q(psi_edge)=0)
             M[junc_rows, _col_right(msing, N)] .= uShootR[msing]
-            if wv !== nothing
+            if rpec
+                # Coil/rpec drive: prescribe the edge displacement u_edge as a Dirichlet unit source
+                # on the value rows (junc_rows[1:N] carry uShootR·c_right|_value, so this couples to
+                # the interior); record u'_edge in the c_edge slack via the derivative rows.
+                # Homogeneously u_edge = 0 (fixed edge, matching the Galerkin identity-edge
+                # convention); each coil column drives one edge poloidal mode to 1 → resonant
+                # response is the delta_coil block. Driving the value rows here (rather than the
+                # -I(N) coupling that leaves c_edge a decoupled slack) is what makes delta_coil nonzero.
+                M[junc_rows[N+1:end], col_edge] .= -I(N)
+                edge_drive_rows = collect(junc_rows[1:N])
+            elseif wv !== nothing
                 M[junc_rows[1:N],     col_edge] .= -I(N)
                 M[junc_rows[N+1:end], col_edge] .= wv .* psio^2
             else
@@ -620,12 +649,13 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
     for j in 1:msing
         ipert_j = ipert_all[j]
         row_offset += 1
-        M[row_offset, _col_left(j, N)[ipert_j]]  = 1
+        M[row_offset, _col_left(j, N)[ipert_j]]  = 1 #(ViaW edit these)
         row_offset += 1
-        M[row_offset, _col_right(j, N)[ipert_j]] = 1
+        M[row_offset, _col_right(j, N)[ipert_j]] = 1 #(ViaW edit these)
     end
     @assert row_offset == nMat "Row count mismatch: expected $nMat, got $row_offset"
-    return M, nMat, col_edge
+    # edge_drive_rows = collect((row_offset - 2N + 1):(row_offset - 2N + N))
+    return M, nMat, col_edge, edge_drive_rows
 end
 
 # Fallback BVP assembly with FM-based axis BC (used when no Riccati S matrices are available).
@@ -636,6 +666,7 @@ function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}},
                                T_left_inv::Vector{Matrix{ComplexF64}},
                                T_right_inv::Vector{Matrix{ComplexF64}}, has_ua::Bool,
                                wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64)
+    @warn "compute_delta_prime_matrix!: FM-axis fallback"
     nMat = (2 + 4 * msing) * N
     col_axis = 1:N
     col_edge = (N + 4N*msing + 1):nMat
@@ -650,7 +681,7 @@ function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}},
         cl = _col_left(j, N)
         cr = _col_right(j, N)
         row_cont = 2N + (4N-2)*(j-1)
-        for i in 1:2N
+        for i in 1:2N #(ViaW: non resonant continuity conditions)
             if i != ipert_j && i != ipert_j + N
                 row_cont += 1
                 M[row_cont, cl[i]] =  1
@@ -678,7 +709,7 @@ function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}},
             M[row_drive_base + 2j,   cr[ipert_j]] = 1
         end
     end
-    return M, nMat, col_edge
+    return M, nMat, col_edge, edge_drive_rows = Int[]
 end
 
 # Solve the BVP for each driving configuration and apply the PEST3 four-term combination.
@@ -687,19 +718,20 @@ end
 # precision lets the imaginary part drift 2–5× on DIIID-class equilibria.
 function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, nMat::Int,
                                       use_S_axis::Bool, ipert_all::Vector{Int}, col_edge,
-                                      ctrl, debug::Bool)
-    s2 = 2 * msing
+                                      ctrl, debug::Bool;
+                                      compute_coil::Bool = false,
+                                      edge_drive_rows=nothing)
+    s2 = 2 * msing #does N=1 in the edge case, but not in the 2N=4 case? (ViaW)
     Tc = (ctrl === nothing || ctrl.extended_precision_bvp) ? Complex{Double64} : ComplexF64
     M_solve = Tc.(M)
-
     M_lu = lu(M_solve; check=false)
     use_lu = issuccess(M_lu)
     M_pinv = use_lu ? nothing : pinv(M_solve)
     if !use_lu
         @warn "Δ' BVP: LU factorization singular (rank $(rank(M))/$nMat), using pseudo-inverse fallback"
     end
-
     dp_raw = zeros(Tc, s2, s2)
+    #one of these 2msing will need a +N (ViaW)
     b = zeros(Tc, nMat)
     for jsing in 1:msing, side in 1:2
         dRow = 2jsing - (2 - side)
@@ -718,6 +750,19 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int,
         end
     end
 
+    delta_coil = zeros(Tc, s2, N)
+    if compute_coil && !isempty(edge_drive_rows)
+        for kmode in 1:N
+            fill!(b,0)
+            b[edge_drive_rows[kmode]] = 1
+            x = use_lu ? (M_lu \ b) : (M_pinv * b)
+            for ksing in 1:msing
+                ipert_k = ipert_all[ksing]
+                delta_coil[2ksing-1, kmode] = x[_col_left(ksing, N)[ipert_k+N]]
+                delta_coil[2ksing,   kmode] = x[_col_right(ksing, N)[ipert_k+N]]
+            end
+        end    
+    end
     # PEST3 four-term combination [Chance PPPL-2527; Glasser-Kolemen 2018 PoP 25, 032501 Eq. 31].
     # Δ'[i,j] = (NW − NE − SW + SE) on each 2×2 block of dp_raw, in extended precision.
     deltap_ext = zeros(Tc, msing, msing)
@@ -727,7 +772,7 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int,
     deltap = ComplexF64.(deltap_ext)
 
     debug && _log_bvp_pest3(dp_raw, deltap, s2, msing, Tc)
-    return deltap
+    return deltap, ComplexF64.(dp_raw), ComplexF64.(delta_coil)
 end
 
 # Logging helpers for `compute_delta_prime_matrix!`. Called only when debug=true.
diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl
index 16661c228..3ecd63628 100644
--- a/src/ForceFreeStates/Sing.jl
+++ b/src/ForceFreeStates/Sing.jl
@@ -117,6 +117,42 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     end
 end
 
+"""
+    sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
+
+Set the lower integration bound `intr.psilow`. Port of Fortran RDCON `sing_min` (sing.f):
+when `qlow > qmin`, the q < qlow core (including any q ≤ 1 sawtooth/internal-kink surfaces) must be
+excluded from the outer-region Galerkin domain — otherwise the Hermite FEM integrates through those
+ideal singularities without imposing the ideal constraint, contaminating Δ′ at the innermost kept
+surface. A Newton iteration locates ψ where q = qlow; scanning starts from the edge inward for
+robustness in reverse-shear cores. When `qlow ≤ qmin` the axis value is kept.
+"""
+function sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
+    profiles = equil.profiles
+    intr.psilow = profiles.xs[1]   # default: equilibrium axis-side bound
+    ctrl.qlow > equil.params.qmin || return intr.psilow
+
+    # Scan from the edge inward for the first node with q < qlow (robust for reverse-shear q).
+    qy = profiles.q_spline.y
+    jpsi = 1
+    for j in (length(profiles.xs)-1):-1:1
+        if qy[j] < ctrl.qlow
+            jpsi = j
+            break
+        end
+    end
+
+    hint = Ref(jpsi)
+    intr.psilow = find_zero(
+        (psi -> profiles.q_spline(psi; hint=hint) - ctrl.qlow,
+            psi -> profiles.q_deriv(psi; hint=hint)),
+        profiles.xs[jpsi], Roots.Newton()
+    )
+    @info "sing_min: qlow=$(@sprintf("%.3f", ctrl.qlow)) > qmin=$(@sprintf("%.3f", equil.params.qmin)); " *
+          "raising psilow from $(@sprintf("%.5f", profiles.xs[1])) to $(@sprintf("%.5f", intr.psilow)) (excludes q<qlow core)"
+    return intr.psilow
+end
+
 """
     compute_sing_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
@@ -694,6 +730,105 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64)
     return ua
 end
 
+"""
+    sing_get_dua(sing_asymp::SingAsymptotics, dpsi::Float64) -> dua
+
+Compute the derivative of the asymptotic series solution with respect to the positive distance
+`dpsi = |ψ − ψ_res|`, consistent with `sing_get_ua`. Port of Fortran `sing_get_dua`
+(sing.f / sing1.f). Shape `(numpert_total, 2*numpert_total, 2)`. Used by the
+outer-region Galerkin solver (`gal_extension`, `sing_matvec`).
+
+Direction is carried by the `SingAsymptotics` (left vs right vmat built with sig=∓1), exactly as in
+`sing_get_ua`, so `dpsi` is always positive and the arithmetic is real (no analytic continuation).
+The term-by-term derivative is built via a `power` array tracking each term's total exponent (in
+half-powers of `dpsi`: the 2*sing_order offset, the ∓1 shearing on resonant rows, and the ∓2α on
+resonant columns), then the same shearing/`pfac` restoration and the chain-rule factor `1/(2·dpsi)`
+are applied. The physical d/dψ sign for the left side (ψ = ψ_res − dpsi) is applied by the caller.
+
+### Arguments
+
+  - `sing_asymp::SingAsymptotics`: Pre-computed asymptotic data (must be left- or right-specific)
+  - `dpsi::Float64`: Positive distance from singular surface = |ψ − ψ_res|
+"""
+function sing_get_dua(sing_asymp::SingAsymptotics, dpsi::Float64)
+
+    r1 = sing_asymp.r1
+    r2 = sing_asymp.r2
+    order = sing_asymp.sing_order
+    sqrtfac = sqrt(dpsi)
+    N = size(sing_asymp.vmat, 1)
+
+    # Per-term total exponent (×2, in half-powers of dpsi). See Fortran sing_get_dua.
+    power = fill(ComplexF64(2 * order), N, 2 * N, 2)
+    power[r1, :, 1] .-= 1
+    power[r1, :, 2] .+= 1
+    for i in eachindex(r1)
+        power[:, r2[2*i-1], :] .-= 2 * sing_asymp.alpha[i]
+        power[:, r2[2*i], :] .+= 2 * sing_asymp.alpha[i]
+    end
+
+    # Power series derivative by Horner's method (Fortran sing_get_dua)
+    dua = sing_asymp.vmat[:, :, :, 2*order+1] .* power
+    for iorder in (2*order-1):-1:0
+        power .-= 1
+        dua .= dua .* sqrtfac .+ sing_asymp.vmat[:, :, :, iorder+1] .* power
+    end
+
+    # Restore shearing and resonant powers, then the chain-rule factor 1/(2·dpsi) (Fortran sing_get_dua)
+    for i in eachindex(r1)
+        pfac = dpsi^sing_asymp.alpha[i]
+        dua[:, r2[2*i-1], :] ./= pfac
+        dua[:, r2[2*i], :] .*= pfac
+        dua[r1[i], :, 1] ./= sqrtfac
+        dua[r1[i], :, 2] .*= sqrtfac
+    end
+    dua ./= (2 * dpsi)
+
+    return dua
+end
+
+"""
+    sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Float64, q::Float64, ua, dua) -> matvec
+
+Apply the Euler-Lagrange residual operator `L u = -(F u' + K u)' + (K† u' + G u)` to the asymptotic
+solutions. Port of Fortran `sing_matvec` (sing.f). Returns `matvec`, shape
+`(numpert_total, size(ua,2))`.
+
+Uses the reduced (Schur-complemented) `ffit.kmats` (= K̄) and `ffit.gmats` (= Ḡ) directly, with the
+**direct** singular factor `singfac = m - n q` applied to `u' = dua[:,:,1]`. The second component
+`ua[:,:,2]` is the canonical momentum `F u' + K u`, so `-dua[:,:,2] = -(F u' + K u)'`.
+
+### Arguments
+
+  - `psi`: flux coordinate; `q`: safety factor at `psi` (passed in to avoid re-evaluating the spline)
+  - `ua`, `dua`: asymptotic solution and its ψ-derivative from `sing_get_ua`/`sing_get_dua`
+"""
+function sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Float64, q::Float64,
+    ua::Array{ComplexF64,3}, dua::Array{ComplexF64,3})
+
+    N = intr.numpert_total
+    msol = size(ua, 2)
+
+    # Direct singular factor (m - n q), NOT 1/(m - n q)
+    sfvec = vec((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)')
+
+    kmat = Matrix{ComplexF64}(undef, N, N)
+    gmat = Matrix{ComplexF64}(undef, N, N)
+    ffit.kmats(vec(kmat), psi; hint=ffit._hint)
+    ffit.gmats(vec(gmat), psi; hint=ffit._hint)
+    kdag = adjoint(kmat)
+
+    matvec = zeros(ComplexF64, N, msol)
+    d1 = Vector{ComplexF64}(undef, N)
+    for isol in 1:msol
+        @views d1 .= dua[:, isol, 1] .* sfvec          # u' · singfac
+        @views matvec[:, isol] .= kdag * d1            # K̄† (u' · singfac)
+        @views matvec[:, isol] .+= gmat * ua[:, isol, 1]  # + Ḡ u
+        @views matvec[:, isol] .-= dua[:, isol, 2]     # - (F u' + K u)'
+    end
+    return matvec
+end
+
 """
     sing_get_ca(u::Array{ComplexF64,3}, ua::Array{ComplexF64,3}, intr::ForceFreeStatesInternal)
 
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index b6988f0c9..0e179883a 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -13,14 +13,14 @@ include("Vacuum/Vacuum.jl")
 import .Vacuum as Vacuum
 export Vacuum
 
-include("ForceFreeStates/ForceFreeStates.jl")
-import .ForceFreeStates as ForceFreeStates
-export ForceFreeStates
-
 include("InnerLayer/InnerLayer.jl")
 import .InnerLayer as InnerLayer
 export InnerLayer
 
+include("ForceFreeStates/ForceFreeStates.jl")
+import .ForceFreeStates as ForceFreeStates
+export ForceFreeStates
+
 include("ForcingTerms/ForcingTerms.jl")
 import .ForcingTerms as ForcingTerms
 export ForcingTerms
@@ -49,11 +49,12 @@ import AdaptiveArrayPools: @with_pool
 
 # Import ForceFreeStates types and functions needed for main
 using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, VacuumData, OdeState, FourFitVars
-using .ForceFreeStates: sing_lim!, sing_find!
+using .ForceFreeStates: sing_lim!, sing_min!, sing_find!
 using .ForceFreeStates: mercier_scan!, compute_ballooning_stability!
 using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix
 using .ForceFreeStates: find_kinetic_singular_surfaces!
 using .ForceFreeStates: eulerlagrange_integration, free_run!
+using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate
 
 const _BANNER = "="^60
 const _SECTION = "-"^40
@@ -221,6 +222,14 @@ function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothin
         end
     end
 
+    # For the outer-region Galerkin solve, exclude the q < qlow core (incl. any q≤1 sawtooth
+    # surfaces) by raising psilow to where q = qlow (RDCON sing_min). Without this, the gal FEM
+    # integrates through the unhandled q≤1 ideal singularity and contaminates Δ′ at the innermost
+    # kept surface when q0 < qlow. No-op (keeps the axis bound) when qlow ≤ qmin.
+    if ctrl.gal_flag
+        sing_min!(intr, ctrl, equil)
+    end
+
     # Determine poloidal mode numbers
     if ctrl.delta_mlow < 0 || ctrl.delta_mhigh < 0
         error("Negative delta_mlow or delta_mhigh not allowed")
@@ -364,8 +373,15 @@ function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothin
         end
     end
 
+    # Outer-region resistive Δ′ matrix via the singular Galerkin method (RDCON gal_solve)
+    gal_data = nothing
+    if ctrl.gal_flag
+        ctrl.mat_flag || error("gal_flag=true requires mat_flag=true (needs the F/G/K matrix splines)")
+        gal_data = galerkin_solve(ctrl, equil, ffit, intr; vac_data=ctrl.vac_flag ? vac_data : nothing)
+    end
+
     if ctrl.write_outputs_to_HDF5
-        write_outputs_to_HDF5(ctrl, equil, intr, odet, ctrl.vac_flag ? vac_data : nothing, ffit, git_version)
+        write_outputs_to_HDF5(ctrl, equil, intr, odet, ctrl.vac_flag ? vac_data : nothing, ffit, git_version, gal_data)
         @info "Results written to $(ctrl.HDF5_filename)"
     end
 
@@ -405,10 +421,20 @@ function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothin
         )
         pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=intr.dir_path)
 
+        # DRIVEN (RPEC): feed the coil-matched gal solution to PE instead of the shooting solution.
+        # The matched OdeState is in the identity-at-edge basis; build_flux_matrix rederives the edge BC
+        # from u_store[:,:,1,step], so PE consumes it unchanged. The shooting odet is left untouched for
+        # the Force-Free States HDF5 output.
+        pe_odet = odet
+        if ctrl.gal_flag && ctrl.gal_match_flag && gal_data !== nothing && gal_data.match !== nothing
+            @info "PerturbedEquilibrium: using the RPEC-matched gal solution"
+            pe_odet = gal_matched_odestate(gal_data, ffit, intr)
+        end
+
         # Run perturbed equilibrium calculations
         # Pass vac_data and intr for response matrix calculations
         pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(
-            equil, odet, ctrl.vac_flag ? vac_data : nothing, intr, ft_ctrl, pe_ctrl, pe_intr,
+            equil, pe_odet, ctrl.vac_flag ? vac_data : nothing, intr, ft_ctrl, pe_ctrl, pe_intr,
             metric, ffit
         )
 
@@ -480,7 +506,8 @@ function write_outputs_to_HDF5(
     odet::OdeState,
     vac_data::Union{VacuumData,Nothing},
     ffit::Union{FourFitVars,Nothing}=nothing,
-    git_version::String="unknown"
+    git_version::String="unknown",
+    gal_data::Union{GalerkinResult,Nothing}=nothing
 )
 
     h5open(joinpath(intr.dir_path, ctrl.HDF5_filename), "w") do out_h5
@@ -488,6 +515,11 @@ function write_outputs_to_HDF5(
         # Store git version for reproducibility
         out_h5["info/git_version"] = git_version
 
+        # Outer-region Galerkin Δ′ matrix (RDCON), if computed
+        if gal_data !== nothing
+            write_galerkin!(out_h5, gal_data)
+        end
+
         # Store input parameters
         for (key, val) in zip(fieldnames(ForceFreeStatesControl), getfield.(Ref(ctrl), fieldnames(ForceFreeStatesControl)))
             out_h5["input/ForceFreeStates/$key"] = val
@@ -613,6 +645,13 @@ function write_outputs_to_HDF5(
         # Shape: [msing × msing] — PEST3-convention deltap (STRIDE BVP with vacuum coupling).
         if intr.msing > 0 && !isempty(intr.delta_prime_matrix)
             out_h5["singular/delta_prime_matrix"] = intr.delta_prime_matrix
+            out_h5["singular/delta_prime_raw"] = intr.delta_raw_matrix #(ViaW edit)
+        end
+        if intr.msing > 0 && !isempty(intr.delta_coil_matrix)
+            out_h5["singular/delta_coil_matrix"] = intr.delta_coil_matrix
+            out_h5["singular/delta_coil_abs"] = abs.(intr.delta_coil_matrix)
+            out_h5["singular/delta_coil_real"] = real.(intr.delta_coil_matrix)
+            out_h5["singular/delta_coil_imag"] = imag.(intr.delta_coil_matrix) #(ViaW edit)
         end
 
         # Write kinetic singular surface data (det(F̄) near-zeros) and the cond(F̄) scan
diff --git a/src/InnerLayer/GGJ/Galerkin.jl b/src/InnerLayer/GGJ/Galerkin.jl
index 93f889018..43b4dbbaf 100644
--- a/src/InnerLayer/GGJ/Galerkin.jl
+++ b/src/InnerLayer/GGJ/Galerkin.jl
@@ -708,7 +708,7 @@ function solve_inner(::GGJModel{:galerkin}, params::GGJParameters, γ::Number;
     emap1 = res_cell.emap[1]
     Δ_raw = SVector{2,ComplexF64}(ws.sol[emap1, 1], ws.sol[emap1, 2])
 
-    # Apply deltac.f's swap convention (line 194-196)
+    # Apply deltac.f's swap convention (deltac_solve)
     Δ_swapped = SVector{2,ComplexF64}(Δ_raw[2], Δ_raw[1])
 
     return rescale_delta(Δ_swapped, params)