Fix COMSOL ferromagnet model, and make a new example for it

I'm a little skeptical of it as a model, though.

examples/ferromagnet.rs

@@ -0,0 +1,43 @@
+use coremem::SimState;
+use coremem::render::ColorTermRenderer as Renderer;
+use coremem::consts;
+use std::{thread, time};
+
+fn main() {
+    let width = 201;
+    let mut state = SimState::new(width, 101, 1e-3 /* feature size */);
+
+    for y in 0..100 {
+        for x in 50..60 {
+            state.get_mut(x, y).mat_mut().conductivity = 10.0;
+        }
+    }
+    for y in 40..60 {
+        for x in 62..70 {
+            state.get_mut(x, y).mat_mut().xi = 150.0;
+            state.get_mut(x, y).mat_mut().ms = 2.0;
+        }
+    }
+
+    let mut step = 0u64;
+    loop {
+        step += 1;
+        let imp = if step < 50 { 
+            25.0 * ((step as f64)*0.02*std::f64::consts::PI).sin()
+        } else {
+            0.0
+        };
+        // state.impulse_ex(50, 50, imp);
+        // state.impulse_ey(50, 50, imp);
+        // state.impulse_bz(20, 20, (imp / 3.0e8) as _);
+        // state.impulse_bz(80, 20, (imp / 3.0e8) as _);
+        for y in 0..100 {
+            for x in 52..58 {
+                state.impulse_ey(x, y, imp as _);
+            }
+        }
+        Renderer.render(&state);
+        state.step();
+        thread::sleep(time::Duration::from_millis(67));
+    }
+}

src/lib.rs

@@ -317,11 +317,13 @@ impl Material for GenericMaterial {
         // let unsaturated_new_m = self.mz + self.xi * unsaturated_delta_h;
         let unsaturated_new_m = if self.mz == self.ms {
             // positively saturated
-            let delta_h = expected_new_h.max(0.0) - context.hz().max(0.0);
+            let delta_h = expected_new_h.min(0.0) - context.hz().min(0.0);
+            // assert!(self.ms == 0.0, "pos sat");
             self.mz + self.xi * delta_h
         } else if self.mz == -self.ms {
             // negatively saturated
-            let delta_h = expected_new_h.min(0.0) - context.hz().min(0.0);
+            let delta_h = expected_new_h.max(0.0) - context.hz().max(0.0);
+            //assert!(self.ms == 0.0, "neg sat");
             self.mz + self.xi * delta_h
         } else {
             // not saturated

src/render.rs

@@ -51,13 +51,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*50.0);
+                let r = norm_color(cell.mat.mz);
                 let b = (55.0*cell.mat().conductivity).min(255.0) as u8;
                 //let b = 0;
                 //let g = norm_color(cell.ex());
                 //let b = norm_color(cell.ey());
                 //let g = norm_color(curl(cell.ex(), cell.ey()));
-                let g = norm_color((cell.bz() * 3.0e8).into());
+                //let g = norm_color((cell.bz() * 3.0e8).into());
+                let g = norm_color(cell.ey().into());
                 write!(&mut buf, "{}", RGB(r, g, b).paint(square));
             }
             write!(&mut buf, "\n");