Implement time stepping of both the B and E field

examples/coremem.rs

@@ -1,3 +1,5 @@
+use coremem::SimState;
 fn main() {
-    println!("loaded");
+    let mut state = SimState::default();
+    state.step();
 }

src/lib.rs

@@ -1,7 +1,102 @@
-#[cfg(test)]
-mod tests {
-    #[test]
-    fn it_works() {
-        assert_eq!(2 + 2, 4);
+//! Magnetic core memory simulator.
+//! Built using a Finite-Difference Time-Domain electromagnetics simulator.
+//! The exhaustive guide for implementing a FDTD simulator by John B. Schneider [1] gives some
+//! overview of the theory.
+//!
+//! [1] https://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf
+
+// Some things to keep in mind:
+// B = mu_r*H + M
+// For a vacuum, B = H
+//
+// mu_0 = vacuum permeability = 1.25663706212(19)×10−6 H/m
+// c_0 = speed of light in vacuum = 299792458
+
+mod consts {
+    /// Speed of light in a vacuum; m/s.
+    /// Also equal to 1/sqrt(epsilon_0 mu_0)
+    pub const C:f32 = 299792458f32;
+}
+
+#[derive(Default)]
+pub struct SimState {
+    cells: Vec<Cell>,
+}
+
+impl SimState {
+    pub fn step(&mut self) {
+        let mut working_cells = self.cells.clone();
+        // first advance all the magnetic fields
+        for (i, left_cell) in self.cells.iter().enumerate() {
+            let right_cell = match self.cells.get(i+1) {
+                Some(&cell) => cell,
+                _ => Cell::default(),
+            };
+            working_cells[i] = left_cell.step_b(right_cell);
+        }
+
+        for (i, right_cell) in working_cells.iter().enumerate() {
+            let left_cell = match i {
+                0 => Cell::default(),
+                _ => working_cells[i-1],
+            };
+            self.cells[i] = right_cell.step_e(left_cell);
+        }
+    }
+    pub fn get(&self) -> &[Cell] {
+        &*self.cells
+    }
+}
+
+/// Conceptually, one cell looks like this:
+///
+/// +-------------+
+/// |             |
+/// .Ez    .By    .  next cell
+/// |             |
+/// +-------------+
+///
+/// Where +By points up and +Ez points into the page, and the cell has a unit length of 1.
+#[derive(Copy, Clone, Default)]
+pub struct Cell {
+    /// electric field
+    ez: f32,
+    /// magnetic field
+    by: f32,
+}
+
+impl Cell {
+    fn step_b(self, right: Cell) -> Self {
+        // Maxwell's equation: del x E = -dB/dt
+        // Expands: dB_y/dt = dE_z/dx
+        // Discretize: (delta B_y) / (delta t) = (delta E_z)/(delta x)
+        // Rearrange: delta B_y = (delta t)/(delta x) * (delta E_z)
+        //
+        // light travels C meters per second, therefore (delta x)/(delta t) = c if we use SI units.
+        // XXX: this differs from [1], which says we should use Z_0 = 1/(epsilon c) here instead of c.
+        let delta_e = right.ez - self.ez;  //< delta E_z
+        let delta_b = delta_e / consts::C;  //< delta B_y
+        Cell {
+            ez: self.ez,
+            by: self.by + delta_b,
+        }
+    }
+
+    fn step_e(self, left: Cell) -> Self {
+        // Maxwell's equation: del x B = mu_0 eps_0 dE/dt
+        // Expands: -dB_z/dx = mu_0 eps_0 dE_y/dt
+        // Rearrange: dE_y/dt = -1/(mu_0 eps_0) dB_z/dx
+        // Discretize: (delta E_y)/(delta t) = -1/(mu_0 eps_0) (delta dB_z)/(delta x)
+        // Rearrange: delta E_y = (delta t)/(delta x) 1/(mu_0 eps_0) (-delta B_z)
+        // Substitute c as in step_b: delta E_y = (mu_0 eps_0)/c (-delta B_z)
+        // Note that c = 1/sqrt(mu_0 eps_0), so this becomes:
+        // delta E_y = c (-delta B_z)
+        // XXX once again this diffes from [1]
+        let delta_b = self.by - left.by;
+        let delta_e = (-delta_b) * consts::C;
+        Cell {
+            ez: self.ez + delta_e,
+            by: self.by,
+        }
     }
 }