Fix the math errors so that 2d simulations appear to now conserve energy.

examples/coremem.rs

@@ -4,13 +4,13 @@ use coremem::consts;
 use std::{thread, time};
 
 fn main() {
-    let mut state = SimState::new(16, 4);
+    let mut state = SimState::new(16, 5);
     state.impulse_bz(0, 0, 1.0/consts::C);
-    //state.impulse_ey(0, 0, 1.0);
+    state.impulse_ey(0, 0, 1.0);
     state.impulse_bz(7, 0, 1.0/consts::C);
-    //state.impulse_ey(7, 0, 1.0);
+    state.impulse_ey(7, 2, 1.0);
     state.impulse_bz(15, 0, 1.0/consts::C);
-    //state.impulse_ey(15, 0, 1.0);
+    state.impulse_ey(15, 0, 1.0);
     loop {
         NumericTermRenderer.render(&state);
         state.step();

src/lib.rs

@@ -21,6 +21,10 @@ pub mod consts {
     /// Also equal to 1/sqrt(epsilon_0 mu_0)
     pub const C: f32 = 299792458f32;
     // pub const Z0: f32 = 376.73031366857f32;
+    // Vacuum Permeability
+    // pub const Mu0: f32 = 1.2566370621219e-6; // H/m
+    // Vacuum Permittivity
+    // pub const Eps0: f32 = 8.854187812813e-12 // F⋅m−1
 }
 
 #[derive(Default)]
@@ -36,8 +40,7 @@ impl SimState {
     }
 
     pub fn step(&mut self) {
-        let mut working_cells = self.cells.clone();
+        let mut working_cells = Array2::default((self.height(), self.width()));
-        //let mut working_cells = vec![Cell::default(); self.cells.len()];
         // first advance all the magnetic fields
         for down_y in 1..self.height() {
             for right_x in 1..self.width() {
@@ -94,7 +97,7 @@ impl SimState {
 /// |               |
 /// +-------.-------+
 ///
-/// Where the right hand rule indicates that positive Bz is pointing out of the page, towards the
+/// Where the right hand rule indicates that positive Bz is pointing into the page, away from the
 /// reader.
 ///
 /// The dot on bottom is Ex of the cell at (x, y+1) and the dot on the right is the Ey of the cell at
@@ -123,11 +126,12 @@ impl Cell {
         // Rearrange: dB_z/dt = dE_x/dy - dE_y/dx
         // Discretize: (delta B_z)/(delta t) = (delta E_x)/(delta y) - (delta E_y)/(delta x)
         // 
-        // light travels C meters per second, therefore (delta x)/(delta t) = (delta y)/(delta t) = c if we use SI units.
+        // light travels C meters per second, and it we're only stepping half a timestep here (the other half when we advance E),
-        // Simplify: delta B_z = (delta E_x)/c - (delta E_y)/c
+        //   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_ey = right.ey - self.ey;
-        let delta_bz = (delta_ex - delta_ey) / consts::C;
+        let delta_bz = (delta_ex - delta_ey) / (2f32*consts::C);
         Cell {
             ex: self.ex,
             ey: self.ey,
@@ -142,15 +146,15 @@ impl Cell {
         // Expand: dE_x/dt = c^2 dB_z/dy   (1);  dE_y/dt = -c^2 dB_z/dx   (2)
         //
         // Discretize (1): (delta E_x)/(delta t) = c^2 (delta B_z)/(delta y)
-        // Recall: (delta y)/(delta t) = c, as from step_b
+        // Recall: (delta y)/(delta t) = 2c, as from step_b
-        // Substitute: (delta E_x) = c (delta B_z,y)
+        // Substitute: (delta E_x) = c/2 (delta B_z,y)
         //
         // Discretize (2): (delta E_y)/(delta t) = -c^2 (delta B_z)/(delta x)
-        // Substitute c: (delta E_y) = -c (delta B_z,x)
+        // Substitute c: (delta E_y) = -c/2 (delta B_z,x)
         let delta_bz_y = self.bz - up.bz;
-        let delta_ex = consts::C * delta_bz_y;
+        let delta_ex = (0.5f32*consts::C) * delta_bz_y;
         let delta_bz_x = self.bz - left.bz;
-        let delta_ey = consts::C * delta_bz_x;
+        let delta_ey = -(0.5f32*consts::C) * delta_bz_x;
         Cell {
             ex: self.ex + delta_ex,
             ey: self.ey + delta_ey,