Implement a BHCurve

Can't implement this style of ferromagnet until ew rework step_e to use
the macroscopic version of Maxwell's equation though, I think.

Cargo.toml

@@ -11,3 +11,4 @@ ansi_term = "0.12"
 decorum = "0.3"
 enum_dispatch = "0.3"
 ndarray = "0.13"
+piecewise-linear = "0.1"

src/lib.rs

@@ -271,6 +271,7 @@ impl CellState {
 pub mod mat {
     use super::{CellState, consts};
     use enum_dispatch::enum_dispatch;
+    use piecewise_linear::PiecewiseLinearFunction;
 
     #[enum_dispatch]
     pub trait Material {
@@ -350,6 +351,70 @@ pub mod mat {
         }
     }
 
+    pub struct BHCurve {
+        /// describes the curve as B rises from 0.
+        rising: PiecewiseLinearFunction<f64>,
+        /// describes the curve as B falls towards 0.
+        falling: PiecewiseLinearFunction<f64>,
+    }
+
+    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()
+                }
+                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)
+                }
+                _ => 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 {
+        /// Instantaneous magnetization in the z direction.
+        pub mz: f64,
+        bh_curve: BHCurve,
+    }
+
+    impl Material for BHFerromagnet {
+        fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
+            unimplemented!()
+        }
+    }
+
     #[enum_dispatch(Material)]
     #[derive(Clone)]
     pub enum GenericMaterial {