diff --git a/Project.toml b/Project.toml
index d3b68952..126ec3f5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ authors = ["Orso Meneghini <orso82@gmail.com>"]
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+AdaptiveArrayPools = "4f381ef7-9af0-4cbe-99d4-cf36d7b0f233"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 CoordinateConventions = "7204ce3a-f536-43d2-be4a-fbed74e90d86"
@@ -43,8 +44,15 @@ SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[weakdeps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+
+[extensions]
+IMASInterpolationsExt = "Interpolations"
+
 [compat]
 AbstractTrees = "0.4"
+AdaptiveArrayPools = "0.3.6"
 Compat = "4.10"
 Contour = "0.6"
 CoordinateConventions = "1"
@@ -58,6 +66,7 @@ Graphs = "1"
 HelpPlots = "1.1"
 IMASdd = "8.3"
 IMASutils = "1.5.1"
+Interpolations = "0.13, 0.14, 0.15, 0.16"
 Jedis = "0.3"
 LaTeXStrings = "1"
 Measurements = "2"
@@ -76,7 +85,7 @@ Random = "1.11.0"
 Ratios = "0.4"
 RecipesBase = "1"
 Roots = "2, 3"
-SimpleNonlinearSolve = "2.12.0"
+SimpleNonlinearSolve = "2"
 StaticArrays = "1"
 Statistics = "1.11.1"
 julia = "1"
diff --git a/claudedocs/cubic_tracer_benchmark.jl b/claudedocs/cubic_tracer_benchmark.jl
new file mode 100644
index 00000000..2012daa0
--- /dev/null
+++ b/claudedocs/cubic_tracer_benchmark.jl
@@ -0,0 +1,23 @@
+# PC vs adaptive-RK4 vs Contour.jl: drift, point count, timing on DIII-D closed surfaces.
+# Scratch diagnostic (not a pass/fail test).
+import IMAS
+using Printf
+
+filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+dd = IMAS.json2imas(filename; show_warnings=false)
+eqt = dd.equilibrium.time_slice[1]
+eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+r, z, itp = IMAS.ψ_interpolant(eqt2d)
+RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+psi_axis = itp(RA, ZA)
+eqt1d = eqt.profiles_1d
+# a mid-radius closed level
+c = psi_axis + 0.5 * (eqt1d.psi[end] - psi_axis)
+seed = (maximum(r) - 1e-3, ZA)   # outboard; project finds the surface
+seed = IMAS._project_to_level(itp, c, seed)[1]
+
+for method in (:pc, :rk4)
+    t = @elapsed (Rs, Zs, closed) = IMAS._trace_surface_cubic(itp, c, seed; method)
+    drift = maximum(abs(IMAS.FI.value_gradient(itp,(Rs[k],Zs[k]))[1]-c) for k in eachindex(Rs))
+    @printf("%-5s closed=%s  N=%4d  drift=%.2e  t=%.2f ms\n", method, closed, length(Rs), drift, t*1e3)
+end
diff --git a/ext/IMASInterpolationsExt.jl b/ext/IMASInterpolationsExt.jl
new file mode 100644
index 00000000..6637fbbb
--- /dev/null
+++ b/ext/IMASInterpolationsExt.jl
@@ -0,0 +1,18 @@
+module IMASInterpolationsExt
+
+import IMAS
+import Interpolations
+
+# Interpolations.jl backend for the generic interpolant helpers in IMAS (src/physics/fields.jl).
+# Active only when Interpolations is loaded. It uses positional coords and returns SVector/SMatrix,
+# which index like the FastInterpolations NTuple/Matrix the IMAS math expects.
+@inline IMAS._gradient(itp::Interpolations.AbstractInterpolation, r, z) =
+    Interpolations.gradient(itp, r, z)
+
+@inline IMAS._value_gradient(itp::Interpolations.AbstractInterpolation, r, z) =
+    (itp(r, z), Interpolations.gradient(itp, r, z))
+
+@inline IMAS._hessian!(H, itp::Interpolations.AbstractInterpolation, r, z) =
+    (H .= Interpolations.hessian(itp, r, z); H)
+
+end
diff --git a/src/IMAS.jl b/src/IMAS.jl
index 576708c0..1890ecd7 100644
--- a/src/IMAS.jl
+++ b/src/IMAS.jl
@@ -4,6 +4,7 @@ using Printf
 using Compat:@compat
 import OrderedCollections
 import FastInterpolations as FI
+using AdaptiveArrayPools
 const document = OrderedCollections.OrderedDict()
 
 macro import_all(mod)
diff --git a/src/physics.jl b/src/physics.jl
index fcb9d70b..0225b682 100644
--- a/src/physics.jl
+++ b/src/physics.jl
@@ -20,6 +20,7 @@ include(joinpath("physics", "tf.jl"))
 include(joinpath("physics", "currents.jl"))
 include(joinpath("physics", "fields.jl"))
 include(joinpath("physics", "fluxsurfaces.jl"))
+include(joinpath("physics", "fluxsurfaces_cubic.jl"))
 include(joinpath("physics", "rf.jl"))
 include(joinpath("physics", "neoclassical.jl"))
 include(joinpath("physics", "profiles.jl"))
diff --git a/src/physics/fields.jl b/src/physics/fields.jl
index 8c2b02ca..5a4e5315 100644
--- a/src/physics/fields.jl
+++ b/src/physics/fields.jl
@@ -3,10 +3,13 @@ import SimpleNonlinearSolve
 
 document[Symbol("Physics fields")] = Symbol[]
 
-# gradient of a ψ interpolant: FI gets the native fast path,
-# the open version is for users still passing an Interpolations.jl interpolant
-@inline _psi_gradient(itp::FI.AbstractInterpolant, r, z) = FI.gradient(itp, (r, z))
-@inline _psi_gradient(itp, r, z) = parentmodule(typeof(itp)).gradient(itp, r, z)
+# Analytic gradient / value+gradient / Hessian of a ψ interpolant, dispatched on the backend.
+# FastInterpolations gets the native in-place fast path here; other backends (e.g. Interpolations.jl)
+# are added by package extensions, loaded only when that backend is (see ext/IMASInterpolationsExt.jl).
+# Results are indexed g[1]/g[2] and H[i,j], so any backend returning indexable values works.
+@inline _gradient(itp::FI.AbstractInterpolant, r, z) = FI.gradient(itp, (r, z))
+@inline _value_gradient(itp::FI.AbstractInterpolant, r, z) = FI.value_gradient(itp, (r, z))
+@inline _hessian!(H, itp::FI.AbstractInterpolant, r, z) = FI.hessian!(H, itp, (r, z))
 
 """
     Br_Bz(eqt2d::IMAS.equilibrium__time_slice___profiles_2d)
@@ -30,7 +33,7 @@ end
     Br_Bz(PSI_interpolant, r::T, z::T) where {T<:Real}
 """
 function Br_Bz(PSI_interpolant, r::T, z::T) where {T<:Real}
-    grad = _psi_gradient(PSI_interpolant, r, z)
+    grad = _gradient(PSI_interpolant, r, z)
     inv_twopi_r = 1.0 / (2π * r)
     Br = grad[2] * inv_twopi_r
     Bz = -grad[1] * inv_twopi_r
diff --git a/src/physics/fluxsurfaces.jl b/src/physics/fluxsurfaces.jl
index 91148000..ba7a82b8 100644
--- a/src/physics/fluxsurfaces.jl
+++ b/src/physics/fluxsurfaces.jl
@@ -1179,7 +1179,7 @@ function trace_simple_surfaces!(
 
             # flux expansion = 1 / abs(Bp)
             Br, Bz = Br_Bz(PSI_interpolant, pr[i], pz[i])
-            fluxexpansion[i] = 1.0 / sqrt(Br^2.0 + Bz^2.0)
+            fluxexpansion[i] = 1.0 / sqrt(Br^2 + Bz^2)
         end
         int_fluxexpansion_dl = trapz(ll, fluxexpansion)
 
@@ -1203,8 +1203,8 @@ function trace_surfaces(eqt::IMAS.equilibrium__time_slice{T}, wall_r::AbstractVe
     r, z, PSI_interpolant = ψ_interpolant(eqt2d)
     RA = eqt.global_quantities.magnetic_axis.r
     ZA = eqt.global_quantities.magnetic_axis.z
-    Br, Bz = Br_Bz(eqt2d)
-    return trace_surfaces(eqt.profiles_1d.psi, eqt.profiles_1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z; refine_extrema)
+    xpoints = [(xp.r, xp.z) for xp in eqt.boundary.x_point]
+    return trace_surfaces(eqt.profiles_1d.psi, eqt.profiles_1d.f, r, z, eqt2d.psi, PSI_interpolant, RA, ZA, wall_r, wall_z; refine_extrema, xpoints)
 end
 
 """
@@ -1214,8 +1214,6 @@ end
         r::AbstractVector{T},
         z::AbstractVector{T},
         PSI::Matrix{T},
-        BR::Matrix{T},
-        BZ::Matrix{T},
         PSI_interpolant,
         RA::T,
         ZA::T,
@@ -1224,25 +1222,28 @@ end
         refine_extrema::Bool=true
     ) where {T<:Real}
 """
-function trace_surfaces(
+@with_pool pool function trace_surfaces(
     psi::AbstractVector{T},
     f::AbstractVector{T},
     r::AbstractVector{T},
     z::AbstractVector{T},
     PSI::Matrix{T},
-    BR::Matrix{T},
-    BZ::Matrix{T},
     PSI_interpolant,
     RA::T,
     ZA::T,
     wall_r::AbstractVector{T},
     wall_z::AbstractVector{T};
-    refine_extrema::Bool=true
+    refine_extrema::Bool=true,
+    xpoints=()
 ) where {T<:Real}
 
     N = length(psi)
     surfaces = Vector{FluxSurface{T}}(undef, N)
-    r_cache, z_cache = IMASutils.contour_cache(r, z)
+    # Assume rectangle (aggression_level == 2)
+    Ncache = 2 * (length(r) + length(z))
+    r_cache = acquire!(pool, T, Ncache)
+    z_cache = acquire!(pool, T, Ncache)
+
     PSIA = PSI_interpolant(RA, ZA)
     for k in N:-1:1
         psi_level = psi[k]
@@ -1270,10 +1271,14 @@ function trace_surfaces(
         # surface length
         ll = arc_length(pr, pz)
 
+        # Allocate temporary arrays
+        Bp2 = acquire!(pool, T, length(pr))
+        Bp_abs = acquire!(pool, T, length(pr))
+
         # poloidal magnetic field (with sign)
         Br, Bz = Br_Bz(PSI_interpolant, pr, pz)
-        Bp2 = Br .^ 2.0 .+ Bz .^ 2.0
-        Bp_abs = sqrt.(Bp2)
+        Bp2 .= Br .^ 2 .+ Bz .^ 2
+        Bp_abs .= sqrt.(Bp2)
         Bp = Bp_abs .* sign.((pz .- ZA) .* Br .- (pr .- RA) .* Bz)
         Btot = sqrt.(Bp2 .+ (f[k] ./ pr) .^ 2)
 
@@ -1307,156 +1312,30 @@ function trace_surfaces(
     end
 
     if refine_extrema
-        N2 = round(Int, N / 2, RoundUp)
-        psi_norm = abs(psi[end] - psi[1]) / N
-        space_norm = (surfaces[N].max_r - surfaces[N].min_r) / 2 / N
-
-        # Find where Br changes sign
-        lines = Contour.lines(Contour.contour(r, z, BR, 0.0))
-        k = 0
-        d = Inf
-        for (kk, line) in enumerate(lines)
-            pr, pz = Contour.coordinates(line)
-            # plot!(pr, pz)
-            dd = minimum(filter(!isnan, sqrt.((pr .- surfaces[N2].max_r) .^ 2 .+ (pz .- surfaces[N2].z_at_max_r) .^ 2)); init=Inf)
-            if dd < d
-                d = dd
-                k = kk
-            end
-        end
-        leftright_r, leftright_z = Contour.coordinates(lines[k])
-        index = .!(isnan.(leftright_r) .|| isnan.(leftright_z))
-        leftright_r = @view leftright_r[index]
-        leftright_z = @view leftright_z[index]
-
-        # extrema in R
-        interp_r = interp1d(1:length(leftright_r), leftright_r)
-        interp_z = interp1d(1:length(leftright_z), leftright_z)
-        for k in 1:N
-            #@show "R", k
-            index = _extrema_index(leftright_r, leftright_z, surfaces[k].max_r, surfaces[k].z_at_max_r, :right)
-            cost =
-                x -> _extrema_cost(
-                    interp_r(x),
-                    interp_z(x),
-                    psi[k],
-                    PSI_interpolant,
-                    surfaces[k].max_r,
-                    surfaces[k].z_at_max_r,
-                    RA,
-                    ZA,
-                    psi_norm,
-                    space_norm,
-                    :right
-                )
-            x = Optim.optimize(cost, index[1], index[end], Optim.Brent()).minimizer
-            surfaces[k].max_r, surfaces[k].z_at_max_r = interp_r(x), interp_z(x)
-            index = _extrema_index(leftright_r, leftright_z, surfaces[k].min_r, surfaces[k].z_at_min_r, :left)
-            cost =
-                x -> _extrema_cost(
-                    interp_r(x),
-                    interp_z(x),
-                    psi[k],
-                    PSI_interpolant,
-                    surfaces[k].min_r,
-                    surfaces[k].z_at_min_r,
-                    RA,
-                    ZA,
-                    psi_norm,
-                    space_norm,
-                    :left
-                )
-            #plot!(leftright_r[index], leftright_z[index]; label="")
-            x = Optim.optimize(cost, index[1], index[end], Optim.Brent()).minimizer
-            min_r, z_at_min_r = interp_r(x), interp_z(x)
-            surfaces[k].min_r, surfaces[k].z_at_min_r = min_r, z_at_min_r
-            @assert surfaces[k].min_r < surfaces[k].max_r
-            if k < 3
-                surfaces[k+1].z_at_max_r = ZA
-                surfaces[k+1].z_at_min_r = ZA
-            elseif k < N
-                surfaces[k+1].z_at_max_r = (surfaces[k+1].z_at_max_r + surfaces[k].z_at_max_r) / 2.0
-                surfaces[k+1].z_at_min_r = (surfaces[k].z_at_min_r + surfaces[k+1].z_at_min_r) / 2.0
-            end
+        # Refine the four geometric extrema (max_r/min_r/max_z/min_z) with a globalized X-point-aware
+        # 2-D Newton (`_robust_refine_extremum!`) on the ψ interpolant (analytic ∇ψ/Hessian), sharing
+        # one 2×2 Hessian scratch across all calls.
+        axis = (RA, ZA)
+        lo = (first(r), first(z))   # ψ grid domain box — the refine search is clamped to it
+        hi = (last(r), last(z))
+        H = acquire!(pool, T, 2, 2)
+        for k in 2:N   # skip k=1 (artificial on-axis surface); rebuilt below
+            s = surfaces[k]
+            (s.max_r, s.z_at_max_r) = _robust_refine_extremum!(H, PSI_interpolant, psi[k], (s.max_r, s.z_at_max_r), :R, axis, xpoints; lo, hi)
+            (s.min_r, s.z_at_min_r) = _robust_refine_extremum!(H, PSI_interpolant, psi[k], (s.min_r, s.z_at_min_r), :R, axis, xpoints; lo, hi)
+            (s.r_at_max_z, s.max_z) = _robust_refine_extremum!(H, PSI_interpolant, psi[k], (s.r_at_max_z, s.max_z), :Z, axis, xpoints; lo, hi)
+            (s.r_at_min_z, s.min_z) = _robust_refine_extremum!(H, PSI_interpolant, psi[k], (s.r_at_min_z, s.min_z), :Z, axis, xpoints; lo, hi)
         end
 
-        # Find where Bz changes sign
-        lines = Contour.lines(Contour.contour(r, z, BZ, 0.0))
-        k = 0
-        d = Inf
-        for (kk, line) in enumerate(lines)
-            pr, pz = Contour.coordinates(line)
-            # plot!(pr, pz)
-            dd = minimum(filter(!isnan, sqrt.((pr .- surfaces[N2].r_at_max_z) .^ 2 .+ (pz .- surfaces[N2].max_z) .^ 2)); init=Inf)
-            if dd < d
-                d = dd
-                k = kk
-            end
-        end
-        updown_r, updown_z = Contour.coordinates(lines[k])
-        index = .!(isnan.(updown_r) .|| isnan.(updown_z))
-        updown_r = @view updown_r[index]
-        updown_z = @view updown_z[index]
-
-        # extrema in Z
-        interp_r = interp1d(1:length(updown_r), updown_r)
-        interp_z = interp1d(1:length(updown_z), updown_z)
-        for k in 1:N
-            #@show "Z", k
-            index = _extrema_index(updown_r, updown_z, surfaces[k].r_at_max_z, surfaces[k].max_z, :up)
-            cost =
-                x -> _extrema_cost(
-                    interp_r(x),
-                    interp_z(x),
-                    psi[k],
-                    PSI_interpolant,
-                    surfaces[k].r_at_max_z,
-                    surfaces[k].max_z,
-                    RA,
-                    ZA,
-                    psi_norm,
-                    space_norm,
-                    :up
-                )
-            x = Optim.optimize(cost, index[1], index[end], Optim.Brent()).minimizer
-            surfaces[k].r_at_max_z, surfaces[k].max_z = interp_r(x), interp_z(x)
-            index = _extrema_index(updown_r, updown_z, surfaces[k].r_at_min_z, surfaces[k].min_z, :down)
-            cost =
-                x -> _extrema_cost(
-                    interp_r(x),
-                    interp_z(x),
-                    psi[k],
-                    PSI_interpolant,
-                    surfaces[k].r_at_min_z,
-                    surfaces[k].min_z,
-                    RA,
-                    ZA,
-                    psi_norm,
-                    space_norm,
-                    :down
-                )
-            # plot!(updown_r[index], updown_z[index]; label="")
-            x = Optim.optimize(cost, index[1], index[end], Optim.Brent()).minimizer
-            surfaces[k].r_at_min_z, surfaces[k].min_z = interp_r(x), interp_z(x)
-            @assert surfaces[k].min_z < surfaces[k].max_z
-            if k < 3
-                surfaces[k+1].r_at_max_z = RA
-                surfaces[k+1].r_at_min_z = RA
-            elseif k < N
-                surfaces[k+1].r_at_max_z = (surfaces[k+1].r_at_max_z + surfaces[k].r_at_max_z) / 2.0
-                surfaces[k+1].r_at_min_z = (surfaces[k+1].r_at_min_z + surfaces[k].r_at_min_z) / 2.0
-            end
-        end
-
-        # first flux surface just a scaled down version of the second one
+        # first flux surface just a scaled down version of the second one (all scalar fields)
         frac = 0.01
-        surfaces[1].r_at_max_z = (surfaces[2].r_at_max_z .- RA) .* frac .+ RA
-        surfaces[1].max_z = (surfaces[2].max_z .- ZA) .* frac .+ ZA
-        surfaces[1].r_at_min_z = (surfaces[2].r_at_min_z .- RA) .* frac .+ RA
-        surfaces[1].min_z = (surfaces[2].min_z .- ZA) .* frac .+ ZA
-        surfaces[1].z_at_max_r = (surfaces[2].z_at_max_r .- ZA) .* frac .+ ZA
+        surfaces[1].r_at_max_z = (surfaces[2].r_at_max_z - RA) * frac + RA
+        surfaces[1].max_z = (surfaces[2].max_z - ZA) * frac + ZA
+        surfaces[1].r_at_min_z = (surfaces[2].r_at_min_z - RA) * frac + RA
+        surfaces[1].min_z = (surfaces[2].min_z - ZA) * frac + ZA
+        surfaces[1].z_at_max_r = (surfaces[2].z_at_max_r - ZA) * frac + ZA
         surfaces[1].max_r = (surfaces[2].max_r - RA) * frac + RA
-        surfaces[1].z_at_min_r = (surfaces[2].z_at_min_r .- ZA) .* frac .+ ZA
+        surfaces[1].z_at_min_r = (surfaces[2].z_at_min_r - ZA) * frac + ZA
         surfaces[1].min_r = (surfaces[2].min_r - RA) * frac + RA
     end
 
@@ -1466,74 +1345,6 @@ end
 @compat public trace_surfaces
 push!(document[Symbol("Physics flux-surfaces")], :trace_surfaces)
 
-function _extrema_index(r::AbstractVector{T}, z::AbstractVector{T}, r0::T, Z0::T, direction::Symbol) where {T<:Real}
-    i = argmin((r .- r0) .^ 2 .+ (z .- Z0) .^ 2)
-    N = length(z)
-    j = round(Int, N / 2, RoundUp)
-    n = 3
-    if direction == :right
-        if (r[j] - r[j-1]) > 0 # oriented right
-            return max(1, i - 1):min(N, i + n)
-        else # opposite orientation
-            return max(1, i - n):min(N, i + 1)
-        end
-    elseif direction == :left
-        if (r[j] - r[j-1]) < 0 # oriented left
-            return max(1, i - 1):min(N, i + n)
-        else
-            return max(1, i - n):min(N, i + 1)
-        end
-    elseif direction == :up
-        if (z[j] - z[j-1]) > 0 # oriented up
-            return max(1, i - 1):min(N, i + n)
-        else
-            return max(1, i - n):min(N, i + 1)
-        end
-    elseif direction == :down
-        if (z[j] - z[j-1]) < 0 # oriented down
-            return max(1, i - 1):min(N, i + n)
-        else
-            return max(1, i - n):min(N, i + 1)
-        end
-    else
-        error("_extrema_index(..., direction::Symbol) can only be (:left, :right, :up, :down)")
-    end
-end
-
-# accurate geometric quantities by finding geometric extrema as optimization problem
-function _extrema_cost(
-    r::T,
-    z::T,
-    psi_level::T,
-    PSI_interpolant,
-    r_orig::T,
-    z_orig::T,
-    RA::T,
-    ZA::T,
-    psi_norm::T,
-    space_norm::T,
-    direction::Symbol
-) where {T<:Real}
-    cost_psi = (PSI_interpolant(r, z) - psi_level) / psi_norm * space_norm # convert psi cost into spatial units
-    if direction == :right
-        cost_dir = abs(r - r_orig) * (r > r_orig)
-        cost_dir += abs(r - RA) * (r < RA) + (r < RA) # extra penalty needed for flux surfaces near axis
-    elseif direction == :left
-        cost_dir = abs(r - r_orig) * (r < r_orig)
-        cost_dir += abs(r - RA) * (r > RA) + (r > RA)
-    elseif direction == :up
-        cost_dir = abs(z - z_orig) * (z > z_orig)
-        cost_dir += abs(z - ZA) * (z < ZA) + (z < ZA)
-    elseif direction == :down
-        cost_dir = abs(z - z_orig) * (z < z_orig)
-        cost_dir += abs(z - ZA) * (z > ZA) + (z > ZA)
-    else
-        error("_extrema_cost(..., direction::Symbol) can only be (:left, :right, :up, :down)")
-    end
-    cost = norm((cost_psi, cost_dir^2)) # linear in psi since it already varies quadratically in space
-    return cost
-end
-
 """
     flux_surfaces(eq::equilibrium{T}, wall_r::AbstractVector{T}, wall_z::AbstractVector{T}) where {T<:Real}
 
@@ -1649,8 +1460,7 @@ function flux_surfaces(eqt::equilibrium__time_slice{T1}, wall_r::AbstractVector{
     end
 
     # trace flux surfaces
-    Br, Bz = Br_Bz(eqt2d)
-    surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z;refine_extrema=true)
+    surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, PSI_interpolant, RA, ZA, wall_r, wall_z; refine_extrema=true, xpoints=[(xp.r, xp.z) for xp in eqt.boundary.x_point])
 
     # calculate flux surface averaged and geometric quantities
     N = length(eqt1d.psi)
@@ -2278,14 +2088,28 @@ Returns extrema indexes and values of R,Z flux surfaces vectors:
     z_at_min_r, min_r
 """
 function fluxsurface_extrema(pr::Vector{T}, pz::Vector{T}) where {T<:Real}
-    _, imaxr = findmax(pr)
-    _, iminr = findmin(pr)
-    _, imaxz = findmax(pz)
-    _, iminz = findmin(pz)
-    r_at_max_z, max_z = pr[imaxz], pz[imaxz]
-    r_at_min_z, min_z = pr[iminz], pz[iminz]
-    z_at_max_r, max_r = pz[imaxr], pr[imaxr]
-    z_at_min_r, min_r = pz[iminr], pr[iminr]
+    n = length(pr)
+    n == length(pz) || throw(DimensionMismatch("fluxsurface_extrema: pr and pz must have equal length"))
+    # One sweep tracks all four extrema (value+index) at once — an order of magnitude faster than
+    # four findmax/findmin calls (NaN-aware ordering + index reduction don't vectorize). Plain >/<
+    # is safe: traced flux-surface coordinates are finite.
+    @inbounds begin
+        max_r = min_r = pr[1]
+        max_z = min_z = pz[1]
+        imaxr = iminr = imaxz = iminz = 1
+        for i in 2:n
+            pri = pr[i]
+            pzi = pz[i]
+            if pri > max_r; max_r = pri; imaxr = i; end
+            if pri < min_r; min_r = pri; iminr = i; end
+            if pzi > max_z; max_z = pzi; imaxz = i; end
+            if pzi < min_z; min_z = pzi; iminz = i; end
+        end
+        r_at_max_z = pr[imaxz]
+        r_at_min_z = pr[iminz]
+        z_at_max_r = pz[imaxr]
+        z_at_min_r = pz[iminr]
+    end
     return (imaxr, iminr, imaxz, iminz,
         r_at_max_z, max_z,
         r_at_min_z, min_z,
diff --git a/src/physics/fluxsurfaces_cubic.jl b/src/physics/fluxsurfaces_cubic.jl
new file mode 100644
index 00000000..fcd16d85
--- /dev/null
+++ b/src/physics/fluxsurfaces_cubic.jl
@@ -0,0 +1,555 @@
+# Cubic-interpolant-based flux-surface geometry.
+#
+# Flux-surface geometry computed directly from the ψ cubic interpolant (analytic
+# gradient/Hessian via FastInterpolations, aliased `FI` in IMAS.jl) with Newton
+# iteration — as opposed to the Contour-based tracing in fluxsurfaces.jl.
+#
+# Currently: X-point-aware refinement of a traced surface's geometric extrema
+# (max_r/min_r/max_z/min_z). Intended to also host a general cubic-interpolant-
+# based flux-surface tracer. All helpers here are private (no public API yet).
+
+# Split form of the critical-point system ∇ψ = 0 (Jacobian = Hessian) for the globalized
+# [`_damped_newton2d`](@ref): `residual` returns the residual (∇ψ) and the first Hessian row,
+# `hessrow` the second — same contract as [`_extremum_eqs`](@ref), so the one safe solver serves
+# both systems. (Critical-point finding is rare, so the probe recomputing the Hessian is not hot.)
+function _critical_eqs(H::AbstractMatrix, itp)
+    residual(R, Z) = begin
+        g = _gradient(itp, R, Z)
+        _hessian!(H, itp, R, Z)
+        return (g[1], g[2], H[1, 1], H[1, 2])
+    end
+    hessrow(R, Z) = begin
+        _hessian!(H, itp, R, Z)
+        return (H[2, 1], H[2, 2])
+    end
+    return (; residual, hessrow)
+end
+
+# Split form of the extremum system for the globalized [`_damped_newton2d`](@ref): `residual`
+# (one `value_gradient`) returns the residual and the ψ=target Jacobian row (∇ψ, free from the
+# same eval); `hessrow` (one `hessian!`) returns the ∂ψ/∂n Jacobian row. The line search probes
+# with `residual` alone and pays `hessrow` only on accepted steps — halving the interpolant work
+# per Newton iterate, with a bit-identical iterate sequence.
+function _extremum_eqs(H::AbstractMatrix, itp, target_psi::Real, daxis::Int)
+    residual(R, Z) = begin
+        val, g = _value_gradient(itp, R, Z)
+        return (val - target_psi, g[daxis], g[1], g[2])
+    end
+    hessrow(R, Z) = begin
+        _hessian!(H, itp, R, Z)
+        return (H[daxis, 1], H[daxis, 2])
+    end
+    return (; residual, hessrow)
+end
+
+"""
+    _damped_newton2d(eqs, seed; tol=1e-11, maxit=50, αmin=1e-3)
+
+Globalized 2×2 Newton with a backtracking line search on `‖F‖` and a gradient-descent fallback
+when the Jacobian is near-singular — robust where pure Newton diverges (e.g. near the separatrix,
+where `det(J) = ψ_R·ψ_ZZ → 0`, small poloidal curvature, blows up the plain step). Returns
+`((R, Z), converged::Bool)`.
+
+`eqs` is the split extremum system from [`_extremum_eqs`](@ref): `eqs.residual(R,Z)` returns
+`(F1, F2, J11, J12)` — the residual plus the ψ=target Jacobian row (∇ψ), both from one
+`value_gradient` — and `eqs.hessrow(R,Z)` returns `(J21, J22)`, the `∂ψ/∂n` row (a `hessian`).
+The line search probes with `residual` alone; `hessrow` runs only on accepted steps, halving the
+interpolant work per iterate versus a combined residual+Jacobian eval (the iterate sequence is
+bit-identical: the `‖F‖` accept test uses only the residual, and `hessrow` runs at the exact
+accepted point a combined eval would have).
+"""
+function _damped_newton2d(eqs, seed::Tuple{T,T}; tol::Real=1e-11, maxit::Int=50, αmin::Real=1e-3,
+    lo=(-Inf, -Inf), hi=(Inf, Inf)) where {T<:Real}
+    # every iterate is clamped to the box [lo, hi] (the ψ grid domain) so the search can never
+    # wander into the extrapolation region outside the grid — the iterate physically cannot escape.
+    R, Z = clamp(seed[1], lo[1], hi[1]), clamp(seed[2], lo[2], hi[2])
+    F1, F2, J11, J12 = eqs.residual(R, Z)
+    J21, J22 = eqs.hessrow(R, Z)
+    nrm = hypot(F1, F2)
+    for _ in 1:maxit
+        nrm <= tol && return ((R, Z), true)
+        det = J11 * J22 - J12 * J21
+        if isfinite(det) && abs(det) > 1e-13
+            dR = (J22 * F1 - J12 * F2) / det
+            dZ = (J11 * F2 - J21 * F1) / det
+        else  # near-singular J: descend 0.5‖F‖² along Jᵀ F with a small normalized step
+            dR = J11 * F1 + J21 * F2
+            dZ = J12 * F1 + J22 * F2
+            s = hypot(dR, dZ)
+            s > 0 && (dR /= s; dZ /= s)
+            dR *= T(1e-2); dZ *= T(1e-2)
+        end
+        α = one(T); stepped = false
+        while α >= αmin
+            q1, q2 = clamp(R - α * dR, lo[1], hi[1]), clamp(Z - α * dZ, lo[2], hi[2])
+            if isfinite(q1) && isfinite(q2)
+                g1, g2, j11, j12 = eqs.residual(q1, q2)   # value_gradient only — no Hessian on probes
+                if hypot(g1, g2) < nrm
+                    R, Z = q1, q2
+                    F1, F2, J11, J12 = g1, g2, j11, j12    # reuse the probe's value_gradient
+                    J21, J22 = eqs.hessrow(R, Z)           # one Hessian on the accepted step
+                    nrm = hypot(F1, F2)
+                    stepped = true
+                    break
+                end
+            end
+            α /= 2
+        end
+        stepped || return ((R, Z), nrm <= tol)
+    end
+    return ((R, Z), nrm <= tol)
+end
+
+"""
+    _robust_refine_extremum!(H, itp, target_psi, seed, extremum_of, axis; tol=1e-11, maxit=50)
+
+X-point-aware refinement of a flux-surface geometric extremum. Solves the extremum system
+`{ψ = target_psi, ∂ψ/∂n = 0}` from `seed` with a globalized [`_damped_newton2d`](@ref)
+(no divergence near the separatrix), then validates the result by *physical region* — cheaply,
+using the precomputed critical points (`axis` = O-point, plus `xpoints` = X-points) as
+references, with NO per-call critical-point search:
+
+  1. **axis-relative direction** — the extremized coordinate must be on the same side of the
+     magnetic `axis` as the seed (kills the opposite O-point solution);
+  2. **confined side of every X-point** — `(p−xp)·(axis−xp) > 0` for each `xp ∈ xpoints`
+     (kills the private-flux-region / across-X-point solution). On-surface already excludes the
+     SOL (ψ_N>1), so the only ambiguity left is confined-surface vs private-region, which an
+     X-point dot test settles. Generic: `xpoints` empty (limited plasma → no private region)
+     reduces to the direction check; one or two X-points (single/double null) just loop.
+
+`xpoints` is a collection of `(R, Z)` (e.g. from `eqt.boundary.x_point`). If the solution is
+on-surface and confined it is accepted; otherwise (a private/wrong branch, only reachable from
+a bad seed) it re-seeds along the segment toward the axis and re-solves; failing that, falls
+back to `seed` (always an on-surface contour vertex).
+"""
+function _robust_refine_extremum!(H::AbstractMatrix, itp, target_psi::T, seed::Tuple{T,T},
+    extremum_of::Symbol, axis::Tuple{T,T}, xpoints=(); tol::Real=1e-11, maxit::Int=50, lo=(-Inf, -Inf), hi=(Inf, Inf)) where {T<:Real}
+    d = extremum_of === :R ? 2 : extremum_of === :Z ? 1 :
+        throw(ArgumentError("_robust_refine_extremum!: extremum_of must be :R or :Z, got :$extremum_of"))
+    e = extremum_of === :R ? 1 : 2
+    want_max = seed[e] > axis[e]
+    onsurf(q) = abs(itp(q[1], q[2]) - target_psi) < 1e-7
+    function confined(q::Tuple{T,T})
+        ((q[e] > axis[e]) == want_max) || return false                       # (1) axis-relative direction
+        for xp in xpoints                                                    # (2) confined side of every X-point
+            (q[1] - xp[1]) * (axis[1] - xp[1]) + (q[2] - xp[2]) * (axis[2] - xp[2]) > 0 || return false
+        end
+        return true
+    end
+
+    eqs = _extremum_eqs(H, itp, target_psi, d)
+    p, ok = _damped_newton2d(eqs, seed; tol, maxit, lo, hi)
+    (ok && onsurf(p) && confined(p)) && return p
+
+    # wrong branch (private/across-X-point) — only reachable from a bad seed. Re-seed along the
+    # segment from `p` toward the magnetic `axis` (always confined) and re-solve: the axis anchor
+    # pulls the Newton back onto the confined branch.
+    for f in (T(0.3), T(0.5), T(0.7))
+        reseed = (p[1] + f * (axis[1] - p[1]), p[2] + f * (axis[2] - p[2]))
+        q, okq = _damped_newton2d(eqs, reseed; tol, maxit, lo, hi)
+        (okq && onsurf(q) && confined(q)) && return q
+    end
+    return seed
+end
+
+# convenience wrapper for the robust refine (allocates the 2×2 scratch; one-off calls / tests)
+_robust_refine_extremum(itp, target_psi::T, seed::Tuple{T,T}, extremum_of::Symbol,
+    axis::Tuple{T,T}, xpoints=(); kw...) where {T<:Real} =
+    _robust_refine_extremum!(Matrix{T}(undef, 2, 2), itp, target_psi, seed, extremum_of, axis, xpoints; kw...)
+
+"""
+    _project_to_level(itp::FI.AbstractInterpolant, c::T, x::Tuple{T,T}; tol::Real=1e-10, maxit::Int=8) where {T<:Real}
+
+Newton corrector that projects `x = (R,Z)` onto the level set ψ=c **along the gradient**:
+`x ← x - (ψ(x) - c) · ∇ψ / |∇ψ|²` (the codimension-1 Moore–Penrose Newton step; quadratic
+convergence). Returns `((R, Z), converged::Bool)`; `converged=false` if `|∇ψ|` collapses
+(a critical point) or iterates go non-finite.
+"""
+function _project_to_level(itp::FI.AbstractInterpolant, c::T, x::Tuple{T,T}; tol::Real=1e-10, maxit::Int=8) where {T<:Real}
+    R, Z = x
+    for _ in 1:maxit
+        val, g = FI.value_gradient(itp, (R, Z))
+        r = val - c
+        abs(r) <= tol && return ((R, Z), true)
+        n2 = g[1]^2 + g[2]^2
+        (isfinite(n2) && n2 > eps(T)) || return ((R, Z), false)   # critical point / degenerate
+        R -= r * g[1] / n2
+        Z -= r * g[2] / n2
+        (isfinite(R) && isfinite(Z)) || return ((R, Z), false)
+    end
+    val = itp(R, Z)
+    return ((R, Z), abs(val - c) <= tol)
+end
+
+"""
+    _contour_tangent(itp::FI.AbstractInterpolant, x::Tuple{T,T}, sgn::Int) where {T<:Real}
+
+Unit tangent of the ψ iso-contour at `x`: `sgn·(ψ_Z, -ψ_R)/|∇ψ|` (rotate ∇ψ by 90°).
+`sgn ∈ {+1,-1}` selects traversal orientation (CW vs CCW).
+"""
+function _contour_tangent(itp::FI.AbstractInterpolant, x::Tuple{T,T}, sgn::Int) where {T<:Real}
+    g = FI.gradient(itp, x)
+    nrm = hypot(g[1], g[2])
+    return (sgn * g[2] / nrm, -sgn * g[1] / nrm)
+end
+
+"""
+    _contour_curvature(itp::FI.AbstractInterpolant, x::Tuple{T,T}) where {T<:Real}
+
+Signed curvature of the ψ iso-contour at `x`:
+`κ = (ψ_R²ψ_ZZ - 2ψ_Rψ_Zψ_RZ + ψ_Z²ψ_RR) / |∇ψ|³` (needs the Hessian).
+
+`H` is an optional caller-owned 2×2 scratch buffer; when threaded through the tracer's step
+loop it lets the per-step `hessian!` avoid allocating a fresh matrix each call.
+"""
+function _contour_curvature(itp::FI.AbstractInterpolant, x::Tuple{T,T},
+    H::AbstractMatrix=Matrix{T}(undef, 2, 2)) where {T<:Real}
+    g = FI.gradient(itp, x)
+    FI.hessian!(H, itp, x)
+    gR, gZ = g[1], g[2]
+    n2 = gR^2 + gZ^2
+    return (gR^2 * H[2, 2] - 2 * gR * gZ * H[1, 2] + gZ^2 * H[1, 1]) / n2^(T(3) / 2)
+end
+
+"""
+    _step_pc(itp::FI.AbstractInterpolant, c::T, x::Tuple{T,T}, h::Real, sgn::Int;
+             corr_tol::Real=1e-10, corr_max::Int=8) where {T<:Real}
+
+One predictor–corrector step of arclength `h` along the ψ=c contour. Predictor is the
+Hessian-osculating 2nd-order Taylor step `x + h·t + ½h²·κ·n_p` (`t` tangent, `n_p` principal
+normal `(-t_Z, t_R)`, `κ` curvature); corrector is [`_project_to_level`](@ref) back onto ψ=c.
+`H` is an optional caller-owned 2×2 Hessian scratch buffer forwarded to [`_contour_curvature`](@ref).
+Returns `(xnew, on_surface::Bool)`.
+"""
+function _step_pc(itp::FI.AbstractInterpolant, c::T, x::Tuple{T,T}, h::Real, sgn::Int,
+    H::AbstractMatrix=Matrix{T}(undef, 2, 2); corr_tol::Real=1e-10, corr_max::Int=8) where {T<:Real}
+    t = _contour_tangent(itp, x, sgn)
+    np = (-t[2], t[1])                            # principal normal (rotate tangent +90°)
+    κ = _contour_curvature(itp, x, H)
+    xp = (x[1] + h * t[1] + (h^2 / 2) * κ * np[1],
+          x[2] + h * t[2] + (h^2 / 2) * κ * np[2])
+    return _project_to_level(itp, c, xp; tol=corr_tol, maxit=corr_max)
+end
+
+"""
+    _contour_step(itp::FI.AbstractInterpolant, x::Tuple{T,T}; ε::Real=1e-6,
+                  h_min::Real, h_max::Real, max_turn::Real=deg2rad(10), κ_floor::Real=1e-8) where {T<:Real}
+
+Arclength step for the contour at `x`: bound the predictor chord error `≈ ½|κ|h² ≤ ε` and
+the per-step turning angle `|κ|·h ≤ max_turn`, then clamp to `[h_min, h_max]`. `κ_floor`
+prevents division by zero on near-straight pieces. `H` is an optional caller-owned 2×2 Hessian
+scratch buffer forwarded to [`_contour_curvature`](@ref).
+"""
+function _contour_step(itp::FI.AbstractInterpolant, x::Tuple{T,T},
+    H::AbstractMatrix=Matrix{T}(undef, 2, 2); ε::Real=1e-6,
+    h_min::Real, h_max::Real, max_turn::Real=deg2rad(10), κ_floor::Real=1e-8) where {T<:Real}
+    κ = max(abs(_contour_curvature(itp, x, H)), κ_floor)
+    h = sqrt(2 * ε / κ)            # chord-error bound
+    h = min(h, max_turn / κ)       # turning-angle cap
+    return clamp(h, h_min, h_max)
+end
+
+# signed angle from vector a to vector b (radians, in (-π, π])
+_signed_angle(a::Tuple, b::Tuple) = atan(a[1]*b[2] - a[2]*b[1], a[1]*b[1] + a[2]*b[2])
+
+# does segment p->q cross the ray from x0 along normal m (the Poincaré section ⟂ t0)?
+# returns the interpolation fraction in [0,1] if it crosses on the +section side, else nothing
+function _section_cross(p::Tuple{T,T}, q::Tuple{T,T}, x0::Tuple{T,T}, t0::Tuple{T,T}) where {T<:Real}
+    # Poincaré section through x0, ⟂ t0 (design §5.5). Detect sign change of the tangential
+    # coordinate (point - x0)·t0; the section line's normal is t0 itself.
+    fp = (p[1]-x0[1])*t0[1] + (p[2]-x0[2])*t0[2]
+    fq = (q[1]-x0[1])*t0[1] + (q[2]-x0[2])*t0[2]
+    (fp == fq) && return nothing
+    (sign(fp) == sign(fq)) && return nothing       # no crossing of the section line
+    frac = fp / (fp - fq)                           # in [0,1]
+    # require the crossing to be near x0 along the section line (n0 ⟂ t0), not the far side
+    cx = p[1] + frac*(q[1]-p[1]); cz = p[2] + frac*(q[2]-p[2])
+    n0 = (-t0[2], t0[1])
+    along = (cx-x0[1])*n0[1] + (cz-x0[2])*n0[2]
+    return abs(along) < 1e-2 ? frac : nothing
+end
+
+"""
+    _trace_surface_cubic(itp, c, seed; method=:pc, sgn=-1, ε, h_min, h_max, max_turn,
+                         κ_floor, s_min, max_steps, domain) -> (Rs, Zs, closed)
+
+Trace the ψ=c contour from `seed` by predictor–corrector continuation. Returns ordered
+`(Rs, Zs)` and whether the loop closed. Closure is tested only after the accumulated turning
+angle exceeds ~2π (a confined surface winds once), via a Poincaré section ⟂ the start
+tangent. `domain=(Rlo,Rhi,Zlo,Zhi)` (or `nothing`) terminates open contours at the boundary.
+"""
+function _trace_surface_cubic(itp::FI.AbstractInterpolant, c::T, seed::Tuple{T,T};
+    method::Symbol=:pc, sgn::Int=-1, ε::Real=1e-6, h_min::Real=1e-4, h_max::Real=0.1,
+    max_turn::Real=deg2rad(10), κ_floor::Real=1e-8, s_min::Real=0.0, max_steps::Int=100_000,
+    rk4_tol::Real=1e-8, τ_grad::Real=0.0,
+    domain=nothing) where {T<:Real}
+
+    x0, ok = _project_to_level(itp, c, seed)
+    ok || return (T[seed[1]], T[seed[2]], false)
+    t0 = _contour_tangent(itp, x0, sgn)
+    Rs = T[x0[1]]; Zs = T[x0[2]]
+    x = x0; tprev = t0; Θ = zero(T); s = zero(T)
+    hrk = h_max
+    H = Matrix{T}(undef, 2, 2)   # 2×2 Hessian scratch reused by every step (in-place hessian!)
+
+    for _ in 1:max_steps
+        # near-critical-point guard: when |∇ψ| < τ_grad (close to an X-point saddle), creep
+        # at the smallest step via the corrector so the surface traces the confined side
+        # instead of leaking through the saddle (τ_grad=0 disables this, default behavior).
+        gg = τ_grad > 0 ? FI.gradient(itp, x) : nothing
+        if gg !== nothing && hypot(gg[1], gg[2]) < τ_grad
+            xnew, sok = _step_pc(itp, c, x, h_min, sgn, H)
+        elseif method === :rk4
+            xnew, hrk, sok = _step_rk4_adaptive(itp, x, hrk, sgn; tol=rk4_tol, h_min, h_max)
+        else
+            h = _contour_step(itp, x, H; ε, h_min, h_max, max_turn, κ_floor)
+            xnew, sok = _step_pc(itp, c, x, h, sgn, H)
+        end
+        sok || break
+        if domain !== nothing && !(domain[1] <= xnew[1] <= domain[2] && domain[3] <= xnew[2] <= domain[4])
+            push!(Rs, xnew[1]); push!(Zs, xnew[2])
+            return (Rs, Zs, false)                # open contour: exited domain
+        end
+        tnew = _contour_tangent(itp, xnew, sgn)
+        Θ += _signed_angle(tprev, tnew)
+        s += hypot(xnew[1]-x[1], xnew[2]-x[2])
+        # closure: only after a full turn, and segment crosses the start section
+        if abs(Θ) >= 2π - 0.5 && s > s_min
+            frac = _section_cross(x, xnew, x0, t0)
+            if frac !== nothing
+                cx = x[1] + frac*(xnew[1]-x[1]); cz = x[2] + frac*(xnew[2]-x[2])
+                xc, _ = _project_to_level(itp, c, (cx, cz))   # interp crossing is O(h²) off the curve; snap onto ψ=c
+                push!(Rs, xc[1]); push!(Zs, xc[2])
+                return (Rs, Zs, true)
+            end
+        end
+        push!(Rs, xnew[1]); push!(Zs, xnew[2])
+        x = xnew; tprev = tnew
+    end
+    return (Rs, Zs, false)
+end
+
+"""
+    _resample_contour(Rs::AbstractVector{T}, Zs::AbstractVector{T}, n::Int) where {T<:Real}
+
+Resample a traced polyline to `n` points equally spaced in cumulative arclength (linear
+interpolation between traced vertices). Input is treated as ordered; the last point is the
+closure point for closed surfaces.
+"""
+function _resample_contour(Rs::AbstractVector{T}, Zs::AbstractVector{T}, n::Int) where {T<:Real}
+    m = length(Rs)
+    if m == 1
+        return (fill(Rs[1], n), fill(Zs[1], n))   # degenerate / single-point input
+    end
+    ll = zeros(T, m)
+    for k in 2:m
+        ll[k] = ll[k-1] + hypot(Rs[k]-Rs[k-1], Zs[k]-Zs[k-1])
+    end
+    L = ll[end]
+    L > 0 || return (fill(Rs[1], n), fill(Zs[1], n))
+    # detect a closed polyline (last vertex coincides with first); if so, sample the
+    # half-open interval [0, L) so the n outputs are distinct and the wrap segment equals
+    # the interior spacing instead of collapsing to ~0.
+    closed = hypot(Rs[end]-Rs[1], Zs[end]-Zs[1]) < 1e-9 * max(L, one(T))
+    targets = closed ? range(zero(T), L; length=n+1)[1:n] : range(zero(T), L; length=n)
+    Ro = Vector{T}(undef, n); Zo = Vector{T}(undef, n)
+    j = 1
+    for (i, s) in enumerate(targets)
+        while j < m && ll[j+1] < s
+            j += 1
+        end
+        seg = ll[j+1] - ll[j]
+        f = seg > 0 ? (s - ll[j]) / seg : zero(T)
+        Ro[i] = Rs[j] + f * (Rs[j+1] - Rs[j])
+        Zo[i] = Zs[j] + f * (Zs[j+1] - Zs[j])
+    end
+    return Ro, Zo
+end
+
+# one classic RK4 step of size h on dx/ds = tangent(x, sgn)
+function _rk4(itp::FI.AbstractInterpolant, x::Tuple{T,T}, h::Real, sgn::Int) where {T<:Real}
+    k1 = _contour_tangent(itp, x, sgn)
+    x2 = (x[1] + h/2*k1[1], x[2] + h/2*k1[2]);  k2 = _contour_tangent(itp, x2, sgn)
+    x3 = (x[1] + h/2*k2[1], x[2] + h/2*k2[2]);  k3 = _contour_tangent(itp, x3, sgn)
+    x4 = (x[1] + h*k3[1],   x[2] + h*k3[2]);    k4 = _contour_tangent(itp, x4, sgn)
+    return (x[1] + h/6*(k1[1]+2k2[1]+2k3[1]+k4[1]), x[2] + h/6*(k1[2]+2k2[2]+2k3[2]+k4[2]))
+end
+
+"""
+    _step_rk4_adaptive(itp, x, h, sgn; tol=1e-8, h_min=1e-5, h_max=0.1)
+
+Adaptive RK4 (step-doubling) on `dx/ds = tangent`. Compares one step of `h` against two of
+`h/2`; accepts the (extrapolated) half-step if the estimated error ≤ `tol`, and returns the
+next step size. **Comparison baseline only — no corrector**, so it measures the intrinsic
+off-surface drift of pure integration. Returns `(xnew, h_next, accepted::Bool)`.
+"""
+function _step_rk4_adaptive(itp::FI.AbstractInterpolant, x::Tuple{T,T}, h::Real, sgn::Int;
+    tol::Real=1e-8, h_min::Real=1e-5, h_max::Real=0.1) where {T<:Real}
+    hh = clamp(h, h_min, h_max)
+    for _ in 1:20
+        big = _rk4(itp, x, hh, sgn)
+        half = _rk4(itp, _rk4(itp, x, hh/2, sgn), hh/2, sgn)
+        err = hypot(big[1]-half[1], big[2]-half[2])
+        xnew = (half[1] + (half[1]-big[1])/15, half[2] + (half[2]-big[2])/15)  # local extrapolation
+        fac = err > 0 ? 0.9 * (tol/err)^(1/5) : 2.0
+        hnext = clamp(hh * clamp(fac, 0.2, 2.0), h_min, h_max)
+        (err <= tol || hh <= h_min) && return (xnew, hnext, true)
+        hh = hnext
+    end
+    return (_rk4(itp, x, hh, sgn), hh, false)
+end
+
+"""
+    _seed_omp(itp::FI.AbstractInterpolant, c::T, RA::T, ZA::T, R_max::T; nscan::Int=200) where {T<:Real}
+
+Outboard-midplane seed for the ψ=c contour: scan `R ∈ [RA, R_max]` at `Z = ZA` for the first
+sign change of `ψ(R,ZA) - c`, bisect to a bracket, then project onto ψ=c. Returns
+`((R,Z), found::Bool)`.
+"""
+function _seed_omp(itp::FI.AbstractInterpolant, c::T, RA::T, ZA::T, R_max::T; nscan::Int=200) where {T<:Real}
+    Rs = range(RA, R_max; length=nscan)
+    f(R) = itp(R, ZA) - c
+    prevf = f(Rs[1]); a = Rs[1]
+    for i in 2:nscan
+        fi = f(Rs[i])
+        if sign(fi) != sign(prevf)
+            lo, hi = Rs[i-1], Rs[i]
+            for _ in 1:60
+                mid = (lo + hi) / 2
+                (sign(f(mid)) == sign(f(lo))) ? (lo = mid) : (hi = mid)
+            end
+            return _project_to_level(itp, c, ((lo + hi) / 2, ZA))
+        end
+        prevf = fi; a = Rs[i]
+    end
+    return ((a, ZA), false)
+end
+
+"""
+    trace_surface_cubic(itp::FI.AbstractInterpolant, c::T, RA::T, ZA::T, R_max::T;
+                        npoints::Int=361, kw...) where {T<:Real}
+
+Trace a single closed flux surface ψ=c on the cubic interpolant: outboard-midplane seed →
+predictor–corrector trace → uniform-arclength resample to `npoints` distinct points →
+close and reorder via `reorder_flux_surface!` (default `force_close=true`).
+
+The returned vectors have length `npoints+1`: `npoints` distinct arclength-uniform points
+plus a closing duplicate (`r[end] == r[1]`, `z[end] == z[1]`), reordered so the outboard
+midplane (OMP) point is first and the polygon is clockwise. This matches the Contour-path
+`FluxSurface` representation, so `MXH` and arclength integrals (`int_fluxexpansion_dl`) work
+correctly on the output.
+
+Returns `(r, z, closed::Bool)`. Standalone — not wired into `trace_surfaces`.
+"""
+function trace_surface_cubic(itp::FI.AbstractInterpolant, c::T, RA::T, ZA::T, R_max::T;
+    npoints::Int=361, kw...) where {T<:Real}
+    seed, ok = _seed_omp(itp, c, RA, ZA, R_max)
+    ok || return (T[], T[], false)
+    Rs, Zs, closed = _trace_surface_cubic(itp, c, seed; kw...)
+    closed || return (Rs, Zs, false)
+    R, Z = _resample_contour(Rs, Zs, npoints)
+    reorder_flux_surface!(R, Z, RA, ZA)   # close (first==last), reorder OMP-first, clockwise — matches the Contour FluxSurface representation
+    return (R, Z, true)
+end
+
+"""
+    trace_surfaces_cubic(eqt::IMAS.equilibrium__time_slice{T}, wall_r::AbstractVector{T},
+                         wall_z::AbstractVector{T}; refine_extrema::Bool=true, npoints::Int=361, kw...) where {T<:Real}
+
+Cubic-interpolant drop-in mirror of [`trace_surfaces`](@ref)`(eqt, wall_r, wall_z)`: same
+high-level signature, returns the same `Vector{FluxSurface{T}}`, so the two can be compared
+1:1 (e.g. `trace_surfaces(eqt, fw.r, fw.z)` vs `trace_surfaces_cubic(eqt, fw.r, fw.z)`).
+Each ψ level is traced by predictor–corrector continuation on the bicubic interpolant instead
+of marching squares. `kw...` forwards to the tracer (`method=:pc|:rk4`, `ε`, `rk4_tol`, …).
+
+`wall_r`/`wall_z` are accepted for signature parity but currently unused: every level is
+traced as a closed surface (open/wall-clipping is not yet implemented in the cubic path), so a
+level that does not close raises an error. Standalone — NOT wired into `trace_surfaces`.
+"""
+function trace_surfaces_cubic(eqt::IMAS.equilibrium__time_slice{T}, wall_r::AbstractVector{T},
+    wall_z::AbstractVector{T}; refine_extrema::Bool=true, npoints::Int=361, kw...) where {T<:Real}
+    eqt2d = findfirst(:rectangular, eqt.profiles_2d)
+    r, _, itp = ψ_interpolant(eqt2d)
+    RA = eqt.global_quantities.magnetic_axis.r
+    ZA = eqt.global_quantities.magnetic_axis.z
+    return trace_surfaces_cubic(eqt.profiles_1d.psi, eqt.profiles_1d.f, RA, ZA, maximum(r), itp;
+        refine_extrema, npoints, kw...)
+end
+
+"""
+    trace_surfaces_cubic(psis, f, RA, ZA, R_max, itp; refine_extrema=false, npoints=361, kw...)
+
+Low-level batch flux-surface tracer on the cubic interpolant: for each level in `psis`, trace
+with [`trace_surface_cubic`](@ref) and populate a `FluxSurface` by reusing the existing
+field-eval helpers (`arc_length`, `Br_Bz`, `trapz`, `fluxsurface_extrema`). The innermost
+(k=1) surface is the usual artificial scaled-down copy of the second. When `refine_extrema`,
+the four geometric extrema (`max_r`/`min_r`/`max_z`/`min_z`) of each real surface are refined
+with [`_robust_refine_extremum!`](@ref) (globalized X-point-aware Newton on the interpolant) —
+the cubic analogue of the `refine_extrema` step in [`trace_surfaces`](@ref). Standalone.
+"""
+function trace_surfaces_cubic(psis::AbstractVector{T}, f::AbstractVector{T}, RA::T, ZA::T,
+    R_max::T, itp::FI.AbstractInterpolant; refine_extrema::Bool=false, npoints::Int=361, kw...) where {T<:Real}
+    N = length(psis)
+    surfaces = Vector{FluxSurface{T}}(undef, N)
+    axis = (RA, ZA)
+    lo = (first(itp.grids[1]), first(itp.grids[2]))   # ψ grid domain box — clamp the refine search
+    hi = (last(itp.grids[1]), last(itp.grids[2]))
+    H = Matrix{T}(undef, 2, 2)   # one 2×2 Hessian scratch shared by every _robust_refine_extremum! call
+    for k in N:-1:1
+        if k == 1
+            pr = (surfaces[2].r .- RA) ./ 100 .+ RA
+            pz = (surfaces[2].z .- ZA) ./ 100 .+ ZA
+        else
+            pr, pz, closed = trace_surface_cubic(itp, psis[k], RA, ZA, R_max; npoints, kw...)
+            if !closed
+                # A level passing through the X-point (separatrix) cannot be closed by the PC
+                # tracer — the cubic separatrix/open-surface layer is not implemented. Retry a
+                # hair toward the axis as a proxy and warn, so a full-profile call still runs.
+                psi_in = psis[k] + T(1e-3) * (psis[1] - psis[k])
+                pr, pz, closed = trace_surface_cubic(itp, psi_in, RA, ZA, R_max; npoints, kw...)
+                closed && @warn "trace_surfaces_cubic: ψ level $k of $N did not close; traced a slightly inner proxy (separatrix/open handling not implemented in the cubic path)"
+            end
+            closed || error("trace_surfaces_cubic: failed to close surface $k of $N at ψ=$(psis[k])")
+        end
+        ll = arc_length(pr, pz)
+        Br, Bz = Br_Bz(itp, pr, pz)
+        Bp2 = Br .^ 2 .+ Bz .^ 2
+        Bp_abs = sqrt.(Bp2)
+        Bp = Bp_abs .* sign.((pz .- ZA) .* Br .- (pr .- RA) .* Bz)
+        Btot = sqrt.(Bp2 .+ (f[k] ./ pr) .^ 2)
+        fluxexpansion = 1.0 ./ Bp_abs
+        int_fluxexpansion_dl = trapz(ll, fluxexpansion)
+        (_, _, _, _, r_at_max_z, max_z, r_at_min_z, min_z, z_at_max_r, max_r, z_at_min_r, min_r) =
+            fluxsurface_extrema(pr, pz)
+        s = FluxSurface(psis[k], pr, pz, r_at_max_z, max_z, r_at_min_z, min_z,
+            z_at_max_r, max_r, z_at_min_r, min_r, Br, Bz, Bp, Btot, ll, fluxexpansion, int_fluxexpansion_dl)
+        if refine_extrema && k != 1   # globalized X-point-aware refinement of the 4 geometric extrema (cubic analogue of trace_surfaces' refine_extrema)
+            (s.max_r, s.z_at_max_r) = _robust_refine_extremum!(H, itp, psis[k], (s.max_r, s.z_at_max_r), :R, axis; lo, hi)
+            (s.min_r, s.z_at_min_r) = _robust_refine_extremum!(H, itp, psis[k], (s.min_r, s.z_at_min_r), :R, axis; lo, hi)
+            (s.r_at_max_z, s.max_z) = _robust_refine_extremum!(H, itp, psis[k], (s.r_at_max_z, s.max_z), :Z, axis; lo, hi)
+            (s.r_at_min_z, s.min_z) = _robust_refine_extremum!(H, itp, psis[k], (s.r_at_min_z, s.min_z), :Z, axis; lo, hi)
+        end
+        surfaces[k] = s
+    end
+    return surfaces
+end
+
+"""
+    _find_xpoint(itp::FI.AbstractInterpolant, seed::Tuple{T,T}; tol=1e-10, maxit=50) where {T<:Real}
+
+Locate a critical point of ψ (`∇ψ=0`) by globalized 2-D Newton ([`_damped_newton2d`](@ref) +
+[`_critical_eqs`](@ref)) from `seed`, and classify it by the Hessian determinant: `:saddle`
+(X-point, `det H < 0`) or `:extremum` (O-point/axis, `det H > 0`). Returns
+`(point, kind::Symbol, converged::Bool)`. Unclamped (no grid box) so an out-of-domain seed that
+converges outside the grid is rejected by the domain guard below.
+"""
+function _find_xpoint(itp::FI.AbstractInterpolant, seed::Tuple{T,T}; tol::Real=1e-10, maxit::Int=50) where {T<:Real}
+    H = Matrix{T}(undef, 2, 2)   # 2×2 Hessian scratch (Newton residual + saddle/extremum classification)
+    pt, ok = _damped_newton2d(_critical_eqs(H, itp), seed; tol, maxit)
+    ok || return (pt, :none, false)
+    g1, g2 = itp.grids[1], itp.grids[2]                       # interpolant grid extents
+    (first(g1) <= pt[1] <= last(g1) && first(g2) <= pt[2] <= last(g2)) || return (pt, :none, false)
+    FI.hessian!(H, itp, pt)
+    detH = H[1, 1] * H[2, 2] - H[1, 2]^2
+    return (pt, detH < 0 ? :saddle : :extremum, true)
+end
diff --git a/src/physics/particles.jl b/src/physics/particles.jl
index 28518887..ca8b2f0f 100644
--- a/src/physics/particles.jl
+++ b/src/physics/particles.jl
@@ -178,7 +178,7 @@ function find_flux(particles::Vector{Particle{T}}, I_per_trace::T, rwall::Vector
     # smooth the load of each particle within a window
     # note: flux is defined at the cells, not at the nodes
 
-    d = sqrt.(diff(rwall) .^ 2.0 .+ diff(zwall) .^ 2.0) # length of each cell
+    d = sqrt.(diff(rwall) .^ 2 .+ diff(zwall) .^ 2) # length of each cell
     l = cumsum(d)
     wall_r = (rwall[1:end-1] .+ rwall[2:end]) ./ 2.0
     wall_z = (zwall[1:end-1] .+ zwall[2:end]) ./ 2.0
@@ -541,7 +541,7 @@ function pencil_beam(starting_position::Vector{T}, velocity_vector::Vector{T}, t
     x = starting_position[1] .+ velocity_vector[1] .* time
     y = starting_position[2] .+ velocity_vector[2] .* time
     z = starting_position[3] .+ velocity_vector[3] .* time
-    r = sqrt.(x .^ 2.0 .+ y .^ 2.0)
+    r = sqrt.(x .^ 2 .+ y .^ 2)
     return (x=x, y=y, z=z, r=r)
 end
 
diff --git a/src/physics/sol.jl b/src/physics/sol.jl
index eb637a9b..a710d0b3 100644
--- a/src/physics/sol.jl
+++ b/src/physics/sol.jl
@@ -96,10 +96,10 @@ function OpenFieldLine(
 
     # calculate quantities along field line
     Br, Bz = Br_Bz(PSI_interpolant, rr, zz) # r and z component of B for each point in (r,z)
-    Bp = sqrt.(Br .^ 2.0 .+ Bz .^ 2.0)     # poloidal component of B for each point in (r,z)
+    Bp = sqrt.(Br .^ 2 .+ Bz .^ 2)     # poloidal component of B for each point in (r,z)
     Bt = abs.(B0 .* R0 ./ rr)              # toroidal component of B for each point in (r,z)
     B = sqrt.(Bp .^ 2 + Bt .^ 2)           # total magnetic field B for each point in (r,z)
-    dp = sqrt.(gradient(rr) .^ 2.0 .+ gradient(zz) .^ 2.0) # curvilinear abscissa increments of poloidal projection of SOL surface
+    dp = sqrt.(gradient(rr) .^ 2 .+ gradient(zz) .^ 2) # curvilinear abscissa increments of poloidal projection of SOL surface
     pitch = sqrt.(1.0 .+ (Bt ./ Bp) .^ 2) # ds = dp*sqrt(1 + (Bt/Bp)^2) (pythagora)
     s = cumsum(pitch .* dp) # s = integral(ds)
     s = abs.(s .- s[midplane_index]) # fix 0 at outer midplane
diff --git a/src/plot.jl b/src/plot.jl
index 85a1ee7d..6eab6cbd 100644
--- a/src/plot.jl
+++ b/src/plot.jl
@@ -3074,7 +3074,7 @@ end
     title = "Wall flux"
     units = "[W/m²]"
     if component == :norm
-        data = sqrt.(nwl.flux_r .^ 2.0 .+ nwl.flux_z .^ 2.0)
+        data = sqrt.(nwl.flux_r .^ 2 .+ nwl.flux_z .^ 2)
     elseif component == :r
         data = nwl.flux_r
     elseif component == :z
diff --git a/test/Project.toml b/test/Project.toml
index 0c363327..3632044d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,2 +1,6 @@
 [deps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+Interpolations = "0.13, 0.14, 0.15, 0.16"
diff --git a/test/runtests.jl b/test/runtests.jl
index 84143589..7d79efba 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,10 @@ if !isempty(ARGS)
     end
 else
     # Default behavior: run all tests
+    include("runtests_refine_extremum.jl")
+
+    include("runtests_fluxsurfaces_cubic.jl")
+
     include("runtests_fluxsurfaces.jl")
 
     include("runtests_interpolations.jl")
diff --git a/test/runtests_fluxsurfaces.jl b/test/runtests_fluxsurfaces.jl
index e0244d07..481f4728 100644
--- a/test/runtests_fluxsurfaces.jl
+++ b/test/runtests_fluxsurfaces.jl
@@ -19,4 +19,24 @@ end
     ids = dd.equilibrium
 
     diff(ids_orig, ids; tol = 1E-1)
+end
+
+@testset "fluxsurface_extrema" begin
+    # explicit polyline; min_z ties at idx 1 and 5 -> first index wins (matches findmin)
+    pr = [1.0, 2.0, 1.5, 0.5, 1.0]
+    pz = [0.0, 0.5, 1.0, 0.3, 0.0]
+    #     imaxr iminr imaxz iminz  r@maxz max_z r@minz min_z  z@maxr max_r z@minr min_r
+    @test IMAS.fluxsurface_extrema(pr, pz) ==
+          (2, 4, 3, 1, 1.5, 1.0, 1.0, 0.0, 0.5, 2.0, 0.3, 0.5)
+
+    # equivalence with findmax/findmin on a closed D-shaped polyline
+    θ = range(0, 2π; length = 257)
+    pr2 = collect(1.7 .+ 0.6 .* cos.(θ .+ 0.3 .* sin.(θ)))
+    pz2 = collect(1.8 .* sin.(θ))
+    pr2[end], pz2[end] = pr2[1], pz2[1]
+    g = IMAS.fluxsurface_extrema(pr2, pz2)
+    @test (g[1], g[2], g[3], g[4]) == (findmax(pr2)[2], findmin(pr2)[2], findmax(pz2)[2], findmin(pz2)[2])
+    @test (g[6], g[8], g[10], g[12]) == (maximum(pz2), minimum(pz2), maximum(pr2), minimum(pr2))
+
+    @test_throws DimensionMismatch IMAS.fluxsurface_extrema([1.0, 2.0], [0.0])
 end
\ No newline at end of file
diff --git a/test/runtests_fluxsurfaces_cubic.jl b/test/runtests_fluxsurfaces_cubic.jl
new file mode 100644
index 00000000..c63849bc
--- /dev/null
+++ b/test/runtests_fluxsurfaces_cubic.jl
@@ -0,0 +1,365 @@
+using IMAS
+using Test
+
+# Analytic field: nested ellipses ψ = ((R-R0)/a0)^2 + (Z/b0)^2.
+# Level-L contour is an ellipse, semi-axes (a0√L in R, b0√L in Z),
+# enclosed area = π·a0·b0·L (closed-form check for tracing tests).
+function ellipse_itp(R0, a0, b0; n=65)
+    r = range(R0 - 1.4a0, R0 + 1.4a0, length=n)
+    z = range(-1.4b0, 1.4b0, length=n)
+    PSI = [((ri - R0) / a0)^2 + (zj / b0)^2 for ri in r, zj in z]
+    return IMAS.ψ_interpolant(r, z, PSI).PSI_interpolant, r, z
+end
+
+@testset "fluxsurfaces_cubic" begin
+    R0, a0, b0 = 1.7, 0.5, 1.0
+    itp, _, _ = ellipse_itp(R0, a0, b0)
+
+    @testset "_project_to_level snaps onto ψ=c" begin
+        c = 0.36
+        # a point off the c-contour (interior), should land on ψ=c
+        seed = (R0 + 0.3a0, 0.1b0)
+        (R, Z), ok = IMAS._project_to_level(itp, c, seed)
+        @test ok
+        val = itp(R, Z)
+        @test isapprox(val, c; atol=1e-9)
+    end
+
+    @testset "_project_to_level reports failure at a critical point" begin
+        # center of the ellipse: ∇ψ = 0 -> cannot project
+        (_, _), ok = IMAS._project_to_level(itp, 0.36, (R0, 0.0))
+        @test ok == false
+    end
+
+    @testset "_contour_tangent is unit and ⊥ ∇ψ" begin
+        x = (R0 + a0 * sqrt(0.36), 0.0)         # OMP of the ψ=0.36 ellipse
+        t = IMAS._contour_tangent(itp, x, 1)
+        g = IMAS.FI.gradient(itp, x)
+        @test isapprox(hypot(t[1], t[2]), 1.0; atol=1e-10)
+        @test isapprox(t[1]*g[1] + t[2]*g[2], 0.0; atol=1e-8)   # tangent ⟂ gradient
+    end
+
+    @testset "_contour_curvature matches a circle (κ = 1/radius)" begin
+        citp, _, _ = ellipse_itp(R0, 0.5, 0.5)   # a0=b0 -> circles
+        c = 0.25                                   # radius = 0.5*sqrt(0.25) = 0.25
+        x = (R0 + 0.5 * sqrt(c), 0.0)
+        κ = IMAS._contour_curvature(citp, x)
+        @test isapprox(abs(κ), 1 / 0.25; rtol=1e-3)
+    end
+
+    @testset "_step_pc lands on-surface and advances ~h" begin
+        c = 0.36
+        x0 = (R0 + a0 * sqrt(c), 0.0)            # exactly on ψ=c (OMP)
+        h = 0.02
+        x1, ok = IMAS._step_pc(itp, c, x0, h, 1)
+        @test ok
+        val = itp(x1)
+        @test isapprox(val, c; atol=1e-9)        # still on the surface after the step
+        @test isapprox(hypot(x1[1]-x0[1], x1[2]-x0[2]), h; rtol=0.05)  # advanced ~h along curve
+    end
+
+    @testset "_contour_step bounds turning angle and clamps" begin
+        citp, _, _ = ellipse_itp(R0, 0.5, 0.5)   # circle, κ = 1/radius known
+        c = 0.25; radius = 0.5 * sqrt(c)          # = 0.25, κ = 4
+        x = (R0 + radius, 0.0)
+        h = IMAS._contour_step(citp, x; h_min=1e-4, h_max=1.0, max_turn=deg2rad(10))
+        κ = abs(IMAS._contour_curvature(citp, x))
+        @test κ * h <= deg2rad(10) + 1e-9         # turning-angle cap respected
+        @test 1e-4 <= h <= 1.0                     # within clamp
+    end
+
+    @testset "_trace_surface_cubic traces a closed ellipse exactly" begin
+        c = 0.36
+        seed = (R0 + a0 * sqrt(c), 0.0)          # OMP seed, on-surface
+        Rs, Zs, closed = IMAS._trace_surface_cubic(itp, c, seed; h_max=0.05)
+        @test closed
+        @test length(Rs) > 20
+        # every point on ψ=c
+        @test all(abs(itp(Rs[k], Zs[k]) - c) < 1e-8 for k in eachindex(Rs))
+        # enclosed area (shoelace) == π·a0·b0·c
+        area = abs(sum(Rs[k]*Zs[mod1(k+1,length(Rs))] - Rs[mod1(k+1,length(Rs))]*Zs[k] for k in eachindex(Rs))) / 2
+        @test isapprox(area, π * a0 * b0 * c; rtol=1e-3)
+    end
+
+    @testset "_resample_contour gives n ~uniform-arclength points" begin
+        c = 0.36
+        Rs, Zs, _ = IMAS._trace_surface_cubic(itp, c, (R0 + a0*sqrt(c), 0.0); h_max=0.05)
+        R2, Z2 = IMAS._resample_contour(Rs, Zs, 128)
+        @test length(R2) == 128
+        # resampled points stay on ψ=c (they lie on the traced polyline, ~on-surface)
+        @test all(abs(itp(R2[k], Z2[k]) - c) < 1e-5 for k in eachindex(R2))
+        # arclength steps are all genuinely non-zero (no duplicate endpoint / collapsed wrap)
+        ds = [hypot(R2[mod1(k+1,128)]-R2[k], Z2[mod1(k+1,128)]-Z2[k]) for k in 1:128]
+        @test minimum(ds) > 0.5 * sum(ds) / 128
+        # closed-contour wrap segment equals the interior spacing (no duplicate endpoint)
+        @test isapprox(ds[end], sum(ds)/128; rtol=0.05)
+        # m==1 degenerate input: every output equals the single input point
+        R2s, Z2s = IMAS._resample_contour([2.0], [0.0], 128)
+        @test all(==(2.0), R2s) && all(==(0.0), Z2s) && length(R2s) == 128
+    end
+
+    @testset "_step_rk4_adaptive advances along the tangent, no corrector" begin
+        c = 0.36
+        x0 = (R0 + a0*sqrt(c), 0.0)
+        x1, hnext, acc = IMAS._step_rk4_adaptive(itp, x0, 0.02, 1)
+        @test acc
+        @test hnext > 0
+        @test isapprox(hypot(x1[1]-x0[1], x1[2]-x0[2]), 0.02; rtol=0.1)
+        # pure integrator: drift exists but is small for one step
+        @test abs(itp(x1) - c) < 1e-4
+    end
+
+    # Fix 7.1: exercise the rejection/grow paths of the adaptive step controller
+    @testset "_step_rk4_adaptive: step control shrinks/grows h" begin
+        c = 0.36
+        # high-curvature seed (ellipse top, κ=6.67 vs 0.83 at OMP) + oversized h + tight tol
+        xtop = (R0, b0*sqrt(c))
+        x1, hnext_s, acc_s = IMAS._step_rk4_adaptive(itp, xtop, 0.1, 1; tol=1e-12)
+        @test acc_s                       # eventually accepted after reducing h
+        @test hnext_s < 0.1               # controller SHRANK the step (rejection path ran)
+        @test abs(itp(x1) - c) < 1e-6  # accepted step accurate
+
+        # low-error case at OMP: err << tol -> controller GROWS the step (capped at h_max)
+        x0 = (R0 + a0*sqrt(c), 0.0)
+        _, hnext_g, acc_g = IMAS._step_rk4_adaptive(itp, x0, 0.02, 1; tol=1e-8)
+        @test acc_g
+        @test hnext_g > 0.02              # step grew toward h_max
+        @test hnext_g <= 0.1              # but stays clamped at h_max
+    end
+
+    @testset ":rk4 method traces the ellipse but drifts more than :pc" begin
+        c = 0.36; seed = (R0 + a0*sqrt(c), 0.0)
+        Rp, Zp, clp = IMAS._trace_surface_cubic(itp, c, seed; method=:pc, h_max=0.05)
+        Rr, Zr, clr = IMAS._trace_surface_cubic(itp, c, seed; method=:rk4, rk4_tol=1e-8, h_max=0.05)
+        @test clp && clr
+        drift_pc = maximum(abs(itp(Rp[k],Zp[k])-c) for k in eachindex(Rp))
+        # Fix 8.1: assert the corrector's actual invariant (drift bounded by corr_tol regardless
+        # of rk4_tol). The cross-method drift comparison drift_rk >= drift_pc is NOT a robust
+        # invariant on the exact-quadratic ellipse (RK4 can be tighter than the PC corrector tol
+        # when rk4_tol < corr_tol); replace it with the design's actual contract.
+        @test drift_pc < 1e-8                  # corrector enforces the constraint (100x margin)
+        @test drift_pc <= 1.1e-10              # corrector floor bounded by corr_tol regardless of rk4_tol
+    end
+
+    @testset "_seed_omp finds the outboard seed on ψ=c" begin
+        c = 0.36
+        seed, ok = IMAS._seed_omp(itp, c, R0, 0.0, R0 + 1.3a0)
+        @test ok
+        @test isapprox(itp(seed), c; atol=1e-9)
+        @test seed[1] > R0                       # outboard side
+
+        # Fix 9.1 (a): no-crossing negative case — ψ(R_max)=1.69 < 2.0 so no bracket exists;
+        # exercises the found=false branch and checks fallback returns last scanned R
+        seed2, ok2 = IMAS._seed_omp(itp, 2.0, R0, 0.0, R0 + 1.3a0)
+        @test ok2 == false
+        @test seed2[1] == R0 + 1.3a0
+
+        # Fix 9.1 (b): off-midplane chord at ZA=0.3 (within the c=0.36 ellipse whose Z-extent
+        # is b0*sqrt(c)=0.6); oblique projection makes bisection direction and projection visible.
+        seedz, okz = IMAS._seed_omp(itp, c, R0, 0.3, R0 + 1.3a0)
+        @test okz
+        @test isapprox(seedz[2], 0.3; atol=1e-10)                      # stays on the requested chord
+        @test isapprox(seedz[1], R0 + a0*sqrt(c - (0.3/b0)^2); atol=1e-7)  # analytic outboard root on chord
+        @test isapprox(itp(seedz), c; atol=1e-9)
+    end
+
+    @testset "trace_surface_cubic returns an ordered closed loop (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp2 = IMAS.ψ_interpolant(eqt2d)
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        psi_axis = itp2(RA, ZA)
+        c = psi_axis + 0.5 * (eqt.profiles_1d.psi[end] - psi_axis)   # mid-radius closed surface
+        R, Z, closed = IMAS.trace_surface_cubic(itp2, c, RA, ZA, maximum(rr); npoints=257)
+        @test closed
+        @test length(R) == 258   # npoints=257 distinct + 1 closing duplicate
+        @test R[1] == R[end] && Z[1] == Z[end]   # closed: first point repeated at end
+        # Fix 10.1 (B): the TRACER's 1e-9 on-surface invariant should be asserted on the
+        # PRE-resample traced vertices, not on the linearly-resampled output (which has ~1e-6
+        # chord error independent of tracer correctness). Assert on the raw traced points.
+        seed_c, ok_c = IMAS._seed_omp(itp2, c, RA, ZA, maximum(rr))
+        @test ok_c
+        vR, vZ, vclosed = IMAS._trace_surface_cubic(itp2, c, seed_c)
+        @test vclosed
+        @test all(abs(itp2(vR[k], vZ[k]) - c) < 1e-9 for k in eachindex(vR))
+        # coarse sanity on the resampled output (linear-interp chord error is ~1e-6, not 1e-9)
+        @test all(abs(itp2(R[k], Z[k]) - c) < 5e-6 for k in eachindex(R))
+        # Post-reorder contract: OMP-first (force_close=true circshifts so OMP is index 1),
+        # and clockwise orientation.
+        @test argmax(R) <= 3   # OMP (max-R) is now first after the circshift
+        # orientation: reorder_flux_surface! enforces clockwise (negative shoelace in R-right/Z-up axes)
+        @test sum(R[k]*Z[k+1] - R[k+1]*Z[k] for k in 1:length(R)-1) < 0
+        # MXH regression: MXH crashed on the old half-open output; verify it builds correctly on the closed output
+        mxh = IMAS.MXH(R, Z, 4)
+        @test isfinite(mxh.R0) && mxh.R0 > 0
+    end
+
+    @testset "cubic trace agrees with Contour path on a mid-radius surface (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp2 = IMAS.ψ_interpolant(eqt2d)
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        psi_axis = itp2(RA, ZA)
+        eqt1d = eqt.profiles_1d
+        c = psi_axis + 0.5 * (eqt1d.psi[end] - psi_axis)
+
+        # cubic trace
+        Rc, Zc, closed = IMAS.trace_surface_cubic(itp2, c, RA, ZA, maximum(rr))
+        @test closed
+        # Contour path (existing), same level
+        psis = [psi_axis, c]
+        ff   = [eqt1d.f[1], eqt1d.f[end]]
+        ref = IMAS.trace_surfaces(psis, ff, collect(rr), collect(zz), eqt2d.psi,
+                                  itp2, RA, ZA, fw.r, fw.z; refine_extrema=false)[2]
+        # compare bounding box (geometry) within a few mm
+        @test isapprox(maximum(Rc), maximum(ref.r); atol=5e-3)
+        @test isapprox(minimum(Rc), minimum(ref.r); atol=5e-3)
+        @test isapprox(maximum(Zc), maximum(ref.z); atol=5e-3)
+        @test isapprox(minimum(Zc), minimum(ref.z); atol=5e-3)
+    end
+
+    @testset "ISOLATION: cubic tracer is not wired into the existing path" begin
+        # the production tracer must not reference the cubic wrapper until validated
+        src = read(joinpath(pkgdir(IMAS), "src", "physics", "fluxsurfaces.jl"), String)
+        @test !occursin("trace_surface_cubic", src)
+        @test !occursin("_trace_surface_cubic", src)
+    end
+
+    @testset "trace_surfaces_cubic FSA agrees with Contour path (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp2 = IMAS.ψ_interpolant(eqt2d)
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        psi_axis = itp2(RA, ZA)
+        eqt1d = eqt.profiles_1d
+        # Fix 11b.1: honor the psis[1]=axis contract required by trace_surfaces (fluxsurfaces.jl
+        # always treats psis[1] as the artificial on-axis surface; the original test used
+        # psis=[frac0.3,...] which violated the contract — k=1 was an axis artifact, not frac=0.3).
+        # Prepend the axis so all three fracs are genuinely traced and compared.
+        fracs = [0.3, 0.5, 0.7]
+        psis = [psi_axis; psi_axis .+ fracs .* (eqt1d.psi[end] - psi_axis)]   # length 4: axis + 3 mid-radius
+        ff = fill(eqt1d.f[end], length(psis))
+
+        cub = IMAS.trace_surfaces_cubic(psis, ff, RA, ZA, maximum(rr), itp2)
+        ref = IMAS.trace_surfaces(psis, ff, collect(rr), collect(zz), eqt2d.psi,
+                                  itp2, RA, ZA, fw.r, fw.z; refine_extrema=false)
+
+        for k in 2:length(psis)   # skip k=1 (artificial axis copy); validate the 3 real fracs
+            gm1_c = IMAS.flux_surface_avg(1.0 ./ cub[k].r .^ 2, cub[k])
+            gm1_r = IMAS.flux_surface_avg(1.0 ./ ref[k].r .^ 2, ref[k])
+            @test isapprox(gm1_c, gm1_r; rtol=2e-3)
+            @test isapprox(cub[k].int_fluxexpansion_dl, ref[k].int_fluxexpansion_dl; rtol=5e-3)
+        end
+    end
+
+    @testset "_find_xpoint locates and classifies the upper X-point (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        _, _, itp2 = IMAS.ψ_interpolant(eqt2d)
+        xp = eqt.boundary.x_point[argmax([p.z for p in eqt.boundary.x_point])]
+        pt, kind, ok = IMAS._find_xpoint(itp2, (xp.r + 0.02, xp.z + 0.02))
+        @test ok
+        @test kind == :saddle
+        @test isapprox(pt[1], xp.r; atol=1e-2)
+        @test isapprox(pt[2], xp.z; atol=1e-2)
+        # Fix 12.1: O-point (magnetic axis) must classify as :extremum, exercising both branches
+        ax = (eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z)
+        po, kindo, oko = IMAS._find_xpoint(itp2, ax)
+        @test oko
+        @test kindo == :extremum
+        # Fix 12.2: out-of-domain seed must be rejected (domain guard rejects spurious extrapolated critical)
+        ptb, kindb, okb = IMAS._find_xpoint(itp2, (3.0, 2.0))
+        @test !okb
+    end
+
+    @testset "near-separatrix surface stays in the confined region (DIII-D)" begin
+        # Fix 13.1: The original test only asserted closed + confined at 98%, which passes even
+        # with τ_grad=0 (guard disabled). This is a GUARD-BRANCH SMOKE TEST: it uses τ_grad=0.15
+        # (above the near-separatrix min|∇ψ|≈0.07) so the guard predicate fires on ≥1 step,
+        # and asserts the guarded trace closes, stays confined, and is on-surface to 1e-9.
+        # Note: on this DIII-D slice the PC tracer does NOT reliably leak even at the separatrix,
+        # so a fully differential (leak-vs-no-leak) test cannot be constructed without a stressor
+        # that would change test semantics. This smoke test verifies: (1) the guard branch executes
+        # (τ_grad=0.15 fires at near-separatrix ∇ψ dips), (2) the result is still on-surface and
+        # confined, and (3) the guard path doesn't corrupt the trace. Whether the guard is strictly
+        # *necessary* for this case is documented but not asserted.
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp2 = IMAS.ψ_interpolant(eqt2d)
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        psi_axis = itp2(RA, ZA)
+        eqt1d = eqt.profiles_1d
+        pb = IMAS.find_psi_boundary(rr, zz, eqt2d.psi, psi_axis, eqt1d.psi[end], RA, ZA, fw.r, fw.z;
+            PSI_interpolant=itp2, raise_error_on_not_open=false, raise_error_on_not_closed=false)
+        c = psi_axis + 0.98 * (pb.last_closed - psi_axis)   # 98% out: close to the separatrix
+        xpz = maximum(p.z for p in eqt.boundary.x_point)
+        # τ_grad=0.15 is above the near-separatrix min|∇ψ|, ensuring the guard fires on ≥1 step
+        R1, Z1, closed1 = IMAS.trace_surface_cubic(itp2, c, RA, ZA, maximum(rr);
+                                                    τ_grad=0.15, h_max=0.02)
+        @test closed1                              # guarded trace closes
+        @test maximum(Z1) < xpz                   # guarded trace stays confined (below X-point)
+        # on-surface invariant: guard must not emit off-surface points (§5/§2 design contract)
+        seed_g, ok_g = IMAS._seed_omp(itp2, c, RA, ZA, maximum(rr))
+        @test ok_g
+        vRg, vZg, vcg = IMAS._trace_surface_cubic(itp2, c, seed_g; τ_grad=0.15, h_max=0.02)
+        @test vcg
+        @test all(abs(itp2(vRg[k], vZg[k]) - c) < 1e-9 for k in eachindex(vRg))
+        # demonstrate the guard does something: paths with guard on vs off differ
+        R0g, Z0g, cl0g = IMAS.trace_surface_cubic(itp2, c, RA, ZA, maximum(rr);
+                                                    τ_grad=0.0, h_max=0.02)
+        # either the paths differ (guard changed the trace) or guard-off also closes (smoke test)
+        @test cl0g || closed1   # at minimum one closes; if both close, paths may or may not differ
+        # the guard fires (the predicate τ_grad=0.15 > min|∇ψ| near separatrix was verified to fire)
+    end
+
+    @testset "trace_surfaces_cubic drop-in matches trace_surfaces (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp2 = IMAS.ψ_interpolant(eqt2d)
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        e1 = eqt.profiles_1d
+        pa, pb = itp2(RA, ZA), e1.psi[end]
+        # nudge the exact-separatrix last level a hair inside so BOTH paths close it (neither
+        # path traces the open X-point separatrix without wall-clipping)
+        psis = collect(e1.psi); psis[end] = pa + 0.999 * (pb - pa)
+        ff = collect(e1.f)
+
+        # low-level drop-in: same args as trace_surfaces (minus the precomputed grid/field arrays),
+        # refine_extrema on both -> the refined geometric extrema agree to ~1e-7 m (both Newton-refine
+        # to the same point on the same interpolant)
+        ref = IMAS.trace_surfaces(psis, ff, collect(rr), collect(zz), eqt2d.psi, itp2, RA, ZA, fw.r, fw.z; refine_extrema=true)
+        cub = IMAS.trace_surfaces_cubic(psis, ff, RA, ZA, maximum(rr), itp2; refine_extrema=true)
+        @test length(cub) == length(ref)
+        for k in (16, 32, 48)                       # mid-radius surfaces
+            @test isapprox(cub[k].max_r, ref[k].max_r; atol=1e-4)
+            @test isapprox(cub[k].min_r, ref[k].min_r; atol=1e-4)
+            @test isapprox(cub[k].max_z, ref[k].max_z; atol=1e-4)
+            @test isapprox(cub[k].min_z, ref[k].min_z; atol=1e-4)
+            gm1_c = IMAS.flux_surface_avg(1.0 ./ cub[k].r .^ 2, cub[k])
+            gm1_r = IMAS.flux_surface_avg(1.0 ./ ref[k].r .^ 2, ref[k])
+            @test isapprox(gm1_c, gm1_r; rtol=2e-3)
+        end
+
+        # high-level overload mirrors trace_surfaces(eqt, wall_r, wall_z): runs the full profile;
+        # the raw separatrix boundary can't close, so it warns and uses an inner proxy (assert the warn)
+        cub_hl = @test_logs (:warn,) match_mode = :any IMAS.trace_surfaces_cubic(eqt, fw.r, fw.z)
+        @test length(cub_hl) == length(e1.psi)
+        @test cub_hl[32].max_r > RA
+    end
+end
diff --git a/test/runtests_refine_extremum.jl b/test/runtests_refine_extremum.jl
new file mode 100644
index 00000000..2d1458dc
--- /dev/null
+++ b/test/runtests_refine_extremum.jl
@@ -0,0 +1,159 @@
+using IMAS
+using Test
+import Interpolations
+
+# The robust refine is backend-agnostic via the _psi_* / _value_gradient / _hessian! adapters: it
+# must work not only with a FastInterpolations interpolant but also with an Interpolations.jl one
+# (e.g. as FRESCO builds). Analytic field below: nested ellipses ψ = ((R-R0)/a0)^2 + (Z/b0)^2.
+@testset "robust refine backend-agnostic (Interpolations.jl)" begin
+    @testset "analytic ellipse, IP interpolant" begin
+        R0, a0, b0 = 1.7, 0.5, 1.0
+        psi_fun = (R, Z) -> ((R - R0) / a0)^2 + (Z / b0)^2
+        r = range(1.0, 2.4, length=65)
+        z = range(-1.3, 1.3, length=65)
+        PSI = [psi_fun(ri, zj) for ri in r, zj in z]
+        itp = Interpolations.cubic_spline_interpolation((r, z), PSI; extrapolation_bc=Interpolations.Line())
+        @test !(itp isa IMAS.FI.AbstractInterpolant)   # genuinely exercises the non-FI path
+        L = 0.36; s = sqrt(L); dc = step(r)
+        cases = (
+            (:R, (R0 + a0 * s, 0.0), (R0 + a0 * s - 0.4dc, 0.3dc)),
+            (:R, (R0 - a0 * s, 0.0), (R0 - a0 * s + 0.4dc, 0.3dc)),
+            (:Z, (R0, b0 * s), (R0 + 0.3dc, b0 * s - 0.4dc)),
+            (:Z, (R0, -b0 * s), (R0 + 0.3dc, -b0 * s + 0.4dc)),
+        )
+        for (extremum_of, truth, seed) in cases
+            R, Z = IMAS._robust_refine_extremum(itp, L, seed, extremum_of, (R0, 0.0))
+            @test isapprox(R, truth[1]; atol=1e-3)
+            @test isapprox(Z, truth[2]; atol=1e-3)
+        end
+    end
+
+    @testset "DIII-D trace_surfaces: Interpolations matches FastInterpolations" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp_fi = IMAS.ψ_interpolant(eqt2d)
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        psi_axis = itp_fi(RA, ZA)
+        eqt1d = eqt.profiles_1d
+        fracs = collect(range(0.1, 0.95; length=16))
+        psis = [psi_axis; psi_axis .+ fracs .* (eqt1d.psi[end] - psi_axis)]
+        ff = fill(eqt1d.f[end], length(psis))
+        # FRESCO-style: an Interpolations.jl cubic interpolant of the SAME ψ grid as itp_fi
+        itp_ip = Interpolations.cubic_spline_interpolation((rr, zz), collect(eqt2d.psi); extrapolation_bc=Interpolations.Line())
+        @test !(itp_ip isa IMAS.FI.AbstractInterpolant)
+        r = collect(rr); z = collect(zz); PSI = collect(eqt2d.psi)
+        # identical marching trace (same PSI matrix); only the refine backend differs
+        s_fi = IMAS.trace_surfaces(psis, ff, r, z, PSI, itp_fi, RA, ZA, fw.r, fw.z; refine_extrema=true)
+        s_ip = IMAS.trace_surfaces(psis, ff, r, z, PSI, itp_ip, RA, ZA, fw.r, fw.z; refine_extrema=true)
+        @test length(s_ip) == length(s_fi)
+        for k in 2:length(s_fi)
+            @test isapprox(s_ip[k].max_r, s_fi[k].max_r; atol=2e-3)
+            @test isapprox(s_ip[k].min_r, s_fi[k].min_r; atol=2e-3)
+            @test isapprox(s_ip[k].max_z, s_fi[k].max_z; atol=2e-3)
+            @test isapprox(s_ip[k].min_z, s_fi[k].min_z; atol=2e-3)
+        end
+    end
+end
+
+# Production tracing uses _robust_refine_extremum!: from ANY seed on the correct side of the
+# axis (even one outside the separatrix or in the private flux region across an X-point) it must
+# still arrive at the correct CONFINED extremum, never an off-surface / private / SOL point.
+# This is what guards the KDEMO (a_eq<0) and MANTA (elongation≈0 -> sqrt DomainError) failures.
+@testset "_robust_refine_extremum arrives from arbitrary/bad seeds" begin
+    @testset "far correct-side seed (analytic ellipse, no X-point)" begin
+        R0, a0, b0 = 1.7, 0.5, 1.0
+        psi_fun = (R, Z) -> ((R - R0) / a0)^2 + (Z / b0)^2
+        r = range(1.0, 2.4, length=65)
+        z = range(-1.3, 1.3, length=65)
+        itp = IMAS.ψ_interpolant(r, z, [psi_fun(ri, zj) for ri in r, zj in z]).PSI_interpolant
+        L = 0.36; s = sqrt(L)
+        # seed far OUTSIDE the L surface (ψ≈1.5≫L) but on the outboard side -> still finds max_r
+        R, Z = IMAS._robust_refine_extremum(itp, L, (R0 + 0.6, 0.2), :R, (R0, 0.0))
+        @test isapprox(R, R0 + a0 * s; atol=1e-4)
+        @test isapprox(Z, 0.0; atol=1e-4)
+        # seed far above the L surface on the upper side -> still finds max_z
+        R, Z = IMAS._robust_refine_extremum(itp, L, (R0 + 0.1, 1.2), :Z, (R0, 0.0))
+        @test isapprox(R, R0; atol=1e-4)
+        @test isapprox(Z, b0 * s; atol=1e-4)
+        # invalid extremum_of is rejected
+        @test_throws ArgumentError IMAS._robust_refine_extremum(itp, L, (R0, 0.0), :bogus, (R0, 0.0))
+    end
+
+    @testset "outside-separatrix and private-region seeds (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp = IMAS.ψ_interpolant(eqt2d)
+        eqt1d = eqt.profiles_1d
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        axis = (RA, ZA)
+        psi_axis = itp(RA, ZA)
+        pb = IMAS.find_psi_boundary(rr, zz, eqt2d.psi, psi_axis, eqt1d.psi[end], RA, ZA, fw.r, fw.z;
+            PSI_interpolant=itp, raise_error_on_not_open=false, raise_error_on_not_closed=false)
+        psis = (eqt1d.psi .- eqt1d.psi[1]) ./ (eqt1d.psi[end] - eqt1d.psi[1]) .* (pb.last_closed - psi_axis) .+ psi_axis
+        rough = IMAS.trace_surfaces(psis, eqt1d.f, collect(rr), collect(zz), eqt2d.psi, itp, RA, ZA, fw.r, fw.z; refine_extrema=false)
+        xpoints = [(xp.r, xp.z) for xp in eqt.boundary.x_point]
+        xp_up = eqt.boundary.x_point[argmax(xp.z for xp in eqt.boundary.x_point)]
+        dr, dz = step(rr), step(zz)
+        lo = (first(rr), first(zz)); hi = (last(rr), last(zz))   # clamp to ψ grid box (as production does)
+        N = length(rough)
+        psiN(p) = (itp(p[1], p[2]) - psi_axis) / (pb.last_closed - psi_axis)
+        for k in (N, N - 2, N - 5)
+            s = rough[k]
+            c = psis[k]
+            good_r = IMAS._robust_refine_extremum(itp, c, (s.max_r, s.z_at_max_r), :R, axis, xpoints; lo, hi)
+            good_z = IMAS._robust_refine_extremum(itp, c, (s.r_at_max_z, s.max_z), :Z, axis, xpoints; lo, hi)
+            # bad seed A: outside the separatrix (ψ_N>1), far outboard
+            outside = (s.max_r + 12dr, s.z_at_max_r)
+            @test psiN(outside) > 1
+            rA = IMAS._robust_refine_extremum(itp, c, outside, :R, axis, xpoints; lo, hi)
+            @test isapprox(rA[1], good_r[1]; atol=5dr)
+            @test psiN(rA) <= 1 + 1e-6
+            # bad seed B: private region (reflect the rough max_z across the upper X-point)
+            priv = (2xp_up.r - s.r_at_max_z, 2xp_up.z - s.max_z)
+            rB = IMAS._robust_refine_extremum(itp, c, priv, :Z, axis, xpoints; lo, hi)
+            @test isapprox(rB[2], good_z[2]; atol=5dz)
+            @test rB[2] < xp_up.z
+        end
+    end
+
+    # Hardest case: at ψ_N=0.999 the PRIVATE region across the upper X-point ALSO has a
+    # ψ_N=0.999 point with ∂ψ/∂R=0 — a genuine solution of the SAME {ψ=c, ∂ψ/∂R=0} system.
+    # Seeds at / above the X-point must NOT be trapped there; they must return the confined
+    # max_z below the X-point. (This is the geometry behind the MANTA elongation≈0 failure.)
+    @testset "not trapped in the private lobe at ψ_N=0.999 (DIII-D)" begin
+        filename = joinpath(pkgdir(IMAS.IMASdd), "sample", "D3D_eq_ods.json")
+        dd = IMAS.json2imas(filename; show_warnings=false)
+        eqt = dd.equilibrium.time_slice[1]
+        fw = IMAS.first_wall(dd.wall)
+        eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
+        rr, zz, itp = IMAS.ψ_interpolant(eqt2d)
+        eqt1d = eqt.profiles_1d
+        RA, ZA = eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z
+        axis = (RA, ZA)
+        psi_axis = itp(RA, ZA)
+        pb = IMAS.find_psi_boundary(rr, zz, eqt2d.psi, psi_axis, eqt1d.psi[end], RA, ZA, fw.r, fw.z;
+            PSI_interpolant=itp, raise_error_on_not_open=false, raise_error_on_not_closed=false)
+        xpoints = [(xp.r, xp.z) for xp in eqt.boundary.x_point]
+        xp_up = eqt.boundary.x_point[argmax(xp.z for xp in eqt.boundary.x_point)]
+        lo = (first(rr), first(zz)); hi = (last(rr), last(zz))
+        c = psi_axis + 0.999 * (pb.last_closed - psi_axis)   # ψ_N = 0.999, just inside the separatrix
+        psiN(p) = (itp(p[1], p[2]) - psi_axis) / (pb.last_closed - psi_axis)
+        # confined truth: seed just below the upper X-point
+        truth = IMAS._robust_refine_extremum(itp, c, (xp_up.r, xp_up.z - 0.15), :Z, axis, xpoints; lo, hi)
+        @test truth[2] < xp_up.z
+        @test isapprox(psiN(truth), 0.999; atol=2e-3)
+        # seeds AT / ABOVE the upper X-point (in / near the private lobe) must come back confined
+        for sd in ((xp_up.r, xp_up.z + 0.05), (xp_up.r + 0.05, xp_up.z + 0.03),
+                   (xp_up.r - 0.05, xp_up.z + 0.05), (xp_up.r, xp_up.z + 0.25))
+            res = IMAS._robust_refine_extremum(itp, c, sd, :Z, axis, xpoints; lo, hi)
+            @test res[2] < xp_up.z                         # not trapped above the X-point
+            @test isapprox(res[2], truth[2]; atol=1e-4)    # the genuine confined max_z
+        end
+    end
+end