diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
index 32aa4c51..979621fa 100644
--- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
+++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
@@ -77,11 +77,13 @@ Read `T_total_*` and `dW_total_*` scalar reference values from the Fortran
 function read_pentrc_reference(pentrc_nc::String)
     NCDatasets.Dataset(pentrc_nc, "r") do ds
         gatt = ds.attrib
+        # Reference runs may enable only a subset of methods (e.g. fgar only); missing attributes → NaN.
+        att(k) = haskey(gatt, k) ? Float64(gatt[k]) : NaN
         return (
-            T_total_fgar=Float64(gatt["T_total_fgar"]),
-            T_total_tgar=Float64(gatt["T_total_tgar"]),
-            dW_total_fgar=Float64(gatt["dW_total_fgar"]),
-            dW_total_tgar=Float64(gatt["dW_total_tgar"])
+            T_total_fgar=att("T_total_fgar"),
+            T_total_tgar=att("T_total_tgar"),
+            dW_total_fgar=att("dW_total_fgar"),
+            dW_total_tgar=att("dW_total_tgar")
         )
     end
 end
@@ -142,7 +144,6 @@ nn_high = 1
 delta_mlow = 8
 delta_mhigh = 8
 mthvac = 512
-thmax0 = 1
 
 kinetic_source = "fixed"
 kinetic_factor = 0.0
@@ -264,15 +265,14 @@ function run_benchmark(fortran_dir::String=DEFAULT_FORTRAN_DIR)
     _p("\n--- Build PE state (Fortran Clebsch → dbob_m/divx_m) ---")
     t1 = time()
 
-    # Construct PerturbedEquilibriumState with Fortran Clebsch data
-    # Fortran xclebsch stores ξ^α directly; PE convention is ξ^α/χ₁
-    chi1 = 2π * equil.psio
+    # Construct PerturbedEquilibriumState with Fortran Clebsch data.
+    # Fortran gpec_xclebsch already writes ξ^α = xms/χ₁, exactly what the KF operator consumes — feed directly.
     pe_state = PE.PerturbedEquilibriumState(;
         psi_grid=psi_grid_f,
         xi_modes=(
             clebsch_psi=clebsch_psi,
             clebsch_psi1=clebsch_psi1,
-            clebsch_alpha=clebsch_alpha ./ chi1  # Store as ξ^α/χ₁ per PE convention
+            clebsch_alpha=clebsch_alpha
         )
     )
 
diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml
index 6693b560..17f6ff58 100644
--- a/regression-harness/cases/diiid_n1.toml
+++ b/regression-harness/cases/diiid_n1.toml
@@ -360,6 +360,17 @@ label = "NTV torque FGAR (Re,Im) [N·m]"
 noise_threshold = 1e-2
 order = 95
 
+# NTV kinetic energy (δW) through the same post-PE FGAR path. This is the
+# quantity most sensitive to the bounce-averaging chain (quadrature scheme,
+# v_par source, bo normalization), so it guards that whole path directly.
+[quantities.ntv_energy_fgar]
+h5path = "kinetic_forces/fgar/total_energy"
+type = "real_vector"
+extract = "all_real"
+label = "NTV energy FGAR (Re,Im) [J]"
+noise_threshold = 1e-3
+order = 96
+
 # Profiles (checksums)
 [quantities.q_profile]
 h5path = "splines/profiles/q"
diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl
index e0153fcb..4d1e7c94 100644
--- a/src/KineticForces/BounceAveraging.jl
+++ b/src/KineticForces/BounceAveraging.jl
@@ -146,10 +146,10 @@ end
 # ============================================================================
 
 """
-    compute_bounce_data(psi, n, l, q, bo, bmax, bmin, ibmax, theta_bmax,
-                        tspl, mfac, chi1, ro, dbob_m_f, divx_m_f,
+    compute_bounce_data(psi, n, l, q, bo, bmax, bmin, theta_bmax,
+                        tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f,
                         divxfac, wdfac, mass, chrg, T_s, method;
-                        nlmda=64, ntheta=128,
+                        nlmda=128, ntheta=128,
                         smat=nothing, tmat=nothing, xmat=nothing,
                         ymat=nothing, zmat=nothing) → BounceData
 
@@ -166,9 +166,10 @@ Ports Fortran torque.F90 lines 530-816 (GAR branch).
 - `q`: Safety factor at this ψ
 - `bo`: On-axis toroidal field [T]
 - `bmax, bmin`: Max/min of B(θ) at this ψ
-- `ibmax`: Index of Bmax in poloidal grid
-- `theta_bmax`: θ location of Bmax
-- `tspl`: Poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ]
+- `theta_bmax`: θ location of Bmax (nodal knot; the passing-transit start)
+- `tspl`: Periodic poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ]
+- `B_extrap`: Endpoint-fit (non-periodic) cubic of B(θ) used for v_par and the
+  bounce-point roots (the Fortran `vspl` equivalent)
 - `mfac`: Poloidal mode numbers [mlow:mhigh]
 - `chi1`: 2π·ψ₀ flux normalization
 - `ro`: Major radius [m]
@@ -181,20 +182,20 @@ Ports Fortran torque.F90 lines 530-816 (GAR branch).
 - `method`: Method string (first char: f/t/p determines λ range)
 
 # Keyword Arguments
-- `nlmda`: Number of pitch angle grid points (default 64)
+- `nlmda`: Number of pitch angle grid points (default 128, matching Fortran pentrc nlmda)
 - `ntheta`: Number of poloidal grid points per bounce (default 128)
 - `smat, tmat, xmat, ymat, zmat`: Geometric matrices (mpert×mpert) for kinetic matrix path
 """
 function compute_bounce_data(
     psi::Float64, n::Int, l::Int, q::Float64,
     bo::Float64, bmax::Float64, bmin::Float64,
-    ibmax::Int, theta_bmax::Float64,
-    tspl, mfac::Vector{Int}, chi1::Float64, ro::Float64,
+    theta_bmax::Float64,
+    tspl, B_extrap, mfac::Vector{Int}, chi1::Float64, ro::Float64,
     dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64},
     divxfac::Float64, wdfac::Float64,
     mass::Float64, chrg::Float64,
     T_s::Float64, method::String;
-    nlmda::Int=64, ntheta::Int=128,
+    nlmda::Int=128, ntheta::Int=128,
     smat::Union{Nothing,Matrix{ComplexF64}}=nothing,
     tmat::Union{Nothing,Matrix{ComplexF64}}=nothing,
     xmat::Union{Nothing,Matrix{ComplexF64}}=nothing,
@@ -240,13 +241,12 @@ function compute_bounce_data(
 
         # Find bounce points and build θ sub-grid
         _, _, tdt_pts, tdt_wts = _find_bounce_points_and_grid(
-            lmda, bo, sigma, tspl, ibmax, theta_bmax,
-            lmdatpb, lmdamax, psi, ntheta, scr)
+            lmda, bo, sigma, B_extrap, theta_bmax, psi, ntheta)
 
         # Bounce integrals over θ (Fortran lines 674-735)
         wbbar, wdbar, dJdJ_val, wmats_lmda = _bounce_integrate(
             tdt_pts, tdt_wts, lmda, lnq, sigma, n, q, bo,
-            tspl, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac,
+            tspl, B_extrap, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac,
             do_matrices, mpert, smat, tmat, xmat, ymat, zmat, scr)
 
         # Physical frequencies (Fortran lines 744-745)
@@ -255,8 +255,8 @@ function compute_bounce_data(
         dJdJ_arr[ilmda] = dJdJ_val
 
         if do_matrices && !isnothing(wmats_lmda)
-            @inbounds for q in 1:length(wmats_lmda)
-                wmats_arr[ilmda, q] = wmats_lmda[q]
+            @inbounds for iq in eachindex(wmats_lmda)
+                wmats_arr[ilmda, iq] = wmats_lmda[iq]
             end
         end
     end
@@ -267,16 +267,17 @@ end
 
 
 """
-    BounceScratch(ntheta, mpert; nfine=256)
+    BounceScratch(ntheta, mpert)
 
-Per-surface scratch buffers for the bounce-averaging inner loops. Allocated once in
-`compute_bounce_data` and reused across every λ (each λ previously reallocated these
-θ- and mode-sized arrays). Sizes are fixed for a flux surface (`ntheta` sub-grid points,
-`mpert` Fourier modes), so a single set of buffers serves the whole λ sweep. Buffers that
-the callers zero-initialize are `fill!`-reset per λ, preserving bit-for-bit results.
+Per-surface scratch for the bounce-averaging inner loops, allocated once in
+`compute_bounce_data` and reused across every λ. Sizes are fixed for a flux surface
+(`ntheta` sub-grid points, `mpert` Fourier modes). Buffers the loops populate only
+partially are `fill!`-reset per λ.
 
 ## Fields
-- `cum_wb_arr::Vector{Float64}`: length `ntheta` — cumulative bounce action
+- `g_wb::Vector{Float64}`: length `ntheta` — bounce-action integrand samples
+- `g_wd::Vector{Float64}`: length `ntheta` — drift integrand samples
+- `cum_wb_arr::Vector{Float64}`: length `ntheta` — cumulative bounce-action integral
 - `jvtheta::Vector{ComplexF64}`: length `ntheta` — action integrand
 - `wmu_mt::Matrix{ComplexF64}`: `mpert × ntheta` — W_μ per θ
 - `wen_mt::Matrix{ComplexF64}`: `mpert × ntheta` — W_E per θ
@@ -286,10 +287,12 @@ the callers zero-initialize are `fill!`-reset per λ, preserving bit-for-bit res
 - `wen_ba::Vector{ComplexF64}`: length `mpert` — bounce-averaged W_E
 - `wmats_lmda::Vector{ComplexF64}`: length `nqty_matrix(mpert)` — packed W outer products
 - `tspl_f::Vector{Float64}`: length 5 — in-place tspl(θ) evaluation
-- `vpar_fine::Vector{Float64}`: length `nfine+1` — v_par(θ) for bounce-point search;
-  `_find_bounce_points_and_grid` derives its fine-grid size from this buffer's length
+- `itp_r::TR`, `itp_c::TC`: real/complex fit-and-integrate interpolants on the unit
+  θ-quadrature grid (`CubicFit` endpoints); see `_refit!`
 """
-struct BounceScratch
+struct BounceScratch{TR, TC}
+    g_wb::Vector{Float64}
+    g_wd::Vector{Float64}
     cum_wb_arr::Vector{Float64}
     jvtheta::Vector{ComplexF64}
     wmu_mt::Matrix{ComplexF64}
@@ -300,11 +303,19 @@ struct BounceScratch
     wen_ba::Vector{ComplexF64}
     wmats_lmda::Vector{ComplexF64}
     tspl_f::Vector{Float64}
-    vpar_fine::Vector{Float64}
+    itp_r::TR
+    itp_c::TC
 end
 
-function BounceScratch(ntheta::Int, mpert::Int; nfine::Int=256)
+function BounceScratch(ntheta::Int, mpert::Int)
+    # Reusable fit-and-integrate interpolants on the unit θ-quadrature grid;
+    # per λ their y contents are overwritten and z refit in place.
+    xs = range(0.0, 1.0, length=ntheta)
+    itp_r = cubic_interp(xs, zeros(Float64, ntheta); bc=CubicFit())
+    itp_c = cubic_interp(xs, zeros(ComplexF64, ntheta); bc=CubicFit())
     return BounceScratch(
+        Vector{Float64}(undef, ntheta),
+        Vector{Float64}(undef, ntheta),
         Vector{Float64}(undef, ntheta),
         Vector{ComplexF64}(undef, ntheta),
         Matrix{ComplexF64}(undef, mpert, ntheta),
@@ -315,10 +326,21 @@ function BounceScratch(ntheta::Int, mpert::Int; nfine::Int=256)
         Vector{ComplexF64}(undef, mpert),
         Vector{ComplexF64}(undef, nqty_matrix(mpert)),
         Vector{Float64}(undef, 5),
-        Vector{Float64}(undef, nfine + 1),
+        itp_r,
+        itp_c,
     )
 end
 
+"""
+Refit `itp` in place to its current `y` contents, reusing the cached tridiagonal
+factorization. Uses the package-internal `_solve_system!` because no exported API
+refits an interpolant without allocating; this keeps the per-λ hot loop allocation-free.
+"""
+@inline function _refit!(itp)
+    FastInterpolations._solve_system!(itp.z, itp.cache, itp.y, itp.bc)
+    return itp
+end
+
 
 # ============================================================================
 # Internal helpers
@@ -358,36 +380,30 @@ function _build_lambda_grid(method_char::Char, lmdatpb::Float64, lmdamax::Float6
 end
 
 
+"""
+Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the endpoint-fit cubic of B
+(`B_extrap`, built where the surface interpolants are constructed), keeping v_par
+consistent with the bounce-point roots as in Fortran's `vspl`.
+"""
+@inline _vpar_from_extrap(B_extrap, lmda::Float64, bo::Float64, θ::Float64) =
+    1.0 - (lmda / bo) * B_extrap(mod(θ, 1.0))
+
+
 """
 Find bounce points for trapped/passing particles and build θ sub-grid.
 Returns (t1, t2, theta_points, theta_weights).
 """
 function _find_bounce_points_and_grid(
     lmda::Float64, bo::Float64, sigma::Int,
-    tspl, ::Int, theta_bmax::Float64,
-    ::Float64, ::Float64, psi::Float64,
-    ntheta::Int, scr::BounceScratch
+    B_extrap, theta_bmax::Float64, psi::Float64,
+    ntheta::Int
 )
     if sigma == 0  # trapped
-        # Build v_par(θ) = 1 - (λ/bo)*B(θ) and find roots on a dense θ grid.
-        # Grid size derives from the scratch buffer so it stays in sync with BounceScratch's nfine.
-        nfine = length(scr.vpar_fine) - 1
-        theta_fine = range(0.0, 1.0, length=nfine+1)
-        vpar_fine = scr.vpar_fine
-        @inbounds for i in 1:(nfine+1)
-            tspl(scr.tspl_f, mod(theta_fine[i], 1.0))
-            vpar_fine[i] = 1.0 - (lmda / bo) * scr.tspl_f[1]
-        end
-
-        # Find zero crossings
-        bpts = Float64[]
-        for i in 1:nfine
-            if vpar_fine[i] * vpar_fine[i+1] < 0
-                # Bisect for better accuracy
-                θ_root = _bisect_vpar(tspl, scr.tspl_f, lmda, bo, theta_fine[i], theta_fine[i+1])
-                push!(bpts, θ_root)
-            end
-        end
+        # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_extrap(θ) in (0,1),
+        # sorted descending — the same order as Fortran spline_roots, which the
+        # marginally-trapped and deepest-well wrap logic below assume.
+        vpar_fn = θ -> _vpar_from_extrap(B_extrap, lmda, bo, θ)
+        bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0); rev=true)
 
         nbpts = length(bpts)
         if nbpts < 1
@@ -399,8 +415,7 @@ function _find_bounce_points_and_grid(
             t1 = bpts[1]
             t2 = bpts[1] + 1.0
         else
-            # Find deepest potential well (Fortran lines 616-639)
-            t1, t2 = _find_deepest_well(bpts, tspl, scr.tspl_f, lmda, bo)
+            t1, t2 = _find_deepest_well(bpts, B_extrap, lmda, bo)
         end
 
         # Power-law grid refined near bounce points
@@ -417,57 +432,10 @@ end
 
 
 """
-    _bisect_vpar(tspl, lmda, bo, θa, θb; tol=1e-12, maxiter=50) → θ
-
-Bisect on `θ ∈ [θa, θb]` to find a bounce point — the angle where the parallel
-velocity vanishes, `v_par(θ) = 1 - (λ/bo) · B(θ) = 0`. The caller is responsible
-for supplying a bracket where `v_par(θa)` and `v_par(θb)` have opposite signs;
-this routine does not check, it just halves toward the sign change.
-
-# Arguments
-- `tspl`: 1D θ-interpolant returning at least `[B(θ), …]` at first index;
-  `tspl(mod(θ, 1.0))[1]` extracts the local field magnitude.
-- `lmda`: pitch-angle parameter λ = μ·bo / E.
-- `bo`: on-axis toroidal field used to normalise λ.
-- `θa, θb`: bracket endpoints (normalised θ ∈ [0,1)); typically straddling a B-peak.
-
-# Keyword arguments
-- `tol`: termination tolerance — exits when either `|v_par(θ_mid)| < tol` or the
-  bracket width `θb − θa < tol`. Default `1e-12` is tight enough that the
-  residual `v_par` is dominated by the spline-evaluation roundoff of `tspl`.
-- `maxiter`: hard iteration cap. Default `50` halves the bracket by ~10⁻¹⁵, well
-  past `tol` on a unit-scale bracket; the cap is a runaway guard rather than the
-  expected exit. On exhaustion the midpoint of the final bracket is returned.
-
-# Returns
-- `θ::Float64`: the converged (or capped) bounce-point angle.
+Find the deepest potential well (largest midpoint v_par) among bounce-point pairs,
+handling pairs that wrap through θ = 0/1.
 """
-function _bisect_vpar(tspl, tspl_f::Vector{Float64}, lmda::Float64, bo::Float64, θa::Float64, θb::Float64; tol=1e-12, maxiter=50)
-    tspl(tspl_f, mod(θa, 1.0))
-    va = 1.0 - (lmda / bo) * tspl_f[1]
-    for _ in 1:maxiter
-        θm = 0.5 * (θa + θb)
-        tspl(tspl_f, mod(θm, 1.0))
-        vm = 1.0 - (lmda / bo) * tspl_f[1]
-        if abs(vm) < tol || (θb - θa) < tol
-            return θm
-        end
-        if va * vm < 0
-            θb = θm
-        else
-            θa = θm
-            va = vm
-        end
-    end
-    return 0.5 * (θa + θb)
-end
-
-
-"""
-Find deepest potential well among bounce point pairs.
-Ports Fortran lines 616-639.
-"""
-function _find_deepest_well(bpts::Vector{Float64}, tspl, tspl_f::Vector{Float64}, lmda::Float64, bo::Float64)
+function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo::Float64)
     nbpts = length(bpts)
     best_vpar = 0.0
     best_t1 = 0.0
@@ -481,8 +449,7 @@ function _find_deepest_well(bpts::Vector{Float64}, tspl, tspl_f::Vector{Float64}
         else
             θmid = 0.5 * (bpts[i] + bpts[j])
         end
-        tspl(tspl_f, mod(θmid, 1.0))
-        vpar_mid = 1.0 - (lmda / bo) * tspl_f[1]
+        vpar_mid = _vpar_from_extrap(B_extrap, lmda, bo, θmid)
         if vpar_mid > best_vpar
             best_t1 = bpts[i]
             best_t2 = bpts[j]
@@ -511,7 +478,7 @@ Ports Fortran torque.F90 lines 674-793.
 function _bounce_integrate(
     tdt_pts::Vector{Float64}, tdt_wts::Vector{Float64},
     lmda::Float64, lnq::Float64, sigma::Int, n::Int, q::Float64, bo::Float64,
-    tspl, chi1::Float64, ro::Float64,
+    tspl, B_extrap, chi1::Float64, ro::Float64,
     mfac::Vector{Int}, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64},
     divxfac::Float64, wdfac::Float64,
     do_matrices::Bool, mpert::Int,
@@ -520,15 +487,12 @@ function _bounce_integrate(
     ntheta = length(tdt_pts)
     theta0 = tdt_pts[1]
 
-    # Cumulative bounce integrals
-    cum_wb = 0.0
-    cum_wd = 0.0
-
-    # θ-scratch arrays reused across λ from the per-surface BounceScratch. Zero-reset
-    # the ones the loop populates only partially (i in 2:ntheta-1, with continue/break),
-    # so unwritten entries stay 0 exactly as the previous `zeros(...)` did.
-    cum_wb_arr = scr.cum_wb_arr
-    fill!(cum_wb_arr, 0.0)
+    # θ-sample integrands (zero-reset per λ: the loop populates 2:ntheta-1 with
+    # continue/break paths that rely on unwritten entries staying 0)
+    g_wb = scr.g_wb    # J·B/√v_par · dθ/dx
+    g_wd = scr.g_wd    # drift integrand · dθ/dx
+    fill!(g_wb, 0.0)
+    fill!(g_wd, 0.0)
 
     # Action integrand
     jvtheta = scr.jvtheta
@@ -558,39 +522,40 @@ function _bounce_integrate(
         jac = tspl_f[4]
         djdpsi = tspl_f[5]
 
-        vpar = 1.0 - (lmda / bo) * B_val
+        # v_par from the endpoint-fit cubic (consistent with the bounce points);
+        # the periodic tspl B_val remains the numerator field in the integrands.
+        vpar = 1.0 - (lmda / bo) * B_extrap(θmod)
 
         if vpar <= 0
-            # Zero crossing near bounce points — handle like Fortran lines 678-697
+            # Negative v_par near a bounce point: same fill rules as the Fortran bounce loop.
             if i < ntheta ÷ 2
-                # Before midpoint: zero out previous entries
+                # Before midpoint: restart — zero everything up to this sample.
+                fill!(view(g_wb, 1:i), 0.0)
+                fill!(view(g_wd, 1:i), 0.0)
+                fill!(view(jvtheta, 1:i), ComplexF64(0.0))
                 continue
             else
-                # After midpoint: clamp remaining entries
+                # After midpoint: hold the previous sample to the end of the grid.
+                # The wd slot is deliberately held from the wb integrand (g_wb),
+                # reproducing the Fortran behavior for parity.
+                fill!(view(g_wb, i:ntheta), g_wb[i-1])
+                fill!(view(g_wd, i:ntheta), g_wb[i-1])
+                fill!(view(jvtheta, i:ntheta), jvtheta[i-1])
                 break
             end
         end
 
         sqrt_vpar = sqrt(vpar)
 
-        # Bounce integrands (Fortran lines 698-701)
-        wb_integrand = dt * jac * B_val / sqrt_vpar
-        wd_integrand = dt * jac * dBdpsi * (1.0 - 1.5 * lmda * B_val / bo) / sqrt_vpar +
-                       dt * djdpsi * B_val * sqrt_vpar
+        # Bounce integrands
+        g_wb[i] = dt * jac * B_val / sqrt_vpar
+        g_wd[i] = dt * jac * dBdpsi * (1.0 - 1.5 * lmda * B_val / bo) / sqrt_vpar +
+                  dt * djdpsi * B_val * sqrt_vpar
 
-        cum_wb += wb_integrand
-        cum_wd += wd_integrand
-        # Trapezoidal cumulative (matches Fortran spline_int semantics on linear fn):
-        # bspl%fsi(j)/Δx = g_1 + ... + g_{j-1} + g_j/2, so subtract half the current sample.
-        cum_wb_arr[i] = cum_wb - wb_integrand / 2
-
-        # Fourier modes at this θ (Fortran lines 702-708) — write into pre-allocated
-        # expm buffer using the ORIGINAL expression order to preserve bit-level parity.
+        # Fourier modes at this θ
         @inbounds for mi in 1:mpert
             expm[mi] = cis(twopi * mfac[mi] * θ)
         end
-        # Replaces `sum(dbob_m_f .* expm)` / `sum(divx_m_f .* expm) * divxfac`
-        # with direct accumulators; same evaluation order as the broadcast + sum.
         dbob = ComplexF64(0.0)
         divx = ComplexF64(0.0)
         @inbounds for mi in 1:mpert
@@ -599,15 +564,13 @@ function _bounce_integrate(
         end
         divx *= divxfac
 
-        # Action integrand (Fortran line 706-708)
+        # Action integrand
         phase = cis(-twopi * n * q * (θ - theta0))
         jvtheta[i] = dt * jac * B_val *
             (divx * sqrt_vpar + dbob * (1.0 - 1.5 * lmda * B_val / bo) / sqrt_vpar) *
             phase
 
-        # W vectors for matrix path (Fortran lines 722-727). Element-by-element
-        # write preserving original broadcast evaluation order exactly to keep
-        # bit-level parity (matters because downstream quadrature is tolerance-sensitive).
+        # W vectors for matrix path
         if do_matrices
             wmu_pre = dt * (lmda / bo)
             wen_pre = dt
@@ -616,78 +579,80 @@ function _bounce_integrate(
                 wen_mt[mi, i] = wen_pre * expm[mi] / (B_val * sqrt_vpar) * phase / (2 * chi1)
             end
         end
+
+        # Smooth backfill for points zeroed before a restart. Exact equality with
+        # 0.0 is safe: the restart branch set these entries with fill!(…, 0.0).
+        if i >= 3 && g_wb[i-1] == 0.0
+            fill!(view(g_wb, 3:i-1), g_wb[i])
+            fill!(view(g_wd, 3:i-1), g_wd[i])
+            fill!(view(jvtheta, 2:i-1), jvtheta[i])
+        end
     end
 
-    # Total bounce integrals — Fortran splines over bspl%xs = linspace(0,1,ntheta)
-    # and integrates via spline_int, which is ≈ (1/(ntheta-1)) × Σ f_i. The tdt(2,i)
-    # weights contain dθ/dx so Σ tdt·f is raw Riemann; divide by (ntheta-1) to get
-    # the integral over the unit linear space [0,1] that Fortran produces.
-    nrm = 1.0 / (ntheta - 1)
-    total_wb = cum_wb * nrm
-    total_wd = cum_wd * nrm
+    # Total bounce integrals: fit an endpoint-fit cubic on the unit x-grid (the tdt
+    # weights contain dθ/dx) and integrate it exactly. With the 1/√v_par endpoint
+    # singularities the quadrature scheme is a leading-order effect, so this matches
+    # the exact-cubic integration of the Fortran spline package (not a trapezoid sum).
+    itp_r = scr.itp_r
+    copyto!(itp_r.y, g_wb)
+    _refit!(itp_r)
+    fsi_wb = scr.cum_wb_arr  # per-surface scratch; fully overwritten
+    FastInterpolations.cumulative_integrate!(fsi_wb, itp_r)
+    total_wb = fsi_wb[ntheta]
+    copyto!(itp_r.y, g_wd)
+    _refit!(itp_r)
+    total_wd = FastInterpolations.integrate(itp_r)
 
     if total_wb ≈ 0.0
         # Degenerate case — return zeros
         return 0.0, 0.0, 0.0, nothing
     end
 
-    # Bounce-averaged frequencies (Fortran lines 740-741)
+    # Bounce-averaged frequencies
     wbbar = ro * twopi / ((2 - sigma) * total_wb)
     wdbar = ro^2 * bo * wdfac * wbbar * 2 * (2 - sigma) * total_wd
 
-    # Phase factor (Fortran line 750). Ratio cum_wb_arr[i]/total_wb is dimensionless —
-    # (ntheta-1) cancels, no scaling. For do_matrices we keep `pl` as a Vector because
-    # it is referenced below; otherwise we fuse pl + bjspl → bj_integral in a single
-    # pass, avoiding two Vector{ComplexF64}(ntheta) allocations.
-    # Trapezoidal quadrature: boundary samples weighted by 0.5. jvtheta is zero at
-    # i=1 and i=ntheta (loop above runs 2:ntheta-1) so the boundary terms contribute
-    # nothing in practice, but writing the weights explicitly keeps the integration
-    # self-correct if the boundary handling ever changes.
+    # Phase factor pl_i = exp(-2πi·lnq·fsi_wb(θ_i)/((2-σ)·total_wb)), using the
+    # cumulative spline integral of the bounce action.
     pl_denom = (2 - sigma) * total_wb
     one_minus_sigma = 1 - sigma
-    bj_integral = ComplexF64(0.0)
-    if do_matrices
-        pl = scr.pl
-        @inbounds for i in 1:ntheta
-            pl[i] = cis(-twopi * lnq * cum_wb_arr[i] * nrm / pl_denom)
-        end
-        @inbounds for i in 1:ntheta
-            w = (i == 1 || i == ntheta) ? 0.5 : 1.0
-            bj_integral += w * conj(jvtheta[i]) * (pl[i] + one_minus_sigma / (pl[i] + SINGULAR_EPS))
-        end
-    else
-        pl = nothing
-        @inbounds for i in 1:ntheta
-            pli = cis(-twopi * lnq * cum_wb_arr[i] * nrm / pl_denom)
-            w = (i == 1 || i == ntheta) ? 0.5 : 1.0
-            bj_integral += w * conj(jvtheta[i]) * (pli + one_minus_sigma / (pli + SINGULAR_EPS))
-        end
+    pl = scr.pl  # per-surface scratch; fully overwritten
+    @inbounds for i in 1:ntheta
+        pl[i] = cis(-twopi * lnq * fsi_wb[i] / pl_denom)
+    end
+
+    # Action bounce integral (exact-cubic quadrature, as for the totals above).
+    itp_c = scr.itp_c
+    bj_samples = itp_c.y
+    @inbounds for i in 1:ntheta
+        bj_samples[i] = conj(jvtheta[i]) * (pl[i] + one_minus_sigma / pl[i])
     end
-    bj_integral *= nrm
+    _refit!(itp_c)
+    bj_integral = FastInterpolations.integrate(itp_c)
 
-    # |δJ|² (Fortran line 756) — division by 2 corrects quadratic form
+    # |δJ|² — division by 2 corrects the quadratic form
     dJdJ_val = wbbar * abs(bj_integral)^2 / 2.0 / ro^2
 
-    # Matrix path: bounce-average W vectors and form outer products (Fortran lines 759-793)
+    # Matrix path: bounce-average W vectors and form outer products
     wmats_lmda = nothing
     if do_matrices
-        # Bounce-average W_μ and W_E vectors (Fortran lines 762-767).
-        # Trapezoidal quadrature: boundary samples weighted by 0.5 (wmu_mt and wen_mt
-        # are zero at i=1 and i=ntheta from the 2:ntheta-1 population loop above).
+        # Per mode, exact-cubic quadrature of conj(W_m(θ))·(pl + (1-σ)/pl),
+        # matching the action bounce integral above.
         wmu_ba = scr.wmu_ba
         wen_ba = scr.wen_ba
-        fill!(wmu_ba, ComplexF64(0.0))
-        fill!(wen_ba, ComplexF64(0.0))
-        @inbounds for i in 1:ntheta
-            w = (i == 1 || i == ntheta) ? 0.5 : 1.0
-            factor = w * (pl[i] + one_minus_sigma / (pl[i] + SINGULAR_EPS))
-            for mi in 1:mpert
-                wmu_ba[mi] += conj(wmu_mt[mi, i]) * factor
-                wen_ba[mi] += conj(wen_mt[mi, i]) * factor
+        wsamp = itp_c.y
+        @inbounds for mi in 1:mpert
+            for i in 1:ntheta
+                wsamp[i] = conj(wmu_mt[mi, i]) * (pl[i] + one_minus_sigma / pl[i])
+            end
+            _refit!(itp_c)
+            wmu_ba[mi] = FastInterpolations.integrate(itp_c)
+            for i in 1:ntheta
+                wsamp[i] = conj(wen_mt[mi, i]) * (pl[i] + one_minus_sigma / pl[i])
             end
+            _refit!(itp_c)
+            wen_ba[mi] = FastInterpolations.integrate(itp_c)
         end
-        wmu_ba .*= nrm
-        wen_ba .*= nrm
 
         # Reshape as 1×mpert for matrix multiply (Fortran lines 771-772)
         wmmt = reshape(wmu_ba, 1, mpert)
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 5e235812..7c3eb016 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -254,9 +254,14 @@ the equilibrium geometry parameters needed for NTV calculations.
 function KineticForcesInternal(equil; verbose::Bool=false)
     mthsurf = length(equil.rzphi_ys) - 1
     nth = mthsurf + 1
+    # Toroidal field at the magnetic axis, which normalizes λ = μ·bo/E: F_spline
+    # stores 2πF (F = R·B_t), so |F_spline(0)|/(2π·ro) = F(0)/ro, the same
+    # normalization as Fortran PENTRC. params.b0 is the edge-extrapolated vacuum
+    # field and differs by ~2%.
+    bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro)
     KineticForcesInternal(;
         ro      = equil.ro,
-        bo      = equil.params.b0,
+        bo      = bo_axis,
         chi1    = 2π * equil.psio,
         mthsurf,
         tpsi_xs       = collect(range(0.0, 1.0, length=nth)),
@@ -309,16 +314,15 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in
     psi_grid = pe_state.psi_grid
     npsi = length(psi_grid)
     mpert = ffs_intr.mpert
-    chi1 = kf_intr.chi1
 
     # Build xs_m: 3 CubicSeriesInterpolants from Clebsch displacement matrices
     # xs_m[1] = ξ^ψ (unregularized), xs_m[2] = ∂ξ^ψ/∂ψ (regularized), xs_m[3] = ξ^α
-    # Note: clebsch_alpha is stored as ξ^α/χ₁, multiply by chi1 to get ξ^α
     itp_opts = (; extrap=ExtendExtrap())
+    # clebsch_alpha is already the physical ξ^α = xms/χ₁ (as Fortran gpout_xclebsch writes it); use directly.
+    clebsch_alpha_mat = xi_modes.clebsch_alpha
     xs_m_1 = cubic_interp(psi_grid, Series(xi_modes.clebsch_psi); itp_opts...)
     xs_m_2 = cubic_interp(psi_grid, Series(xi_modes.clebsch_psi1); itp_opts...)
-    clebsch_alpha_raw = xi_modes.clebsch_alpha .* chi1
-    xs_m_3 = cubic_interp(psi_grid, Series(clebsch_alpha_raw); itp_opts...)
+    xs_m_3 = cubic_interp(psi_grid, Series(clebsch_alpha_mat); itp_opts...)
     kf_intr.xs_m = [xs_m_1, xs_m_2, xs_m_3]
 
     # Build geometric matrices (S,T,X,Y,Z) for JBB deweighting
@@ -355,7 +359,7 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in
         # Get Clebsch displacement vectors at this ψ
         xsp  = view(xi_modes.clebsch_psi,  ipsi, :)       # ξ^ψ [mpert]
         xmp1 = view(xi_modes.clebsch_psi1, ipsi, :)       # ∂ξ^ψ/∂ψ [mpert]
-        xms  = view(clebsch_alpha_raw, ipsi, :)            # ξ^α [mpert]
+        xms  = view(clebsch_alpha_mat, ipsi, :)            # ξ^α [mpert]
 
         # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert
         geom_mats.smats(smat_flat, psi; hint=hint_s)
diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl
index 240358b1..c5d85dcd 100644
--- a/src/KineticForces/Torque.jl
+++ b/src/KineticForces/Torque.jl
@@ -94,6 +94,10 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
 
     # Create periodic interpolant for poloidal quantities
     tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC())
+    # v_par and the bounce points use a separate endpoint-fit (non-periodic) cubic
+    # of B, like Fortran's vspl. The endpoint fit of 1−(λ/bo)B equals 1−(λ/bo)
+    # times the fit of B, so one B_extrap per surface serves every λ.
+    B_extrap = cubic_interp(xs, B_vals; bc=CubicFit())
 
     bmax = maximum(B_vals)
     ibmax = argmax(B_vals)
@@ -138,7 +142,8 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
         end
         if Bθ > bmax
             bmax = Bθ
-            theta_bmax = θn
+            # theta_bmax stays at the nodal knot xs[ibmax]: the transit starts at
+            # the knot (as in Fortran), not at the refined extremum.
         end
     end
 
@@ -225,7 +230,8 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
                                    tspl, dbob_m_f, divx_m_f, divxfac, wdfac,
                                    method, op_wmats;
                                    chi1=intr.chi1, ro=intr.ro, mfac=intr.mfac,
-                                   mpert=intr.mpert, ibmax=ibmax, theta_bmax=theta_bmax,
+                                   mpert=intr.mpert, theta_bmax=theta_bmax,
+                                   B_extrap=B_extrap,
                                    smat=smat_f, tmat=tmat_f, xmat=xmat_f,
                                    ymat=ymat_f, zmat=zmat_f,
                                    energy_atol=atol_xlmda, energy_rtol=rtol_xlmda,
@@ -400,10 +406,10 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
                        bmax, bmin, n_s::Float64, T_s::Float64, mass, chrg, tspl,
                        dbob_m_f, divx_m_f, divxfac, wdfac, method, op_wmats;
                        chi1::Float64, ro::Float64, mfac::Vector{Int}, mpert::Int,
-                       ibmax::Int, theta_bmax::Float64,
+                       theta_bmax::Float64, B_extrap,
                        smat=nothing, tmat=nothing, xmat=nothing,
                        ymat=nothing, zmat=nothing,
-                       nlmda::Int=64, ntheta::Int=128,
+                       nlmda::Int=128, ntheta::Int=128,
                        nutype::String="harmonic", f0type::String="maxwellian",
                        nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false,
                        energy_atol::Float64=1e-7, energy_rtol::Float64=1e-5,
@@ -413,8 +419,8 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
 
     # 1. Compute bounce-averaged quantities
     bounce = compute_bounce_data(
-        psi, n, l, q, bo, bmax, bmin, ibmax, theta_bmax,
-        tspl, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac,
+        psi, n, l, q, bo, bmax, bmin, theta_bmax,
+        tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac,
         mass, chrg, T_s, method;
         nlmda, ntheta, smat, tmat, xmat, ymat, zmat)
 
@@ -458,8 +464,9 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
         end
     end
 
-    # Build CubicSeriesInterpolant on normalized data (single build)
-    fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=ZeroCurvBC())
+    # Build CubicSeriesInterpolant on normalized data (single build); CubicFit
+    # endpoint BC matches the Fortran "extrap" fit.
+    fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=CubicFit())
 
     # 3. Set rex/imx multipliers (Fortran lines 839-847)
     method_suffix = length(method) >= 4 ? method[2:4] : ""
@@ -617,6 +624,8 @@ function _setup_surface_state(
     end
 
     tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC())
+    # Endpoint-fit (non-periodic) cubic of B for v_par and bounce points (Fortran vspl equivalent).
+    B_extrap = cubic_interp(xs, B_vals; bc=CubicFit())
 
     bmax = maximum(B_vals)
     ibmax = argmax(B_vals)
@@ -650,7 +659,8 @@ function _setup_surface_state(
             end
             if Bθ > bmax
                 bmax = Bθ
-                theta_bmax = θn
+                # theta_bmax stays at the nodal knot xs[ibmax]: the transit starts at
+                # the knot (as in Fortran), not at the refined extremum.
             end
         end
     end
@@ -687,7 +697,7 @@ function _setup_surface_state(
 
     return (;
         chrg, mass,
-        tspl, bmax, bmin, ibmax, theta_bmax,
+        tspl, B_extrap, bmax, bmin, theta_bmax,
         q, n_s, T_s, welec,
         wdian, wdiat, wtran, wgyro, nuk,
         epsr,
@@ -751,8 +761,8 @@ function kinetic_energy_matrices_for_euler_lagrange!(
     divx_m_f = zeros(ComplexF64, mpert)
 
     bounce = compute_bounce_data(
-        psi, n, l, state.q, bo, state.bmax, state.bmin, state.ibmax, state.theta_bmax,
-        state.tspl, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac,
+        psi, n, l, state.q, bo, state.bmax, state.bmin, state.theta_bmax,
+        state.tspl, state.B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac,
         state.mass, state.chrg, state.T_s, "fwmm";
         nlmda, ntheta, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f)
 
@@ -789,7 +799,8 @@ function kinetic_energy_matrices_for_euler_lagrange!(
         end
     end
 
-    fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=ZeroCurvBC())
+    # CubicFit endpoint BC matches Fortran cspline_fit(fbnce,'extrap').
+    fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=CubicFit())
 
     # Dual-output pitch integration: one energy sweep per (λ, E) produces both
     # halves. Packed return is [wmm | tmm], each length nqty.