Store magnetic state as H_z instead of B_z

This is consistent with the Macroscopic version of Maxwell's equations,
and will allow for storing Magnetic state with less trouble than if we
used B.

src/lib.rs

@@ -79,7 +79,7 @@ impl SimState {
                 let cell = self.get(right_x-1, down_y-1);
                 let right_cell = self.get(right_x, down_y-1);
                 let down_cell = self.get(right_x-1, down_y);
-                working_cells[[down_y-1, right_x-1]] = cell.step_b(right_cell, down_cell, half_time_step.into());
+                working_cells[[down_y-1, right_x-1]] = cell.step_h(right_cell, down_cell, half_time_step.into());
             }
         }
         std::mem::swap(&mut working_cells, &mut self.cells);
@@ -98,13 +98,13 @@ impl SimState {
     }
 
     pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f64) {
-        self.cells[[y, x]].ex += ex;
+        self.cells[[y, x]].state.ex += ex;
     }
     pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f64) {
-        self.cells[[y, x]].ey += ey;
+        self.cells[[y, x]].state.ey += ey;
     }
     pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f64) {
-        self.cells[[y, x]].bz += bz;
+        self.cells[[y, x]].impulse_bz(bz);
     }
 
     pub fn width(&self) -> usize {
@@ -145,21 +145,19 @@ impl SimState {
 /// the pluses.
 #[derive(Copy, Clone, Default)]
 pub struct Cell<M> {
-    ex: R64,
+    state: CellState,
-    ey: R64,
-    bz: R64,
     mat: M,
 }
 
 impl<M> Cell<M> {
     pub fn ex(&self) -> f64 {
-        self.ex.into()
+        self.state.ex()
     }
     pub fn ey(&self) -> f64 {
-        self.ey.into()
+        self.state.ey()
     }
-    pub fn bz(&self) -> f64 {
+    pub fn hz(&self) -> f64 {
-        self.bz.into()
+        self.state.hz()
     }
     pub fn mat(&self) -> &M {
         &self.mat
@@ -170,7 +168,15 @@ impl<M> Cell<M> {
 }
 
 impl<M: Material + Clone> Cell<M> {
-    fn step_b(self, right: Self, down: Self, _delta_t: R64) -> Self {
+    pub fn bz(&self) -> f64 {
+        consts::MU0 * (self.hz() + self.mat.mz())
+    }
+
+    fn impulse_bz(&mut self, delta_bz: f64) {
+        self.state.hz = self.mat.step_h(&self.state, delta_bz).into();
+    }
+
+    fn step_h(mut self, right: Self, down: Self, _delta_t: R64) -> Self {
         // ```tex
         // Maxwell's equation: $\nabla x E = -dB/dt$
         // Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$
@@ -181,13 +187,15 @@ impl<M: Material + Clone> Cell<M> {
         //   therefore $(\Delta x)/(2 \Delta t) = (\Delta y)/(2 \Delta t) = c$ if we use SI units.
         // Simplify: $\Delta B_z = (\Delta E_x)/(2c) - (\Delta E_y)/(2c)$
         // ```
-        let delta_ex = down.ex - self.ex;
+        let delta_ex = down.ex() - self.ex();
-        let delta_ey = right.ey - self.ey;
+        let delta_ey = right.ey() - self.ey();
         let delta_bz = (delta_ex - delta_ey) / (2.0*consts::C);
+        let hz = self.mat.step_h(&self.state, delta_bz);
         Cell {
-            ex: self.ex,
+            state: CellState {
-            ey: self.ey,
+                hz: hz.into(),
-            bz: self.bz + delta_bz,
+                ..self.state
+            },
             mat: self.mat,
         }
     }
@@ -222,29 +230,58 @@ impl<M: Material + Clone> Cell<M> {
 
         let sigma = self.mat.conductivity();
 
-        let delta_bz_y = self.bz - up.bz;
+        let delta_bz_y = self.bz() - up.bz();
-        let ex_rhs = self.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size;
+        let ex_rhs = self.state.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size;
         let ex_next = ex_rhs / (ONE() + C2()*MU0()*sigma*delta_t);
 
-        let delta_bz_x = self.bz - left.bz;
+        let delta_bz_x = self.bz() - left.bz();
-        let ey_rhs = self.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size;
+        let ey_rhs = self.state.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size;
         let ey_next = ey_rhs / (ONE() + C2()*MU0()*sigma*delta_t);
 
         Cell {
-            ex: ex_next,
+            state: CellState {
-            ey: ey_next,
+                ex: ex_next,
-            bz: self.bz,
+                ey: ey_next,
+                hz: self.state.hz,
+            },
             mat: self.mat,
         }
     }
 }
 
+#[derive(Copy, Clone, Default)]
+pub struct CellState {
+    ex: R64,
+    ey: R64,
+    hz: R64,
+}
+
+impl CellState {
+    pub fn ex(&self) -> f64 {
+        self.ex.into()
+    }
+    pub fn ey(&self) -> f64 {
+        self.ey.into()
+    }
+    pub fn hz(&self) -> f64 {
+        self.hz.into()
+    }
+}
 
 pub trait Material {
     /// Return \sigma, the electrical conductivity.
     /// For a vacuum, this is zero. For a perfect conductor, \inf.
     fn conductivity(&self) -> f64 {
-        0.0.into()
+        0.0
+    }
+    /// Return the magnetization in the z direction (M_z).
+    fn mz(&self) -> f64 {
+        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 {
+        // B = mu0*(H+M), and in a default material M=0.
+        context.hz() + delta_bz / consts::MU0
     }
 }