From f50dbf7d4ed9c15c7edbd37a18fb15ca450a7eed Mon Sep 17 00:00:00 2001 From: Colin Date: Sun, 13 Sep 2020 22:38:00 -0700 Subject: [PATCH] Make the float depth compile-time configurable --- examples/em_reflection.rs | 12 ++----- examples/toroid2.rs | 1 + src/driver.rs | 12 +++---- src/geom.rs | 72 +++++++++++++++++++-------------------- src/lib.rs | 50 ++++++++++++++++++--------- src/mat.rs | 55 +++++++++++++++++------------- src/meas.rs | 24 ++++++++++++- src/render.rs | 30 ++++++++-------- src/sim.rs | 52 ++++++++++++++-------------- 9 files changed, 175 insertions(+), 133 deletions(-) diff --git a/examples/em_reflection.rs b/examples/em_reflection.rs index 9356a77..b11e079 100644 --- a/examples/em_reflection.rs +++ b/examples/em_reflection.rs @@ -5,17 +5,9 @@ fn main() { let width = 201; let mut driver = Driver::new(width, 101, 1e-3 /* feature size */); driver.add_y4m_renderer("em_reflection.y4m"); - driver.add_term_renderer(); + // driver.add_term_renderer(); - for inset in 0..20 { - for x in 0..width { - *driver.state.get_mut(x, inset).mat_mut() = mat::Static::conductor(0.1*(20.0 - inset as f64)).into(); - } - for y in 0..100 { - *driver.state.get_mut(inset, y).mat_mut() = mat::Static::conductor(0.1*(20.0 - inset as f64)).into(); - *driver.state.get_mut(width-1 - inset, y).mat_mut() = mat::Static::conductor(0.1*(20.0 - inset as f64)).into(); - } - } + driver.add_boundary(20, 0.1); for y in 75..100 { for x in 0..width { // from https://www.thoughtco.com/table-of-electrical-resistivity-conductivity-608499 diff --git a/examples/toroid2.rs b/examples/toroid2.rs index bffc381..e1471e7 100644 --- a/examples/toroid2.rs +++ b/examples/toroid2.rs @@ -16,6 +16,7 @@ fn main() { //driver.set_steps_per_frame(40); driver.add_y4m_renderer("toroid2.y4m"); // driver.add_term_renderer(); + driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current))); driver.add_measurement(meas::Current(width / 2 + (conductor_inner_rad + conductor_outer_rad) / 2, width / 2)); driver.add_measurement(meas::Current(width / 2 + conductor_inner_rad + 2, width / 2)); driver.add_measurement(meas::Magnetization(width / 2 + ferro_inner_rad + 2, width / 2)); diff --git a/src/driver.rs b/src/driver.rs index cd14d7b..268cdd7 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -1,4 +1,4 @@ -use crate::mat; +use crate::{flt, mat}; use crate::meas::{self, AbstractMeasurement}; use crate::render::{self, MultiRenderer, Renderer}; use crate::sim::SimState; @@ -17,11 +17,11 @@ pub struct Driver { } impl Driver { - pub fn new(width: u32, height: u32, feature_size: f64) -> Self { + pub fn new(width: u32, height: u32, feature_size: flt::Pub) -> Self { Driver { state: SimState::new(width, height, feature_size), steps_per_frame: 1, - measurements: vec![Box::new(meas::Time), Box::new(meas::Energy)], + measurements: vec![Box::new(meas::Time), Box::new(meas::Meta), Box::new(meas::Energy)], ..Default::default() } } @@ -49,10 +49,10 @@ impl Driver { self.add_renderer(render::ColorTermRenderer, "terminal"); } - pub fn add_boundary(&mut self, thickness: u32, base_conductivity: f64) { + pub fn add_boundary(&mut self, thickness: u32, base_conductivity: flt::Pub) { for inset in 0..thickness { let depth = thickness - inset; - let conductivity = base_conductivity * (depth*depth) as f64; + let conductivity = base_conductivity * (depth*depth) as flt::Pub; for x in inset..self.state.width() - inset { // left *self.state.get_mut(x, inset).mat_mut() = mat::Static::conductor(conductivity).into(); @@ -76,7 +76,7 @@ impl Driver { self.state.step(); self.time_spent_stepping += start_time.elapsed().unwrap(); if self.state.step_no() % (10*self.steps_per_frame) == 0 { - info!("fps: {:6.2}", (self.state.step_no() as f64) / self.time_spent_stepping.as_secs_f64()); + info!("fps: {:6.2}", (self.state.step_no() as flt::Pub) / (self.time_spent_stepping.as_secs_f64() as flt::Pub)); } } } diff --git a/src/geom.rs b/src/geom.rs index 2922f0c..40d3312 100644 --- a/src/geom.rs +++ b/src/geom.rs @@ -1,18 +1,18 @@ -use crate::F64; +use crate::flt; use std::ops::{Add, AddAssign, Mul, Neg, Sub}; -fn real(f: f64) -> F64 { - F64::from_inner(f) +fn real(f: flt::Pub) -> flt::Priv { + flt::Priv::from_inner(f) } -fn round(f: F64) -> F64 { - F64::from_inner(f.into_inner().round()) +fn round(f: flt::Priv) -> flt::Priv { + flt::Priv::from_inner(f.into_inner().round()) } #[derive(Copy, Clone, Debug, Default)] pub struct Point { - pub x: F64, - pub y: F64, + pub x: flt::Priv, + pub y: flt::Priv, } impl Add for Point { @@ -48,9 +48,9 @@ impl Sub for Point { } } -impl Mul for Point { +impl Mul for Point { type Output = Point; - fn mul(self, s: f64) -> Point { + fn mul(self, s: flt::Pub) -> Point { Point { x: self.x * s, y: self.y * s, @@ -58,25 +58,25 @@ impl Mul for Point { } } -impl Into<(f64, f64)> for Point { - fn into(self) -> (f64, f64) { +impl Into<(flt::Pub, flt::Pub)> for Point { + fn into(self) -> (flt::Pub, flt::Pub) { (self.x(), self.y()) } } impl Point { - pub fn new(x: f64, y: f64) -> Self { + pub fn new(x: flt::Pub, y: flt::Pub) -> Self { Self { x: x.into(), y: y.into(), } } - pub fn x(&self) -> f64 { + pub fn x(&self) -> flt::Pub { self.x.into() } - pub fn y(&self) -> f64 { + pub fn y(&self) -> flt::Pub { self.y.into() } @@ -87,20 +87,20 @@ impl Point { } } - pub fn mag_sq(&self) -> f64 { + pub fn mag_sq(&self) -> flt::Pub { (self.x*self.x + self.y*self.y).into() } - pub fn mag(&self) -> f64 { + pub fn mag(&self) -> flt::Pub { self.mag_sq().sqrt() } - pub fn with_mag(&self, new_mag: f64) -> Point { + pub fn with_mag(&self, new_mag: flt::Pub) -> Point { if new_mag == 0.0 { // avoid div-by-zero if self.mag() == 0 and new_mag == 0 Point::new(0.0, 0.0) } else { - let scale = F64::from_inner(new_mag) / self.mag(); + let scale = flt::Priv::from_inner(new_mag) / self.mag(); self.clone() * scale.into_inner() } } @@ -130,17 +130,17 @@ impl Line { } } - pub fn x(&self, y: f64) -> f64 { + pub fn x(&self, y: flt::Pub) -> flt::Pub { self.shifted(Point::new(0.0, -y)).x_intercept() } - pub fn y(&self, x: f64) -> f64 { + pub fn y(&self, x: flt::Pub) -> flt::Pub { (self.from.y + (real(x) - self.from.x) * self.slope()).into() } - pub fn at_x(&self, x: f64) -> Point { + pub fn at_x(&self, x: flt::Pub) -> Point { Point::new(x, self.y(x)) } - pub fn x_intercept(&self) -> f64 { + pub fn x_intercept(&self) -> flt::Pub { let t = (-self.from.y) / (self.to.y - self.from.y); (self.from.x + t * (self.to.x - self.from.x)).into() } @@ -149,12 +149,12 @@ impl Line { self.to + (-self.from) } - pub fn length_sq(&self) -> f64 { + pub fn length_sq(&self) -> flt::Pub { let Point { x, y } = self.as_vector(); (x*x + y*y).into() } - pub fn distance_sq(&self, pt: Point) -> f64 { + pub fn distance_sq(&self, pt: Point) -> flt::Pub { // source: https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points let d = self.as_vector(); let twice_area = d.y*pt.x - d.x*pt.y + self.to.x*self.from.y - self.to.y*self.from.x; @@ -170,7 +170,7 @@ impl Line { self.to.x > self.from.x } - pub fn contains_x(&self, x: f64) -> bool { + pub fn contains_x(&self, x: flt::Pub) -> bool { let x = real(x); if self.is_ascending_x() { x >= self.from.x && x < self.to.x @@ -179,7 +179,7 @@ impl Line { } } - pub fn contains_y(&self, y: f64) -> bool { + pub fn contains_y(&self, y: flt::Pub) -> bool { let y = real(y); if self.is_ascending_x() { y >= self.from.y && y < self.to.y @@ -188,15 +188,15 @@ impl Line { } } - pub fn min_x(&self) -> f64 { + pub fn min_x(&self) -> flt::Pub { self.from.x.min(self.to.x).into() } - pub fn max_x(&self) -> f64 { + pub fn max_x(&self) -> flt::Pub { self.from.x.max(self.to.x).into() } - pub fn clamp_x(&self, x: f64) -> f64 { + pub fn clamp_x(&self, x: flt::Pub) -> flt::Pub { match x { v if v < self.min_x() => self.min_x(), v if v > self.max_x() => self.max_x(), @@ -204,18 +204,18 @@ impl Line { } } - pub fn clamp_by_x(&self, x: f64) -> Point { + pub fn clamp_by_x(&self, x: flt::Pub) -> Point { self.at_x(self.clamp_x(x)) } - pub fn slope(&self) -> f64 { + pub fn slope(&self) -> flt::Pub { let s = (self.to.y - self.from.y) / (self.to.x - self.from.x); s.into() } /// Move the coordinate along the line toward some desired `y(x)` value. /// The returned value could be OOB: check `self.contains_x` to find out. - pub fn move_toward_y_unclamped(&self, start_x: f64, target_y: f64) -> f64 { + pub fn move_toward_y_unclamped(&self, start_x: flt::Pub, target_y: flt::Pub) -> flt::Pub { let start_y = self.y(start_x); let delta_x = real(target_y - start_y) / self.slope(); (real(start_x) + delta_x).into() @@ -249,19 +249,19 @@ impl Polygon { }) } - pub fn min_x(&self) -> f64 { + pub fn min_x(&self) -> flt::Pub { self.points.iter().map(|p| p.x).min().unwrap().into() } - pub fn max_x(&self) -> f64 { + pub fn max_x(&self) -> flt::Pub { self.points.iter().map(|p| p.x).max().unwrap().into() } - pub fn min_y(&self) -> f64 { + pub fn min_y(&self) -> flt::Pub { self.points.iter().map(|p| p.y).min().unwrap().into() } - pub fn max_y(&self) -> f64 { + pub fn max_y(&self) -> flt::Pub { self.points.iter().map(|p| p.y).max().unwrap().into() } } diff --git a/src/lib.rs b/src/lib.rs index 68f5e81..e45b662 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -25,42 +25,60 @@ pub use sim::*; // mu_0 = vacuum permeability = 1.25663706212(19)×10−6 H/m // c_0 = speed of light in vacuum = 299792458 -#[cfg(debug)] -pub(crate) type F64 = decorum::Real; -#[cfg(not(debug))] -pub(crate) type F64 = decorum::Total; +#[allow(unused)] +pub(crate) mod flt { + use decorum::{R32, R64}; + pub(crate) type T64 = decorum::Total; + pub(crate) type T32 = decorum::Total; + + // Use checked floats for debug and fast floats for release + #[cfg(debug)] + pub(crate) mod defaults { + use super::*; + pub(crate) type Priv = R64; + pub type Pub = f64; + } + #[cfg(not(debug))] + pub(crate) mod defaults { + use super::*; + pub(crate) type Priv = T32; + pub type Pub = f32; + } + pub use defaults::*; +} #[allow(non_snake_case)] pub mod consts { + use super::*; /// Speed of light in a vacuum; m/s. /// Also equal to 1/sqrt(epsilon_0 mu_0) - pub const C: f64 = 299792458.0; - // pub const Z0: f64 = 376.73031366857f32; + pub const C: flt::Pub = 299792458.0; + // pub const Z0: flt::Pub = 376.73031366857f32; /// Vacuum Permeability - pub const MU0: f64 = 1.2566370621219e-6; // H/m + pub const MU0: flt::Pub = 1.2566370621219e-6; // H/m // Vacuum Permittivity - pub const EPS0: f64 = 8.854187812813e-12; // F⋅m−1 + pub const EPS0: flt::Pub = 8.854187812813e-12; // F⋅m−1 pub mod real { - use crate::F64; - pub fn C() -> F64 { + use super::*; + pub(crate) fn C() -> flt::Priv { super::C.into() } - pub fn C2() -> F64 { + pub(crate) fn C2() -> flt::Priv { C() * C() } - pub fn EPS0() -> F64 { + pub(crate) fn EPS0() -> flt::Priv { super::EPS0.into() } - pub fn MU0() -> F64 { + pub(crate) fn MU0() -> flt::Priv { super::MU0.into() } - pub fn ONE() -> F64 { + pub(crate) fn ONE() -> flt::Priv { 1.0.into() } - pub fn TWO() -> F64 { + pub(crate) fn TWO() -> flt::Priv { 2.0.into() } - pub fn HALF() -> F64 { + pub(crate) fn HALF() -> flt::Priv { 0.5.into() } } diff --git a/src/mat.rs b/src/mat.rs index 5b472f1..a048724 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -1,4 +1,4 @@ -use crate::{CellState, consts}; +use crate::{CellState, consts, flt}; use crate::geom::{Line, Point, Polygon}; use enum_dispatch::enum_dispatch; use std::cmp::Ordering; @@ -7,15 +7,15 @@ use std::cmp::Ordering; pub trait Material { /// Return \sigma, the electrical conductivity. /// For a vacuum, this is zero. For a perfect conductor, \inf. - fn conductivity(&self) -> f64 { + fn conductivity(&self) -> flt::Pub { 0.0 } /// Return the magnetization in the z direction (M_z). - fn mz(&self) -> f64 { + fn mz(&self) -> flt::Pub { 0.0 } /// Return the new h_z, and optionally change any internal state (e.g. magnetization). - fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + fn step_h(&mut self, context: &CellState, delta_bz: flt::Pub) -> flt::Pub { // B = mu0*(H+M), and in a default material M=0. context.hz() + delta_bz / consts::MU0 } @@ -23,12 +23,12 @@ pub trait Material { #[derive(Clone, Default)] pub struct Static { - pub conductivity: f64, - pub mz: f64, + pub conductivity: flt::Pub, + pub mz: flt::Pub, } impl Static { - pub fn conductor(conductivity: f64) -> Self { + pub fn conductor(conductivity: flt::Pub) -> Self { Self { conductivity, ..Default::default() @@ -37,10 +37,10 @@ impl Static { } impl Material for Static { - fn conductivity(&self) -> f64 { + fn conductivity(&self) -> flt::Pub { self.conductivity } - fn mz(&self) -> f64 { + fn mz(&self) -> flt::Pub { self.mz } } @@ -48,19 +48,19 @@ impl Material for Static { #[derive(Clone, Default)] pub struct ComsolFerromagnet { /// Instantaneous magnetization in the z direction. - pub mz: f64, + pub mz: flt::Pub, // Ferromagnetic parameters: /// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms). - pub ms: f64, + pub ms: flt::Pub, /// Rate at which M changes w.r.t. H. - pub xi: f64, + pub xi: flt::Pub, } impl Material for ComsolFerromagnet { - fn mz(&self) -> f64 { + fn mz(&self) -> flt::Pub { self.mz } - fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + fn step_h(&mut self, context: &CellState, delta_bz: flt::Pub) -> flt::Pub { // COMSOL ferromagnetic model: https://www.comsol.com/blogs/accessing-external-material-models-for-magnetic-simulations/ // delta(M) = xi * delta*(H) # delta(M) = M(t+1) - M(t) // new_M = clamp(M+delta(M), -Ms, Ms)) # i.e. clamp to saturation @@ -113,7 +113,7 @@ impl MHCurve { /// Moves (h, m) towards some location in the MH curve where H + M = target_hm. /// Returns `Ok((h, m))` if complete; `Err((h, m))` if there's more work to be done (call it /// again). - fn step_toward(&self, h: f64, m: f64, target_hm: f64) -> Result<(f64, f64), (f64, f64)> { + fn step_toward(&self, h: flt::Pub, m: flt::Pub, target_hm: flt::Pub) -> Result<(flt::Pub, flt::Pub), (flt::Pub, flt::Pub)> { let is_ascending = match target_hm.partial_cmp(&(h + m)).unwrap() { Ordering::Greater => true, Ordering::Less => false, @@ -153,7 +153,7 @@ impl MHCurve { Err(active_segment.clamp_by_x(new_h).into()) } } - fn move_to(&self, mut h: f64, mut m: f64, target_hm: f64) -> (f64, f64) { + fn move_to(&self, mut h: flt::Pub, mut m: flt::Pub, target_hm: flt::Pub) -> (flt::Pub, flt::Pub) { loop { match self.step_toward(h, m, target_hm) { Ok((x, y)) => break (x, y), @@ -169,7 +169,8 @@ impl MHCurve { #[derive(Clone)] pub struct PiecewiseLinearFerromagnet { /// Instantaneous magnetization in the z direction. - pub mz: f64, + pub mz: flt::Pub, + pub conductivity: flt::Pub, mh_curve: MHCurve, } @@ -177,23 +178,27 @@ impl PiecewiseLinearFerromagnet { /// Construct a B(H) curve from a sweep from H = 0 to Hs and back down to H = 0. /// The curve below B = 0 is derived by symmetry. /// Points are (H, B). - pub fn from_bh(points: &[(f64, f64)]) -> Self { + pub fn from_bh(points: &[(flt::Pub, flt::Pub)]) -> Self { let mh_points: Vec = points.iter().cloned().map(|(h, b)| { Point::new(h, b / consts::MU0 - h) }).collect(); Self { mz: 0.0, + conductivity: 0.0, mh_curve: MHCurve::new(&*mh_points), } } } impl Material for PiecewiseLinearFerromagnet { - fn mz(&self) -> f64 { + fn mz(&self) -> flt::Pub { self.mz } - fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + fn conductivity(&self) -> flt::Pub { + self.conductivity + } + fn step_h(&mut self, context: &CellState, delta_bz: flt::Pub) -> flt::Pub { let (h, m) = (context.hz(), self.mz); let target_hm = h + m + delta_bz/consts::MU0; @@ -223,7 +228,7 @@ pub mod db { use super::*; /// https://www.ferroxcube.com/upload/media/product/file/MDS/3r1.pdf pub fn ferroxcube_3r1() -> GenericMaterial { - PiecewiseLinearFerromagnet::from_bh(&[ + let mut m = PiecewiseLinearFerromagnet::from_bh(&[ ( 35.0, 0.0), ( 50.0, 0.250), ( 100.0, 0.325), @@ -234,7 +239,9 @@ pub mod db { ( 100.0, 0.345), ( 50.0, 0.340), ( 0.0, 0.325), - ]).into() + ]); + m.conductivity = 1e3; + m.into() } } @@ -258,7 +265,7 @@ mod test { ]) } - fn assert_step_toward_symmetric(h: f64, m: f64, target_mh: f64, target: Result<(f64, f64), (f64, f64)>) { + fn assert_step_toward_symmetric(h: flt::Pub, m: flt::Pub, target_mh: flt::Pub, target: Result<(flt::Pub, flt::Pub), (flt::Pub, flt::Pub)>) { let curve = mh_curve_for_test(); let neg_target = match target { Ok((a, b)) => Ok((-a, -b)), @@ -268,7 +275,7 @@ mod test { assert_eq!(curve.step_toward(-h, -m, -target_mh), neg_target); } - fn assert_move_to_symmetric(h: f64, m: f64, target_mh: f64, target: (f64, f64)) { + fn assert_move_to_symmetric(h: flt::Pub, m: flt::Pub, target_mh: flt::Pub, target: (flt::Pub, flt::Pub)) { let curve = mh_curve_for_test(); assert_eq!(curve.move_to(h, m, target_mh), target); assert_eq!(curve.move_to(-h, -m, -target_mh), (-target.0, -target.1)); diff --git a/src/meas.rs b/src/meas.rs index 8ce6154..9b7346c 100644 --- a/src/meas.rs +++ b/src/meas.rs @@ -9,7 +9,29 @@ pub struct Time; impl AbstractMeasurement for Time { fn eval(&self, state: &SimSnapshot) -> String { - format!("time: {:.3e} (step {})", state.time(), state.step_no()) + format!("{:.3e}s (step {})", state.time(), state.step_no()) + } +} + +pub struct Meta; + +impl AbstractMeasurement for Meta { + fn eval(&self, state: &SimSnapshot) -> String { + format!("{}x{} feat: {:.1e}m", state.width(), state.height(), state.feature_size()) + } +} + +pub struct Label(pub String); + +impl Label { + pub fn new>(s: S) -> Self { + Self(s.into()) + } +} + +impl AbstractMeasurement for Label { + fn eval(&self, _state: &SimSnapshot) -> String { + self.0.clone() } } diff --git a/src/render.rs b/src/render.rs index 40b37f9..f2f236b 100644 --- a/src/render.rs +++ b/src/render.rs @@ -1,10 +1,10 @@ use ansi_term::Color::RGB; use crate::geom::Point; -use crate::{F64, Material as _, SimSnapshot, SimState}; +use crate::{flt, Material as _, SimSnapshot, SimState}; use crate::mat; use crate::sim::Cell; use crate::meas::AbstractMeasurement; -use font8x8::{BASIC_FONTS, UnicodeFonts as _}; +use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _}; use image::{RgbImage, Rgb}; use imageproc::{pixelops, drawing}; use std::fs::File; @@ -14,7 +14,7 @@ use y4m; /// Accept a value from (-\inf, \inf) and return a value in (-1, 1). /// If the input is equal to `typical`, it will be mapped to 0.5. /// If the input is equal to -`typical`, it will be mapped to -0.5. -fn scale_signed(x: f64, typical: f64) -> f64 { +fn scale_signed(x: flt::Pub, typical: flt::Pub) -> flt::Pub { if x >= 0.0 { scale_unsigned(x, typical) } else { @@ -24,29 +24,29 @@ fn scale_signed(x: f64, typical: f64) -> f64 { /// Accept a value from [0, \inf) and return a value in [0, 1). /// If the input is equal to `typical`, it will be mapped to 0.5. -fn scale_unsigned(x: f64, typical: f64) -> f64 { +fn scale_unsigned(x: flt::Pub, typical: flt::Pub) -> flt::Pub { // f(0) => 0 // f(1) => 0.5 // f(\inf) => 1 // f(x) = 1 - 1/(x+1) - let x = F64::from_inner(x); - let typ = F64::from_inner(typical); - let y = F64::from_inner(1.0) - F64::from_inner(1.0)/(x/typ + 1.0); + let x = flt::Priv::from_inner(x); + let typ = flt::Priv::from_inner(typical); + let y = flt::Priv::from_inner(1.0) - flt::Priv::from_inner(1.0)/(x/typ + 1.0); y.into() } -fn scale_signed_to_u8(x: f64, typ: f64) -> u8 { +fn scale_signed_to_u8(x: flt::Pub, typ: flt::Pub) -> u8 { let norm = 128.0 + 128.0*scale_signed(x, typ); norm as _ } -fn scale_unsigned_to_u8(x: f64, typ: f64) -> u8 { +fn scale_unsigned_to_u8(x: flt::Pub, typ: flt::Pub) -> u8 { let norm = 256.0*scale_unsigned(x, typ); norm as _ } /// Scale a vector to have magnitude between [0, 1). -fn scale_vector(x: Point, typical_mag: f64) -> Point { +fn scale_vector(x: Point, typical_mag: flt::Pub) -> Point { let new_mag = scale_unsigned(x.mag(), typical_mag); x.with_mag(new_mag) } @@ -88,7 +88,7 @@ impl<'a> RenderSteps<'a> { } } } - fn render_vector_field) -> Point>(&mut self, color: Rgb, typical: f64, measure: F) { + fn render_vector_field) -> Point>(&mut self, color: Rgb, typical: flt::Pub, measure: F) { let w = self.sim.width(); let h = self.sim.height(); let vec_spacing = 10; @@ -98,7 +98,7 @@ impl<'a> RenderSteps<'a> { let vec = self.field_vector(x, y, vec_spacing, &measure); let norm_vec = scale_vector(vec, typical); let alpha = 0.7*scale_unsigned(vec.mag_sq(), typical * 5.0); - let vec = norm_vec * (vec_spacing as f64); + let vec = norm_vec * (vec_spacing as flt::Pub); let center = Point::new(x as _, y as _) + Point::new(vec_spacing as _, vec_spacing as _)*0.5; self.im.draw_field_arrow(center, vec, color, alpha as f32); } @@ -120,7 +120,9 @@ impl<'a> RenderSteps<'a> { for (meas_no, m) in self.meas.iter().enumerate() { let meas_string = m.eval(self.sim); for (i, c) in meas_string.chars().enumerate() { - let glyph = BASIC_FONTS.get(c).unwrap(); + let glyph = BASIC_FONTS.get(c) + .or_else(|| GREEK_FONTS.get(c)) + .unwrap_or_else(|| BASIC_FONTS.get('?').unwrap()); for (y, bmp) in glyph.iter().enumerate() { for x in 0..8 { if (bmp & 1 << x) != 0 { @@ -155,7 +157,7 @@ impl<'a> RenderSteps<'a> { // avoid division by zero Point::new(0.0, 0.0) } else { - field * (1.0 / ((xw*yw) as f64)) + field * (1.0 / ((xw*yw) as flt::Pub)) } } } diff --git a/src/sim.rs b/src/sim.rs index 94395d1..1228e9d 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,4 +1,4 @@ -use crate::{F64, consts}; +use crate::{flt, consts}; use crate::geom::Point; use crate::mat::{self, GenericMaterial, Material}; @@ -12,12 +12,12 @@ pub type SimSnapshot = SimState; pub struct SimState { cells: Array2>, scratch: Array2>, - feature_size: F64, + feature_size: flt::Priv, step_no: u64, } impl SimState { - pub fn new(width: u32, height: u32, feature_size: f64) -> Self { + pub fn new(width: u32, height: u32, feature_size: flt::Pub) -> Self { Self { cells: Array2::default((height as _, width as _)), scratch: Array2::default((height as _, width as _)), @@ -64,7 +64,7 @@ impl SimState { pub fn snapshot(&self, dec: u32) -> SimSnapshot { let rows = self.height() / dec; let cols = self.width() / dec; - let sample_area = (dec * dec) as f64; + let sample_area = (dec * dec) as flt::Pub; // XXX Zip::par_apply_collect requires the Cell type to be Copy, // so we assign into a pre-allocated type instead. @@ -102,28 +102,28 @@ impl SimState { } impl SimState { - pub fn impulse_bz(&mut self, x: u32, y: u32, bz: f64) { + pub fn impulse_bz(&mut self, x: u32, y: u32, bz: flt::Pub) { self.get_mut(x, y).impulse_bz(bz); } } impl SimState { - pub fn time(&self) -> f64 { - (self.timestep() * self.step_no as f64).into() + pub fn time(&self) -> flt::Pub { + (self.timestep() * self.step_no as flt::Pub).into() } pub fn step_no(&self) -> u64 { self.step_no } - pub fn feature_size(&self) -> f64 { + pub fn feature_size(&self) -> flt::Pub { self.feature_size.into() } - pub fn impulse_ex(&mut self, x: u32, y: u32, ex: f64) { + pub fn impulse_ex(&mut self, x: u32, y: u32, ex: flt::Pub) { self.get_mut(x, y).state.ex += ex; } - pub fn impulse_ey(&mut self, x: u32, y: u32, ey: f64) { + pub fn impulse_ey(&mut self, x: u32, y: u32, ey: flt::Pub) { self.get_mut(x, y).state.ey += ey; } @@ -140,7 +140,7 @@ impl SimState { &mut self.cells[[y as _, x as _]] } - fn timestep(&self) -> F64 { + fn timestep(&self) -> flt::Priv { self.feature_size / consts::real::C() } } @@ -170,16 +170,16 @@ pub struct Cell { } impl Cell { - pub fn ex(&self) -> f64 { + pub fn ex(&self) -> flt::Pub { self.state.ex() } - pub fn ey(&self) -> f64 { + pub fn ey(&self) -> flt::Pub { self.state.ey() } pub fn e(&self) -> Point { Point::new(self.ex(), self.ey()) } - pub fn hz(&self) -> f64 { + pub fn hz(&self) -> flt::Pub { self.state.hz() } pub fn mat(&self) -> &M { @@ -191,7 +191,7 @@ impl Cell { } impl Cell { - pub fn bz(&self) -> f64 { + pub fn bz(&self) -> flt::Pub { consts::MU0 * (self.hz() + self.mat.mz()) } @@ -200,16 +200,16 @@ impl Cell { Point::new(self.ex()*conductivity, self.ey()*conductivity) } - fn bz_to_hz(&self, bz: f64) -> f64 { + fn bz_to_hz(&self, bz: flt::Pub) -> flt::Pub { // B = mu0*(H + M) => H = B/mu0 - M bz/consts::MU0 - self.mat.mz() } - fn impulse_bz(&mut self, delta_bz: f64) { + fn impulse_bz(&mut self, delta_bz: flt::Pub) { self.state.hz = self.mat.step_h(&self.state, delta_bz).into(); } - fn step_h(mut self, right: &Self, down: &Self, _delta_t: F64) -> Self { + fn step_h(mut self, right: &Self, down: &Self, _delta_t: flt::Priv) -> Self { // ```tex // Maxwell-Faraday equation: $\nabla x E = -dB/dt$ // Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$ @@ -235,7 +235,7 @@ impl Cell { /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation /// feature_size = how many units apart is the center of each adjacent cell on the grid. - fn step_e(self, left: &Self, up: &Self, delta_t: F64, feature_size: F64) -> Self { + fn step_e(self, left: &Self, up: &Self, delta_t: flt::Priv, feature_size: flt::Priv) -> Self { // ```tex // Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"): // $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity") @@ -262,7 +262,7 @@ impl Cell { use consts::real::{EPS0, ONE, TWO}; - let sigma: F64 = self.mat.conductivity().into(); + let sigma: flt::Priv = self.mat.conductivity().into(); // XXX not obvious that bz_to_hz is sensible // let delta_hz_y = self.hz() - self.bz_to_hz(up.bz()); @@ -288,19 +288,19 @@ impl Cell { #[derive(Copy, Clone, Default)] pub struct CellState { - ex: F64, - ey: F64, - hz: F64, + ex: flt::Priv, + ey: flt::Priv, + hz: flt::Priv, } impl CellState { - pub fn ex(&self) -> f64 { + pub fn ex(&self) -> flt::Pub { self.ex.into() } - pub fn ey(&self) -> f64 { + pub fn ey(&self) -> flt::Pub { self.ey.into() } - pub fn hz(&self) -> f64 { + pub fn hz(&self) -> flt::Pub { self.hz.into() } }