diff --git a/crates/coremem/src/geom/units.rs b/crates/coremem/src/geom/units.rs
index ebc82d8..f59cbf3 100644
--- a/crates/coremem/src/geom/units.rs
+++ b/crates/coremem/src/geom/units.rs
@@ -188,6 +188,17 @@ impl Index {
     }
 }
 
+impl Into<Vec3u> for Index {
+    fn into(self) -> Vec3u {
+        self.0
+    }
+}
+impl From<Vec3u> for Index {
+    fn from(v: Vec3u) -> Self {
+        Self(v)
+    }
+}
+
 impl Coord for Index {
     fn to_meters(&self, feature_size: f32) -> Meters {
         Meters(Vec3::from(self.0) * feature_size)
diff --git a/crates/coremem/src/sim/legacy/mod.rs b/crates/coremem/src/sim/legacy/mod.rs
index da84350..93c9524 100644
--- a/crates/coremem/src/sim/legacy/mod.rs
+++ b/crates/coremem/src/sim/legacy/mod.rs
@@ -304,8 +304,8 @@ impl<R: Real, M: Material<R> + cross::mat::Material<R> + Send + Sync> SimState<R
         trace!("apply_stimulus begin");
         Zip::from(ndarray::indices_of(&self.e)).and(&mut self.e).and(&mut self.h).par_for_each(
             |(z, y, x), e, h| {
-                let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size);
-                let densities = stim.at(t_sec, pos_meters);
+                let idx = Index::new(x as u32, y as u32, z as u32);
+                let densities = stim.at(t_sec, feature_size, idx);
                 let value_e = densities.e.cast::<R>() * timestep.cast::<R>();
                 let value_h = densities.h.cast::<R>() * timestep.cast::<R>();
                 *e += value_e;
@@ -792,7 +792,7 @@ impl<R: Real> CellState<R> {
 #[cfg(test)]
 mod test {
     use super::*;
-    use crate::geom::{Cube, Meters, Region, WorldRegion};
+    use crate::geom::{Cube, Region, WorldRegion};
     use crate::sim::legacy::mat::{AnisomorphicConductor, IsomorphicConductor, Pml, Static};
     use crate::stim::Fields;
     use crate::real::{R64, ToFloat as _};
@@ -1170,13 +1170,12 @@ mod test {
 
     struct ConstStim {
         loc: Index,
-        feature_size: f32,
         e: Vec3<f32>,
         h: Vec3<f32>,
     }
     impl Stimulus for ConstStim {
-        fn at(&self, _t_sec: f32, loc: Meters) -> Fields {
-            if loc.to_index(self.feature_size) == self.loc {
+        fn at(&self, _t_sec: f32, _feat_size: f32, loc: Index) -> Fields {
+            if loc == self.loc {
                 Fields { e: self.e, h: self.h }
             } else {
                 Default::default()
@@ -1194,7 +1193,6 @@ mod test {
             //println!("{}", g.exp());
             state.step_multiple(1, &ConstStim {
                 loc: Index::new(0, 0, 1000),
-                feature_size: 1e-6,
                 e: Vec3::new_z(gate*sig/state.timestep()),
                 h: Vec3::zero(),
             });
@@ -1213,7 +1211,6 @@ mod test {
             let gate = 1e9 * 0.5*(1.0 - angle.cos());
             state.step_multiple(1, &ConstStim {
                 loc: Index::new(10, 10, 10),
-                feature_size: 1e-6,
                 e: Vec3::new_z(gate*sig/state.timestep()),
                 h: Vec3::zero(),
             });
@@ -1236,7 +1233,6 @@ mod test {
             let gate = (progress*f32::pi()).sin();
             state.step_multiple(1, &ConstStim {
                 loc: Index::new(100, 0, 0),
-                feature_size: 1e-6,
                 e: Vec3::new_z(gate*sig/state.timestep()),
                 h: Vec3::zero(),
             });
@@ -1312,14 +1308,13 @@ mod test {
         /// Maps location -> (field direction, frequency)
         f: F,
         t_end: f32,
-        feat_size: f32,
         sqrt_vol: f32,
     }
     impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> Stimulus for PmlStim<F> {
-        fn at(&self, t_sec: f32, pos: Meters) -> Fields {
+        fn at(&self, t_sec: f32, _feat_size: f32, loc: Index) -> Fields {
             let angle = t_sec/self.t_end * f32::two_pi();
             let gate = 0.5*(1.0 - angle.cos());
-            let (e, hz) = (self.f)(pos.to_index(self.feat_size));
+            let (e, hz) = (self.f)(loc);
             let sig_angle = angle*hz;
             let sig = sig_angle.sin();
             // TODO: test H component?
@@ -1337,7 +1332,6 @@ mod test {
         let stim = PmlStim {
             f,
             t_end: 50.0 * state.timestep(),
-            feat_size: state.feature_size(),
             sqrt_vol: state.volume().sqrt(),
         };
         for _t in 0..50 {
diff --git a/crates/coremem/src/sim/spirv/mod.rs b/crates/coremem/src/sim/spirv/mod.rs
index 0308a3c..51006cf 100644
--- a/crates/coremem/src/sim/spirv/mod.rs
+++ b/crates/coremem/src/sim/spirv/mod.rs
@@ -2,7 +2,7 @@ use serde::{Deserialize, Serialize};
 use log::{info, trace, warn};
 
 use crate::diagnostics::SyncDiagnostics;
-use crate::geom::{Coord, Index};
+use crate::geom::Index;
 use crate::real::Real;
 use crate::sim::{AbstractSim, Fields, GenericSim, StaticSim};
 use crate::stim::Stimulus;
@@ -269,6 +269,7 @@ where
     }
 }
 
+// used for test
 fn iterate_stim<R, S, F>(
     stim: &S, dim: Index, feature_size: f32, t_sec: f32, timestep: R, mut f: F
 )
@@ -282,8 +283,7 @@ where
         for y in 0..dim.y() {
             for x in 0..dim.x() {
                 let pos_idx = Index::new(x, y, z);
-                let pos_meters = pos_idx.to_meters(feature_size);
-                let densities = stim.at(t_sec, pos_meters);
+                let densities = stim.at(t_sec, feature_size, pos_idx);
                 let (value_e, value_h) = (densities.e.cast::<R>() * timestep, densities.h.cast::<R>() * timestep);
                 f(pos_idx, value_e, value_h)
             }
diff --git a/crates/coremem/src/stim.rs b/crates/coremem/src/stim.rs
index ccb29ce..40ed994 100644
--- a/crates/coremem/src/stim.rs
+++ b/crates/coremem/src/stim.rs
@@ -3,6 +3,7 @@ use crate::cross::vec::Vec3;
 use crate::geom::{Coord as _, HasCrossSection, Index, Meters, Region};
 use coremem_cross::compound::list::{Visit, Visitor};
 use coremem_cross::dim::OffsetDimSlice;
+use coremem_cross::vec::Vec3u;
 use rand;
 
 /// field densities
@@ -70,32 +71,48 @@ impl std::ops::Mul<f32> for FieldMags {
 }
 
 pub trait Stimulus: Sync {
-    // TODO: using floats for time and position is not great.
     /// Return the (E, H) field which should be added PER-SECOND to the provided position/time.
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields;
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
+        let mut dat = OffsetDimSlice::new(loc.into(), Vec3u::new(1, 1, 1), [Fields::default()]);
+        self.eval_into(t_sec, feat_size, 1.0, dat.as_mut());
+        let [fields] = dat.into_inner();
+        fields
+    }
 
     // TODO: could remove the `scale` param if we parameterized this array over some `F: FnMut(Fields)` instead `Fields`.
     // that would also allow easier unzipping in the SpirvSim code.
     /// bulk version of at. evaluate several positions at once and populate the output region.
     fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, into: OffsetDimSlice<&mut [Fields]>) {
         for (idx, out) in into.enumerated() {
-            *out += self.at(t_sec, Index(idx).to_meters(feat_size)) * scale;
+            *out += self.at(t_sec, feat_size, Index::from(idx)) * scale;
         }
     }
 }
 
 /// used as a MapVisitor in order to evaluate each Stimulus in a List at a specific time/place.
-struct StimulusEvaluator {
-    fields: Fields,
-    t_sec: f32,
-    pos: Meters,
-}
-
-impl<S: Stimulus> Visitor<&S> for &mut StimulusEvaluator {
-    fn visit(&mut self, next: &S) {
-        self.fields += next.at(self.t_sec, self.pos);
-    }
-}
+// struct StimulusEvaluator {
+//     fields: Fields,
+//     t_sec: f32,
+//     feat_size: f32,
+//     loc: Index,
+// }
+//
+// impl<S: Stimulus> Visitor<&S> for &mut StimulusEvaluator {
+//     fn visit(&mut self, next: &S) {
+//         self.fields += next.at(self.t_sec, self.feat_size, self.loc);
+//     }
+// }
+//
+// impl<L: Sync> Stimulus for L
+// where
+//     for<'a, 'b> &'a L: Visit<&'b mut StimulusEvaluator>,
+// {
+//     fn at(&self, t_sec: f32, pos: Meters) -> Fields {
+//         let mut ev = StimulusEvaluator { t_sec, pos, fields: Fields::default()};
+//         self.visit(&mut ev);
+//         ev.fields
+//     }
+// }
 
 struct StimEvalInto<'a> {
     into: OffsetDimSlice<&'a mut [Fields]>,
@@ -110,26 +127,10 @@ impl<'a, S: Stimulus> Visitor<&S> for StimEvalInto<'a> {
     }
 }
 
-// impl<L: Sync> Stimulus for L
-// where
-//     for<'a, 'b> &'a L: Visit<&'b mut StimulusEvaluator>,
-// {
-//     fn at(&self, t_sec: f32, pos: Meters) -> Fields {
-//         let mut ev = StimulusEvaluator { t_sec, pos, fields: Fields::default()};
-//         self.visit(&mut ev);
-//         ev.fields
-//     }
-// }
-
 impl<L: Sync> Stimulus for L
 where
     for<'a, 'b> &'a L: Visit<StimEvalInto<'b>>,
 {
-    fn at(&self, _t_sec: f32, _pos: Meters) -> Fields {
-        // TODO: if we replace `_pos` with `feat_size`, `Index`, we can replace this with a generic
-        // impl that calls eval_into one a one-len slice.
-        unimplemented!();
-    }
     fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, into: OffsetDimSlice<&mut [Fields]>) {
         let ev = StimEvalInto { t_sec, feat_size, scale, into };
         self.visit(ev);
@@ -159,8 +160,8 @@ impl DynStimuli {
     }
 }
 impl Stimulus for DynStimuli {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
-        self.0.iter().map(|i| i.at(t_sec, pos))
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
+        self.0.iter().map(|i| i.at(t_sec, feat_size, loc))
             .fold(Fields::default(), core::ops::Add::add)
     }
     fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, mut into: OffsetDimSlice<&mut [Fields]>) {
@@ -178,7 +179,7 @@ impl Stimulus for DynStimuli {
 
 pub struct NoopStimulus;
 impl Stimulus for NoopStimulus {
-    fn at(&self, _t_sec: f32, _pos: Meters) -> Fields {
+    fn at(&self, _t_sec: f32, _feat_size: f32, _loc: Index) -> Fields {
         Fields::default()
     }
 }
@@ -197,7 +198,7 @@ impl UniformStimulus {
 }
 
 impl Stimulus for UniformStimulus {
-    fn at(&self, _t_sec: f32, _pos: Meters) -> Fields {
+    fn at(&self, _t_sec: f32, _feat_size: f32, _loc: Index) -> Fields {
         self.fields
     }
 }
@@ -233,10 +234,10 @@ impl RngStimulus {
 }
 
 impl Stimulus for RngStimulus {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
         Fields {
-            e: self.gen(t_sec, pos, self.e_scale, 0),
-            h: self.gen(t_sec, pos, self.h_scale, 0x7de3)
+            e: self.gen(t_sec, loc.to_meters(feat_size), self.e_scale, 0),
+            h: self.gen(t_sec, loc.to_meters(feat_size), self.h_scale, 0x7de3)
         }
     }
 }
@@ -257,9 +258,9 @@ impl<R, S> RegionGated<R, S> {
 }
 
 impl<R: Region + Sync, S: Stimulus + Sync> Stimulus for RegionGated<R, S> {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
-        if self.region.contains(pos) {
-            self.stim.at(t_sec, pos)
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
+        if self.region.contains(loc.to_meters(feat_size)) {
+            self.stim.at(t_sec, feat_size, loc)
         } else {
             Fields::default()
         }
@@ -281,7 +282,8 @@ impl<R, T> CurlStimulus<R, T> {
 }
 
 impl<R: Region + HasCrossSection, T: TimeVarying + Sync> Stimulus for CurlStimulus<R, T> {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
+        let pos = loc.to_meters(feat_size);
         if self.region.contains(pos) {
             let FieldMags { e, h } = self.stim.at(t_sec);
             let rotational = self.region.cross_section_normal(pos);
@@ -395,9 +397,9 @@ impl<A: TimeVarying> TimeVarying for Exp<A> {
 }
 
 impl<A: Stimulus> Stimulus for Exp<A> {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
         let scale = (t_sec * -self.tau).exp();
-        let inner = self.amp.at(t_sec, pos);
+        let inner = self.amp.at(t_sec, feat_size, loc);
         Fields {
             e: inner.e * scale,
             h: inner.h * scale,
@@ -433,9 +435,9 @@ impl<T: TimeVarying> TimeVarying for Gated<T> {
 }
 
 impl<T: Stimulus> Stimulus for Gated<T> {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
         if (self.start..self.end).contains(&t_sec) {
-            self.inner.at(t_sec, pos)
+            self.inner.at(t_sec, feat_size, loc)
         } else {
             Default::default()
         }
@@ -466,8 +468,8 @@ impl<T: TimeVarying> TimeVarying for Shifted<T> {
 }
 
 impl<T: Stimulus> Stimulus for Shifted<T> {
-    fn at(&self, t_sec: f32, pos: Meters) -> Fields {
-        self.inner.at(t_sec - self.start_at, pos)
+    fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
+        self.inner.at(t_sec - self.start_at, feat_size, loc)
     }
     fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, into: OffsetDimSlice<&mut [Fields]>) {
         self.inner.eval_into(t_sec - self.start_at, feat_size, scale, into)
@@ -529,7 +531,7 @@ mod test {
         };
         // stim acts as e: 1.0, h: 0.0
         let stim = CurlStimulus::new(region, 1.0);
-        assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
+        assert_eq!(stim.at(0.0, 1.0, Index::new(0, 0, 0)), Fields {
             e: Vec3::new(1.0, 0.0, 0.0),
             h: Vec3::new(0.0, 0.0, 0.0),
         });
@@ -544,7 +546,7 @@ mod test {
         };
         // stim acts as e: -3.0, h: 0.0
         let stim = CurlStimulus::new(region, -3.0);
-        assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
+        assert_eq!(stim.at(0.0, 1.0, Index::new(0, 0, 0)), Fields {
             e: Vec3::new(-3.0, 0.0, 0.0),
             h: Vec3::new(0.0, 0.0, 0.0),
         });
@@ -557,7 +559,7 @@ mod test {
         };
         // stim acts as e: 2.0, h: 0.0
         let stim = CurlStimulus::new(region, 2.0);
-        let Fields { e, h: _ } = stim.at(0.0, Meters::new(0.0, 0.0, 0.0));
+        let Fields { e, h: _ } = stim.at(0.0, 1.0, Index::new(0, 0, 0));
         assert!(e.distance(
                 Vec3::new(0.0, -1.0, 1.0).with_mag(2.0).unwrap()
             ) < 1e-6