diff --git a/CHANGELOG b/CHANGELOG.md
similarity index 94%
rename from CHANGELOG
rename to CHANGELOG.md
index ad3f236..d730bd3 100644
--- a/CHANGELOG
+++ b/CHANGELOG.md
@@ -26,6 +26,15 @@ Everything below is merged to `main` but not yet tagged/released.
   every dimensionality and included in the corresponding `Strategy*Enum` types.
 
 ### Changed
+- **Breaking:** `find_nearest_index` is renamed to `locate_lower_index` and, along with
+  the other grid/index search helpers (`step_index` -> `locate_step_index`,
+  `uniform_lower_index` -> `locate_lower_index_uniform`, `exact_index`,
+  `check_uniform_grid`), moves from `strategy::traits` to a new `strategy::utils`
+  module — `traits` now holds only the `Strategy1D`/`2D`/`3D`/`ND` trait definitions.
+  No deprecation shim, matching the other breaking renames in this release.
+  `locate_lower_index` also now clamps out-of-range points to `[0, len - 2]` itself,
+  rather than relying on each `Linear` call site to inline the same clamp before
+  calling it.
 - **Breaking:** `LeftNearest` and `RightNearest` are removed. Migrate to
   `Step::from(StepDirection::Lower)` / `Step::from(StepDirection::Upper)`, or the
   leaner `StepLower` / `StepUpper` markers.
diff --git a/LICENSE b/LICENSE.md
similarity index 96%
rename from LICENSE
rename to LICENSE.md
index fde28a4..d36f3b3 100644
--- a/LICENSE
+++ b/LICENSE.md
@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2024, Alliance for Sustainable Energy, LLC
+Copyright (c) 2024, Alliance for Energy Innovation, LLC
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
diff --git a/src/interpolator/n/strategies.rs b/src/interpolator/n/strategies.rs
index 8d8c046..b55ce4d 100644
--- a/src/interpolator/n/strategies.rs
+++ b/src/interpolator/n/strategies.rs
@@ -27,15 +27,7 @@ where
             if grid[dim].is_empty() {
                 continue;
             }
-            // Binary search for an exact match: find_nearest_index returns the lower bracket,
-            // so the point can only be exactly equal to grid[lower] or grid[lower+1].
-            let lower = if &point[dim] < grid[dim].first().unwrap() {
-                0
-            } else if &point[dim] > grid[dim].last().unwrap() {
-                grid[dim].len() - 2
-            } else {
-                find_nearest_index(grid[dim].view(), &point[dim])
-            };
+            let lower = locate_lower_index(grid[dim].view(), &point[dim]);
             let pos = exact_index(grid[dim].view(), lower, &point[dim]);
             if let Some(pos) = pos {
                 point.remove(dim);
@@ -57,13 +49,7 @@ where
         for dim in 0..n {
             // Extrapolation is checked previously in Interpolator::interpolate,
             // meaning by now, point is within grid bounds or extrapolation is enabled
-            let lower_idx = if &point[dim] < grid[dim].first().unwrap() {
-                0
-            } else if &point[dim] > grid[dim].last().unwrap() {
-                grid[dim].len() - 2
-            } else {
-                find_nearest_index(grid[dim].view(), &point[dim])
-            };
+            let lower_idx = locate_lower_index(grid[dim].view(), &point[dim]);
             let interp_diff = (point[dim] - grid[dim][lower_idx])
                 / (grid[dim][lower_idx + 1] - grid[dim][lower_idx]);
             lower_idxs.push(lower_idx);
@@ -123,7 +109,8 @@ where
         let mut interp_diffs = Vec::with_capacity(n);
         for (grid_dim, &point_dim) in data.grid.iter().zip(point.iter()) {
             let step = grid_dim[1] - grid_dim[0];
-            let lower_idx = uniform_lower_index(grid_dim[0], step, grid_dim.len(), point_dim);
+            let lower_idx =
+                locate_lower_index_uniform(grid_dim[0], step, grid_dim.len(), point_dim);
             let diff = (point_dim - grid_dim[lower_idx]) / step;
             lower_idxs.push(lower_idx);
             interp_diffs.push(diff);
@@ -169,7 +156,7 @@ where
         // dimensionality reduction needed — the distance comparison handles exact matches correctly.
         let mut idx = vec![0usize; n];
         for dim in 0..n {
-            let lower_idx = find_nearest_index(data.grid[dim].view(), &point[dim]);
+            let lower_idx = locate_lower_index(data.grid[dim].view(), &point[dim]);
             idx[dim] = if point[dim] - data.grid[dim][lower_idx]
                 < data.grid[dim][lower_idx + 1] - point[dim]
             {
@@ -212,7 +199,7 @@ where
         let n = data.values.ndim();
         let mut idx = vec![0usize; n];
         for dim in 0..n {
-            idx[dim] = step_index(self.dir(dim), data.grid[dim].view(), &point[dim]);
+            idx[dim] = locate_step_index(self.dir(dim), data.grid[dim].view(), &point[dim]);
         }
         Ok(data.values.view()[idx.as_slice()])
     }
@@ -236,7 +223,7 @@ where
         let n = data.values.ndim();
         let mut idx = vec![0usize; n];
         for dim in 0..n {
-            idx[dim] = step_index(StepDirection::Lower, data.grid[dim].view(), &point[dim]);
+            idx[dim] = locate_step_index(StepDirection::Lower, data.grid[dim].view(), &point[dim]);
         }
         Ok(data.values.view()[idx.as_slice()])
     }
@@ -259,7 +246,7 @@ where
         let n = data.values.ndim();
         let mut idx = vec![0usize; n];
         for dim in 0..n {
-            idx[dim] = step_index(StepDirection::Upper, data.grid[dim].view(), &point[dim]);
+            idx[dim] = locate_step_index(StepDirection::Upper, data.grid[dim].view(), &point[dim]);
         }
         Ok(data.values.view()[idx.as_slice()])
     }
diff --git a/src/interpolator/one/strategies.rs b/src/interpolator/one/strategies.rs
index c939105..ab8c258 100644
--- a/src/interpolator/one/strategies.rs
+++ b/src/interpolator/one/strategies.rs
@@ -13,16 +13,13 @@ where
     ) -> Result<D::Elem, InterpolateError> {
         // Extrapolation is checked previously in Interpolator::interpolate,
         // meaning by now, point is within grid bounds or extrapolation is enabled
-        let x_l = if &point[0] < data.grid[0].first().unwrap() {
-            0
-        } else if &point[0] > data.grid[0].last().unwrap() {
-            data.grid[0].len() - 2
-        } else {
-            find_nearest_index(data.grid[0].view(), &point[0])
-        };
-        let x_u = x_l + 1;
-        let x_diff = (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-        Ok(data.values[x_l] * (D::Elem::one() - x_diff) + data.values[x_u] * x_diff)
+        match locate_axis(data.grid[0].view(), &point[0]) {
+            AxisLocation::Exact(i) => Ok(data.values[i]),
+            AxisLocation::Interp { lower, frac } => {
+                let upper = lower + 1;
+                Ok(data.values[lower] * (D::Elem::one() - frac) + data.values[upper] * frac)
+            }
+        }
     }
 
     /// Returns `true`.
@@ -48,7 +45,7 @@ where
     ) -> Result<D::Elem, InterpolateError> {
         let grid = data.grid[0].view();
         let step = grid[1] - grid[0];
-        let x_l = uniform_lower_index(grid[0], step, grid.len(), point[0]);
+        let x_l = locate_lower_index_uniform(grid[0], step, grid.len(), point[0]);
         let x_u = x_l + 1;
         let x_diff = (point[0] - grid[x_l]) / step;
         Ok(data.values[x_l] * (D::Elem::one() - x_diff) + data.values[x_u] * x_diff)
@@ -70,7 +67,7 @@ where
         data: &InterpData1D<D>,
         point: &[D::Elem; 1],
     ) -> Result<D::Elem, InterpolateError> {
-        let x_l = find_nearest_index(data.grid[0].view(), &point[0]);
+        let x_l = locate_lower_index(data.grid[0].view(), &point[0]);
         let x_u = x_l + 1;
         let i = if point[0] - data.grid[0][x_l] < data.grid[0][x_u] - point[0] {
             x_l
@@ -107,7 +104,7 @@ where
         data: &InterpData1D<D>,
         point: &[D::Elem; 1],
     ) -> Result<D::Elem, InterpolateError> {
-        Ok(data.values[step_index(self.dir(0), data.grid[0].view(), &point[0])])
+        Ok(data.values[locate_step_index(self.dir(0), data.grid[0].view(), &point[0])])
     }
 
     fn allow_extrapolate(&self) -> bool {
@@ -125,7 +122,7 @@ where
         data: &InterpData1D<D>,
         point: &[D::Elem; 1],
     ) -> Result<D::Elem, InterpolateError> {
-        Ok(data.values[step_index(StepDirection::Lower, data.grid[0].view(), &point[0])])
+        Ok(data.values[locate_step_index(StepDirection::Lower, data.grid[0].view(), &point[0])])
     }
 
     /// Returns `false`.
@@ -144,7 +141,7 @@ where
         data: &InterpData1D<D>,
         point: &[D::Elem; 1],
     ) -> Result<D::Elem, InterpolateError> {
-        Ok(data.values[step_index(StepDirection::Upper, data.grid[0].view(), &point[0])])
+        Ok(data.values[locate_step_index(StepDirection::Upper, data.grid[0].view(), &point[0])])
     }
 
     /// Returns `false`.
diff --git a/src/interpolator/three/strategies.rs b/src/interpolator/three/strategies.rs
index 6b6b864..04aed65 100644
--- a/src/interpolator/three/strategies.rs
+++ b/src/interpolator/three/strategies.rs
@@ -12,103 +12,144 @@ where
         point: &[D::Elem; 3],
     ) -> Result<D::Elem, InterpolateError> {
         // Extrapolation is checked previously in Interpolator::interpolate,
-        // meaning by now, point is within grid bounds or extrapolation is enabled
-        let lowers: [usize; 3] = std::array::from_fn(|dim| {
-            if &point[dim] < data.grid[dim].first().unwrap() {
-                0
-            } else if &point[dim] > data.grid[dim].last().unwrap() {
-                data.grid[dim].len() - 2
-            } else {
-                find_nearest_index(data.grid[dim].view(), &point[dim])
+        // meaning by now, point is within grid bounds or extrapolation is enabled.
+        //
+        // Short-circuit if the point lies exactly on a grid coordinate in one or more
+        // dimensions, reducing value lookups from 8 down to 4, 2, or 1.
+        match (
+            locate_axis(data.grid[0].view(), &point[0]),
+            locate_axis(data.grid[1].view(), &point[1]),
+            locate_axis(data.grid[2].view(), &point[2]),
+        ) {
+            (AxisLocation::Exact(i), AxisLocation::Exact(j), AxisLocation::Exact(k)) => {
+                Ok(data.values[[i, j, k]])
             }
-        });
-        let x_l = lowers[0];
-        let x_u = x_l + 1;
-        let y_l = lowers[1];
-        let y_u = y_l + 1;
-        let z_l = lowers[2];
-        let z_u = z_l + 1;
-
-        // Short-circuit if the point lies exactly on a grid coordinate in one or more dimensions,
-        // reducing value lookups from 8 down to 4, 2, or 1.
-        let x_exact = exact_index(data.grid[0].view(), x_l, &point[0]);
-        let y_exact = exact_index(data.grid[1].view(), y_l, &point[1]);
-        let z_exact = exact_index(data.grid[2].view(), z_l, &point[2]);
-        match (x_exact, y_exact, z_exact) {
-            (Some(i), Some(j), Some(k)) => return Ok(data.values[[i, j, k]]),
-            (Some(i), Some(j), None) => {
-                let z_diff =
-                    (point[2] - data.grid[2][z_l]) / (data.grid[2][z_u] - data.grid[2][z_l]);
-                return Ok(data.values[[i, j, z_l]] * (D::Elem::one() - z_diff)
-                    + data.values[[i, j, z_u]] * z_diff);
+            (
+                AxisLocation::Exact(i),
+                AxisLocation::Exact(j),
+                AxisLocation::Interp {
+                    lower: z_l,
+                    frac: z_diff,
+                },
+            ) => {
+                let z_u = z_l + 1;
+                Ok(data.values[[i, j, z_l]] * (D::Elem::one() - z_diff)
+                    + data.values[[i, j, z_u]] * z_diff)
             }
-            (Some(i), None, Some(k)) => {
-                let y_diff =
-                    (point[1] - data.grid[1][y_l]) / (data.grid[1][y_u] - data.grid[1][y_l]);
-                return Ok(data.values[[i, y_l, k]] * (D::Elem::one() - y_diff)
-                    + data.values[[i, y_u, k]] * y_diff);
+            (
+                AxisLocation::Exact(i),
+                AxisLocation::Interp {
+                    lower: y_l,
+                    frac: y_diff,
+                },
+                AxisLocation::Exact(k),
+            ) => {
+                let y_u = y_l + 1;
+                Ok(data.values[[i, y_l, k]] * (D::Elem::one() - y_diff)
+                    + data.values[[i, y_u, k]] * y_diff)
             }
-            (None, Some(j), Some(k)) => {
-                let x_diff =
-                    (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-                return Ok(data.values[[x_l, j, k]] * (D::Elem::one() - x_diff)
-                    + data.values[[x_u, j, k]] * x_diff);
+            (
+                AxisLocation::Interp {
+                    lower: x_l,
+                    frac: x_diff,
+                },
+                AxisLocation::Exact(j),
+                AxisLocation::Exact(k),
+            ) => {
+                let x_u = x_l + 1;
+                Ok(data.values[[x_l, j, k]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, j, k]] * x_diff)
             }
-            (Some(i), None, None) => {
-                let y_diff =
-                    (point[1] - data.grid[1][y_l]) / (data.grid[1][y_u] - data.grid[1][y_l]);
-                let z_diff =
-                    (point[2] - data.grid[2][z_l]) / (data.grid[2][z_u] - data.grid[2][z_l]);
+            (
+                AxisLocation::Exact(i),
+                AxisLocation::Interp {
+                    lower: y_l,
+                    frac: y_diff,
+                },
+                AxisLocation::Interp {
+                    lower: z_l,
+                    frac: z_diff,
+                },
+            ) => {
+                let y_u = y_l + 1;
+                let z_u = z_l + 1;
                 let f0 = data.values[[i, y_l, z_l]] * (D::Elem::one() - y_diff)
                     + data.values[[i, y_u, z_l]] * y_diff;
                 let f1 = data.values[[i, y_l, z_u]] * (D::Elem::one() - y_diff)
                     + data.values[[i, y_u, z_u]] * y_diff;
-                return Ok(f0 * (D::Elem::one() - z_diff) + f1 * z_diff);
+                Ok(f0 * (D::Elem::one() - z_diff) + f1 * z_diff)
             }
-            (None, Some(j), None) => {
-                let x_diff =
-                    (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-                let z_diff =
-                    (point[2] - data.grid[2][z_l]) / (data.grid[2][z_u] - data.grid[2][z_l]);
+            (
+                AxisLocation::Interp {
+                    lower: x_l,
+                    frac: x_diff,
+                },
+                AxisLocation::Exact(j),
+                AxisLocation::Interp {
+                    lower: z_l,
+                    frac: z_diff,
+                },
+            ) => {
+                let x_u = x_l + 1;
+                let z_u = z_l + 1;
                 let f0 = data.values[[x_l, j, z_l]] * (D::Elem::one() - x_diff)
                     + data.values[[x_u, j, z_l]] * x_diff;
                 let f1 = data.values[[x_l, j, z_u]] * (D::Elem::one() - x_diff)
                     + data.values[[x_u, j, z_u]] * x_diff;
-                return Ok(f0 * (D::Elem::one() - z_diff) + f1 * z_diff);
+                Ok(f0 * (D::Elem::one() - z_diff) + f1 * z_diff)
             }
-            (None, None, Some(k)) => {
-                let x_diff =
-                    (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-                let y_diff =
-                    (point[1] - data.grid[1][y_l]) / (data.grid[1][y_u] - data.grid[1][y_l]);
+            (
+                AxisLocation::Interp {
+                    lower: x_l,
+                    frac: x_diff,
+                },
+                AxisLocation::Interp {
+                    lower: y_l,
+                    frac: y_diff,
+                },
+                AxisLocation::Exact(k),
+            ) => {
+                let x_u = x_l + 1;
+                let y_u = y_l + 1;
                 let f0 = data.values[[x_l, y_l, k]] * (D::Elem::one() - x_diff)
                     + data.values[[x_u, y_l, k]] * x_diff;
                 let f1 = data.values[[x_l, y_u, k]] * (D::Elem::one() - x_diff)
                     + data.values[[x_u, y_u, k]] * x_diff;
-                return Ok(f0 * (D::Elem::one() - y_diff) + f1 * y_diff);
+                Ok(f0 * (D::Elem::one() - y_diff) + f1 * y_diff)
+            }
+            (
+                AxisLocation::Interp {
+                    lower: x_l,
+                    frac: x_diff,
+                },
+                AxisLocation::Interp {
+                    lower: y_l,
+                    frac: y_diff,
+                },
+                AxisLocation::Interp {
+                    lower: z_l,
+                    frac: z_diff,
+                },
+            ) => {
+                let x_u = x_l + 1;
+                let y_u = y_l + 1;
+                let z_u = z_l + 1;
+                // interpolate in the x-direction
+                let f00 = data.values[[x_l, y_l, z_l]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, y_l, z_l]] * x_diff;
+                let f01 = data.values[[x_l, y_l, z_u]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, y_l, z_u]] * x_diff;
+                let f10 = data.values[[x_l, y_u, z_l]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, y_u, z_l]] * x_diff;
+                let f11 = data.values[[x_l, y_u, z_u]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, y_u, z_u]] * x_diff;
+                // interpolate in the y-direction
+                let f0 = f00 * (D::Elem::one() - y_diff) + f10 * y_diff;
+                let f1 = f01 * (D::Elem::one() - y_diff) + f11 * y_diff;
+                // interpolate in the z-direction
+                Ok(f0 * (D::Elem::one() - z_diff) + f1 * z_diff)
             }
-            (None, None, None) => {}
         }
-
-        let x_diff = (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-        // y
-        let y_diff = (point[1] - data.grid[1][y_l]) / (data.grid[1][y_u] - data.grid[1][y_l]);
-        // z
-        let z_diff = (point[2] - data.grid[2][z_l]) / (data.grid[2][z_u] - data.grid[2][z_l]);
-        // interpolate in the x-direction
-        let f00 = data.values[[x_l, y_l, z_l]] * (D::Elem::one() - x_diff)
-            + data.values[[x_u, y_l, z_l]] * x_diff;
-        let f01 = data.values[[x_l, y_l, z_u]] * (D::Elem::one() - x_diff)
-            + data.values[[x_u, y_l, z_u]] * x_diff;
-        let f10 = data.values[[x_l, y_u, z_l]] * (D::Elem::one() - x_diff)
-            + data.values[[x_u, y_u, z_l]] * x_diff;
-        let f11 = data.values[[x_l, y_u, z_u]] * (D::Elem::one() - x_diff)
-            + data.values[[x_u, y_u, z_u]] * x_diff;
-        // interpolate in the y-direction
-        let f0 = f00 * (D::Elem::one() - y_diff) + f10 * y_diff;
-        let f1 = f01 * (D::Elem::one() - y_diff) + f11 * y_diff;
-        // interpolate in the z-direction
-        Ok(f0 * (D::Elem::one() - z_diff) + f1 * z_diff)
     }
 
     /// Returns `true`.
@@ -137,9 +178,9 @@ where
         let x_step = data.grid[0][1] - data.grid[0][0];
         let y_step = data.grid[1][1] - data.grid[1][0];
         let z_step = data.grid[2][1] - data.grid[2][0];
-        let x_l = uniform_lower_index(data.grid[0][0], x_step, data.grid[0].len(), point[0]);
-        let y_l = uniform_lower_index(data.grid[1][0], y_step, data.grid[1].len(), point[1]);
-        let z_l = uniform_lower_index(data.grid[2][0], z_step, data.grid[2].len(), point[2]);
+        let x_l = locate_lower_index_uniform(data.grid[0][0], x_step, data.grid[0].len(), point[0]);
+        let y_l = locate_lower_index_uniform(data.grid[1][0], y_step, data.grid[1].len(), point[1]);
+        let z_l = locate_lower_index_uniform(data.grid[2][0], z_step, data.grid[2].len(), point[2]);
         let x_u = x_l + 1;
         let y_u = y_l + 1;
         let z_u = z_l + 1;
@@ -176,7 +217,7 @@ where
         point: &[D::Elem; 3],
     ) -> Result<D::Elem, InterpolateError> {
         // x
-        let x_l = find_nearest_index(data.grid[0].view(), &point[0]);
+        let x_l = locate_lower_index(data.grid[0].view(), &point[0]);
         let x_u = x_l + 1;
         let i = if point[0] - data.grid[0][x_l] < data.grid[0][x_u] - point[0] {
             x_l
@@ -184,7 +225,7 @@ where
             x_u
         };
         // y
-        let y_l = find_nearest_index(data.grid[1].view(), &point[1]);
+        let y_l = locate_lower_index(data.grid[1].view(), &point[1]);
         let y_u = y_l + 1;
         let j = if point[1] - data.grid[1][y_l] < data.grid[1][y_u] - point[1] {
             y_l
@@ -192,7 +233,7 @@ where
             y_u
         };
         // z
-        let z_l = find_nearest_index(data.grid[2].view(), &point[2]);
+        let z_l = locate_lower_index(data.grid[2].view(), &point[2]);
         let z_u = z_l + 1;
         let k = if point[2] - data.grid[2][z_l] < data.grid[2][z_u] - point[2] {
             z_l
@@ -230,9 +271,9 @@ where
         data: &InterpData3D<D>,
         point: &[D::Elem; 3],
     ) -> Result<D::Elem, InterpolateError> {
-        let i = step_index(self.dir(0), data.grid[0].view(), &point[0]);
-        let j = step_index(self.dir(1), data.grid[1].view(), &point[1]);
-        let k = step_index(self.dir(2), data.grid[2].view(), &point[2]);
+        let i = locate_step_index(self.dir(0), data.grid[0].view(), &point[0]);
+        let j = locate_step_index(self.dir(1), data.grid[1].view(), &point[1]);
+        let k = locate_step_index(self.dir(2), data.grid[2].view(), &point[2]);
         Ok(data.values[[i, j, k]])
     }
 
@@ -252,9 +293,9 @@ where
         data: &InterpData3D<D>,
         point: &[D::Elem; 3],
     ) -> Result<D::Elem, InterpolateError> {
-        let i = step_index(StepDirection::Lower, data.grid[0].view(), &point[0]);
-        let j = step_index(StepDirection::Lower, data.grid[1].view(), &point[1]);
-        let k = step_index(StepDirection::Lower, data.grid[2].view(), &point[2]);
+        let i = locate_step_index(StepDirection::Lower, data.grid[0].view(), &point[0]);
+        let j = locate_step_index(StepDirection::Lower, data.grid[1].view(), &point[1]);
+        let k = locate_step_index(StepDirection::Lower, data.grid[2].view(), &point[2]);
         Ok(data.values[[i, j, k]])
     }
 
@@ -273,9 +314,9 @@ where
         data: &InterpData3D<D>,
         point: &[D::Elem; 3],
     ) -> Result<D::Elem, InterpolateError> {
-        let i = step_index(StepDirection::Upper, data.grid[0].view(), &point[0]);
-        let j = step_index(StepDirection::Upper, data.grid[1].view(), &point[1]);
-        let k = step_index(StepDirection::Upper, data.grid[2].view(), &point[2]);
+        let i = locate_step_index(StepDirection::Upper, data.grid[0].view(), &point[0]);
+        let j = locate_step_index(StepDirection::Upper, data.grid[1].view(), &point[1]);
+        let k = locate_step_index(StepDirection::Upper, data.grid[2].view(), &point[2]);
         Ok(data.values[[i, j, k]])
     }
 
diff --git a/src/interpolator/two/strategies.rs b/src/interpolator/two/strategies.rs
index dc6aa28..ce5434e 100644
--- a/src/interpolator/two/strategies.rs
+++ b/src/interpolator/two/strategies.rs
@@ -12,54 +12,58 @@ where
         point: &[D::Elem; 2],
     ) -> Result<D::Elem, InterpolateError> {
         // Extrapolation is checked previously in Interpolator::interpolate,
-        // meaning by now, point is within grid bounds or extrapolation is enabled
-        let lowers: [usize; 2] = std::array::from_fn(|dim| {
-            if &point[dim] < data.grid[dim].first().unwrap() {
-                0
-            } else if &point[dim] > data.grid[dim].last().unwrap() {
-                data.grid[dim].len() - 2
-            } else {
-                find_nearest_index(data.grid[dim].view(), &point[dim])
+        // meaning by now, point is within grid bounds or extrapolation is enabled.
+        //
+        // Short-circuit if the point lies exactly on a grid coordinate in one or both
+        // dimensions, reducing value lookups from 4 to 2 or 1.
+        match (
+            locate_axis(data.grid[0].view(), &point[0]),
+            locate_axis(data.grid[1].view(), &point[1]),
+        ) {
+            (AxisLocation::Exact(i), AxisLocation::Exact(j)) => Ok(data.values[[i, j]]),
+            (
+                AxisLocation::Exact(i),
+                AxisLocation::Interp {
+                    lower: y_l,
+                    frac: y_diff,
+                },
+            ) => {
+                let y_u = y_l + 1;
+                Ok(data.values[[i, y_l]] * (D::Elem::one() - y_diff)
+                    + data.values[[i, y_u]] * y_diff)
             }
-        });
-        // x
-        let x_l = lowers[0];
-        let x_u = x_l + 1;
-        // y
-        let y_l = lowers[1];
-        let y_u = y_l + 1;
-
-        // Short-circuit if the point lies exactly on a grid coordinate in one or both dimensions,
-        // reducing value lookups from 4 to 2 or 1. find_nearest_index returns the lower bracket,
-        // so exact matches appear at grid[lower] or grid[lower+1].
-        let x_exact = exact_index(data.grid[0].view(), x_l, &point[0]);
-        let y_exact = exact_index(data.grid[1].view(), y_l, &point[1]);
-        match (x_exact, y_exact) {
-            (Some(i), Some(j)) => return Ok(data.values[[i, j]]),
-            (Some(i), None) => {
-                let y_diff =
-                    (point[1] - data.grid[1][y_l]) / (data.grid[1][y_u] - data.grid[1][y_l]);
-                return Ok(data.values[[i, y_l]] * (D::Elem::one() - y_diff)
-                    + data.values[[i, y_u]] * y_diff);
+            (
+                AxisLocation::Interp {
+                    lower: x_l,
+                    frac: x_diff,
+                },
+                AxisLocation::Exact(j),
+            ) => {
+                let x_u = x_l + 1;
+                Ok(data.values[[x_l, j]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, j]] * x_diff)
             }
-            (None, Some(j)) => {
-                let x_diff =
-                    (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-                return Ok(data.values[[x_l, j]] * (D::Elem::one() - x_diff)
-                    + data.values[[x_u, j]] * x_diff);
+            (
+                AxisLocation::Interp {
+                    lower: x_l,
+                    frac: x_diff,
+                },
+                AxisLocation::Interp {
+                    lower: y_l,
+                    frac: y_diff,
+                },
+            ) => {
+                let x_u = x_l + 1;
+                let y_u = y_l + 1;
+                // interpolate in the x-direction
+                let f0 = data.values[[x_l, y_l]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, y_l]] * x_diff;
+                let f1 = data.values[[x_l, y_u]] * (D::Elem::one() - x_diff)
+                    + data.values[[x_u, y_u]] * x_diff;
+                // interpolate in the y-direction
+                Ok(f0 * (D::Elem::one() - y_diff) + f1 * y_diff)
             }
-            (None, None) => {}
         }
-
-        let x_diff = (point[0] - data.grid[0][x_l]) / (data.grid[0][x_u] - data.grid[0][x_l]);
-        let y_diff = (point[1] - data.grid[1][y_l]) / (data.grid[1][y_u] - data.grid[1][y_l]);
-        // interpolate in the x-direction
-        let f0 =
-            data.values[[x_l, y_l]] * (D::Elem::one() - x_diff) + data.values[[x_u, y_l]] * x_diff;
-        let f1 =
-            data.values[[x_l, y_u]] * (D::Elem::one() - x_diff) + data.values[[x_u, y_u]] * x_diff;
-        // interpolate in the y-direction
-        Ok(f0 * (D::Elem::one() - y_diff) + f1 * y_diff)
     }
 
     /// Returns `true`.
@@ -86,8 +90,8 @@ where
     ) -> Result<D::Elem, InterpolateError> {
         let x_step = data.grid[0][1] - data.grid[0][0];
         let y_step = data.grid[1][1] - data.grid[1][0];
-        let x_l = uniform_lower_index(data.grid[0][0], x_step, data.grid[0].len(), point[0]);
-        let y_l = uniform_lower_index(data.grid[1][0], y_step, data.grid[1].len(), point[1]);
+        let x_l = locate_lower_index_uniform(data.grid[0][0], x_step, data.grid[0].len(), point[0]);
+        let y_l = locate_lower_index_uniform(data.grid[1][0], y_step, data.grid[1].len(), point[1]);
         let x_u = x_l + 1;
         let y_u = y_l + 1;
         let x_diff = (point[0] - data.grid[0][x_l]) / x_step;
@@ -116,7 +120,7 @@ where
         point: &[D::Elem; 2],
     ) -> Result<D::Elem, InterpolateError> {
         // x
-        let x_l = find_nearest_index(data.grid[0].view(), &point[0]);
+        let x_l = locate_lower_index(data.grid[0].view(), &point[0]);
         let x_u = x_l + 1;
         let i = if point[0] - data.grid[0][x_l] < data.grid[0][x_u] - point[0] {
             x_l
@@ -124,7 +128,7 @@ where
             x_u
         };
         // y
-        let y_l = find_nearest_index(data.grid[1].view(), &point[1]);
+        let y_l = locate_lower_index(data.grid[1].view(), &point[1]);
         let y_u = y_l + 1;
         let j = if point[1] - data.grid[1][y_l] < data.grid[1][y_u] - point[1] {
             y_l
@@ -162,8 +166,8 @@ where
         data: &InterpData2D<D>,
         point: &[D::Elem; 2],
     ) -> Result<D::Elem, InterpolateError> {
-        let i = step_index(self.dir(0), data.grid[0].view(), &point[0]);
-        let j = step_index(self.dir(1), data.grid[1].view(), &point[1]);
+        let i = locate_step_index(self.dir(0), data.grid[0].view(), &point[0]);
+        let j = locate_step_index(self.dir(1), data.grid[1].view(), &point[1]);
         Ok(data.values[[i, j]])
     }
 
@@ -183,8 +187,8 @@ where
         data: &InterpData2D<D>,
         point: &[D::Elem; 2],
     ) -> Result<D::Elem, InterpolateError> {
-        let i = step_index(StepDirection::Lower, data.grid[0].view(), &point[0]);
-        let j = step_index(StepDirection::Lower, data.grid[1].view(), &point[1]);
+        let i = locate_step_index(StepDirection::Lower, data.grid[0].view(), &point[0]);
+        let j = locate_step_index(StepDirection::Lower, data.grid[1].view(), &point[1]);
         Ok(data.values[[i, j]])
     }
 
@@ -203,8 +207,8 @@ where
         data: &InterpData2D<D>,
         point: &[D::Elem; 2],
     ) -> Result<D::Elem, InterpolateError> {
-        let i = step_index(StepDirection::Upper, data.grid[0].view(), &point[0]);
-        let j = step_index(StepDirection::Upper, data.grid[1].view(), &point[1]);
+        let i = locate_step_index(StepDirection::Upper, data.grid[0].view(), &point[0]);
+        let j = locate_step_index(StepDirection::Upper, data.grid[1].view(), &point[1]);
         Ok(data.values[[i, j]])
     }
 
diff --git a/src/lib.rs b/src/lib.rs
index 7552e04..264081d 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -44,6 +44,7 @@ pub(crate) use interpolator::data::*;
 
 pub(crate) use error::*;
 pub(crate) use strategy::traits::*;
+pub(crate) use strategy::utils::*;
 
 pub(crate) use std::fmt::Debug;
 
diff --git a/src/strategy/mod.rs b/src/strategy/mod.rs
index dcdddaa..87b4821 100644
--- a/src/strategy/mod.rs
+++ b/src/strategy/mod.rs
@@ -4,6 +4,7 @@ use super::*;
 
 pub mod enums;
 pub mod traits;
+pub mod utils;
 
 /// Linear interpolation: <https://en.wikipedia.org/wiki/Linear_interpolation>
 #[derive(Debug, Clone, PartialEq)]
diff --git a/src/strategy/traits.rs b/src/strategy/traits.rs
index b35b2e7..35eeefa 100644
--- a/src/strategy/traits.rs
+++ b/src/strategy/traits.rs
@@ -2,128 +2,6 @@
 
 use super::*;
 
-/// Find nearest index in `arr` left of `target`
-///
-/// This method contains code from RouteE Compass, another open-source NLR-developed tool
-/// <https://www.nlr.gov/transportation/route-energy-prediction-model.html>
-/// <https://github.com/NatLabRockies/routee-compass/>
-pub fn find_nearest_index<T: PartialOrd>(arr: ArrayView1<T>, target: &T) -> usize {
-    if target == arr.last().unwrap() {
-        return arr.len() - 2;
-    }
-
-    let mut low = 0;
-    let mut high = arr.len() - 1;
-
-    while low < high {
-        let mid = low + (high - low) / 2;
-
-        if &arr[mid] >= target {
-            high = mid;
-        } else {
-            low = mid + 1;
-        }
-    }
-
-    if low > 0 && &arr[low] >= target {
-        low - 1
-    } else {
-        low
-    }
-}
-
-/// Returns the step index for `point` in `grid` using the given [`StepDirection`].
-///
-/// Handles all exact grid-point edge cases that arise from [`find_nearest_index`]'s
-/// interval semantics (returning the lower bracket rather than the exact position).
-pub(crate) fn step_index<T: PartialOrd + Copy>(
-    dir: StepDirection,
-    grid: ArrayView1<T>,
-    point: &T,
-) -> usize {
-    match dir {
-        StepDirection::Lower => {
-            let x_l = find_nearest_index(grid, point);
-            // find_nearest_index returns i where grid[i] < point <= grid[i+1] for interior
-            // matches, so an exact match at grid[i+1] gives i instead of i+1. Correct both:
-            if point == grid.last().unwrap() {
-                grid.len() - 1
-            } else if *point == grid[x_l + 1] {
-                x_l + 1
-            } else {
-                x_l
-            }
-        }
-        StepDirection::Upper => {
-            // find_nearest_index returns 0 when point == grid[0], giving x_l+1 = 1
-            // which would skip values[0]. Handle the first-element case explicitly:
-            if point == grid.first().unwrap() {
-                0
-            } else {
-                find_nearest_index(grid, point) + 1
-            }
-        }
-    }
-}
-
-/// Returns the exact grid index if `point` lies on `grid[lower]` or `grid[lower+1]`, else `None`.
-///
-/// Used to short-circuit interpolation when a query point coincides with a grid coordinate.
-pub(crate) fn exact_index<T: PartialOrd>(
-    grid: ArrayView1<T>,
-    lower: usize,
-    point: &T,
-) -> Option<usize> {
-    if grid[lower] == *point {
-        Some(lower)
-    } else if grid[lower + 1] == *point {
-        Some(lower + 1)
-    } else {
-        None
-    }
-}
-
-/// Computes the lower bracket index for a uniformly-spaced grid in O(1).
-///
-/// Equivalent to [`find_nearest_index`] but replaces binary search with direct arithmetic.
-/// Only valid when the grid spacing is uniform — validate with [`check_uniform_grid`] first.
-pub(crate) fn uniform_lower_index<T: Float>(grid0: T, step: T, n: usize, point: T) -> usize {
-    let t = (point - grid0) / step;
-    if t < T::zero() {
-        0
-    } else {
-        t.floor().to_usize().unwrap_or(0).min(n - 2)
-    }
-}
-
-/// Validates that `grid` is uniformly spaced within floating-point tolerance.
-///
-/// Uses a relative tolerance of 1024 × ε to accommodate accumulated floating-point rounding
-/// error in grids constructed from repeated arithmetic.
-pub(crate) fn check_uniform_grid<T: Float>(
-    grid: ArrayView1<T>,
-    dim: usize,
-) -> Result<(), ValidateError> {
-    let step = grid[1] - grid[0];
-    // 1024 * epsilon via 10 doublings — avoids numeric literal casting
-    let tolerance = {
-        let mut tol = T::epsilon();
-        for _ in 0..10 {
-            tol = tol + tol;
-        }
-        step.abs() * tol
-    };
-    for i in 1..grid.len() - 1 {
-        let gap = grid[i + 1] - grid[i];
-        if (gap - step).abs() > tolerance {
-            return Err(ValidateError::Other(format!(
-                "LinearUniform: grid[{dim}] is not uniformly spaced (gap at index {i})"
-            )));
-        }
-    }
-    Ok(())
-}
-
 /// 1-D interpolation strategy.
 pub trait Strategy1D<D>: Debug + DynClone
 where
diff --git a/src/strategy/utils.rs b/src/strategy/utils.rs
new file mode 100644
index 0000000..bb81d34
--- /dev/null
+++ b/src/strategy/utils.rs
@@ -0,0 +1,153 @@
+//! Single-axis primitives for locating a query point within one grid dimension.
+//! Strategies compose these per axis to get 1D/2D/3D/ND behavior; none of them
+//! iterate over dimensions themselves.
+
+use super::*;
+
+/// Returns the lower bracket index for `point` in `grid`: the largest `i` such that
+/// `grid[i] <= point`, clamped to `[0, grid.len() - 2]` for out-of-range points.
+///
+/// This method contains code from RouteE Compass, another open-source NLR-developed tool
+/// <https://www.nlr.gov/transportation/route-energy-prediction-model.html>
+/// <https://github.com/NatLabRockies/routee-compass/>
+pub fn locate_lower_index<T: PartialOrd>(grid: ArrayView1<T>, point: &T) -> usize {
+    if point < grid.first().unwrap() {
+        return 0;
+    }
+    if point >= grid.last().unwrap() {
+        return grid.len() - 2;
+    }
+
+    let mut low = 0;
+    let mut high = grid.len() - 1;
+
+    while low < high {
+        let mid = low + (high - low) / 2;
+
+        if &grid[mid] >= point {
+            high = mid;
+        } else {
+            low = mid + 1;
+        }
+    }
+
+    if low > 0 && &grid[low] >= point {
+        low - 1
+    } else {
+        low
+    }
+}
+
+/// Per-axis locate for linear-family strategies: either an exact grid hit,
+/// or an interior interpolation position.
+pub(crate) enum AxisLocation<T> {
+    Exact(usize),
+    Interp { lower: usize, frac: T },
+}
+
+/// Locates `point` along `grid`, resolving to an exact grid hit or an interpolation
+/// position. Combines [`locate_lower_index`] (search + extrapolation clamp) with
+/// [`exact_index`] (exact-match short-circuit) into the single call linear-family
+/// strategies need per axis.
+pub(crate) fn locate_axis<T: Float>(grid: ArrayView1<T>, point: &T) -> AxisLocation<T> {
+    let lower = locate_lower_index(grid, point);
+    match exact_index(grid, lower, point) {
+        Some(idx) => AxisLocation::Exact(idx),
+        None => {
+            let frac = (*point - grid[lower]) / (grid[lower + 1] - grid[lower]);
+            AxisLocation::Interp { lower, frac }
+        }
+    }
+}
+
+/// Returns the step index for `point` in `grid` using the given [`StepDirection`].
+///
+/// Handles all exact grid-point edge cases that arise from [`locate_lower_index`]'s
+/// interval semantics (returning the lower bracket rather than the exact position).
+pub(crate) fn locate_step_index<T: PartialOrd + Copy>(
+    dir: StepDirection,
+    grid: ArrayView1<T>,
+    point: &T,
+) -> usize {
+    match dir {
+        StepDirection::Lower => {
+            let x_l = locate_lower_index(grid, point);
+            // locate_lower_index returns i where grid[i] < point <= grid[i+1] for interior
+            // matches, so an exact match at grid[i+1] gives i instead of i+1. Correct both:
+            if point == grid.last().unwrap() {
+                grid.len() - 1
+            } else if *point == grid[x_l + 1] {
+                x_l + 1
+            } else {
+                x_l
+            }
+        }
+        StepDirection::Upper => {
+            // locate_lower_index returns 0 when point == grid[0], giving x_l+1 = 1
+            // which would skip values[0]. Handle the first-element case explicitly:
+            if point == grid.first().unwrap() {
+                0
+            } else {
+                locate_lower_index(grid, point) + 1
+            }
+        }
+    }
+}
+
+/// Returns the exact grid index if `point` lies on `grid[lower]` or `grid[lower+1]`, else `None`.
+///
+/// Used to short-circuit interpolation when a query point coincides with a grid coordinate.
+pub(crate) fn exact_index<T: PartialOrd>(
+    grid: ArrayView1<T>,
+    lower: usize,
+    point: &T,
+) -> Option<usize> {
+    if grid[lower] == *point {
+        Some(lower)
+    } else if grid[lower + 1] == *point {
+        Some(lower + 1)
+    } else {
+        None
+    }
+}
+
+/// Computes the lower bracket index for a uniformly-spaced grid in O(1).
+///
+/// Equivalent to [`locate_lower_index`] but replaces binary search with direct arithmetic.
+/// Only valid when the grid spacing is uniform — validate with [`check_uniform_grid`] first.
+pub(crate) fn locate_lower_index_uniform<T: Float>(grid0: T, step: T, n: usize, point: T) -> usize {
+    let t = (point - grid0) / step;
+    if t < T::zero() {
+        0
+    } else {
+        t.floor().to_usize().unwrap_or(0).min(n - 2)
+    }
+}
+
+/// Validates that `grid` is uniformly spaced within floating-point tolerance.
+///
+/// Uses a relative tolerance of 1024 × ε to accommodate accumulated floating-point rounding
+/// error in grids constructed from repeated arithmetic.
+pub(crate) fn check_uniform_grid<T: Float>(
+    grid: ArrayView1<T>,
+    dim: usize,
+) -> Result<(), ValidateError> {
+    let step = grid[1] - grid[0];
+    // 1024 * epsilon via 10 doublings — avoids numeric literal casting
+    let tolerance = {
+        let mut tol = T::epsilon();
+        for _ in 0..10 {
+            tol = tol + tol;
+        }
+        step.abs() * tol
+    };
+    for i in 1..grid.len() - 1 {
+        let gap = grid[i + 1] - grid[i];
+        if (gap - step).abs() > tolerance {
+            return Err(ValidateError::Other(format!(
+                "LinearUniform: grid[{dim}] is not uniformly spaced (gap at index {i})"
+            )));
+        }
+    }
+    Ok(())
+}