diff --git a/docs/make.jl b/docs/make.jl
index c5f66fde..5e2b9ad3 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -24,12 +24,11 @@ makedocs(;
         "Home" => "index.md",
         "Setup" => "set_up.md",
         "API Reference" => [
-            "Splines" => "splines.md",
             "Vacuum" => "vacuum.md",
             "Equilibrium" => "equilibrium.md",
             "Utilities" => "utilities.md",
             "Perturbed Equilibrium" => "perturbed_equilibrium.md"
-        ],
+        ]
     ],
     checkdocs=:exports
 )
diff --git a/docs/src/splines.md b/docs/src/splines.md
deleted file mode 100644
index aa83867b..00000000
--- a/docs/src/splines.md
+++ /dev/null
@@ -1,128 +0,0 @@
-# Splines Module
-
-The Splines module provides helpers and utilities for working with FastInterpolations.jl cubic splines, used for smooth representation of MHD equilibria.
-
-## Overview
-
-JPEC uses [FastInterpolations.jl](https://github.com/ThummeTo/FastInterpolations.jl) for high-performance cubic spline interpolation. The SplinesMod provides:
-- Helper functions for exact spline integration (`cumulative_integral`, `total_integral`)
-- Re-exports of FastInterpolations boundary condition types
-- Re-exports of FastInterpolations search strategies for optimized evaluation
-
-## API Reference
-
-```@autodocs
-Modules = [JPEC.SplinesMod]
-```
-
-## Boundary Condition Types
-
-The module re-exports FastInterpolations boundary condition types:
-- `NaturalBC()`: Natural boundary condition (f''(x) = 0 at endpoints)
-- `PeriodicBC()`: Periodic boundary condition (requires data[end] == data[1])
-- `CubicFit()`: Automatic endpoint derivative estimation using 4-point fit
-- `BCPair(bc_left, bc_right)`: Different BCs at left and right boundaries
-- `Deriv1(value)`: First derivative boundary condition
-- `Deriv2(value)`: Second derivative boundary condition
-
-## Search Strategies
-
-For optimized evaluation (especially with monotonic access patterns):
-- `LinearBinary()`: Linear search + binary fallback (~4 ns/point for monotonic access)
-- `Binary()`: Binary search
-- `HintedBinary()`: Binary search with hint
-
-## Integration Functions
-
-### cumulative_integral
-```@docs
-JPEC.SplinesMod.cumulative_integral
-```
-
-### total_integral
-```@docs
-JPEC.SplinesMod.total_integral
-```
-
-### integrate_spline
-```@docs
-JPEC.SplinesMod.integrate_spline
-```
-
-## Example Usage
-
-### 1D Cubic Spline
-```julia
-using JPEC
-using FastInterpolations: cubic_interp, deriv1
-
-# Create data points
-xs = collect(range(0.0, 2π, length=21))
-fs = sin.(xs)
-
-# Create spline with natural boundary conditions
-spline = cubic_interp(xs, fs; bc=JPEC.SplinesMod.NaturalBC())
-
-# Evaluate at a single point
-x = 1.0
-f = spline(x)
-
-# Evaluate derivative at a point
-d1_spline = deriv1(spline)
-f_deriv = d1_spline(x)
-
-# Evaluate at vector of new points
-xs_fine = collect(range(0.0, 2π, length=100))
-fs_fine = spline.(xs_fine)
-
-# Compute cumulative integral
-fs_integral = JPEC.SplinesMod.cumulative_integral(xs, fs; bc=JPEC.SplinesMod.NaturalBC())
-```
-
-### 2D Cubic Spline (CubicInterpolantND)
-```julia
-using FastInterpolations: cubic_interp
-
-# Create 2D grid
-xs = collect(range(0.0, 0.99, length=20))
-ys = collect(range(0.0, 2π, length=20))
-
-# Create 2D function data
-fs = zeros(20, 20)
-for i in 1:20, j in 1:20
-    fs[i, j] = sqrt(xs[i]) * cos(ys[j])
-end
-
-# Set up 2D cubic interpolant
-# Use PeriodicBC on second dimension (θ direction)
-spline_2d = cubic_interp((xs, ys), fs; 
-    bc=(JPEC.SplinesMod.CubicFit(), JPEC.SplinesMod.PeriodicBC()),
-    extrap=(:extension, :wrap))
-
-# Evaluate spline at a point
-x_eval, y_eval = 0.5, π/4
-f = spline_2d((x_eval, y_eval))
-
-# Evaluate with derivative (specify which dimension)
-fx = spline_2d((x_eval, y_eval); deriv=(1, 0))  # ∂f/∂x
-fy = spline_2d((x_eval, y_eval); deriv=(0, 1))  # ∂f/∂y
-```
-
-### Multi-series 1D Spline (CubicSeriesInterpolant)
-```julia
-using FastInterpolations: cubic_interp
-
-# Create multiple data series
-xs = collect(range(0.0, 1.0, length=50))
-fs = hcat(sin.(2π*xs), cos.(2π*xs), xs.^2)  # 3 series
-
-# Create series interpolant (more efficient than separate splines)
-spline_series = cubic_interp(xs, fs; bc=JPEC.SplinesMod.CubicFit())
-
-# Evaluate all series at once
-x = 0.5
-f_all = spline_series(x)  # Returns vector of length 3
-
-# Compute total integral for all series
-integrals = JPEC.SplinesMod.total_integral(xs, fs; bc=JPEC.SplinesMod.CubicFit())
-```
diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl
index c4370004..ef46c8bc 100644
--- a/src/Equilibrium/AnalyticEquilibrium.jl
+++ b/src/Equilibrium/AnalyticEquilibrium.jl
@@ -208,10 +208,10 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
     # Create separate interpolants for R and Z coordinates
     rz_in_xs = r_nodes
     rz_in_ys = collect(rzphi_y_nodes)
-    rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     # LAR equilibrium has midplane symmetry, so magnetic axis is at Z = 0
     return InverseRunInput(equil_input, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, lar_r0, 0.0, psio)
@@ -288,7 +288,7 @@ function sol_run(equil_inputs::EquilibriumConfig, sol_inputs::SolovevConfig)
     end
     psi_in_xs = r
     psi_in_ys = z
-    psi_in = cubic_interp((psi_in_xs, psi_in_ys), psifs;
+    psi_in = cubic_interp((psi_in_xs, psi_in_ys), psifs; search=LinearBinary(),
         bc=(CubicFit(), CubicFit()), extrap=(:extension, :extension))
 
     # Print out equilibrium info
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 8393c89f..08a0b947 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -191,7 +191,7 @@ function direct_position!(raw_profile::DirectRunInput)
     # Access nodal values from psi_in interpolant: partials[1,:,:] = function values
     new_psi_fs = raw_profile.psi_in.nodal_derivs.partials[1, :, :] .* raw_profile.psio / bfield.psi
     # Because DirectRunInput is a mutable struct, we can update the spline here
-    raw_profile.psi_in = cubic_interp((x_coords, y_coords), new_psi_fs;
+    raw_profile.psi_in = cubic_interp((x_coords, y_coords), new_psi_fs; search=LinearBinary(),
         bc=(CubicFit(), CubicFit()), extrap=(:extension, :extension))
 
     # Helper function for robust Newton-Raphson search with restarts
@@ -467,7 +467,7 @@ function equilibrium_solver(raw_profile::DirectRunInput)
         ff_fs_nodes[end, :] .= ff_fs_nodes[1, :]
 
         # Create series interpolant for all columns
-        ff_interp = cubic_interp(ff_x_nodes, ff_fs_nodes; bc=Spl.PeriodicBC())
+        ff_interp = cubic_interp(ff_x_nodes, ff_fs_nodes; bc=PeriodicBC())
         ff_deriv = deriv1(ff_interp)
 
         # Interpolate `ff` onto the uniform `theta` grid for `rzphi`
@@ -524,14 +524,14 @@ function equilibrium_solver(raw_profile::DirectRunInput)
     # theta_nodes includes both 0 and 1 (closed periodic grid).
     rzphi_xs = psi_nodes
     rzphi_ys = collect(theta_nodes)
-    rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 3];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 4];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 3]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 4]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     # Calculate physics quantities (B-field, metric components, etc.) in 2D spline `eqfun`
     # for use in stability and transport codes
@@ -596,12 +596,12 @@ function equilibrium_solver(raw_profile::DirectRunInput)
         end
     end
     # Create 2D interpolants for physics quantities (eqfun)
-    eqfun_B = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    eqfun_metric1 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    eqfun_metric2 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 3];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_B = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_metric1 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_metric2 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 3]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     return PlasmaEquilibrium(raw_profile.config, EquilibriumParameters(), profiles,
         rzphi_xs, rzphi_ys,
diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl
index 40b5789b..370f3858 100644
--- a/src/Equilibrium/Equilibrium.jl
+++ b/src/Equilibrium/Equilibrium.jl
@@ -1,12 +1,11 @@
 module Equilibrium
 
 # --- Module-level Dependencies ---
-import ..Spl
 
 using Printf, OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra, HDF5
 using TOML
 import FastInterpolations
-using FastInterpolations: cubic_interp, deriv1, deriv2, deriv3, LinearBinary, CubicFit
+using FastInterpolations: cubic_interp, deriv1, deriv2, deriv3, LinearBinary, CubicFit, PeriodicBC
 import StaticArrays: @MMatrix, SVector
 
 # --- Internal Module Structure ---
@@ -107,11 +106,12 @@ function equilibrium_separatrix_find!(pe::PlasmaEquilibrium)
 
     for iside in 1:2
         it = 0
+        hint2d = (Ref(1), Ref(1))  # Shared 2D hint for Newton iteration
         while true
             it += 1
             # Evaluate offset and its derivative
-            offset = pe.rzphi_offset((psi_edge, theta))
-            offset_y = pe.rzphi_offset((psi_edge, theta); deriv=(0,1))
+            offset = pe.rzphi_offset((psi_edge, theta); hint=hint2d)
+            offset_y = pe.rzphi_offset((psi_edge, theta); deriv=(0, 1), hint=hint2d)
 
             eta = theta + offset - eta0
             eta_theta = 1 + offset_y
@@ -143,15 +143,16 @@ function equilibrium_separatrix_find!(pe::PlasmaEquilibrium)
         z = 0.0
         max_iter = 1000
         iter = 0
+        hint2d = (Ref(1), Ref(1))  # Shared 2D hint for Newton iteration
         while iter < max_iter
             iter += 1
             # Evaluate rcoord and offset with derivatives
-            r2 = pe.rzphi_rsquared((psi_edge, theta))
-            r2y = pe.rzphi_rsquared((psi_edge, theta); deriv=(0,1))
-            r2yy = pe.rzphi_rsquared((psi_edge, theta); deriv=(0,2))
-            eta = pe.rzphi_offset((psi_edge, theta))
-            eta1 = pe.rzphi_offset((psi_edge, theta); deriv=(0,1))
-            eta2 = pe.rzphi_offset((psi_edge, theta); deriv=(0,2))
+            r2 = pe.rzphi_rsquared((psi_edge, theta); hint=hint2d)
+            r2y = pe.rzphi_rsquared((psi_edge, theta); deriv=(0, 1), hint=hint2d)
+            r2yy = pe.rzphi_rsquared((psi_edge, theta); deriv=(0, 2), hint=hint2d)
+            eta = pe.rzphi_offset((psi_edge, theta); hint=hint2d)
+            eta1 = pe.rzphi_offset((psi_edge, theta); deriv=(0, 1), hint=hint2d)
+            eta2 = pe.rzphi_offset((psi_edge, theta); deriv=(0, 2), hint=hint2d)
 
             rfac = sqrt(r2)
             rfac1 = r2y / (2 * rfac)
@@ -441,17 +442,18 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
         end
     end
     # Create flux interpolants for Grad-Shafranov diagnostics
-    flux1 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    flux1 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
     # Compute flux derivatives at all grid points for diagnostics
+    hint2d = (Ref(1), Ref(1))  # Shared 2D hint for hot loop optimization
     for ipsi in 0:mpsi
         for itheta in 0:mtheta
-            flux_fsx[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1,0))
-            flux_fsx[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1,0))
-            flux_fsy[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0,1))
-            flux_fsy[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0,1))
+            flux_fsx[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1, 0), hint=hint2d)
+            flux_fsx[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1, 0), hint=hint2d)
+            flux_fsy[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0, 1), hint=hint2d)
+            flux_fsy[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0, 1), hint=hint2d)
         end
     end
 
@@ -492,8 +494,9 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
         fs_matrix[:, 1] = flux_fsx[ipsi, :, 1]
         fs_matrix[:, 2] = source[ipsi, :]
 
-        # Compute total integral using exact spline integration (only final value needed)
-        term[ipsi, :] .= Spl.total_integral(equil.rzphi_ys, fs_matrix; bc=Spl.PeriodicBC())
+        # Compute total integral using FastInterpolations native integration
+        itp = cubic_interp(equil.rzphi_ys, fs_matrix; bc=PeriodicBC())
+        term[ipsi, :] .= FastInterpolations.integrate(itp)
     end
 
     totali = sum(term; dims=2)
diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl
index e589fd8d..2ab55b7e 100644
--- a/src/Equilibrium/InverseEquilibrium.jl
+++ b/src/Equilibrium/InverseEquilibrium.jl
@@ -118,10 +118,10 @@ function equilibrium_solver(input::InverseRunInput)
     deta[:, end] .= deta[:, 1]
 
     # Create 2D interpolants for r² and dη
-    rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2;
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rz_deta = cubic_interp((rz_in_xs, rz_in_ys), deta;
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rz_deta = cubic_interp((rz_in_xs, rz_in_ys), deta; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     # c-----------------------------------------------------------------------
     # c     prepare new spline type for surface quantities.
@@ -137,7 +137,7 @@ function equilibrium_solver(input::InverseRunInput)
     if grid_type == "ldp"
         sq_xs = psilow .+ (psihigh - psilow) .* (sin.(range(0.0, 1.0; length=mpsi+1) .* (π/2))) .^ 2
         sq_fs = zeros(Float64, mpsi+1, 4)
-        sq = cubic_interp(sq_xs, sq_fs; bc=Spl.CubicFit(), extrap=:extension)
+        sq = cubic_interp(sq_xs, sq_fs; bc=CubicFit(), extrap=:extension)
     else
         error("Only 'ldp' grid_type is implemented for now.")
     end
@@ -206,8 +206,8 @@ function equilibrium_solver(input::InverseRunInput)
         # (Numerical operations may have broken exact periodicity)
         spl_fs[end, :] .= spl_fs[1, :]
 
-        spl = cubic_interp(spl_xs, spl_fs; bc=Spl.PeriodicBC())
-        spl_fsi = Spl.cumulative_integral(spl_xs, spl_fs; bc=Spl.PeriodicBC())
+        spl = cubic_interp(spl_xs, spl_fs; bc=PeriodicBC())
+        spl_fsi = FastInterpolations.cumulative_integrate(spl)
 
         spl_xs = spl_fsi[:, 5] ./ spl_fsi[mtheta+1, 5]
         spl_fs[:, 2] .+= rzphi_ys .- spl_xs
@@ -226,7 +226,7 @@ function equilibrium_solver(input::InverseRunInput)
         sq_fs[ipsi+1, 4] = spl_fsi[mtheta+1, 4] * sq_fs[ipsi+1, 1] / (2 * twopi * psio) # q-profile
     end
 
-    sq = cubic_interp(sq_xs, sq_fs; bc=Spl.CubicFit(), extrap=:extension)
+    sq = cubic_interp(sq_xs, sq_fs; bc=CubicFit(), extrap=:extension)
 
     # Evaluate sq and its derivative at all grid points
     f_sq = zeros(Float64, mpsi+1, 4)
@@ -251,18 +251,18 @@ function equilibrium_solver(input::InverseRunInput)
             sq_fs[ipsi+1, 4] *= ffac
             rzphi_fs[ipsi+1, :, 3] *= ffac
         end
-        sq = cubic_interp(sq_xs, sq_fs; bc=Spl.CubicFit(), extrap=:extension)
+        sq = cubic_interp(sq_xs, sq_fs; bc=CubicFit(), extrap=:extension)
     end
     qa = f_sq[mpsi+1, 4] + f1_sq[mpsi+1, 4] * (1 - sq_xs[mpsi+1])
     # Create 2D interpolants for geometric quantities (rzphi)
-    rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 3];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 4];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 3]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 4]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     for ipsi in 0:mpsi
         f_sq = sq(sq_xs[ipsi+1])
@@ -311,12 +311,12 @@ function equilibrium_solver(input::InverseRunInput)
         end
     end
     # Create 2D interpolants for physics quantities (eqfun)
-    eqfun_B = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    eqfun_metric1 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    eqfun_metric2 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 3];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_B = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_metric1 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_metric2 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 3]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     # Create ProfileSplines from sq interpolant
     # sq_fs columns: [F*2π, P, dV/dψ, q]
diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl
index 1ce8a974..db45255f 100644
--- a/src/Equilibrium/ReadEquilibrium.jl
+++ b/src/Equilibrium/ReadEquilibrium.jl
@@ -115,7 +115,7 @@ function read_efit(config::EquilibriumConfig)
 
     psi_in_xs = collect(r_grid)
     psi_in_ys = collect(z_grid)
-    psi_in = cubic_interp((psi_in_xs, psi_in_ys), psi_proc;
+    psi_in = cubic_interp((psi_in_xs, psi_in_ys), psi_proc; search=LinearBinary(),
         bc=(CubicFit(), CubicFit()), extrap=(:extension, :extension))
 
     # --- Bundle everything for the solver ---
@@ -171,8 +171,9 @@ function read_chease_binary(config::EquilibriumConfig)
         fs[:, 2] .= zcppr
         fs[:, 3] .= zq
 
-        # Compute cumulative integral of pressure column for normalization
-        fsi_pressure = Spl.cumulative_integral(xs, fs[:, 2]; bc=CubicFit())
+        # Compute cumulative integral of pressure column for normalization using FastInterpolations
+        itp_pressure = cubic_interp(xs, fs[:, 2]; bc=CubicFit())
+        fsi_pressure = FastInterpolations.cumulative_integrate(itp_pressure)
         # Make a writable copy and normalize pressure integral
         fs_copy = copy(fs)
         fs_copy[:, 2] .= (fsi_pressure .- fsi_pressure[ma]) .* psio
@@ -209,10 +210,10 @@ function read_chease_binary(config::EquilibriumConfig)
         # Create separate interpolants for R and Z coordinates
         rz_in_xs = xs
         rz_in_ys = range(0, 2π; length=mtau) |> collect
-        rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1];
-            bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-        rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2];
-            bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+        rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1]; search=LinearBinary(),
+            bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+        rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2]; search=LinearBinary(),
+            bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
         println("--> Finished reading CHEASE equilibrium (Binary).")
         return InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio)
@@ -351,8 +352,9 @@ function read_chease_ascii(config::EquilibriumConfig)
     fs[:, 2] .= zcppr # normalized Pressure
     fs[:, 3] .= zq # q profile
     # Fit spline with extrapolation boundary condition (bctype = 3)
-    # Compute cumulative integral of pressure column for normalization
-    fsi_pressure = Spl.cumulative_integral(xs, fs[:, 2]; bc=CubicFit())
+    # Compute cumulative integral of pressure column for normalization using FastInterpolations
+    itp_pressure = cubic_interp(xs, fs[:, 2]; bc=CubicFit())
+    fsi_pressure = FastInterpolations.cumulative_integrate(itp_pressure)
     # Make a writable copy and normalize pressure integral
     fs_copy = copy(fs)
     fs_copy[:, 2] .= (fsi_pressure .- fsi_pressure[ma]) .* psio
@@ -372,11 +374,11 @@ function read_chease_ascii(config::EquilibriumConfig)
 
     # Create separate interpolants for R and Z coordinates
     rz_in_xs = xs
-    rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_data;
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_data;
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_data; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_data; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
     println("--> Finished reading CHEASE equilibrium.")
     println("    Magnetic axis at (ro=$ro, zo=$zo), psio=$psio")
     return InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio)
-end
\ No newline at end of file
+end
diff --git a/src/ForceFreeStates/Bal.jl b/src/ForceFreeStates/Bal.jl
index 1f858319..52f81b9a 100644
--- a/src/ForceFreeStates/Bal.jl
+++ b/src/ForceFreeStates/Bal.jl
@@ -310,7 +310,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
 
     # initialize bg values
     spl0_bg_fs = -pressure_gradient .* q_derivative .* two_pi_f ./ (bsq .* chi_prime^2)
-    spl0_bg_fsi = Spl.cumulative_integral(theta_grid, spl0_bg_fs; bc=PeriodicBC())
+    itp_bg = cubic_interp(theta_grid, spl0_bg_fs; bc=PeriodicBC())
+    spl0_bg_fsi = FastInterpolations.cumulative_integrate(itp_bg)
     bg_fs = zeros(mtheta + 1, 5)
     bg_fs[:, 5] = spl0_bg_fsi .- spl0_bg_fsi[end]
 
@@ -328,7 +329,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
         spl1_fs[itheta, 4] = -spl1_fs[itheta, 1]
     end
 
-    spl1_totals = Spl.total_integral(theta_grid, spl1_fs; bc=PeriodicBC())
+    itp_spl1 = cubic_interp(theta_grid, spl1_fs; bc=PeriodicBC())
+    spl1_totals = FastInterpolations.integrate(itp_spl1)
 
     d0bar = [spl1_totals[1] spl1_totals[2]; spl1_totals[3] spl1_totals[4]]
 
@@ -361,7 +363,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
         spl2_fs[itheta, 4] = spl1_fs[itheta, 3] * v0[1, 2] + (spl1_fs[itheta, 4] + alpha) * v0[2, 2]
     end
 
-    spl2_fsi = Spl.cumulative_integral(theta_grid, spl2_fs; bc=PeriodicBC())
+    itp_spl2 = cubic_interp(theta_grid, spl2_fs; bc=PeriodicBC())
+    spl2_fsi = FastInterpolations.cumulative_integrate(itp_spl2)
     # CRITICAL: Use integrated values as the new fs (matching Fortran's spl2%fs=spl2%fsi)
     spl2_fs_new = copy(spl2_fsi)
 
@@ -400,7 +403,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
                              d1_21 * v0[1, 2] + d1_22 * v0[2, 2]
     end
 
-    spl3_totals = Spl.total_integral(theta_grid, spl3_fs; bc=PeriodicBC())
+    itp_spl3 = cubic_interp(theta_grid, spl3_fs; bc=PeriodicBC())
+    spl3_totals = FastInterpolations.integrate(itp_spl3)
 
     # Compute first-order constants for both eigenfunctions
     # First eigenfunction (with -alpha correction)
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
index 5d5adfb6..06a48fe4 100644
--- a/src/ForceFreeStates/ForceFreeStates.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -9,9 +9,9 @@ using OrdinaryDiffEq
 using HDF5
 using JLD2
 using FastInterpolations
+using FastInterpolations: cubic_interp, deriv1, PeriodicBC, LinearBinary
 
 import ..Equilibrium
-import ..Spl
 import ..Utilities
 import ..Vacuum
 using Printf
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index f4090503..4124454b 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -466,6 +466,9 @@ and a small set of temporary matrices and factors used to compute singular-layer
     spline_hint::Base.RefValue{Int} = Ref(1)
     # Separate hint for wvmat splines (different grid size than equilibrium profiles)
     wv_hint::Base.RefValue{Int} = Ref(1)
+    # Shared 2D hint for CubicInterpolantND (rzphi splines) during ODE integration
+    # Tuple of (psi_hint, theta_hint) for O(1) interval lookups in 2D bicubic splines
+    rzphi_hint::Tuple{Base.RefValue{Int},Base.RefValue{Int}} = (Ref(1), Ref(1))
 end
 
 # Initialize function for OdeState with relevant parameters for array initialization
diff --git a/src/ForceFreeStates/Mercier.jl b/src/ForceFreeStates/Mercier.jl
index 79f6621e..b23be520 100644
--- a/src/ForceFreeStates/Mercier.jl
+++ b/src/ForceFreeStates/Mercier.jl
@@ -60,8 +60,9 @@ function mercier_scan!(locstab_fs::Matrix{Float64}, plasma_eq::Equilibrium.Plasm
             @views ff_fs[itheta, :] .*= jac / v1
         end
 
-        # Integrate quantities with respect to theta using exact spline integral
-        avg = Spl.total_integral(plasma_eq.rzphi_ys, ff_fs; bc=Spl.PeriodicBC())
+        # Integrate quantities with respect to theta using FastInterpolations
+        itp = cubic_interp(plasma_eq.rzphi_ys, ff_fs; bc=PeriodicBC())
+        avg = FastInterpolations.integrate(itp)
 
         # Evaluate Mercier criterion and related quantities
         term = twopif * p1 * v1 / (q1 * chi1^3) * avg[2]
diff --git a/src/JPEC.jl b/src/JPEC.jl
index 26a94f0f..309f8779 100755
--- a/src/JPEC.jl
+++ b/src/JPEC.jl
@@ -1,10 +1,6 @@
 # JPEC.jl
 module JPEC
 
-include("Splines/Splines.jl")
-import .SplinesMod as Spl
-export SplinesMod, Spl
-
 include("Utilities/Utilities.jl")
 import .Utilities as Utilities
 export Utilities
@@ -95,7 +91,7 @@ function main(args::Vector{String}=String[])
     else
         intr.debug_settings = DebugSettings()
     end
-    
+
     # Set up variables
     # TODO: parallel threads logic
     ctrl.delta_mhigh *= 2 # for consistency with Fortran DCON TODO: why is this present in the Fortran?
@@ -185,10 +181,14 @@ function main(args::Vector{String}=String[])
     if ctrl.mat_flag || ctrl.ode_flag
         if ctrl.verbose
             println("Run parameters:")
-            println("   q0 = $(@sprintf("%.3f", equil.params.q0)), qmin = $(@sprintf("%.3f", equil.params.qmin)), qmax = $(@sprintf("%.3f", equil.params.qmax)), q95 = $(@sprintf("%.3f", equil.params.q95))")
+            println(
+                "   q0 = $(@sprintf("%.3f", equil.params.q0)), qmin = $(@sprintf("%.3f", equil.params.qmin)), qmax = $(@sprintf("%.3f", equil.params.qmax)), q95 = $(@sprintf("%.3f", equil.params.q95))"
+            )
             println("   qlim = $(@sprintf("%.5f", intr.qlim)), psilim = $(@sprintf("%.9f", intr.psilim))")
             println("   betat = $(@sprintf("%.3f", equil.params.betat)), betan = $(@sprintf("%.3f", equil.params.betan)), betap1 = $(@sprintf("%.3f", equil.params.betap1))")
-            println("   mlow = $(@sprintf("%4i", intr.mlow)), mhigh = $(@sprintf("%4i", intr.mhigh)), mpert = $(@sprintf("%4i", intr.mpert)), mband = $(@sprintf("%4i", intr.mband))")
+            println(
+                "   mlow = $(@sprintf("%4i", intr.mlow)), mhigh = $(@sprintf("%4i", intr.mhigh)), mpert = $(@sprintf("%4i", intr.mpert)), mband = $(@sprintf("%4i", intr.mband))"
+            )
             println("   nlow = $(@sprintf("%4i", intr.nlow)), nhigh = $(@sprintf("%4i", intr.nhigh)), npert = $(@sprintf("%4i", intr.npert))")
         end
 
@@ -314,9 +314,8 @@ vacuum data if `vac_flag` is true.
 ### TODOs
 
 Combine spline unpacking if possible, too many extra lines
-
 """
-function write_outputs_to_HDF5(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState, vac::Union{VacuumData, Nothing})
+function write_outputs_to_HDF5(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState, vac::Union{VacuumData,Nothing})
 
     h5open(joinpath(intr.dir_path, ctrl.HDF5_filename), "w") do out_h5
 
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
index e38f1cfa..90d08b6b 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
@@ -7,7 +7,6 @@ using LinearAlgebra
 using Statistics
 
 # Import parent modules
-import ..Spl
 import ..Equilibrium
 import ..ForceFreeStates
 import ..ForceFreeStates: OdeState, VacuumData, ForceFreeStatesInternal
@@ -53,6 +52,7 @@ Computes plasma response to external forcing and calculates singular layer
 coupling metrics.
 
 ## Arguments
+
   - `equil`: Equilibrium solution from Equilibrium module
   - `ForceFreeStates_results`: Stability calculation results from ForceFreeStates module
   - `vac_data`: Vacuum response data from ForceFreeStates free boundary calculation
@@ -62,18 +62,20 @@ coupling metrics.
   - `intr`: Internal state variables
 
 ## Returns
+
   - `PerturbedEquilibriumState`: Calculation results
 
 ## Workflow
-1. Load forcing data from file
-2. Compute plasma response (if enabled)
-3. Calculate singular coupling metrics (if enabled)
-4. Output eigenmode fields (if enabled)
+
+ 1. Load forcing data from file
+ 2. Compute plasma response (if enabled)
+ 3. Calculate singular coupling metrics (if enabled)
+ 4. Output eigenmode fields (if enabled)
 """
 function compute_perturbed_equilibrium(
     equil::Equilibrium.PlasmaEquilibrium,
     ForceFreeStates_results::OdeState,
-    vac_data::Union{VacuumData, Nothing},
+    vac_data::Union{VacuumData,Nothing},
     ffs_intr::ForceFreeStates.ForceFreeStatesInternal,
     ft_ctrl::ForcingTerms.ForcingTermsControl,
     ctrl::PerturbedEquilibriumControl,
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index 924b8ba2..b88e1dd8 100644
--- a/src/PerturbedEquilibrium/SingularCoupling.jl
+++ b/src/PerturbedEquilibrium/SingularCoupling.jl
@@ -423,17 +423,18 @@ function compute_current_density(
     last_delpsi = 0.0
     last_jcfun = 0.0
 
+    hint2d = (Ref(1), Ref(1))  # Shared 2D hint for hot loop optimization
     for itheta in 0:mthsurf
         # Theta coordinate normalized to [0, 1]
         theta = itheta / mthsurf
 
         # Evaluate bicubic splines with derivatives at (psi, theta)
         # New API uses separate interpolants for each component
-        r2 = equil.rzphi_rsquared((psi, theta))           # rfac²
-        deta = equil.rzphi_offset((psi, theta))           # angle offset
-        jac = equil.rzphi_jac((psi, theta))               # Jacobian
-        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1))  # ∂(rfac²)/∂theta
-        deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1))  # ∂(deta)/∂theta
+        r2 = equil.rzphi_rsquared((psi, theta); hint=hint2d)           # rfac²
+        deta = equil.rzphi_offset((psi, theta); hint=hint2d)           # angle offset
+        jac = equil.rzphi_jac((psi, theta); hint=hint2d)               # Jacobian
+        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1), hint=hint2d)  # ∂(rfac²)/∂theta
+        deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1), hint=hint2d)  # ∂(deta)/∂theta
 
         rfac = sqrt(abs(r2))
         fy_rfac2 = r2_y
@@ -767,16 +768,17 @@ function compute_surface_area(
     last_jac = 0.0
     last_delpsi = 0.0
 
+    hint2d = (Ref(1), Ref(1))  # Shared 2D hint for hot loop optimization
     for itheta in 0:mthsurf
         # Theta coordinate normalized to [0, 1]
         theta = itheta / mthsurf
 
         # Evaluate bicubic splines with derivatives at (psi, theta)
-        r2 = equil.rzphi_rsquared((psi, theta))
-        jac = equil.rzphi_jac((psi, theta))
-        deta = equil.rzphi_offset((psi, theta))
-        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1))
-        deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1))
+        r2 = equil.rzphi_rsquared((psi, theta); hint=hint2d)
+        jac = equil.rzphi_jac((psi, theta); hint=hint2d)
+        deta = equil.rzphi_offset((psi, theta); hint=hint2d)
+        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1), hint=hint2d)
+        deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1), hint=hint2d)
 
         # Compute rfac
         rfac = sqrt(abs(r2))
diff --git a/src/Splines/FastInterpolationsAdaptor.jl b/src/Splines/FastInterpolationsAdaptor.jl
deleted file mode 100644
index 82bec8fb..00000000
--- a/src/Splines/FastInterpolationsAdaptor.jl
+++ /dev/null
@@ -1,223 +0,0 @@
-"""
-FastInterpolationsAdaptor - Helpers and types for FastInterpolations.jl in JPEC
-
-This module provides:
-- cumulative_integral: Exact spline integration (uses cubic spline coefficients)
-- total_integral: Optimized final-value-only integration
-- Re-exports of FastInterpolations boundary condition types
-
-## Boundary Conditions
-
-Use FastInterpolations' native BC types:
-- `NaturalBC()`: f''(x) = 0 at endpoints
-- `PeriodicBC()`: Periodic (requires fs[end] == fs[1])
-- `CubicFit()`: Automatic endpoint derivative estimation (4-point fit)
-
-## Performance Notes
-
-For monotonic access patterns (typical in ODE integration), use `hint=Ref(1)` with
-`search=LinearBinary()` for optimal performance (~4 ns/pt).
-"""
-
-using FastInterpolations
-using LinearAlgebra
-
-# Re-export FastInterpolations BC types for convenience
-const NaturalBC = FastInterpolations.NaturalBC
-const PeriodicBC = FastInterpolations.PeriodicBC
-const CubicFit = FastInterpolations.CubicFit
-const BCPair = FastInterpolations.BCPair
-const Deriv1 = FastInterpolations.Deriv1
-const Deriv2 = FastInterpolations.Deriv2
-
-# =============================================================================
-# Exact Spline Integration
-# =============================================================================
-
-"""
-    cumulative_integral(xs, fs; bc=NaturalBC()) -> Vector/Matrix
-
-Compute exact cumulative integral of cubic spline fitted to (xs, fs).
-Returns array where result[1] = 0 and result[i+1] = result[i] + ∫_{xs[i]}^{xs[i+1]} S(x) dx,
-where S(x) is the cubic spline interpolant.
-
-For FastInterpolations' moment formulation, the exact integral over interval [x_i, x_{i+1}] is:
-∫S(x)dx = h/2 * (y_i + y_{i+1}) - h³/24 * (z_i + z_{i+1})
-where h is the interval width, y are function values, and z are second derivatives (moments).
-
-This is more accurate than trapezoidal rule, matching the Fortran GPEC behavior.
-
-# Arguments
-
-  - `xs`: Grid points (sorted ascending)
-  - `fs`: Function values (vector or matrix with columns as quantities)
-  - `bc`: Boundary condition for spline fitting (default: NaturalBC())
-"""
-function cumulative_integral(xs::AbstractVector{Float64}, fs::AbstractVector{T}; bc=NaturalBC()) where {T}
-    npts = length(xs)
-    fsi = zeros(T, npts)
-
-    # Create spline to get second derivative coefficients
-    itp = cubic_interp(xs, Vector{Float64}(fs); bc=bc)
-    y = itp.y
-    z = itp.z
-    h = itp.cache.spacing.h  # interval widths: h[i] = xs[i+1] - xs[i]
-
-    @inbounds for i in 1:(npts-1)
-        hi = h[i]
-        # Exact integral: h/2*(y_i + y_{i+1}) - h³/24*(z_i + z_{i+1})
-        fsi[i+1] = fsi[i] + hi/2 * (y[i] + y[i+1]) - hi^3/24 * (z[i] + z[i+1])
-    end
-    return fsi
-end
-
-function cumulative_integral(xs::AbstractVector{Float64}, fs::AbstractMatrix{T}; bc=NaturalBC()) where {T}
-    npts, nqty = size(fs)
-    fsi = zeros(T, npts, nqty)
-
-    @inbounds for k in 1:nqty
-        # Create spline for this column
-        itp = cubic_interp(xs, Vector{Float64}(fs[:, k]); bc=bc)
-        y = itp.y
-        z = itp.z
-        h = itp.cache.spacing.h
-
-        for i in 1:(npts-1)
-            hi = h[i]
-            fsi[i+1, k] = fsi[i, k] + hi/2 * (y[i] + y[i+1]) - hi^3/24 * (z[i] + z[i+1])
-        end
-    end
-    return fsi
-end
-
-"""
-    integrate_spline(itp::CubicInterpolant) -> Float64
-
-Compute exact integral of cubic spline over its entire domain.
-"""
-function integrate_spline(itp::FastInterpolations.CubicInterpolant{T}) where {T}
-    y = itp.y
-    z = itp.z
-    h = itp.cache.spacing.h
-    n = length(y)
-
-    total = zero(T)
-    @inbounds for i in 1:(n-1)
-        hi = h[i]
-        total += hi/2 * (y[i] + y[i+1]) - hi^3/24 * (z[i] + z[i+1])
-    end
-    return total
-end
-
-# =============================================================================
-# Total Integral (Optimized for final value only)
-# =============================================================================
-
-"""
-    total_integral(xs, fs; bc=NaturalBC()) -> Float64 or Vector{Float64}
-
-Compute the total integral of a cubic spline over the entire domain [xs[1], xs[end]].
-This is an optimized version of `cumulative_integral(xs, fs)[end, :]` that avoids
-allocating the full cumulative array.
-
-For the exact spline integration formula over interval [x_i, x_{i+1}]:
-integral = h/2 * (y_i + y_{i+1}) - h^3/24 * (z_i + z_{i+1})
-where h is the interval width, y are function values, and z are second derivatives (moments).
-
-# Arguments
-
-  - `xs::AbstractVector{Float64}`: Grid points (sorted ascending)
-  - `fs`: Function values - vector (returns scalar) or matrix (returns vector, one value per column)
-  - `bc`: Boundary condition for spline fitting (default: NaturalBC())
-
-# Returns
-
-  - For vector `fs`: Returns a scalar Float64
-  - For matrix `fs`: Returns a Vector{Float64} of length `size(fs, 2)`
-
-# Performance
-
-Avoids allocating the full (npts,) or (npts, nqty) cumulative integral array.
-For matrix input, uses CubicSeriesInterpolant internally for efficiency.
-
-# Example
-
-```julia
-xs = collect(0.0:0.01:1.0)
-fs = xs .^ 2
-total = total_integral(xs, fs)  # Returns approx 1/3
-
-# Matrix version
-Y = hcat(xs .^ 2, sin.(xs))
-totals = total_integral(xs, Y)  # Returns [1/3, 1-cos(1)] approximately    # Create spline to get second derivative coefficients
-```
-"""
-function total_integral(xs::AbstractVector{Float64}, fs::AbstractVector{Float64}; bc=NaturalBC())
-    # Create spline to get second derivative coefficients
-    itp = cubic_interp(xs, fs; bc=bc)
-    return integrate_spline(itp)
-end
-
-"""
-    total_integral(xs, fs::AbstractMatrix; bc=NaturalBC()) -> Vector{Float64}
-
-Matrix version: compute total integral for each column of fs.
-Returns a vector of length `size(fs, 2)`.
-"""
-function total_integral(xs::AbstractVector{Float64}, fs::AbstractMatrix{Float64}; bc=NaturalBC())
-    npts, nqty = size(fs)
-
-    # Use CubicSeriesInterpolant for efficient multi-quantity handling
-    itp = cubic_interp(xs, fs; bc=bc)
-    y = itp.y   # (npts, nqty)
-    z = itp.z   # (npts, nqty)
-    h = itp.cache.spacing.h  # (npts-1,)
-
-    # Allocate only the output vector, not the full (npts, nqty) cumulative array
-    result = zeros(Float64, nqty)
-
-    @inbounds for i in 1:(npts-1)
-        hi = h[i]
-        hi_half = hi / 2
-        hi3_24 = hi^3 / 24
-        for k in 1:nqty
-            result[k] += hi_half * (y[i, k] + y[i+1, k]) - hi3_24 * (z[i, k] + z[i+1, k])
-        end
-    end
-
-    return result
-end
-
-"""
-    total_integral!(result::AbstractVector{Float64}, xs, fs::AbstractMatrix; bc=NaturalBC()) -> result
-
-In-place matrix version: compute total integral for each column of fs into pre-allocated result.
-Zero allocations beyond the spline creation.
-"""
-function total_integral!(result::AbstractVector{Float64}, xs::AbstractVector{Float64},
-    fs::AbstractMatrix{Float64}; bc=NaturalBC())
-    npts, nqty = size(fs)
-    @assert length(result) >= nqty "result vector must have at least $nqty elements"
-
-    # Use CubicSeriesInterpolant for efficient multi-quantity handling
-    itp = cubic_interp(xs, fs; bc=bc)
-    y = itp.y   # (npts, nqty)
-    z = itp.z   # (npts, nqty)
-    h = itp.cache.spacing.h  # (npts-1,)
-
-    # Zero the result
-    @inbounds for k in 1:nqty
-        result[k] = 0.0
-    end
-
-    @inbounds for i in 1:(npts-1)
-        hi = h[i]
-        hi_half = hi / 2
-        hi3_24 = hi^3 / 24
-        for k in 1:nqty
-            result[k] += hi_half * (y[i, k] + y[i+1, k]) - hi3_24 * (z[i, k] + z[i+1, k])
-        end
-    end
-
-    return result
-end
diff --git a/src/Splines/Splines.jl b/src/Splines/Splines.jl
deleted file mode 100644
index 5a96f8f0..00000000
--- a/src/Splines/Splines.jl
+++ /dev/null
@@ -1,18 +0,0 @@
-module SplinesMod
-
-# Pure Julia spline implementations
-# FastInterpolationsAdaptor.jl provides helpers for FastInterpolations CubicInterpolant/CubicSeriesInterpolant
-# The codebase uses native interpolants from FastInterpolations directly
-include("FastInterpolationsAdaptor.jl")
-
-# Re-export FastInterpolations search strategies for use with CubicInterpolant hint/search kwargs
-using FastInterpolations: LinearBinary, Binary, HintedBinary
-
-# Exports
-export evaluate!, deriv1!, integrate!, cumulative_integral, integrate_spline, total_integral, total_integral!
-export LinearBinary, Binary, HintedBinary  # Search strategies for hint-based evaluation
-
-# Export boundary condition types
-export NaturalBC, PeriodicBC, CubicFit, BCPair, Deriv1, Deriv2
-
-end
diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl
index 1da6bc1d..ee3c06d2 100644
--- a/src/Vacuum/Kernel2D.jl
+++ b/src/Vacuum/Kernel2D.jl
@@ -109,8 +109,8 @@ function kernel!(
     simpson_weights = dtheta / 3 .* [(k == 1 || k == nsrc) ? 1 : (iseven(k) ? 4 : 2) for k in 1:nsrc]
 
     # Set up periodic splines used for off-grid Gaussian quadrature points
-    spline_x = cubic_interp(theta_grid, source.x; bc=PeriodicBC(; endpoint=:exclusive))
-    spline_z = cubic_interp(theta_grid, source.z; bc=PeriodicBC(; endpoint=:exclusive))
+    spline_x = cubic_interp(theta_grid, source.x; bc=PeriodicBC(; endpoint=:exclusive, period=2π))
+    spline_z = cubic_interp(theta_grid, source.z; bc=PeriodicBC(; endpoint=:exclusive, period=2π))
     d1_spline_x = deriv1(spline_x)
     d1_spline_z = deriv1(spline_z)
 
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 1b848b53..c24acefc 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -4,7 +4,6 @@ using TOML, SpecialFunctions, LinearAlgebra, Printf
 using FastInterpolations: cubic_interp, deriv1, PeriodicBC, NaturalBC
 
 # Import parent modules
-import ..Spl
 import ..Equilibrium
 
 # Import FourierTransforms utility for coefficient calculation and transforms
diff --git a/src/Vacuum/VacuumFromEquilibrium.jl b/src/Vacuum/VacuumFromEquilibrium.jl
index 0b0481dd..93970390 100644
--- a/src/Vacuum/VacuumFromEquilibrium.jl
+++ b/src/Vacuum/VacuumFromEquilibrium.jl
@@ -48,11 +48,12 @@ function extract_plasma_surface_at_psi(equil::Equilibrium.PlasmaEquilibrium, ψ:
     ν = zeros(mtheta)
 
     # Evaluate equilibrium around the flux surface
+    hint2d = (Ref(1), Ref(1))
     for (i, θ_sfl) in enumerate(equil.rzphi_ys)
         # Get minor radius, geometric poloidal angle, and toroidal angle offset
-        r_minor[i] = sqrt(equil.rzphi_rsquared((ψ, θ_sfl)))
-        θ_cyl[i] = 2π * (θ_sfl + equil.rzphi_offset((ψ, θ_sfl)))
-        ν[i] = equil.rzphi_nu((ψ, θ_sfl))
+        r_minor[i] = sqrt(equil.rzphi_rsquared((ψ, θ_sfl); hint=hint2d))
+        θ_cyl[i] = 2π * (θ_sfl + equil.rzphi_offset((ψ, θ_sfl); hint=hint2d))
+        ν[i] = equil.rzphi_nu((ψ, θ_sfl); hint=hint2d)
     end
 
     # Compute R and Z using the computed cylindrical coordinates
diff --git a/test/runtests_spline.jl b/test/runtests_spline.jl
deleted file mode 100644
index b2de24d3..00000000
--- a/test/runtests_spline.jl
+++ /dev/null
@@ -1,107 +0,0 @@
-@testset "Test Spline Module" begin
-    @info "Testing pure Julia spline implementations"
-
-    @testset "MultiQuantityProfile - Polynomial" begin
-        # Make x^3 spline - cubic splines can exactly represent cubics
-        xs = collect(range(1.0; stop=2.0, length=21))
-        fx3 = xs .^ 3
-        fs = reshape(fx3, :, 1)  # Make it a matrix
-        spline = JPEC.Spl.MultiQuantityProfile(xs, fs)
-
-        # Test value interpolation accuracy (primary use case)
-        xs_fine = collect(range(1.1; stop=1.9, length=20))
-        for x in xs_fine
-            f = JPEC.Spl.evaluate!(spline, x)
-            @test abs(f[1] - x^3) < 1e-4  # Value interpolation
-        end
-
-        # Test first derivative (commonly used)
-        f1 = JPEC.Spl.deriv1!(spline, 1.5)
-        f = JPEC.Spl.evaluate!(spline, 1.5)
-        @test abs(f[1] - 1.5^3) < 1e-4
-        @test abs(f1[1] - 3 * 1.5^2) < 0.1
-    end
-
-    @testset "MultiQuantityProfile - Sine" begin
-        # Make sine spline
-        xs = collect(range(Float64(pi) / 8; stop=3 * Float64(pi) / 8, length=50))
-        fs = reshape(sin.(xs), :, 1)
-        spline = JPEC.Spl.MultiQuantityProfile(xs, fs)
-
-        # Test value interpolation
-        xs_fine = collect(range(Float64(pi) / 6; stop=Float64(pi) / 3, length=20))
-        for x in xs_fine
-            f = JPEC.Spl.evaluate!(spline, x)
-            @test abs(f[1] - sin(x)) < 1e-6
-        end
-
-        # Test first derivative
-        f1 = JPEC.Spl.deriv1!(spline, Float64(pi) / 4)
-        f = JPEC.Spl.evaluate!(spline, Float64(pi) / 4)
-        @test abs(f[1] - sin(Float64(pi) / 4)) < 1e-6
-        @test abs(f1[1] - cos(Float64(pi) / 4)) < 1e-4
-    end
-
-    # NOTE: BicubicSpline has been removed - use native FastInterpolations.cubic_interp instead
-    # @testset "BicubicSpline" begin
-    #     ... (tests removed)
-    # end
-end
-
-# NOTE: BicubicSpline has been removed - use native FastInterpolations.cubic_interp instead
-# @testset "BicubicSpline - PeriodicBC" begin
-#     ... (tests removed)
-# end
-
-@testset "Exact Spline Integration" begin
-    # Test that cumulative_integral uses exact spline integration formula
-    # For periodic functions with PeriodicBC, the spline can match well
-
-    # Test with periodic function (sin) and PeriodicBC - spline fits well
-    xs_periodic = collect(range(0.0, 2π; length=65))  # Closed grid
-    fs_periodic = sin.(xs_periodic)
-    fs_periodic[end] = fs_periodic[1]  # Ensure closed
-
-    fsi_periodic = JPEC.Spl.cumulative_integral(xs_periodic, fs_periodic; bc=JPEC.Spl.PeriodicBC())
-
-    # ∫sin(x)dx = -cos(x), so cumulative from 0 should be -cos(x) + cos(0) = 1 - cos(x)
-    exact_periodic = 1.0 .- cos.(xs_periodic)
-
-    # Spline integration should be very accurate for smooth periodic functions
-    @test maximum(abs.(fsi_periodic .- exact_periodic)) < 1e-6
-
-    # Compare accuracy: spline vs trapezoidal for periodic case
-    fsi_trap = zeros(length(xs_periodic))
-    for i in 1:(length(xs_periodic)-1)
-        h = xs_periodic[i+1] - xs_periodic[i]
-        fsi_trap[i+1] = fsi_trap[i] + h * (fs_periodic[i] + fs_periodic[i+1]) / 2
-    end
-    trap_error = maximum(abs.(fsi_trap .- exact_periodic))
-    spline_error = maximum(abs.(fsi_periodic .- exact_periodic))
-
-    # Exact spline integration should be more accurate than trapezoidal
-    @test spline_error < trap_error
-
-    # Test integrate_spline function for total integral
-    using FastInterpolations: cubic_interp
-    itp = cubic_interp(xs_periodic, fs_periodic; bc=JPEC.Spl.PeriodicBC())
-    total_integral = JPEC.Spl.integrate_spline(itp)
-
-    # ∫₀^{2π} sin(x)dx = 0
-    @test abs(total_integral) < 1e-10
-end
-
-@testset "Empty Spline Constructors" begin
-    @info "Testing empty spline constructors for type stability"
-
-    # Test empty MultiQuantityProfile
-    empty_mqp = JPEC.Spl.empty_MultiQuantityProfile()
-    @test length(empty_mqp.xs) >= 4
-    @test typeof(empty_mqp) <: JPEC.Spl.MultiQuantityProfile
-
-    # Test empty BicubicSpline
-    empty_bcs = JPEC.Spl.empty_BicubicSpline()
-    @test length(empty_bcs.xs) >= 2
-    @test length(empty_bcs.ys) >= 2
-    @test typeof(empty_bcs) <: JPEC.Spl.BicubicSpline
-end