diff --git a/examples/em_reflection.rs b/examples/em_reflection.rs index 5c5fdcd..11be6a4 100644 --- a/examples/em_reflection.rs +++ b/examples/em_reflection.rs @@ -29,7 +29,7 @@ fn main() { loop { step += 1; let imp = if step < 50 { - 2.5 * ((step as f64)*0.04*std::f64::consts::PI).sin() + 250000.0 * ((step as f64)*0.04*std::f64::consts::PI).sin() } else { 0.0 }; diff --git a/examples/ferromagnet.rs b/examples/ferromagnet.rs index 75ffb8b..6a5e23b 100644 --- a/examples/ferromagnet.rs +++ b/examples/ferromagnet.rs @@ -13,11 +13,23 @@ fn main() { } for y in 40..60 { for x in 62..70 { - *state.get_mut(x, y).mat_mut() = mat::ComsolFerromagnet { - xi: 150.0, - ms: 2.0, - ..mat::ComsolFerromagnet::default() - }.into(); + *state.get_mut(x, y).mat_mut() = mat::PiecewiseLinearFerromagnet::from_bh(&[ + ( 35.0, 0.0), + ( 50.0, 0.250), + ( 100.0, 0.325), + ( 200.0, 0.350), + (1000.0, 0.390), + // Falling + ( 200.0, 0.360), + ( 100.0, 0.345), + ( 50.0, 0.340), + ( 0.0, 0.325), + ]).into(); + //*state.get_mut(x, y).mat_mut() = mat::ComsolFerromagnet { + // xi: 150.0, + // ms: 2.0, + // ..mat::ComsolFerromagnet::default() + //}.into(); } } @@ -25,7 +37,7 @@ fn main() { loop { step += 1; let imp = if step < 50 { - 25.0 * ((step as f64)*0.02*std::f64::consts::PI).sin() + 250000.0 * ((step as f64)*0.02*std::f64::consts::PI).sin() } else { 0.0 }; diff --git a/src/geom.rs b/src/geom.rs new file mode 100644 index 0000000..a2a1a60 --- /dev/null +++ b/src/geom.rs @@ -0,0 +1,174 @@ +use decorum::R64; +use std::ops::{Add, Neg}; + +#[derive(Copy, Clone, Debug, Default)] +pub struct Point { + pub x: f64, + pub y: f64, +} + +impl Add for Point { + type Output = Point; + fn add(self, other: Point) -> Point { + Point { + x: self.x + other.x, + y: self.y + other.y, + } + } +} + +impl Neg for Point { + type Output = Point; + fn neg(self) -> Point { + Point { + x: -self.x, + y: -self.y, + } + } +} + +impl Point { + pub fn new(x: f64, y: f64) -> Self { + Self { + x, + y, + } + } +} + +#[derive(Copy, Clone, Debug, Default)] +pub struct Line { + from: Point, + to: Point, +} + +impl Add for Line { + type Output = Line; + fn add(self, other: Line) -> Line { + Line { + from: self.from + Point::new(0.0, other.y(self.from.x)), + to: self.to + Point::new(0.0, other.y(self.to.x)), + } + } +} + +impl Line { + pub fn new(from: Point, to: Point) -> Self { + Self { + from, + to, + } + } + + pub fn x(&self, y: f64) -> f64 { + self.shifted(Point::new(0.0, -y)).x_intercept() + } + pub fn y(&self, x: f64) -> f64 { + self.from.y + (x - self.from.x) * self.slope() + } + pub fn at_x(&self, x: f64) -> Point { + Point::new(x, self.y(x)) + } + + pub fn x_intercept(&self) -> f64 { + let t = (0.0 - self.from.y) / (self.to.y - self.from.y); + self.from.x + t * (self.to.x - self.from.x) + } + + pub fn is_ascending(&self) -> bool { + self.to.y > self.from.y + } + + pub fn is_ascending_x(&self) -> bool { + self.to.x > self.from.x + } + + pub fn contains_x(&self, x: f64) -> bool { + if self.is_ascending_x() { + x >= self.from.x && x < self.to.x + } else { + x >= self.to.x && x < self.from.x + } + } + + pub fn contains_y(&self, y: f64) -> bool { + if self.is_ascending_x() { + y >= self.from.y && y < self.to.y + } else { + y >= self.to.y && y < self.from.y + } + } + + pub fn clamp_x(&self, x: f64) -> f64 { + match x { + v if v < self.from.x => self.from.x, + v if v > self.to.x => self.to.x, + v => v, + } + } + + pub fn clamp_by_x(&self, x: f64) -> Point { + self.at_x(self.clamp_x(x)) + } + + pub fn slope(&self) -> f64 { + (self.to.y - self.from.y) / (self.to.x - self.from.x) + } + + /// 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 { + let start_y = self.y(start_x); + let delta_x = (target_y - start_y) / self.slope(); + start_x + delta_x + } + + pub fn move_toward_y(&self, start_x: f64, target_y: f64) -> f64 { + let x = self.move_toward_y(start_x, target_y); + self.clamp_x(x) + } + + pub fn shifted(&self, p: Point) -> Line { + Line { + from: self.from + p, + to: self.to + p + } + } +} + +#[derive(Clone, Debug)] +pub struct Polygon { + points: Vec, +} + +impl Polygon { + pub fn new(points: Vec) -> Self { + Self { + points + } + } + + pub fn segments<'a>(&'a self) -> impl Iterator + 'a { + (0..self.points.len()).into_iter().map(move |i| { + let from = self.points[i]; + let to = *self.points.get(i+1).unwrap_or_else(|| &self.points[0]); + Line { from, to } + }) + } + + pub fn min_x(&self) -> f64 { + self.points.iter().map(|p| R64::from_inner(p.x)).min().unwrap().into() + } + + pub fn max_x(&self) -> f64 { + self.points.iter().map(|p| R64::from_inner(p.x)).max().unwrap().into() + } + + pub fn min_y(&self) -> f64 { + self.points.iter().map(|p| R64::from_inner(p.y)).min().unwrap().into() + } + + pub fn max_y(&self) -> f64 { + self.points.iter().map(|p| R64::from_inner(p.y)).max().unwrap().into() + } +} diff --git a/src/lib.rs b/src/lib.rs index 532dfeb..de1d1d4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -7,6 +7,7 @@ use decorum::R64; +pub mod geom; pub mod mat; pub mod render; pub mod sim; diff --git a/src/mat.rs b/src/mat.rs index 1860f2a..6260e37 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -1,6 +1,6 @@ use crate::{CellState, consts}; +use crate::geom::{Line, Point, Polygon}; use enum_dispatch::enum_dispatch; -use piecewise_linear::PiecewiseLinearFunction; #[enum_dispatch] pub trait Material { @@ -80,67 +80,108 @@ impl Material for ComsolFerromagnet { } } -pub struct BHCurve { - /// describes the curve as B rises from 0. - rising: PiecewiseLinearFunction, - /// describes the curve as B falls towards 0. - falling: PiecewiseLinearFunction, +/// M as a function of H +#[derive(Clone)] +struct MHCurve { + geom: Polygon, } -impl BHCurve { - fn b_bounds_for_h(&self, h: f64) -> (f64, f64) { - let (rise_start, rise_end) = self.rising.domain(); - let bottom = match h { - x if (..-rise_end).contains(&x) => { - // negative saturation - -self.falling.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end) - }, - x if (-rise_end..rise_start).contains(&x) => { - // negative - -self.falling.y_at_x(-h).unwrap() +impl MHCurve { + /// Construct a M(H) curve from a sweep from M = 0 to Ms and back down to M = 0. + /// The curve below M = 0 is derived by symmetry. + pub fn new(points: &[Point]) -> Self { + let full_pts: Vec = points.iter().cloned().chain(points.iter().cloned().map(|p| -p)).collect(); + + Self { + geom: Polygon::new(full_pts) + } + } + /// 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 { + let is_ascending = target_hm > h + m; + if (is_ascending && m == self.geom.max_y()) || (!is_ascending && m == -self.geom.min_y()) { + // Fully saturated. m is fixed, while h moves freely + return Ok(Point::new(target_hm - m, m)); + } + // Locate the segment which would contain the current point + let mut segments = self.geom.segments(); + let active_segment = loop { + let line = segments.next().unwrap_or_else(|| { + panic!("failed to find segment for h:{}, m:{}, {:?}", h, m, self.geom.segments().collect::>()); + }); + if line.contains_y(m) && line.is_ascending() == is_ascending { + if line.contains_x(h) { + break line; + } else { + // need to move the point toward this line + let h_intercept = line.x(m); + break Line::new(Point::new(h, m), Point::new(h_intercept, m)); + } } - x if (rise_start..rise_end).contains(&x) => { - // positive - self.rising.y_at_x(h).unwrap() - }, - x if (rise_end..).contains(&x) => { - // positive saturation - self.rising.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end) + }; + + // Find some m(h) on the active_segment such that sum(h) = h + m(h) = target_hm + let sum_h = active_segment + Line::new(Point::new(0.0, 0.0), Point::new(1.0, 1.0)); + let new_h = sum_h.move_toward_y_unclamped(h, target_hm); + + if sum_h.contains_x(new_h) { + // the segment contains a point with the target H+M + Ok(active_segment.at_x(new_h)) + } else { + // the segment doesn't contain the desired point: clamp and try the next segment + Err(active_segment.clamp_by_x(new_h)) + } + } + fn move_to(&self, mut h: f64, mut m: f64, target_hm: f64) -> (f64, f64) { + loop { + match self.step_toward(h, m, target_hm) { + Ok(p) => break (p.x, p.y), + Err(Point{ x: better_h, y: better_m}) => { + h = better_h; + m = better_m; + }, } - _ => unreachable!("invalid h: {}", h), - }; - let top = match h { - x if (..-rise_end).contains(&x) => { - // negative saturation - -self.rising.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end) - }, - x if (-rise_end..-rise_start).contains(&x) => { - // negative - -self.rising.y_at_x(-h).unwrap() - }, - x if (-rise_start..rise_end).contains(&x) => { - // positive - self.falling.y_at_x(h).unwrap() - }, - x if (rise_end..).contains(&x) => { - // positive saturation - self.falling.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end) - }, - _ => unreachable!("invalid h: {}", h), - }; - (bottom, top) + } } } -pub struct BHFerromagnet { +#[derive(Clone)] +pub struct PiecewiseLinearFerromagnet { /// Instantaneous magnetization in the z direction. pub mz: f64, - bh_curve: BHCurve, + mh_curve: MHCurve, } -impl Material for BHFerromagnet { +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 { + let mh_points: Vec = points.iter().cloned().map(|(h, b)| { + Point::new(h, b / consts::MU0 - h) + }).collect(); + + Self { + mz: 0.0, + mh_curve: MHCurve::new(&*mh_points), + } + } +} + +impl Material for PiecewiseLinearFerromagnet { + fn mz(&self) -> f64 { + self.mz + } fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { - unimplemented!() + let (h, m) = (context.hz(), self.mz); + let target_hm = h + m + delta_bz/consts::MU0; + + let (h, m) = self.mh_curve.move_to(h, m, target_hm); + + self.mz = m; + h } } @@ -148,7 +189,8 @@ impl Material for BHFerromagnet { #[derive(Clone)] pub enum GenericMaterial { Conductor(Conductor), - ComsolFerromagnet(ComsolFerromagnet), + // ComsolFerromagnet(ComsolFerromagnet), + PiecewiseLinearFerromagnet(PiecewiseLinearFerromagnet) } impl Default for GenericMaterial { diff --git a/src/render.rs b/src/render.rs index 6c2c62c..7c5e277 100644 --- a/src/render.rs +++ b/src/render.rs @@ -50,14 +50,14 @@ impl ColorTermRenderer { let cell = state.get(x, y); //let r = norm_color(cell.bz() * consts::C); //let r = 0; - let r = norm_color(cell.mat().mz()); + let r = norm_color(cell.mat().mz()*1.0e-2); let b = (55.0*cell.mat().conductivity()).min(255.0) as u8; //let b = 0; //let b = norm_color(cell.ey()); //let g = 0; //let g = norm_color(cell.ex()); //let g = norm_color(curl(cell.ex(), cell.ey())); - let g = norm_color((cell.bz() * 3.0e8).into()); + let g = norm_color((cell.bz() * 1.0e4).into()); //let g = norm_color(cell.ey().into()); write!(&mut buf, "{}", RGB(r, g, b).paint(square)); }