diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 961a0f7f5..586d81001 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -227,6 +227,32 @@ function direct_position!(raw_profile::DirectRunInput)
     return ro, zo, rs1, rs2
 end
 
+"""
+    eta_at_sfl_angle(sol, y_out, x, total_x) -> Float64
+
+Integration angle η at which the normalised straight-fieldline angle ∫jac·dl/Bp reaches `x`.
+
+`y_out[:, 5]` is monotone in η, so it brackets the root to one solver step and Brent converges in
+a handful of dense-output evaluations. Used to sample every flux surface at the *same* SFL angles
+instead of resampling each surface's own solver steps (issue #376).
+"""
+function eta_at_sfl_angle(sol, y_out::Matrix{Float64}, x::Float64, total_x::Float64)
+    x <= 0 && return y_out[1, 1]
+    x >= 1 && return y_out[end, 1]
+    target = x * total_x
+    hi = searchsortedfirst(view(y_out, :, 5), target)
+    hi = clamp(hi, 2, size(y_out, 1))
+    lo = hi - 1
+    eta_lo, eta_hi = y_out[lo, 1], y_out[hi, 1]
+    f(eta) = sol(eta)[4] - target
+    flo, fhi = f(eta_lo), f(eta_hi)
+    # Degenerate bracket (repeated η, or the root sitting exactly on a step) needs no solve.
+    flo == 0 && return eta_lo
+    fhi == 0 && return eta_hi
+    (flo * fhi > 0 || eta_hi <= eta_lo) && return eta_lo + (eta_hi - eta_lo) * (target - y_out[lo, 5]) / max(y_out[hi, 5] - y_out[lo, 5], eps())
+    return find_zero(f, (eta_lo, eta_hi), Roots.Brent())
+end
+
 """
     direct_fieldline_int(psifac, raw_profile, ro, zo, rs2)
 
@@ -250,9 +276,12 @@ from 1:5 rather than 0:4 as in Fortran.
         - `y_out[:, 4]`: ∫(dl/(R²Bp))
         - `y_out[:, 5]`: ∫(jac*dl/Bp)
 
+  - `sol`: the dense ODE solution, so callers can evaluate the trace at prescribed SFL angles
+    rather than resampling this surface's own solver steps (`nothing` for tracers without it).
+
   - `bfield`: A `DirectBField` object with values at the integration start point.
 """
-function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64)::Tuple{Matrix{Float64},DirectBField}
+function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64)
 
     # Find the starting point on the flux surface (outboard midplane)
     psi0_guess = raw_profile.psio * (1.0 - psifac)
@@ -291,10 +320,12 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
     callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false))
 
     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
-    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=false)
+    # Dense output lets the caller evaluate the trace at the SFL angles it actually wants, instead
+    # of splining this surface's solver-chosen steps and resampling (issue #376).
+    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=true)
 
     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
-    return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
+    return hcat(sol.t::Vector{Float64}, sol_matrix), bfield, sol
 end
 
 """
@@ -495,18 +526,39 @@ robustness.
     ff_deriv_val = zeros!(pool, Float64, 4)
 
     for ipsi in (mpsi+1):-1:1  # outermost to innermost
-        y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2)
+        y_out, bfield, sol = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2)
         checkpoint!(pool, Float64)
 
-        # Fit data into temporary straight fieldline poloidal angle splines
-        ff_x_nodes = acquire!(pool, Float64, size(y_out, 1))
-        @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5]
-
-        ff_fs_nodes = acquire!(pool, Float64, size(y_out, 1), 4)
-        @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2
-        @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes
-        @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4])
-        @. ff_fs_nodes[:, 4] = @view(y_out[:, 2]) / y_out[end, 2] - ff_x_nodes
+        # Straight-fieldline angle x = normalised ∫jac·dl/Bp, monotone in the integration angle η.
+        #
+        # Sampling x at this surface's own solver steps and resampling onto theta_nodes leaves a
+        # resample error that is uncorrelated between neighbouring surfaces, i.e. white noise in ψ
+        # that grid refinement then amplifies (issue #376). With dense output we instead solve for
+        # the η where x hits each target node and evaluate there, so every surface is sampled at
+        # the same abscissae and the resample error at the output nodes is zero.
+        nff = sol === nothing ? size(y_out, 1) : mtheta + 1
+        ff_x_nodes = acquire!(pool, Float64, nff)
+        ff_fs_nodes = acquire!(pool, Float64, nff, 4)
+
+        if sol === nothing
+            @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5]
+            @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2
+            @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes
+            @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4])
+            @. ff_fs_nodes[:, 4] = @view(y_out[:, 2]) / y_out[end, 2] - ff_x_nodes
+        else
+            total_x = y_out[end, 5]
+            for itheta in 1:(mtheta+1)
+                x = theta_nodes[itheta]
+                eta = eta_at_sfl_angle(sol, y_out, x, total_x)
+                u = sol(eta)
+                ff_x_nodes[itheta] = x
+                ff_fs_nodes[itheta, 1] = u[2]^2
+                ff_fs_nodes[itheta, 2] = eta / (2π) - x
+                ff_fs_nodes[itheta, 3] = bfield.f * (u[3] - x * y_out[end, 4])
+                ff_fs_nodes[itheta, 4] = u[1] / y_out[end, 2] - x
+            end
+        end
 
         ff_fs_nodes[end, :] .= ff_fs_nodes[1, :]  # enforce periodic endpoint
 
diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
index 32b93a0ac..ec4c661d3 100644
--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
+++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
@@ -75,7 +75,7 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
 """
 @with_pool pool function arclength_fieldline_int(
     psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64
-)::Tuple{Matrix{Float64},DirectBField}
+)::Tuple{Matrix{Float64},DirectBField,Nothing}
 
     psi0_guess = raw_profile.psio * (1.0 - psifac)
     r = ro + sqrt(psifac) * (rs2 - ro)
@@ -142,6 +142,6 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     end
 
     # bfield at the starting point carries F and P for the surface-averaged quantities
-    return y_out, bfield
+    return y_out, bfield, nothing
 end
 
diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
index f5d188abe..0f8095e3c 100644
--- a/src/ForceFreeStates/Fourfit.jl
+++ b/src/ForceFreeStates/Fourfit.jl
@@ -311,7 +311,7 @@ function build_kinetic_metric_matrices(equil::Equilibrium.PlasmaEquilibrium,
             v[3, 3] = 2π * rs
 
             # Raw g^ij (Fortran dcon_interface.f:933-937)
-            g12 = v[1, 1]*v[2, 1] + v[1, 2]*v[2, 2] + v[1, 3]*v[2, 3]
+            g12 = v[1, 1] * v[2, 1] + v[1, 2] * v[2, 2] + v[1, 3] * v[2, 3]
             g13 = v[3, 3] * v[1, 3]
             g22 = v[2, 1]^2 + v[2, 2]^2 + v[2, 3]^2
             g23 = v[2, 3] * v[3, 3]
@@ -428,6 +428,33 @@ function build_kinetic_metric_matrices(equil::Equilibrium.PlasmaEquilibrium,
 end
 
 
+"""
+    core_capped_knots(xs, rationals) -> Vector{Int}
+
+Indices of `xs` the EL coefficient splines keep, decoupling their knots from the equilibrium grid
+in the packed core. A cubic spline's third-derivative jumps scale as (node error)/Δψ³, so
+equilibrium-grade core packing amplifies tolerance-level node error into jumps that slave the
+Euler-Lagrange step size. Near the axis every component is a Frobenius power law in ψ, and power
+laws are scale-free, so log-uniform sampling (Δψ ≥ 0.05·ψ) resolves them at constant relative
+accuracy — cubic interpolation of ψ^p on that grid errs by ~(0.05p)⁴/384, resolving even the
+steepest component (p = m_max/2) to ~2e-4, well below where the physics responds. The capped region
+ends at the innermost rational (or ψ = 0.1, whichever is smaller) and never drops a knot inside a
+rational's resolution window, preserving the Δ′-stencil structure the equilibrium grid encodes
+(`Equilibrium.RATIONAL_RES_RADIUS`).
+"""
+function core_capped_knots(xs::Vector{Float64}, rationals::Vector{Float64})::Vector{Int}
+    cap_edge = 0.1
+    isempty(rationals) || (cap_edge = min(cap_edge, minimum(rationals) - Equilibrium.RATIONAL_RES_RADIUS))
+    in_rational_window(x) = any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rationals)
+    keep = Int[1]
+    for i in 2:(length(xs)-1)
+        x = xs[i]
+        (x >= cap_edge || in_rational_window(x) || (x - xs[keep[end]]) >= 0.05 * x) && push!(keep, i)
+    end
+    push!(keep, length(xs))
+    return keep
+end
+
 """
     build_matrix_splines(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData) -> MatrixSplines
 
@@ -607,19 +634,23 @@ function build_matrix_splines(equil::Equilibrium.PlasmaEquilibrium, intr::ForceF
     # Create complex series interpolants with per-column extrap BC
     # TODO: set powers. Do we need this yet? Only called if power_flag = true
     itp_opts = (; extrap=ExtendExtrap())
+    keep = core_capped_knots(metric.xs, [s.psifac for s in intr.sing])
+    mxs = metric.xs[keep]
+    length(mxs) < length(metric.xs) &&
+        @info "EL coefficient-spline grid: $(length(metric.xs)) -> $(length(mxs)) knots after core density cap"
     ideal = IdealMatrices(;
-        A_spline=cubic_interp(metric.xs, Series(amats_flat); itp_opts...),
-        B_spline=cubic_interp(metric.xs, Series(bmats_flat); itp_opts...),
-        C_spline=cubic_interp(metric.xs, Series(cmats_flat); itp_opts...),
-        D_spline_prim=cubic_interp(metric.xs, Series(dmats_flat); itp_opts...),
-        E_spline_prim=cubic_interp(metric.xs, Series(emats_flat); itp_opts...),
-        H_spline=cubic_interp(metric.xs, Series(hmats_flat); itp_opts...),
-        F_spline_lower=cubic_interp(metric.xs, Series(fmats_lower_flat); itp_opts...),
-        F_spline_prim=cubic_interp(metric.xs, Series(fmats_prim_flat); itp_opts...),
-        F_spline_gal=cubic_interp(metric.xs, Series(fmats_gal_flat); itp_opts...),
-        G_spline=cubic_interp(metric.xs, Series(gmats_flat); itp_opts...),
-        K_spline=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...),
+        A_spline=cubic_interp(mxs, Series(amats_flat[keep, :]); itp_opts...),
+        B_spline=cubic_interp(mxs, Series(bmats_flat[keep, :]); itp_opts...),
+        C_spline=cubic_interp(mxs, Series(cmats_flat[keep, :]); itp_opts...),
+        D_spline_prim=cubic_interp(mxs, Series(dmats_flat[keep, :]); itp_opts...),
+        E_spline_prim=cubic_interp(mxs, Series(emats_flat[keep, :]); itp_opts...),
+        H_spline=cubic_interp(mxs, Series(hmats_flat[keep, :]); itp_opts...),
+        F_spline_lower=cubic_interp(mxs, Series(fmats_lower_flat[keep, :]); itp_opts...),
+        F_spline_prim=cubic_interp(mxs, Series(fmats_prim_flat[keep, :]); itp_opts...),
+        F_spline_gal=cubic_interp(mxs, Series(fmats_gal_flat[keep, :]); itp_opts...),
+        G_spline=cubic_interp(mxs, Series(gmats_flat[keep, :]); itp_opts...),
+        K_spline=cubic_interp(mxs, Series(kmats_flat[keep, :]); itp_opts...),
         # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl
-        J_spline=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...))
+        J_spline=cubic_interp(mxs, Series(jmats_flat[keep, :]); itp_opts...))
     return MatrixSplines(; ideal)
 end