diff --git a/examples/minimal_torus.rs b/examples/minimal_torus.rs
index d01e40f..58e1f88 100644
--- a/examples/minimal_torus.rs
+++ b/examples/minimal_torus.rs
@@ -30,14 +30,14 @@ fn main() {
     driver.set_steps_per_frame(400);
     driver.set_steps_per_stim(1);
     let base = "minimal_torus-6";
-    let ferro_region = Torus::new_xy(Meters((half_width, half_width, half_depth).into()), ferro_major, ferro_minor);
-    let drive_region = Torus::new_xz(Meters((half_width - ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
-    let sense_region = Torus::new_xz(Meters((half_width + ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
+    let ferro_region = Torus::new_xy(Meters::new(half_width, half_width, half_depth), ferro_major, ferro_minor);
+    let drive_region = Torus::new_xz(Meters::new(half_width - ferro_major, half_width, half_depth), wire_major, wire_minor);
+    let sense_region = Torus::new_xz(Meters::new(half_width + ferro_major, half_width, half_depth), wire_major, wire_minor);
 
     let boundary_xy = half_width - ferro_major - ferro_minor - buffer;
     let boundary_z = half_depth - wire_major - wire_minor - buffer;
-    let boundary_lower = Meters((boundary_xy, boundary_xy, boundary_z).into());
-    let boundary_upper = Meters((width - boundary_xy, width - boundary_xy, depth - boundary_z).into());
+    let boundary_lower = Meters::new(boundary_xy, boundary_xy, boundary_z);
+    let boundary_upper = Meters::new(width - boundary_xy, width - boundary_xy, depth - boundary_z);
     println!("boundary: {}um; {}um", m_to_um(boundary_xy), m_to_um(boundary_z));
     let boundary_region = InvertedRegion::new(Cube::new(boundary_lower, boundary_upper));
 
@@ -52,7 +52,7 @@ fn main() {
     driver.fill_region(&drive_region, mat::db::conductor(conductivity).into());
     driver.fill_region(&sense_region, mat::db::conductor(conductivity).into());
     // driver.add_pml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
-    driver.add_classical_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
+    driver.add_classical_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z));
 
     // J=\sigma E
     // dJ/dt = \sigma dE/dT
diff --git a/examples/toroid25d.rs b/examples/toroid25d.rs
index f22821c..34ffe1f 100644
--- a/examples/toroid25d.rs
+++ b/examples/toroid25d.rs
@@ -64,31 +64,31 @@ fn main() {
         driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current)));
         //driver.add_measurement(meas::Current(conductor_region.clone()));
         driver.add_measurement(meas::MagnetizationAt(
-            Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into())
+            Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
         ));
         driver.add_measurement(meas::MagneticFluxAt(
-            Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into())
+            Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
         ));
         driver.add_measurement(meas::MagneticStrengthAt(
-            Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into())
+            Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
         ));
         driver.add_measurement(meas::MagnetizationAt(
-            Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into())
+            Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
         ));
         driver.add_measurement(meas::MagneticFluxAt(
-            Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into())
+            Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
         ));
         driver.add_measurement(meas::MagneticStrengthAt(
-            Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into())
+            Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
         ));
         driver.add_measurement(meas::MagnetizationAt(
-            Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into())
+            Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
         ));
         driver.add_measurement(meas::MagneticFluxAt(
-            Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into())
+            Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
         ));
         driver.add_measurement(meas::MagneticStrengthAt(
-            Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into())
+            Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
         ));
 
         let center = Vec2::new(half_width as Flt, half_width as Flt);
diff --git a/examples/wrapped_torus.rs b/examples/wrapped_torus.rs
index 2ca0187..dfe562d 100644
--- a/examples/wrapped_torus.rs
+++ b/examples/wrapped_torus.rs
@@ -17,7 +17,8 @@ fn main() {
     let peak_current = 5e6;
     let current_duration = 0.1e-9; // half-wavelength of the sine wave
     let current_break = 0.1e-9; // time between 'set' pulse and 'clear' pulse
-    let conductivity = 5.0e7;
+    let drive_conductivity = 5.0;
+    let sense_conductivity = 5.0;
 
     let from_m = |m: Flt| (m/feat_size).round() as u32;
     let m_to_um = |m: Flt| (m * 1e6).round() as u32;
@@ -34,28 +35,28 @@ fn main() {
     let mut driver: Driver<mat::GenericMaterial> = Driver::new(size_px, feat_size);
     driver.set_steps_per_frame(500);
     // driver.set_steps_per_stim(10);
-    let base = "wrapped_torus-19";
+    let base = "wrapped_torus-21-high-input-resistance";
 
-    let ferro1_region = Torus::new_xy(Meters((q1_width, half_height, half_depth).into()), ferro_major, ferro_minor);
-    let drive1_region = Torus::new_xz(Meters((q1_width - ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
-    let sense1_region = Torus::new_xz(Meters((q1_width + ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
+    let ferro1_region = Torus::new_xy(Meters::new(q1_width, half_height, half_depth), ferro_major, ferro_minor);
+    let drive1_region = Torus::new_xz(Meters::new(q1_width - ferro_major, half_height, half_depth), wire_major, wire_minor);
+    let sense1_region = Torus::new_xz(Meters::new(q1_width + ferro_major, half_height, half_depth), wire_major, wire_minor);
 
     driver.fill_region(&ferro1_region, mat::db::linear_iron().into());
-    driver.fill_region(&drive1_region, mat::db::conductor(conductivity).into());
-    driver.fill_region(&sense1_region, mat::db::conductor(conductivity).into());
+    driver.fill_region(&drive1_region, mat::db::conductor(drive_conductivity).into());
+    driver.fill_region(&sense1_region, mat::db::conductor(sense_conductivity).into());
 
-    let ferro2_region = Torus::new_xy(Meters((q3_width, half_height, half_depth).into()), ferro_major, ferro_minor);
-    let drive2_region = Torus::new_xz(Meters((q3_width - ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
-    let sense2_region = Torus::new_xz(Meters((q3_width + ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
+    let ferro2_region = Torus::new_xy(Meters::new(q3_width, half_height, half_depth), ferro_major, ferro_minor);
+    let drive2_region = Torus::new_xz(Meters::new(q3_width - ferro_major, half_height, half_depth), wire_major, wire_minor);
+    let sense2_region = Torus::new_xz(Meters::new(q3_width + ferro_major, half_height, half_depth), wire_major, wire_minor);
 
     driver.fill_region(&ferro2_region, mat::db::ferroxcube_3r1().into());
-    driver.fill_region(&drive2_region, mat::db::conductor(conductivity).into());
-    driver.fill_region(&sense2_region, mat::db::conductor(conductivity).into());
+    driver.fill_region(&drive2_region, mat::db::conductor(drive_conductivity).into());
+    driver.fill_region(&sense2_region, mat::db::conductor(sense_conductivity).into());
 
     let boundary_xy = q1_width - ferro_major - ferro_minor - buffer;
     let boundary_z = half_depth - wire_major - wire_minor - buffer;
     println!("boundary: {}um; {}um", m_to_um(boundary_xy), m_to_um(boundary_z));
-    driver.add_pml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
+    driver.add_pml_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z));
 
     // J=\sigma E
     // dJ/dt = \sigma dE/dT
@@ -63,7 +64,7 @@ fn main() {
     // dE/dt = dI/dt / (A*\sigma)
     // if I = k*sin(w t) then dE/dt = k*w sin(w t) / (A*\sigma)
     // i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength)
-    let peak_stim = peak_current/current_duration / (drive1_region.cross_section() * conductivity);
+    let peak_stim = peak_current/current_duration / (drive1_region.cross_section() * drive_conductivity);
     let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0)
         .half_cycle();
 
diff --git a/src/bin/pml_tuning.rs b/src/bin/pml_tuning.rs
index 42c238e..8bad120 100644
--- a/src/bin/pml_tuning.rs
+++ b/src/bin/pml_tuning.rs
@@ -1,13 +1,15 @@
 use coremem::*;
 use coremem::geom::*;
+use coremem::real::Real as _;
 use coremem::stim::AbstractStimulus;
-    use rand::rngs::StdRng;
-    use rand::{Rng as _, SeedableRng as _};
+use rand::rngs::StdRng;
+use rand::{Rng as _, SeedableRng as _};
 
 fn energy<R: Region>(s: &dyn GenericSim, reg: &R) -> Flt {
-    0.5 * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
+    let e = Real::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
         cell.e().mag_sq()
-    })
+    });
+    Flt::from_primitive(e)
 }
 
 fn energy_now_and_then<R: Region>(state: &mut StaticSim, reg: &R, frames: u32) -> (Flt, Flt) {
@@ -31,7 +33,7 @@ impl<F: Fn(Index) -> (Vec3, Flt) + Sync> AbstractStimulus for PmlStim<F> {
         let (e, hz) = (self.f)(pos.to_index(self.feat_size));
         let sig_angle = angle*hz;
         let sig = sig_angle.sin();
-        e*gate*sig
+        e*Real::from_primitive(gate*sig)
     }
 }
 
@@ -104,8 +106,8 @@ fn state_for_pml(size: Index, boundary: Index, feat_size: Flt, sc_coeff: Flt, co
     let mut state = StaticSim::new(size, feat_size);
     let timestep = state.timestep();
     state.fill_boundary_using(boundary, |boundary_ness| {
-        let b = boundary_ness.elem_pow(pow);
-        let coord_stretch = b * (sc_coeff / timestep);
+        let b = boundary_ness.elem_pow(Real::from_primitive(pow));
+        let coord_stretch = b * Real::from_primitive(sc_coeff / timestep);
         let conductivity = Vec3::unit() * (b.mag() * cond_coeff / timestep);
         unimplemented!()
         // Static {
diff --git a/src/driver.rs b/src/driver.rs
index 25564ec..b1dbfed 100644
--- a/src/driver.rs
+++ b/src/driver.rs
@@ -1,7 +1,8 @@
-use crate::{flt::Flt, mat};
+use crate::{flt, flt::Flt, mat};
 use crate::geom::{Coord, Index, Meters, Region, Vec3, Vec3u};
 use crate::mat::{GenericMaterial, Material, Pml};
 use crate::meas::{self, AbstractMeasurement};
+use crate::real::Real as _;
 use crate::render::{self, MultiRenderer, Renderer};
 use crate::sim::{GenericSim, SimState};
 use crate::stim::AbstractStimulus;
@@ -193,8 +194,8 @@ impl<M: Material + From<mat::Pml>> Driver<M> {
     pub fn add_pml_boundary<C: Coord>(&mut self, thickness: C) {
         let timestep = self.state.timestep();
         self.state.fill_boundary_using(thickness, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let conductivity = b * (0.5 / timestep);
+            let b = boundary_ness.elem_pow(flt::Real::three());
+            let conductivity = b * (flt::Real::half() / timestep);
             Pml::new(conductivity).into()
         });
     }
@@ -204,11 +205,11 @@ impl<M: Material + From<mat::Conductor>> Driver<M> {
     pub fn add_classical_boundary<C: Coord>(&mut self, thickness: C) {
         let timestep = self.state.timestep();
         self.state.fill_boundary_using(thickness, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let cond = b * (0.5 / timestep);
+            let b = boundary_ness.elem_pow(flt::Real::three());
+            let cond = b * (flt::Real::half() / timestep);
 
             let iso_cond = cond.x() + cond.y() + cond.z();
-            let iso_conductor = mat::db::conductor(iso_cond);
+            let iso_conductor = mat::db::conductor(Flt::from_primitive(iso_cond));
             iso_conductor.into()
         });
     }
@@ -223,7 +224,7 @@ struct StimuliAdapter {
 
 impl AbstractStimulus for StimuliAdapter {
     fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
-        self.stim.at(t_sec, pos) * (self.frame_interval as Flt)
+        self.stim.at(t_sec, pos) * flt::Real::from_f64(self.frame_interval as f64)
     }
 }
 
diff --git a/src/geom/region.rs b/src/geom/region.rs
index bbf2e10..2886381 100644
--- a/src/geom/region.rs
+++ b/src/geom/region.rs
@@ -1,5 +1,6 @@
 use crate::flt::{Flt, Real};
 use crate::geom::{Meters, Vec2, Vec3};
+use crate::real::Real as _;
 
 use dyn_clone::{self, DynClone};
 use serde::{Serialize, Deserialize};
@@ -98,7 +99,7 @@ impl Region for Torus {
             // and they all give the same answer
             Vec3::new(self.major_rad.into(), 0.0, 0.0)
         } else {
-            p_on_plane.with_mag(self.major_rad.into_inner())
+            p_on_plane.with_mag(self.major_rad)
         };
         let distance_to_circle = (rel_p - q).mag();
         distance_to_circle < self.minor_rad.into_inner()
@@ -140,22 +141,22 @@ impl Cube {
         Self { lower, upper }
     }
     pub fn x_range(&self) -> Range<Flt> {
-        self.lower.x()..self.upper.x()
+        Flt::from_primitive(self.lower.x())..Flt::from_primitive(self.upper.x())
     }
     pub fn y_range(&self) -> Range<Flt> {
-        self.lower.y()..self.upper.y()
+        Flt::from_primitive(self.lower.y())..Flt::from_primitive(self.upper.y())
     }
     pub fn z_range(&self) -> Range<Flt> {
-        self.lower.z()..self.upper.z()
+        Flt::from_primitive(self.lower.z())..Flt::from_primitive(self.upper.z())
     }
 }
 
 #[typetag::serde]
 impl Region for Cube {
     fn contains(&self, p: Meters) -> bool {
-        self.x_range().contains(&p.x()) &&
-            self.y_range().contains(&p.y()) &&
-            self.z_range().contains(&p.z())
+        self.x_range().contains(&Flt::from_primitive(p.x())) &&
+            self.y_range().contains(&Flt::from_primitive(p.y())) &&
+            self.z_range().contains(&Flt::from_primitive(p.z()))
     }
 }
 
@@ -214,103 +215,103 @@ mod test {
     use super::*;
     #[test]
     fn torus_as_sphere() {
-        let tor = Torus::new_xy(Meters((0.0, 0.0, 0.0).into()), 0.0, 10.0);
-        assert!(tor.contains(Meters((0.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((5.0, 0.0, 0.0).into())));
+        let tor = Torus::new_xy(Meters::new(0.0, 0.0, 0.0), 0.0, 10.0);
+        assert!(tor.contains(Meters::new(0.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(5.0, 0.0, 0.0)));
 
-        assert!(tor.contains(Meters((9.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((-9.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((0.0, 9.0, 0.0).into())));
-        assert!(tor.contains(Meters((0.0, -9.0, 0.0).into())));
-        assert!(tor.contains(Meters((0.0, 0.0, 9.0).into())));
-        assert!(tor.contains(Meters((0.0, 0.0, -9.0).into())));
+        assert!(tor.contains(Meters::new(9.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(-9.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(0.0, 9.0, 0.0)));
+        assert!(tor.contains(Meters::new(0.0, -9.0, 0.0)));
+        assert!(tor.contains(Meters::new(0.0, 0.0, 9.0)));
+        assert!(tor.contains(Meters::new(0.0, 0.0, -9.0)));
 
-        assert!(tor.contains(Meters((5.0, 5.0, 5.0).into())));
-        assert!(tor.contains(Meters((-5.0, -5.0, -5.0).into())));
+        assert!(tor.contains(Meters::new(5.0, 5.0, 5.0)));
+        assert!(tor.contains(Meters::new(-5.0, -5.0, -5.0)));
 
         // boundary
-        assert!(!tor.contains(Meters((10.0, 0.0, 0.0).into())));
-        assert!(!tor.contains(Meters((-10.0, 0.0, 0.0).into())));
-        assert!(!tor.contains(Meters((0.0, 10.0, 0.0).into())));
-        assert!(!tor.contains(Meters((0.0, -10.0, 0.0).into())));
-        assert!(!tor.contains(Meters((0.0, 0.0, 10.0).into())));
-        assert!(!tor.contains(Meters((0.0, 0.0, -10.0).into())));
+        assert!(!tor.contains(Meters::new(10.0, 0.0, 0.0)));
+        assert!(!tor.contains(Meters::new(-10.0, 0.0, 0.0)));
+        assert!(!tor.contains(Meters::new(0.0, 10.0, 0.0)));
+        assert!(!tor.contains(Meters::new(0.0, -10.0, 0.0)));
+        assert!(!tor.contains(Meters::new(0.0, 0.0, 10.0)));
+        assert!(!tor.contains(Meters::new(0.0, 0.0, -10.0)));
 
-        assert!(!tor.contains(Meters((11.0, 0.0, 0.0).into())));
-        assert!(!tor.contains(Meters((9.0, 9.0, 9.0).into())));
-        assert!(!tor.contains(Meters((-9.0, -9.0, -9.0).into())));
-        assert!(!tor.contains(Meters((-9.0, 9.0, 0.0).into())));
+        assert!(!tor.contains(Meters::new(11.0, 0.0, 0.0)));
+        assert!(!tor.contains(Meters::new(9.0, 9.0, 9.0)));
+        assert!(!tor.contains(Meters::new(-9.0, -9.0, -9.0)));
+        assert!(!tor.contains(Meters::new(-9.0, 9.0, 0.0)));
     }
 
     #[test]
     fn torus_xy() {
-        let tor = Torus::new_xy(Meters((0.0, 0.0, 0.0).into()), 10.0, 2.0);
-        assert!(!tor.contains(Meters((0.0, 0.0, 0.0).into())));
+        let tor = Torus::new_xy(Meters::new(0.0, 0.0, 0.0), 10.0, 2.0);
+        assert!(!tor.contains(Meters::new(0.0, 0.0, 0.0)));
 
-        assert!(tor.contains(Meters((10.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((-10.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((0.0, 10.0, 0.0).into())));
-        assert!(tor.contains(Meters((0.0, -10.0, 0.0).into())));
-        assert!(!tor.contains(Meters((0.0, 0.0, 10.0).into())));
-        assert!(!tor.contains(Meters((0.0, 0.0, -10.0).into())));
+        assert!(tor.contains(Meters::new(10.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(-10.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(0.0, 10.0, 0.0)));
+        assert!(tor.contains(Meters::new(0.0, -10.0, 0.0)));
+        assert!(!tor.contains(Meters::new(0.0, 0.0, 10.0)));
+        assert!(!tor.contains(Meters::new(0.0, 0.0, -10.0)));
 
-        assert!(tor.contains(Meters((11.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((-11.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((11.0, 2.0, 0.0).into())));
-        assert!(tor.contains(Meters((11.0, 0.0, 1.0).into())));
-        assert!(tor.contains(Meters((-11.0, 0.0, 1.0).into())));
-        assert!(tor.contains(Meters((-7.0, 7.0, 1.0).into())));
-        assert!(tor.contains(Meters((-7.0, -7.0, -1.0).into())));
+        assert!(tor.contains(Meters::new(11.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(-11.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(11.0, 2.0, 0.0)));
+        assert!(tor.contains(Meters::new(11.0, 0.0, 1.0)));
+        assert!(tor.contains(Meters::new(-11.0, 0.0, 1.0)));
+        assert!(tor.contains(Meters::new(-7.0, 7.0, 1.0)));
+        assert!(tor.contains(Meters::new(-7.0, -7.0, -1.0)));
 
         // boundary
-        assert!(!tor.contains(Meters((12.0, 0.0, 0.0).into())));
-        assert!(!tor.contains(Meters((10.0, 0.0, 2.0).into())));
-        assert!(!tor.contains(Meters((-10.0, 0.0, 2.0).into())));
+        assert!(!tor.contains(Meters::new(12.0, 0.0, 0.0)));
+        assert!(!tor.contains(Meters::new(10.0, 0.0, 2.0)));
+        assert!(!tor.contains(Meters::new(-10.0, 0.0, 2.0)));
     }
 
     #[test]
     fn torus_xz() {
-        let tor = Torus::new_xz(Meters((0.0, 0.0, 0.0).into()), 10.0, 2.0);
-        assert!(!tor.contains(Meters((0.0, 0.0, 0.0).into())));
+        let tor = Torus::new_xz(Meters::new(0.0, 0.0, 0.0), 10.0, 2.0);
+        assert!(!tor.contains(Meters::new(0.0, 0.0, 0.0)));
 
-        assert!(tor.contains(Meters((10.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((-10.0, 0.0, 0.0).into())));
-        assert!(!tor.contains(Meters((0.0, 10.0, 0.0).into())));
-        assert!(!tor.contains(Meters((0.0, -10.0, 0.0).into())));
-        assert!(tor.contains(Meters((0.0, 0.0, 10.0).into())));
-        assert!(tor.contains(Meters((0.0, 0.0, -10.0).into())));
+        assert!(tor.contains(Meters::new(10.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(-10.0, 0.0, 0.0)));
+        assert!(!tor.contains(Meters::new(0.0, 10.0, 0.0)));
+        assert!(!tor.contains(Meters::new(0.0, -10.0, 0.0)));
+        assert!(tor.contains(Meters::new(0.0, 0.0, 10.0)));
+        assert!(tor.contains(Meters::new(0.0, 0.0, -10.0)));
 
-        assert!(tor.contains(Meters((11.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((-11.0, 0.0, 0.0).into())));
-        assert!(tor.contains(Meters((11.0, 1.0, 0.0).into())));
-        assert!(tor.contains(Meters((11.0, 0.0, 2.0).into())));
-        assert!(tor.contains(Meters((-11.0, 1.0, 0.0).into())));
-        assert!(tor.contains(Meters((-7.0, 1.0, 7.0).into())));
-        assert!(tor.contains(Meters((-7.0, -1.0, -7.0).into())));
+        assert!(tor.contains(Meters::new(11.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(-11.0, 0.0, 0.0)));
+        assert!(tor.contains(Meters::new(11.0, 1.0, 0.0)));
+        assert!(tor.contains(Meters::new(11.0, 0.0, 2.0)));
+        assert!(tor.contains(Meters::new(-11.0, 1.0, 0.0)));
+        assert!(tor.contains(Meters::new(-7.0, 1.0, 7.0)));
+        assert!(tor.contains(Meters::new(-7.0, -1.0, -7.0)));
 
         // boundary
-        assert!(!tor.contains(Meters((12.0, 0.0, 0.0).into())));
-        assert!(!tor.contains(Meters((10.0, 2.0, 0.0).into())));
-        assert!(!tor.contains(Meters((-10.0, 2.0, 0.0).into())));
+        assert!(!tor.contains(Meters::new(12.0, 0.0, 0.0)));
+        assert!(!tor.contains(Meters::new(10.0, 2.0, 0.0)));
+        assert!(!tor.contains(Meters::new(-10.0, 2.0, 0.0)));
     }
 
     #[test]
     fn torus_with_center() {
-        let tor = Torus::new_xy(Meters((1.0, 2.0, 3.0).into()), 0.0, 10.0);
-        assert!(tor.contains(Meters((1.0, 2.0, 3.0).into())));
-        assert!(tor.contains(Meters((1.0, 2.0, 12.0).into())));
-        assert!(!tor.contains(Meters((1.0, 2.0, 13.0).into())));
+        let tor = Torus::new_xy(Meters::new(1.0, 2.0, 3.0), 0.0, 10.0);
+        assert!(tor.contains(Meters::new(1.0, 2.0, 3.0)));
+        assert!(tor.contains(Meters::new(1.0, 2.0, 12.0)));
+        assert!(!tor.contains(Meters::new(1.0, 2.0, 13.0)));
 
-        let tor = Torus::new_xy(Meters((1.0, 2.0, 3.0).into()), 10.0, 3.0);
-        assert!(!tor.contains(Meters((1.0, 2.0, 3.0).into())));
-        assert!(tor.contains(Meters((1.0, 14.0, 3.0).into())));
-        assert!(!tor.contains(Meters((1.0, 15.0, 3.0).into())));
+        let tor = Torus::new_xy(Meters::new(1.0, 2.0, 3.0), 10.0, 3.0);
+        assert!(!tor.contains(Meters::new(1.0, 2.0, 3.0)));
+        assert!(tor.contains(Meters::new(1.0, 14.0, 3.0)));
+        assert!(!tor.contains(Meters::new(1.0, 15.0, 3.0)));
 
-        let tor = Torus::new_xz(Meters((1.0, 2.0, 3.0).into()), 5.0, 5.0);
-        assert!(!tor.contains(Meters((1.0, 2.0, 3.0).into())));
-        assert!(tor.contains(Meters((2.0, 2.0, 3.0).into())));
-        assert!(tor.contains(Meters((1.0, 2.0, 12.0).into())));
-        assert!(!tor.contains(Meters((1.0, 2.0, 13.0).into())));
+        let tor = Torus::new_xz(Meters::new(1.0, 2.0, 3.0), 5.0, 5.0);
+        assert!(!tor.contains(Meters::new(1.0, 2.0, 3.0)));
+        assert!(tor.contains(Meters::new(2.0, 2.0, 3.0)));
+        assert!(tor.contains(Meters::new(1.0, 2.0, 12.0)));
+        assert!(!tor.contains(Meters::new(1.0, 2.0, 13.0)));
     }
 }
 
diff --git a/src/geom/units.rs b/src/geom/units.rs
index 503cd96..9966549 100644
--- a/src/geom/units.rs
+++ b/src/geom/units.rs
@@ -1,4 +1,5 @@
-use crate::flt::Flt;
+use crate::flt::{Flt, Real};
+use crate::real::{Real as _, ToFloat};
 use serde::{Serialize, Deserialize};
 use super::{Vec3, Vec3u};
 use std::fmt::{self, Display};
@@ -15,12 +16,18 @@ pub trait Coord: Copy + Clone {
 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
 pub struct Meters(pub Vec3);
 
+impl Meters {
+    pub fn new<R: ToFloat>(x: R, y: R, z: R) -> Self {
+        Self(Vec3::new(x, y, z).into())
+    }
+}
+
 impl Coord for Meters {
     fn to_meters(&self, _feature_size: Flt) -> Meters {
         *self
     }
     fn to_index(&self, feature_size: Flt) -> Index {
-        Index((self.0 / feature_size).round().into())
+        Index((self.0 / Real::from_f64(feature_size)).round().into())
     }
     fn from_meters(other: Meters, _feature_size: Flt) -> Self {
         other
@@ -74,7 +81,7 @@ impl Index {
 
 impl Coord for Index {
     fn to_meters(&self, feature_size: Flt) -> Meters {
-        Meters(Vec3::from(self.0) * feature_size)
+        Meters(Vec3::from(self.0) * Real::from_f64(feature_size))
     }
     fn to_index(&self, _feature_size: Flt) -> Index {
         *self
diff --git a/src/geom/vec.rs b/src/geom/vec.rs
index 658d474..6d3ddaa 100644
--- a/src/geom/vec.rs
+++ b/src/geom/vec.rs
@@ -1,5 +1,5 @@
 use crate::flt::{self, Flt};
-use crate::real::Real;
+use crate::real::{Real, ToFloat};
 use super::Vec3u;
 
 use serde::{Serialize, Deserialize};
@@ -63,11 +63,11 @@ impl<R> Into<(R, R)> for Vec2<R> {
     }
 }
 
-impl<R> Vec2<R> {
-    pub fn new<R2: Into<R>>(x: R2, y: R2) -> Self {
+impl<R: Real> Vec2<R> {
+    pub fn new<R2: ToFloat>(x: R2, y: R2) -> Self {
         Self {
-            x: x.into(),
-            y: y.into(),
+            x: R::from_primitive(x),
+            y: R::from_primitive(y),
         }
     }
 }
@@ -124,71 +124,113 @@ impl<R: Real> Vec2<R> {
 }
 
 #[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
-pub struct Vec3 {
-    pub(crate) x: flt::Real,
-    pub(crate) y: flt::Real,
-    pub(crate) z: flt::Real,
+pub struct Vec3<R=flt::Real> {
+    pub(crate) x: R,
+    pub(crate) y: R,
+    pub(crate) z: R,
 }
 
-impl Vec3 {
-    pub fn new(x: Flt, y: Flt, z: Flt) -> Self {
-        Self {
-            x: x.into(),
-            y: y.into(),
-            z: z.into()
-        }
-    }
-    pub fn to_tuple(self) -> (Flt, Flt, Flt) {
+impl<R> Vec3<R> {
+    pub fn to_tuple(self) -> (R, R, R) {
         self.into()
     }
-    pub fn unit() -> Self {
-        Self::new(1.0, 1.0, 1.0)
+}
+
+impl<R: Real> Vec3<R> {
+    pub fn new<R2: ToFloat>(x: R2, y: R2, z: R2) -> Self {
+        Self {
+            x: R::from_primitive(x),
+            y: R::from_primitive(y),
+            z: R::from_primitive(z),
+        }
     }
-    pub fn unit_x() -> Self {
-        Self::new(1.0, 0.0, 0.0)
+    pub fn new_x<R2: ToFloat>(x: R2) -> Self {
+        Self {
+            x: R::from_primitive(x),
+            y: R::zero(),
+            z: R::zero(),
+        }
     }
-    pub fn unit_y() -> Self {
-        Self::new(0.0, 1.0, 0.0)
+    pub fn new_y<R2: ToFloat>(y: R2) -> Self {
+        Self {
+            x: R::zero(),
+            y: R::from_primitive(y),
+            z: R::zero(),
+        }
     }
-    pub fn unit_z() -> Self {
-        Self::new(0.0, 0.0, 1.0)
+    pub fn new_z<R2: ToFloat>(z: R2) -> Self {
+        Self {
+            x: R::zero(),
+            y: R::zero(),
+            z: R::from_primitive(z),
+        }
     }
-    pub fn zero() -> Self {
-        Self::new(0.0, 0.0, 0.0)
-    }
-    pub fn x(&self) -> Flt {
-        self.x.into()
-    }
-    pub fn y(&self) -> Flt {
-        self.y.into()
-    }
-    pub fn z(&self) -> Flt {
-        self.z.into()
-    }
-    pub fn set_x(&mut self, x: Flt) {
-        self.x = flt::Real::from_inner(x);
-    }
-    pub fn set_y(&mut self, y: Flt) {
-        self.y = flt::Real::from_inner(y);
-    }
-    pub fn set_z(&mut self, z: Flt) {
-        self.z = flt::Real::from_inner(z);
-    }
-    pub fn xy(&self) -> Vec2 {
-        Vec2::new(self.x(), self.y())
-    }
-    pub fn mag(&self) -> Flt {
-        self.mag_sq().sqrt()
-    }
-    pub fn mag_sq(&self) -> Flt {
-        let Vec3 { x, y, z } = *self;
-        (x*x + y*y + z*z).into_inner()
-    }
-    pub fn component_sum(&self) -> Flt {
-        (self.x + self.y + self.z).into()
+    pub fn uniform<R2: ToFloat>(u: R2) -> Self {
+        let u = R::from_primitive(u);
+        Self {
+            x: u,
+            y: u,
+            z: u,
+        }
     }
 
-    pub fn dot(&self, other: Self) -> Flt {
+    pub fn cast<R2: Real>(&self) -> Vec3<R2> {
+        Vec3::new(
+            R2::from_primitive(self.x),
+            R2::from_primitive(self.y),
+            R2::from_primitive(self.z)
+        )
+    }
+
+    pub fn unit() -> Self {
+        Self::new(R::one(), R::one(), R::one())
+    }
+    pub fn unit_x() -> Self {
+        Self::new(R::one(), R::zero(), R::zero())
+    }
+    pub fn unit_y() -> Self {
+        Self::new(R::zero(), R::one(), R::zero())
+    }
+    pub fn unit_z() -> Self {
+        Self::new(R::zero(), R::zero(), R::one())
+    }
+    pub fn zero() -> Self {
+        Self::new(R::zero(), R::zero(), R::zero())
+    }
+
+    pub fn set_x<R2: ToFloat>(&mut self, x: R2) {
+        self.x = R::from_primitive(x);
+    }
+    pub fn set_y<R2: ToFloat>(&mut self, y: R2) {
+        self.y = R::from_primitive(y);
+    }
+    pub fn set_z<R2: ToFloat>(&mut self, z: R2) {
+        self.z = R::from_primitive(z);
+    }
+    pub fn x(&self) -> R {
+        self.x
+    }
+    pub fn y(&self) -> R {
+        self.y
+    }
+    pub fn z(&self) -> R {
+        self.z
+    }
+    pub fn xy(&self) -> Vec2<R> {
+        Vec2::new(self.x(), self.y())
+    }
+    pub fn mag(&self) -> R {
+        self.mag_sq().sqrt()
+    }
+    pub fn mag_sq(&self) -> R {
+        let Vec3 { x, y, z } = *self;
+        x*x + y*y + z*z
+    }
+    pub fn component_sum(&self) -> R {
+        self.x + self.y + self.z
+    }
+
+    pub fn dot(&self, other: Self) -> R {
         self.elem_mul(other).component_sum()
     }
 
@@ -211,7 +253,7 @@ impl Vec3 {
         }
     }
 
-    pub fn elem_pow(&self, pow: Flt) -> Self {
+    pub fn elem_pow(&self, pow: R) -> Self {
         Self::new(
             self.x().powf(pow),
             self.y().powf(pow),
@@ -230,61 +272,60 @@ impl Vec3 {
 
     /// Return (x^-1, y^-1, z^-1)
     pub fn elem_inv(&self) -> Self {
-        Self::new(1.0 / self.x(), 1.0 / self.y(), 1.0 / self.z())
+        let one = R::one();
+        Self::new(one / self.x(), one / self.y(), one / self.z())
     }
 
     pub fn exp(&self) -> Self {
         Self::new(self.x().exp(), self.y().exp(), self.z().exp())
     }
 
-    pub fn with_mag(&self, new_mag: Flt) -> Vec3 {
-        if new_mag == 0.0 {
+    pub fn with_mag(&self, new_mag: R) -> Self {
+        if new_mag.is_zero() {
             // avoid div-by-zero if self.mag() == 0 and new_mag == 0
             Self::zero()
         } else {
-            let scale = flt::Real::from_inner(new_mag) / self.mag();
-            self.clone() * scale.into_inner()
+            let scale = new_mag / self.mag();
+            *self * scale
         }
     }
 
     /// Normalizes to 1.0 length, unless it's the zero vector, in which case length remains zero.
-    pub fn norm(&self) -> Vec3 {
+    pub fn norm(&self) -> Self {
         if *self == Self::zero() {
             *self
         } else {
-            self.with_mag(1.0)
+            self.with_mag(R::one())
         }
     }
-    pub fn round(&self) -> Vec3 {
+    pub fn round(&self) -> Self {
         Self::new(self.x().round(), self.y().round(), self.z().round())
     }
 }
 
-impl Into<(Flt, Flt, Flt)> for Vec3 {
-    fn into(self) -> (Flt, Flt, Flt) {
-        (self.x(), self.y(), self.z())
+impl<R> Into<(R, R, R)> for Vec3<R> {
+    fn into(self) -> (R, R, R) {
+        (self.x, self.y, self.z)
     }
 }
 
-impl From<(Flt, Flt, Flt)> for Vec3 {
-    fn from((x, y, z): (Flt, Flt, Flt)) -> Self {
-        Self::new(x, y, z)
-    }
-}
-
-impl From<(flt::Real, flt::Real, flt::Real)> for Vec3 {
-    fn from((x, y, z): (flt::Real, flt::Real, flt::Real)) -> Self {
+impl<R> From<(R, R, R)> for Vec3<R> {
+    fn from((x, y, z): (R, R, R)) -> Self {
         Self { x, y, z }
     }
 }
 
-impl From<Vec3u> for Vec3 {
+impl<R: Real> From<Vec3u> for Vec3<R> {
     fn from(v: Vec3u) -> Self {
-        Self::new(v.x() as _, v.y() as _, v.z() as _)
+        Self::new(
+            R::from_primitive(v.x()),
+            R::from_primitive(v.y()),
+            R::from_primitive(v.z())
+        )
     }
 }
 
-impl AddAssign for Vec3 {
+impl<R: Real> AddAssign for Vec3<R> {
     fn add_assign(&mut self, other: Self) {
         self.x += other.x;
         self.y += other.y;
@@ -292,7 +333,7 @@ impl AddAssign for Vec3 {
     }
 }
 
-impl Add for Vec3 {
+impl<R: Real> Add for Vec3<R> {
     type Output = Self;
     fn add(self, other: Self) -> Self {
         let mut work = self.clone();
@@ -301,7 +342,7 @@ impl Add for Vec3 {
     }
 }
 
-impl Sub for Vec3 {
+impl<R: Real> Sub for Vec3<R> {
     type Output = Self;
     fn sub(self, other: Self) -> Self {
         Self {
@@ -312,7 +353,7 @@ impl Sub for Vec3 {
     }
 }
 
-impl Neg for Vec3 {
+impl<R: Real> Neg for Vec3<R> {
     type Output = Self;
     fn neg(self) -> Self {
         Self {
@@ -323,77 +364,45 @@ impl Neg for Vec3 {
     }
 }
 
-impl MulAssign<flt::Real> for Vec3 {
-    fn mul_assign(&mut self, other: flt::Real) {
+impl<R: Real> MulAssign<R> for Vec3<R> {
+    fn mul_assign(&mut self, other: R) {
         self.x *= other;
         self.y *= other;
         self.z *= other;
     }
 }
 
-impl Mul<flt::Real> for Vec3 {
+impl<R: Real> Mul<R> for Vec3<R> {
     type Output = Self;
-    fn mul(self, other: flt::Real) -> Self {
+    fn mul(self, other: R) -> Self {
         let mut work = self.clone();
         work *= other;
         work
     }
 }
 
-impl MulAssign<Flt> for Vec3 {
-    fn mul_assign(&mut self, other: Flt) {
-        self.x *= other;
-        self.y *= other;
-        self.z *= other;
+impl<R: Real> DivAssign<R> for Vec3<R> {
+    fn div_assign(&mut self, other: R) {
+        *self *= R::one() / other;
     }
 }
 
-impl Mul<Flt> for Vec3 {
+impl<R: Real> Div<R> for Vec3<R> {
     type Output = Self;
-    fn mul(self, other: Flt) -> Self {
-        let mut work = self.clone();
-        work *= other;
-        work
-    }
-}
-
-impl DivAssign<flt::Real> for Vec3 {
-    fn div_assign(&mut self, other: flt::Real) {
-        *self *= flt::Real::from_inner(1.0) / other;
-    }
-}
-
-impl Div<flt::Real> for Vec3 {
-    type Output = Self;
-    fn div(self, other: flt::Real) -> Self {
+    fn div(self, other: R) -> Self {
         let mut work = self.clone();
         work /= other;
         work
     }
 }
 
-impl DivAssign<Flt> for Vec3 {
-    fn div_assign(&mut self, other: Flt) {
-        *self *= 1.0 / other;
-    }
-}
-
-impl Div<Flt> for Vec3 {
-    type Output = Self;
-    fn div(self, other: Flt) -> Self {
-        let mut work = self.clone();
-        work /= other;
-        work
-    }
-}
-
-impl Sum for Vec3 {
+impl<R: Real> Sum for Vec3<R> {
     fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
-        iter.fold(Self::default(), |a, b| a + b)
+        iter.fold(Self::zero(), |a, b| a + b)
     }
 }
 
-impl fmt::LowerExp for Vec3 {
+impl<R: Real> fmt::LowerExp for Vec3<R> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
         write!(f, "(")?;
         fmt::LowerExp::fmt(&self.x(), f)?;
@@ -405,7 +414,7 @@ impl fmt::LowerExp for Vec3 {
     }
 }
 
-impl fmt::Display for Vec3 {
+impl<R: Real> fmt::Display for Vec3<R> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
         write!(f, "(")?;
         fmt::Display::fmt(&self.x(), f)?;
diff --git a/src/geom/vecu.rs b/src/geom/vecu.rs
index 9a79ba4..117447a 100644
--- a/src/geom/vecu.rs
+++ b/src/geom/vecu.rs
@@ -1,4 +1,5 @@
-use super::Vec3;
+use crate::geom::Vec3;
+use crate::real::{Real as _, ToFloat as _};
 
 use serde::{Serialize, Deserialize};
 use std::fmt::{self, Display};
@@ -41,7 +42,7 @@ impl From<(u32, u32, u32)> for Vec3u {
 
 impl From<Vec3> for Vec3u {
     fn from(v: Vec3) -> Self {
-        Self::new(v.x() as _, v.y() as _, v.z() as _)
+        Self::new(v.x().to_f64() as _, v.y().to_f64() as _, v.z().to_f64() as _)
     }
 }
 
diff --git a/src/lib.rs b/src/lib.rs
index 2238b60..78d5177 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -33,8 +33,8 @@ pub(crate) mod flt {
     use decorum::{R32, R64};
     // Total enforces zero constraints (allows NaN, INFs).
     // decorum::Real forbids NaN and INFs
-    pub(crate) type T64 = decorum::Total<f64>;
-    pub(crate) type T32 = decorum::Total<f32>;
+    pub type T64 = decorum::Total<f64>;
+    pub type T32 = decorum::Total<f32>;
 
     // Use checked floats for debug and fast floats for release
     #[cfg(debug)]
@@ -42,7 +42,7 @@ pub(crate) mod flt {
         use super::*;
         //pub(crate) type Real = R32;
         //pub type Flt = f32;
-        pub(crate) type Real = R64;
+        pub type Real = R64;
         pub type Flt = f64;
     }
     #[cfg(not(debug))]
@@ -51,12 +51,13 @@ pub(crate) mod flt {
         // pub(crate) type Real = T32;
         // pub type Flt = f32;
         // pub(crate) type Real = T64;
-        pub(crate) type Real = R64;
+        pub type Real = R64;
         pub type Flt = f64;
     }
     pub use defaults::*;
 }
 pub use flt::Flt;
+pub use flt::Real;
 
 #[allow(non_snake_case, unused)]
 pub mod consts {
diff --git a/src/mat.rs b/src/mat.rs
index 2576fa7..7468929 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -1,6 +1,7 @@
 use crate::{CellState, consts};
 use crate::flt::{Flt, Real};
 use crate::geom::{Line2d, Vec2, Vec3, Polygon2d};
+use crate::real::Real as _;
 use crate::sim::{PmlParameters, PmlState, StepParameters, StepParametersMut};
 
 use lazy_static::lazy_static;
@@ -95,7 +96,7 @@ pub struct Conductor {
 
 impl Conductor {
     pub fn new(conductivity: Flt) -> Self {
-        Self::new_anisotropic(Vec3::unit() * conductivity)
+        Self::new_anisotropic(Vec3::unit() * Real::from_f64(conductivity))
     }
     pub fn new_anisotropic(conductivity: Vec3) -> Self {
         Self {
@@ -120,7 +121,7 @@ pub struct LinearMagnet {
 
 impl LinearMagnet {
     pub fn new(relative_permeability: Flt) -> Self {
-        Self::new_anisotropic(Vec3::unit() * relative_permeability)
+        Self::new_anisotropic(Vec3::unit() * Real::from_f64(relative_permeability))
     }
     pub fn new_anisotropic(relative_permeability: Vec3) -> Self {
         Self { relative_permeability, m: Vec3::zero() }
@@ -182,9 +183,21 @@ impl<T: PiecewiseLinearFerromagnet> Material for T {
         let target_hm = h + m + delta_b * consts::real::MU0_INV();
 
         // TODO: this is probably not the best way to generalize a BH curve into 3d.
-        let (_hx, mx) = mh_curve.move_to(h.x(), m.x(), target_hm.x());
-        let (_hy, my) = mh_curve.move_to(h.y(), m.y(), target_hm.y());
-        let (_hz, mz) = mh_curve.move_to(h.z(), m.z(), target_hm.z());
+        let (_hx, mx) = mh_curve.move_to(
+            Flt::from_primitive(h.x()),
+            Flt::from_primitive(m.x()),
+            Flt::from_primitive(target_hm.x())
+        );
+        let (_hy, my) = mh_curve.move_to(
+            Flt::from_primitive(h.y()),
+            Flt::from_primitive(m.y()),
+            Flt::from_primitive(target_hm.y())
+        );
+        let (_hz, mz) = mh_curve.move_to(
+            Flt::from_primitive(h.z()),
+            Flt::from_primitive(m.z()),
+            Flt::from_primitive(target_hm.z())
+        );
 
         *self.m_mut() = Vec3::new(mx, my, mz);
         // let ret = Vec3::new(hx, hy, hz);
diff --git a/src/meas.rs b/src/meas.rs
index ba6bff6..c321009 100644
--- a/src/meas.rs
+++ b/src/meas.rs
@@ -1,6 +1,7 @@
-use crate::flt::Flt;
+use crate::flt::{Flt, Real};
 use crate::geom::{Meters, Region, Torus, Vec3, WorldRegion};
 use crate::mat::Material as _;
+use crate::real::Real as _;
 use crate::sim::GenericSim;
 use common_macros::b_tree_map;
 use dyn_clone::{self, DynClone};
@@ -84,19 +85,19 @@ impl Current {
             region: Box::new(r)
         }
     }
-    fn data(&self, state: &dyn GenericSim) -> (Flt, Vec3) {
+    fn data(&self, state: &dyn GenericSim) -> (Real, Vec3) {
         let FieldSample(volume, current_mag, current_vec) = state.map_sum_over_enumerated(&*self.region, |coord: Meters, _cell| {
             let current = state.current(coord);
             FieldSample(1, current.mag(), current)
         });
-        let mean_current_mag = current_mag / (volume as Flt);
-        let mean_current_vec = current_vec / (volume as Flt);
+        let mean_current_mag = current_mag / (Real::from_u32(volume));
+        let mean_current_vec = current_vec / (Real::from_u32(volume));
         (mean_current_mag, mean_current_vec)
     }
 }
 
 #[derive(Default)]
-struct FieldSample(usize, Flt, Vec3);
+struct FieldSample(u32, Real, Vec3);
 
 impl std::iter::Sum for FieldSample {
     fn sum<I>(iter: I) -> Self
@@ -172,7 +173,7 @@ impl CurrentLoop {
             region: r,
         }
     }
-    fn data(&self, state: &dyn GenericSim) -> Flt {
+    fn data(&self, state: &dyn GenericSim) -> Real {
         let FieldSample(volume, directed_current, _current_vec) = state.map_sum_over_enumerated(&self.region, |coord: Meters, _cell| {
             let normal = self.region.axis();
             let to_coord = *coord - *self.region.center();
@@ -181,7 +182,7 @@ impl CurrentLoop {
             let directed_current = current.dot(tangent);
             FieldSample(1, directed_current, current)
         });
-        let mean_directed_current = directed_current / (volume as Flt);
+        let mean_directed_current = directed_current / Real::from_u32(volume);
         let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size());
         let cross_sectional_current = mean_directed_current * cross_section;
         cross_sectional_current
@@ -217,7 +218,7 @@ impl MagneticLoop {
             region: r,
         }
     }
-    fn data(&self, state: &dyn GenericSim) -> (Flt, Flt, Flt) {
+    fn data(&self, state: &dyn GenericSim) -> (Real, Real, Real) {
         let FieldSamples([
             FieldSample(volume, directed_m, _m_vec),
             FieldSample(_, directed_b, _b_vec),
@@ -241,11 +242,11 @@ impl MagneticLoop {
             ])
         });
         // let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size());
-        let mean_directed_m = directed_m / (volume as Flt);
+        let mean_directed_m = directed_m / Real::from_u32(volume);
         // let cross_sectional_m = mean_directed_m * cross_section;
-        let mean_directed_b = directed_b / (volume as Flt);
+        let mean_directed_b = directed_b / Real::from_u32(volume);
         // let cross_sectional_b = mean_directed_b * cross_section;
-        let mean_directed_h = directed_h / (volume as Flt);
+        let mean_directed_h = directed_h / Real::from_u32(volume);
         // let cross_sectional_h = mean_directed_h * cross_section;
         // format!(
         //     "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}",
@@ -298,7 +299,7 @@ impl MagneticFlux {
             let mag = b.mag();
             FieldSample(1, mag, b)
         });
-        let mean_mag = mag_vec / (volume as Flt);
+        let mean_mag = mag_vec / Real::from_u32(volume);
         mean_mag
     }
 }
@@ -337,7 +338,7 @@ impl Magnetization {
             let mag = m.mag();
             FieldSample(1, mag, m)
         });
-        let mean_mag = mag_vec / (volume as Flt);
+        let mean_mag = mag_vec / Real::from_u32(volume);
         mean_mag
     }
 }
@@ -357,7 +358,7 @@ impl AbstractMeasurement for Magnetization {
 }
 
 fn loc(v: Meters) -> String {
-    format!("{:.0} um", *v * 1e6)
+    format!("{:.0} um", *v * Real::from_u64(1_000_000))
 }
 
 /// M
@@ -447,7 +448,7 @@ impl Energy {
             region: Box::new(region),
         }
     }
-    fn data(&self, state: &dyn GenericSim) -> Flt {
+    fn data(&self, state: &dyn GenericSim) -> Real {
         // Potential energy stored in a E/M field:
         //   https://en.wikipedia.org/wiki/Magnetic_energy
         //   https://en.wikipedia.org/wiki/Electric_potential_energy#Energy_stored_in_an_electrostatic_field_distribution
@@ -456,7 +457,7 @@ impl Energy {
         //   U(E) = 1/2 \int E . D dV
         #[allow(non_snake_case)]
         let dV = state.feature_volume();
-        let e = 0.5 * dV * state.map_sum_over(&*self.region, |cell| {
+        let e = Real::from_primitive(0.5 * dV) * state.map_sum_over(&*self.region, |cell| {
             // E . D = E . (E + P) = E.E since we don't model polarization fields
             cell.h().dot(cell.b()) + cell.e().mag_sq()
         });
@@ -498,11 +499,11 @@ impl Power {
         // Power is P = IV = A*J*V = L^2*J.(LE) = L^3 J.E
         // where L is feature size.
         #[allow(non_snake_case)]
-        let dV = state.feature_volume();
+        let dV = Real::from_primitive(state.feature_volume());
         let power = dV * state.map_sum_over(&*self.region, |cell| {
             cell.current_density().dot(cell.e())
         });
-        power
+        Flt::from_primitive(power)
     }
 }
 
diff --git a/src/real.rs b/src/real.rs
index d8609ff..8c8f682 100644
--- a/src/real.rs
+++ b/src/real.rs
@@ -1,9 +1,10 @@
-use std::ops::AddAssign;
+use std::fmt;
+use std::ops::{AddAssign, MulAssign};
 use decorum::cmp::IntrinsicOrd;
+pub use decorum::Real as decorum_Real;
+pub use num::Zero;
 
-/// This exists to allow configuration over # of bits (f32 v.s. f64) as well as
-/// constraints.
-pub trait Real: decorum::Real + IntrinsicOrd + AddAssign {
+pub trait ToFloat {
     fn to_f32(&self) -> f32 {
         self.to_f64() as _
     }
@@ -12,39 +13,142 @@ pub trait Real: decorum::Real + IntrinsicOrd + AddAssign {
     }
 }
 
-impl Real for f32 {
+/// This exists to allow configuration over # of bits (f32 v.s. f64) as well as
+/// constraints.
+pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign + fmt::LowerExp + fmt::Display {
+    fn from_f32(f: f32) -> Self {
+        Self::from_f64(f as _)
+    }
+    fn from_f64(f: f64) -> Self {
+        Self::from_f32(f as _)
+    }
+    fn from_u64(u: u64) -> Self {
+        Self::from_f64(u as _)
+    }
+    fn from_u32(u: u32) -> Self {
+        Self::from_u64(u as _)
+    }
+
+    fn from_primitive<P: ToFloat>(p: P) -> Self {
+        // TODO: fold with from_<blah>
+        Self::from_f64(p.to_f64())
+    }
+
+    fn tenth() -> Self {
+        Self::one() / Self::ten()
+    }
+    fn third() -> Self {
+        Self::one() / Self::three()
+    }
+    fn half() -> Self {
+        Self::from_f32(0.5)
+    }
+    fn two() -> Self {
+        Self::from_f32(2.0)
+    }
+    fn three() -> Self {
+        Self::from_f32(3.0)
+    }
+    fn ten() -> Self {
+        Self::from_f32(10.0)
+    }
+}
+
+impl ToFloat for f32 {
     fn to_f32(&self) -> f32 {
         *self
     }
 }
 
+impl Real for f32 {
+    fn from_f32(f: f32) -> Self {
+        f
+    }
+    fn from_u64(u: u64) -> Self {
+        Self::from_f32(u as _)
+    }
+}
+
+impl ToFloat for f64 {
+    fn to_f64(&self) -> f64 {
+        *self
+    }
+}
+
 impl Real for f64 {
-    fn to_f64(&self) -> f64 {
-        *self
+    fn from_f64(f: f64) -> Self {
+        f
+    }
+}
+
+impl ToFloat for decorum::R32 {
+    fn to_f32(&self) -> f32 {
+        self.into_inner()
     }
 }
 
 impl Real for decorum::R32 {
-    fn to_f32(&self) -> f32 {
+    fn from_f32(f: f32) -> Self {
+        Self::from_inner(f)
+    }
+    fn from_u64(u: u64) -> Self {
+        Self::from_f32(u as _)
+    }
+}
+
+impl ToFloat for decorum::R64 {
+    fn to_f64(&self) -> f64 {
         self.into_inner()
     }
 }
 
 impl Real for decorum::R64 {
-    fn to_f64(&self) -> f64 {
-        self.into_inner()
+    fn from_f64(f: f64) -> Self {
+        Self::from_inner(f)
     }
 }
 
-
-impl Real for decorum::N32 {
+impl ToFloat for decorum::N32 {
     fn to_f32(&self) -> f32 {
         self.into_inner()
     }
 }
 
-impl Real for decorum::N64 {
+impl Real for decorum::N32 {
+    fn from_f32(f: f32) -> Self {
+        Self::from_inner(f)
+    }
+    fn from_u64(u: u64) -> Self {
+        Self::from_f32(u as _)
+    }
+}
+
+impl ToFloat for decorum::N64 {
     fn to_f64(&self) -> f64 {
         self.into_inner()
     }
 }
+
+impl Real for decorum::N64 {
+    fn from_f64(f: f64) -> Self {
+        Self::from_inner(f)
+    }
+}
+
+impl ToFloat for u32 {
+    fn to_f32(&self) -> f32 {
+        *self as _
+    }
+    fn to_f64(&self) -> f64 {
+        *self as _
+    }
+}
+
+impl ToFloat for u64 {
+    fn to_f32(&self) -> f32 {
+        *self as _
+    }
+    fn to_f64(&self) -> f64 {
+        *self as _
+    }
+}
diff --git a/src/render.rs b/src/render.rs
index 00db479..bc9ffc8 100644
--- a/src/render.rs
+++ b/src/render.rs
@@ -1,6 +1,7 @@
 use crate::geom::{Index, Meters, Vec2, Vec3, Vec3u};
 use crate::{Material as _, MaterialExt as _};
 use crate::mat;
+use crate::real::{Real as _, ToFloat as _};
 use crate::sim::{Cell, GenericSim, StaticSim};
 use crate::meas::AbstractMeasurement;
 use crossterm::{cursor, QueueableCommand as _};
@@ -152,7 +153,7 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
         trace!("rendering at {}x{} with z={}", width, height, z);
         let mut me = Self::new(state, measurements, width, height, z);
         me.render_scalar_field(10.0,  false, 2, |cell| {
-            cell.mat().conductivity().mag() as f32 + if cell.mat().is_vacuum() {
+            cell.mat().conductivity().mag().to_f32() + if cell.mat().is_vacuum() {
                 0.0
             } else {
                 5.0
@@ -193,12 +194,12 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
         let y_prop = y_px as f32 / self.im.height() as f32;
         let y_m = y_prop * (self.sim.height() as f32 * self.sim.feature_size() as f32);
         let z_m = self.z as f32 * self.sim.feature_size() as f32;
-        self.sim.sample(Meters(Vec3::new(x_m as _, y_m as _, z_m as _)))
+        self.sim.sample(Meters(Vec3::new(x_m, y_m, z_m)))
     }
 
     //////////////   Ex/Ey/Bz configuration ////////////
     fn render_b_z_field(&mut self, scale: f32) {
-        self.render_scalar_field(1.0e-4 * scale, true, 1, |cell| cell.b().z() as f32);
+        self.render_scalar_field(1.0e-4 * scale, true, 1, |cell| cell.b().z().to_f32());
     }
     fn render_e_xy_field(&mut self, scale: f32) {
         self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 100.0 * scale, |cell| cell.e().xy().to_f32());
@@ -209,24 +210,24 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
     }
     //////////////   Magnitude configuration /////////////
     fn render_b(&mut self, scale: f32) {
-        self.render_scalar_field(1.0e-3 * scale, false, 1, |cell| cell.b().mag() as f32);
+        self.render_scalar_field(1.0e-3 * scale, false, 1, |cell| cell.b().mag().to_f32());
     }
     fn render_current(&mut self, scale: f32) {
         self.render_scalar_field(1.0e1 * scale, false, 0, |cell| {
-            cell.e().elem_mul(cell.mat().conductivity()).mag() as f32
+            cell.e().elem_mul(cell.mat().conductivity()).mag().to_f32()
         });
     }
 
     //////////////   Bx/By/Ez configuration ////////////
     fn render_e_z_field(&mut self, scale: f32) {
-        self.render_scalar_field(1e4 * scale, true, 1, |cell| cell.e().z() as f32);
+        self.render_scalar_field(1e4 * scale, true, 1, |cell| cell.e().z().to_f32());
     }
     fn render_b_xy_field(&mut self, scale: f32) {
         self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e-9 * scale, |cell| cell.b().xy().to_f32());
     }
 
     fn render_m(&mut self, scale: f32) {
-        self.render_scalar_field(1.0e5 * scale, false, 1, |cell| cell.mat().m().mag() as f32);
+        self.render_scalar_field(1.0e5 * scale, false, 1, |cell| cell.mat().m().mag().to_f32());
         self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e5 * scale, |cell| cell.mat().m().xy().to_f32());
     }
 
@@ -558,7 +559,7 @@ impl<S: GenericSim> Renderer<S> for PlotlyRenderer {
                     yv.push(y);
                     zv.push(z);
                     // opacities.push((cell.e().mag() * 0.1).min(1.0) as f64)
-                    let mat = cell.mat().conductivity().mag() + if cell.mat().is_vacuum() {
+                    let mat = cell.mat().conductivity().mag().to_f32() + if cell.mat().is_vacuum() {
                         0.0
                     } else {
                         5.0
@@ -566,9 +567,9 @@ impl<S: GenericSim> Renderer<S> for PlotlyRenderer {
                     //let g = scale_unsigned_to_u8(mat, 10.0);
                     //let r = scale_unsigned_to_u8(cell.mat().m().mag(), 100.0);
                     //let b = scale_unsigned_to_u8(cell.e().mag(), 1e2);
-                    let r = scale_unsigned_to_u8(cell.mat().m().mag() as f32, 100.0);
-                    let g = scale_unsigned_to_u8(cell.e().mag() as f32, 1e2);
-                    let b = scale_unsigned_to_u8(mat as f32, 10.0);
+                    let r = scale_unsigned_to_u8(cell.mat().m().mag().to_f32(), 100.0);
+                    let g = scale_unsigned_to_u8(cell.e().mag().to_f32(), 1e2);
+                    let b = scale_unsigned_to_u8(mat, 10.0);
                     let alpha = 1.0;
                     colors.push(Rgba::new(r, g, b, alpha));
                 }
@@ -666,6 +667,7 @@ impl<S: GenericSim + Clone + Serialize> Renderer<S> for SerializerRenderer {
     }
     fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], _config: RenderConfig) {
         let frame = SerializedFrame {
+            // N.B.: `to_static` isn't 100% necessary here, but it tends to result in a smaller file.
             state: state.to_static(),
             measurements: measurements.iter().cloned().collect(),
         };
diff --git a/src/sim.rs b/src/sim.rs
index 0420c91..9510dba 100644
--- a/src/sim.rs
+++ b/src/sim.rs
@@ -1,6 +1,7 @@
 use crate::{flt::{Flt, Real}, consts};
 use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u};
 use crate::mat::{self, GenericMaterial, Material, MaterialExt as _};
+use crate::real::{self, Real as _};
 use crate::stim::AbstractStimulus;
 use dyn_clone::{self, DynClone};
 use log::trace;
@@ -207,7 +208,7 @@ pub trait GenericSim: Send + Sync + DynClone {
     }
     fn volume(&self) -> Flt {
         let s = self.size().to_meters(self.feature_size());
-        s.x() * s.y() * s.z()
+        Flt::from_primitive(s.x() * s.y() * s.z())
     }
     fn timestep(&self) -> Flt;
     fn step_no(&self) -> u64;
@@ -309,7 +310,7 @@ impl<'a> dyn GenericSim + 'a {
     }
 
     pub fn current<C: Coord>(&self, c: C) -> Vec3 {
-        self.get(c).current_density() * (self.feature_size() * self.feature_size())
+        self.get(c).current_density() * Real::from_primitive(self.feature_size() * self.feature_size())
     }
 
     pub fn impulse_e<C: Coord>(&mut self, pos: C, amt: Vec3) {
@@ -454,7 +455,7 @@ impl<M: Material + Clone + Default + Send + Sync + 'static> SimState<M> {
         Zip::from(ndarray::indices_of(&self.e)).and(&mut self.e).par_apply(
             |(z, y, x), e| {
                 let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size);
-                let value = stim.at(t_sec, pos_meters) * timestep;
+                let value = stim.at(t_sec, pos_meters) * Real::from_primitive(timestep);
                 *e += value;
             });
         trace!("apply_stimulus end");
@@ -740,26 +741,26 @@ impl<M> Cell<M> {
         self.state.e()
     }
     pub fn ex(&self) -> Flt {
-        self.e().x()
+        Flt::from_primitive(self.e().x())
     }
     pub fn ey(&self) -> Flt {
-        self.e().y()
+        Flt::from_primitive(self.e().y())
     }
     pub fn ez(&self) -> Flt {
-        self.e().z()
+        Flt::from_primitive(self.e().z())
     }
 
     pub fn h(&self) -> Vec3 {
         self.state.h()
     }
     pub fn hx(&self) -> Flt {
-        self.h().x()
+        Flt::from_primitive(self.h().x())
     }
     pub fn hy(&self) -> Flt {
-        self.h().y()
+        Flt::from_primitive(self.h().y())
     }
     pub fn hz(&self) -> Flt {
-        self.h().z()
+        Flt::from_primitive(self.h().z())
     }
 }
 
@@ -768,13 +769,13 @@ impl<M: Material> Cell<M> {
         (self.state.h() + self.mat.m()) * consts::real::MU0()
     }
     pub fn bx(&self) -> Flt {
-        self.b().x()
+        Flt::from_primitive(self.b().x())
     }
     pub fn by(&self) -> Flt {
-        self.b().y()
+        Flt::from_primitive(self.b().y())
     }
     pub fn bz(&self) -> Flt {
-        self.b().z()
+        Flt::from_primitive(self.b().z())
     }
 
     pub fn mat(&self) -> &M {
@@ -795,17 +796,17 @@ impl<M: Material> Cell<M> {
 /// Solve df/dt in:
 /// nabla_g = df_dt_coeff * df/dt + f_coeff * f + f_int_coeff * \int f + f_int_int_coeff * \int \int f
 #[allow(unused)]
-fn solve_step_diff_eq(
-    nabla_g: Vec3,
-    df_dt_coeff: Vec3,
-    f_coeff: Vec3,
-    f_int_coeff: Vec3,
-    f_int_int_coeff: Vec3,
-    delta_t: Real,
-    f_prev: Vec3,
-    f_int_prev: Vec3,
-    f_int_int_prev: Vec3,
-) -> Vec3 {
+fn solve_step_diff_eq<R: real::Real>(
+    nabla_g: Vec3<R>,
+    df_dt_coeff: Vec3<R>,
+    f_coeff: Vec3<R>,
+    f_int_coeff: Vec3<R>,
+    f_int_int_coeff: Vec3<R>,
+    delta_t: R,
+    f_prev: Vec3<R>,
+    f_int_prev: Vec3<R>,
+    f_int_int_prev: Vec3<R>,
+) -> Vec3<R> {
     // ```tex
     // $f(t)$ from $(T-\Delta t)$ to $(T+\Delta t)$ is a line. As such:
     // $f(t') = f_prev + t'*db/dt$ where $t'$ is the time since the last sample of $f$, i.e.  $f_prev$
@@ -813,18 +814,16 @@ fn solve_step_diff_eq(
     // $\int f|t' = f_int_prev + t'*f_prev + 1/2 t'^2 df/dt$
     // $\int \int f|t' = f_int_int_prev + t'*f_int_prev + 1/2 t'^2*f_prev + 1/6 t'^3*df/dt$
     // ```
-    use consts::real::{HALF, THIRD};
-
     // Solve coeffs in:
     // nabla_g = df_dt_coeff * df/dt + f_prev_coeff * f_prev
     //   + f_int_prev_coeff * f_int_prev + f_int_int_prev_coeff * f_int_int_prev
     let f_int_int_prev_coeff = f_int_int_coeff;
     let f_int_prev_coeff = f_int_coeff + f_int_int_coeff*delta_t;
     let f_prev_coeff = f_coeff + (f_int_coeff +
-        f_int_int_coeff*HALF()
+        f_int_int_coeff*R::half()
     )*delta_t;
     let df_dt_coeff = df_dt_coeff + (f_coeff +
-        (f_int_coeff + f_int_int_coeff*THIRD()*delta_t)*HALF()*delta_t
+        (f_int_coeff + f_int_int_coeff*R::third()*delta_t)*R::half()*delta_t
     )*delta_t;
 
     // Solve coeffs in nabla_g = df_dt_coeff + offset;
@@ -1161,6 +1160,7 @@ mod test {
     use super::*;
     use crate::geom::{Cube, Region, WorldRegion};
     use crate::mat::{Conductor, Static};
+    use crate::real::Zero as _;
     use float_eq::assert_float_eq;
     use more_asserts::*;
     use rand::rngs::StdRng;
@@ -1170,7 +1170,7 @@ mod test {
         let mut cf = ConvFlt(Real::from(0.0));
         for (s, e) in signal.iter().zip(expected.iter()) {
             let v = cf.step(a.into(), b.into(), c.into(), s.clone().into(), delta_t.into());
-            assert_float_eq!(v.into_inner(), e, r2nd <= 1e-6);
+            assert_float_eq!(Flt::from_primitive(v), e, r2nd <= 1e-6);
         }
     }
 
@@ -1273,7 +1273,9 @@ mod test {
         do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
     }
 
-    fn assert_vec3_eq(actual: Vec3, expected: Vec3) {
+    fn assert_vec3_eq<R: real::Real>(actual: Vec3<R>, expected: Vec3<R>) {
+        let actual = actual.cast::<f64>();
+        let expected = expected.cast::<f64>();
         assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6);
         assert_float_eq!(actual.y(), expected.y(), r2nd <= 1e-6);
         assert_float_eq!(actual.z(), expected.z(), r2nd <= 1e-6);
@@ -1285,11 +1287,11 @@ mod test {
         // nabla_g = 0.1 df/dt
         let df_dt = solve_step_diff_eq(
             Vec3::new(2.0, 3.0, 4.0), // nabla_g
-            Vec3::unit()*0.1, // df_dt_coeff
+            Vec3::uniform(0.1), // df_dt_coeff
             Vec3::zero(),
             Vec3::zero(),
             Vec3::zero(),
-            consts::real::HALF(), // delta_t
+            Real::half(), // delta_t
             Vec3::zero(),
             Vec3::zero(),
             Vec3::zero(),
@@ -1303,11 +1305,11 @@ mod test {
         // nabla_g = 0.1 df/dt - f - 3 f_int + 0.4 f_int_int
         let df_dt = solve_step_diff_eq(
             Vec3::new(2.0, 3.0, 4.0), // nabla_g
-            Vec3::unit()*0.1, // df_dt_coeff
-            Vec3::unit()*-1.0, // f_coeff
-            Vec3::unit()*-3.0, // f_int_coeff
-            Vec3::unit()*0.4, // f_int_int_coeff
-            consts::real::ZERO(), // delta_t
+            Vec3::uniform(0.1), // df_dt_coeff
+            Vec3::uniform(-1.0), // f_coeff
+            Vec3::uniform(-3.0), // f_int_coeff
+            Vec3::uniform(0.4), // f_int_int_coeff
+            Real::zero(), // delta_t
             Vec3::zero(),
             Vec3::zero(),
             Vec3::zero(),
@@ -1320,16 +1322,16 @@ mod test {
         // nabla_g = 0.1 df/dt + 2 f
         // Let f(t) = sin(2 t)
         // then f'(t) = 2 cos(2 t)
-        let f = Vec3::unit()*(5.0).sin();
-        let df_dt = Vec3::unit()*2.0 * (5.0).cos();
-        let nabla_g = df_dt*0.1 + f*2.0;
+        let f = Vec3::uniform((5.0).sin());
+        let df_dt = Vec3::uniform(2.0 * (5.0).cos());
+        let nabla_g = df_dt*Real::tenth() + f*Real::two();
         let actual_df_dt = solve_step_diff_eq(
             nabla_g,
-            Vec3::unit()*0.1,
-            Vec3::unit()*2.0,
+            Vec3::uniform(0.1),
+            Vec3::uniform(2.0),
             Vec3::zero(),
             Vec3::zero(),
-            consts::real::ZERO(), // delta_t
+            Real::zero(), // delta_t
             f, // f_prev
             Vec3::zero(),
             Vec3::zero(),
@@ -1343,17 +1345,17 @@ mod test {
         // Let f(t) = sin(2 t)
         // then f'(t) = 2 cos(2 t)
         // then \int f(t) = -1/2 cos(2 t)
-        let f = Vec3::unit()*(5.0).sin();
-        let df_dt = Vec3::unit()*2.0 * (5.0).cos();
-        let f_int = Vec3::unit()*-0.5 * (5.0).cos();
-        let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4;
+        let f = Vec3::uniform((5.0).sin());
+        let df_dt = Vec3::uniform(2.0 * (5.0).cos());
+        let f_int = Vec3::uniform(-0.5 * (5.0).cos());
+        let nabla_g = df_dt*Real::tenth() + f*Real::two() + f_int*Real::from_primitive(0.4);
         let actual_df_dt = solve_step_diff_eq(
             nabla_g,
-            Vec3::unit()*0.1,
-            Vec3::unit()*2.0,
-            Vec3::unit()*0.4,
+            Vec3::uniform(0.1),
+            Vec3::uniform(2.0),
+            Vec3::uniform(0.4),
             Vec3::zero(),
-            consts::real::ZERO(), // delta_t
+            Real::zero(), // delta_t
             f, // f_prev
             f_int, // f_int_prev
             Vec3::zero(),
@@ -1375,11 +1377,11 @@ mod test {
         let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4 + f_int_int*0.3;
         let actual_df_dt = solve_step_diff_eq(
             nabla_g,
-            Vec3::unit()*0.1,
-            Vec3::unit()*2.0,
-            Vec3::unit()*0.4,
-            Vec3::unit()*0.3,
-            consts::real::ZERO(), // delta_t
+            Vec3::uniform(0.1),
+            Vec3::uniform(2.0),
+            Vec3::uniform(0.4),
+            Vec3::uniform(0.3),
+            f32::zero(), // delta_t
             f, // f_prev
             f_int, // f_int_prev
             f_int_int, // f_int_int_prev
@@ -1391,16 +1393,16 @@ mod test {
     fn step_diff_eq_1st_order_homogenous() {
         // 0 = f'(t) + f(t)
         // solution: f(t) = k*e^-t
-        let f = |t: Flt| -> Vec3 {
+        let f = |t: f32| -> Vec3<f32> {
             Vec3::unit()*(-t).exp()
         };
-        let df_dt = |t: Flt| -> Vec3 {
+        let df_dt = |t: f32| -> Vec3<f32> {
             Vec3::unit()*-(-t).exp()
         };
-        let f_int = |t: Flt| -> Vec3 {
+        let f_int = |t: f32| -> Vec3<f32> {
             Vec3::unit()*-(-t).exp()
         };
-        let f_int_int = |t: Flt| -> Vec3 {
+        let f_int_int = |t: f32| -> Vec3<f32> {
             Vec3::unit()*(-t).exp()
         };
         let tl = 3.0000;
@@ -1411,7 +1413,7 @@ mod test {
             Vec3::unit(), // f_coeff
             Vec3::zero(),
             Vec3::zero(),
-            (tm - tl).into(), // delta_t
+            tm - tl, // delta_t
             f(tl), // f_prev
             f_int(tl),
             f_int_int(tl),
@@ -1423,16 +1425,16 @@ mod test {
     fn step_diff_eq_1st_order_nonhomogenous() {
         // 1 = f'(t) + f(t)
         // let f(t) = 1 + e^-t
-        let f = |t: Flt| -> Vec3 {
+        let f = |t: f32| -> Vec3<f32> {
             Vec3::unit() + Vec3::unit()*(-t).exp()
         };
-        let df_dt = |t: Flt| -> Vec3 {
+        let df_dt = |t: f32| -> Vec3<f32> {
             Vec3::unit()*-(-t).exp()
         };
-        let f_int = |t: Flt| -> Vec3 {
+        let f_int = |t: f32| -> Vec3<f32> {
             Vec3::unit()*-(-t).exp()
         };
-        let f_int_int = |t: Flt| -> Vec3 {
+        let f_int_int = |t: f32| -> Vec3<f32> {
             Vec3::unit()*(-t).exp()
         };
         let tl = 3.0000;
@@ -1443,7 +1445,7 @@ mod test {
             Vec3::unit(), // f_coeff
             Vec3::zero(),
             Vec3::zero(),
-            (tm - tl).into(), // delta_t
+            tm - tl, // delta_t
             f(tl), // f_prev
             f_int(tl),
             f_int_int(tl),
@@ -1464,16 +1466,16 @@ mod test {
         let b = Vec3::new(0.5, -0.3, 0.0);
         let c = Vec3::new(-5.0, 0.0, -0.1);
         let d = -w.elem_mul(c + w.elem_mul(b + w.elem_mul(a)));
-        let f = |t: Flt| -> Vec3 {
+        let f = |t: f64| -> Vec3<f64> {
             V.elem_mul((w*t).exp())
         };
-        let df_dt = |t: Flt| -> Vec3 {
+        let df_dt = |t: f64| -> Vec3<f64> {
             V.elem_mul((w*t).exp()).elem_mul(w)
         };
-        let f_int = |t: Flt| -> Vec3 {
+        let f_int = |t: f64| -> Vec3<f64> {
             V.elem_mul((w*t).exp()).elem_div(w)
         };
-        let f_int_int = |t: Flt| -> Vec3 {
+        let f_int_int = |t: f64| -> Vec3<f64> {
             V.elem_mul((w*t).exp()).elem_div(w).elem_div(w) + C.elem_div(d)
         };
         let tl = 3.00000000;
@@ -1484,7 +1486,7 @@ mod test {
             b, // f_coeff
             c, // int_f_coeff
             d, // int_int_f_coeff
-            (tm - tl).into(), // delta_t
+            tm - tl, // delta_t
             f(tl), // f_prev
             f_int(tl), // f_int_prev
             f_int_int(tl), // f_int_int_prev
@@ -1495,7 +1497,7 @@ mod test {
     fn energy(s: &dyn GenericSim) -> Flt {
         0.5 * s.feature_volume() * s.map_sum_enumerated(|_pos: Index, cell| {
             // println!("{}: {}", _pos, cell.e());
-            Real::from(cell.e().mag_sq()).into_inner()
+            Flt::from_primitive(cell.e().mag_sq())
         })
     }
 
@@ -1582,7 +1584,7 @@ mod test {
     #[test]
     fn anisotropic_conductor_inactive_x() {
         let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
-            (1e3, 0e3, 0e3).into()
+            Vec3::new_x(1e3)
         ));
         assert_gt!(energy_1, 0.9*energy_0);
         // XXX: I'm not sure why it _gains_ energy :o
@@ -1592,7 +1594,7 @@ mod test {
     #[test]
     fn anisotropic_conductor_inactive_y() {
         let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
-            (0e3, 1e3, 0e3).into()
+            Vec3::new_y(1e3)
         ));
         assert_gt!(energy_1, 0.9*energy_0);
         assert_lt!(energy_1, 1.5*energy_0);
@@ -1601,7 +1603,7 @@ mod test {
     #[test]
     fn anisotropic_conductor_active_z() {
         let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
-            (0e3, 0e3, 1e3).into()
+            Vec3::new_z(1e3)
         ));
         assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6);
     }
@@ -1610,8 +1612,8 @@ mod test {
         let mut state = StaticSim::new(size, 1e-6);
         let timestep = state.timestep();
         state.fill_boundary_using(size/4, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let conductivity = b * (0.5 / timestep);
+            let b = boundary_ness.elem_pow(Real::three());
+            let conductivity = b * (Real::half() / timestep);
             // K scales up to 11; dc_shift scales from 0.05:
             // https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.461.4606&rep=rep1&type=pdf
             // let propagation = Vec3::unit() + b * 10.0;
@@ -1629,8 +1631,8 @@ mod test {
         let mut state = StaticSim::new(size, 1e-6);
         let timestep = state.timestep();
         state.fill_boundary_using(size/4, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let conductivity = Vec3::unit() * (consts::real::EPS0() * b.mag() * 0.5 / timestep);
+            let b = boundary_ness.elem_pow(Real::three());
+            let conductivity = Vec3::unit() * (consts::real::EPS0() * b.mag() * Real::half() / timestep);
             Static {
                 conductivity, ..Default::default()
             }
@@ -1652,7 +1654,7 @@ mod test {
             let (e, hz) = (self.f)(pos.to_index(self.feat_size));
             let sig_angle = angle*hz;
             let sig = sig_angle.sin();
-            e*gate*sig / (self.t_end * self.sqrt_vol)
+            e*Real::from_primitive(gate*sig / (self.t_end * self.sqrt_vol))
         }
     }
 
@@ -1774,14 +1776,14 @@ mod test {
         let mut state = StaticSim::new(size, 1e-6);
         let timestep = state.timestep();
         state.fill_boundary_using(size/4, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let coord_stretch = b * (0.5 / timestep);
+            let b = boundary_ness.elem_pow(Real::three());
+            let coord_stretch = b * (Real::half() / timestep);
             // let coord_stretch = b * 0.5;
             // Force the stretch to be only along one axis
             let coord_stretch = match coord_stretch.to_tuple() {
-                (x, y, z) if x >= y && x >= z => (x, 0.0, 0.0),
-                (x, y, z) if y >= x && y >= z => (0.0, y, 0.0),
-                (x, y, z) if z >= x && z >= y => (0.0, 0.0, z),
+                (x, y, z) if x >= y && x >= z => (x, Real::zero(), Real::zero()),
+                (x, y, z) if y >= x && y >= z => (Real::zero(), y, Real::zero()),
+                (x, y, z) if z >= x && z >= y => (Real::zero(), Real::zero(), z),
                 _ => unreachable!(),
             }.into();
             Static::from_pml(coord_stretch)
@@ -1796,8 +1798,8 @@ mod test {
         let mut state = StaticSim::new(size, 1e-6);
         let timestep = state.timestep();
         state.fill_boundary_using(size/4, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let cs = b * (0.5 / timestep);
+            let b = boundary_ness.elem_pow(Real::three());
+            let cs = b * (Real::half() / timestep);
             // let cs = b * 0.5;
             // permute the stretching so that it shouldn't interfere with the wave
             let coord_stretch = shuffle(cs);
diff --git a/src/stim.rs b/src/stim.rs
index 4b73d95..1f9d8c9 100644
--- a/src/stim.rs
+++ b/src/stim.rs
@@ -1,5 +1,6 @@
 use crate::consts;
-use crate::flt::Flt;
+use crate::flt::{self, Flt};
+use crate::real::*;
 use crate::geom::{Meters, Region, Vec3};
 
 pub trait AbstractStimulus: Sync {
@@ -66,7 +67,7 @@ impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus
             let amount = self.stim.at(t_sec);
             let from_center_to_point = *pos - *self.center;
             let rotational = from_center_to_point.cross(*self.axis);
-            let impulse = rotational.with_mag(amount);
+            let impulse = rotational.with_mag(flt::Real::from_primitive(amount));
             impulse
         } else {
             Vec3::zero()
@@ -140,7 +141,7 @@ impl TimeVarying1 for Sinusoid1 {
 
 impl TimeVarying3 for Sinusoid3 {
     fn at(&self, t_sec: Flt) -> Vec3 {
-        self.amp * (t_sec * self.omega).sin()
+        self.amp * flt::Real::from_primitive(t_sec * self.omega).sin()
     }
 }
 
@@ -216,19 +217,19 @@ mod test {
     
     #[test]
     fn sinusoid3() {
-        let s = Sinusoid3::new((10.0, 1.0, -100.0).into(), 1000.0);
+        let s = Sinusoid3::new(Vec3::new(10.0, 1.0, -100.0), 1000.0);
         assert_eq!(s.at(0.0), Vec3::zero());
-        assert_approx_eq!(s.at(0.00025), (10.0, 1.0, -100.0).into());
+        assert_approx_eq!(s.at(0.00025), Vec3::new(10.0, 1.0, -100.0));
         assert_approx_eq!(s.at(0.00050), Vec3::zero());
-        assert_approx_eq!(s.at(0.00075), (-10.0, -1.0, 100.0).into());
+        assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0));
     }
     
     #[test]
     fn sinusoid3_from_wavelength() {
-        let s = Sinusoid3::from_wavelength((10.0, 1.0, -100.0).into(), 0.001);
+        let s = Sinusoid3::from_wavelength(Vec3::new(10.0, 1.0, -100.0), 0.001);
         assert_eq!(s.at(0.0), Vec3::zero());
-        assert_approx_eq!(s.at(0.00025), (10.0, 1.0, -100.0).into());
+        assert_approx_eq!(s.at(0.00025), Vec3::new(10.0, 1.0, -100.0));
         assert_approx_eq!(s.at(0.00050), Vec3::zero());
-        assert_approx_eq!(s.at(0.00075), (-10.0, -1.0, 100.0).into());
+        assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0));
     }
 }