Add [broken] conductivity support

It doesn't seem to be solved by using f64

examples/coremem.rs

@@ -5,13 +5,18 @@ use std::{thread, time};
 
 fn main() {
     let mut state = SimState::new(101, 101);
+
+    for x in 0..100 {
+        state.get_mut(x, 70).mat_mut().conductivity = 0.00000001;
+    }
+
     let mut step = 0u64;
     loop {
         step += 1;
-        let imp = 50.0 * ((step as f64)*0.05).sin() as f32;
+        let imp = 50.0 * ((step as f64)*0.05).sin();
         // state.impulse_ex(50, 50, imp);
         // state.impulse_ey(50, 50, imp);
-        state.impulse_bz(50, 50, imp / 3e8f32);
+        state.impulse_bz(50, 50, (imp / 3.0e8) as _);
         Renderer.render(&state);
         state.step();
         thread::sleep(time::Duration::from_millis(33));

src/lib.rs

@@ -21,25 +21,27 @@ 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
+    /// Vacuum Permeability
-    // pub const Mu0: f32 = 1.2566370621219e-6; // H/m
+    pub const MU0: f32 = 1.2566370621219e-6; // H/m
     // Vacuum Permittivity
     // pub const Eps0: f32 = 8.854187812813e-12 // F⋅m−1
 }
 
 #[derive(Default)]
 pub struct SimState {
-    cells: Array2<Cell>,
+    cells: Array2<Cell<GenericMaterial>>,
 }
 
 impl SimState {
     pub fn new(width: usize, height: usize) -> Self {
         Self {
-            cells: Array2::default((height, width))
+            cells: Array2::default((height, width)),
         }
     }
 
     pub fn step(&mut self) {
+        // feature size: 1mm.
+        let half_time_step = 0.0005 * consts::C;
         let mut working_cells = Array2::default((self.height(), self.width()));
         // first advance all the magnetic fields
         for down_y in 1..self.height() {
@@ -47,7 +49,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);
+                working_cells[[down_y-1, right_x-1]] = cell.step_b(right_cell, down_cell, half_time_step);
             }
         }
         std::mem::swap(&mut working_cells, &mut self.cells);
@@ -58,7 +60,7 @@ impl SimState {
                 let cell = self.get(x, y);
                 let left_cell = self.get(x-1, y);
                 let up_cell = self.get(x, y-1);
-                working_cells[[y, x]] = cell.step_e(left_cell, up_cell);
+                working_cells[[y, x]] = cell.step_e(left_cell, up_cell, half_time_step);
             }
         }
         std::mem::swap(&mut working_cells, &mut self.cells);
@@ -80,8 +82,11 @@ impl SimState {
     pub fn height(&self) -> usize {
         self.cells.shape()[0]
     }
-    pub fn get(&self, x: usize, y: usize) -> Cell {
+    pub fn get(&self, x: usize, y: usize) -> Cell<GenericMaterial> {
-        self.cells[[y, x]]
+        self.cells[[y, x]].clone()
+    }
+    pub fn get_mut(&mut self, x: usize, y: usize) -> &mut Cell<GenericMaterial> {
+        &mut self.cells[[y, x]]
     }
 }
 
@@ -104,13 +109,14 @@ impl SimState {
 /// (x+1, y). The `+` only indicates the corner of the cell -- nothing of interest is measured at
 /// the pluses.
 #[derive(Copy, Clone, Default)]
-pub struct Cell {
+pub struct Cell<M> {
     ex: f32,
     ey: f32,
     bz: f32,
+    mat: M,
 }
 
-impl Cell {
+impl<M> Cell<M> {
     pub fn ex(&self) -> f32 {
         self.ex
     }
@@ -120,7 +126,16 @@ impl Cell {
     pub fn bz(&self) -> f32 {
         self.bz
     }
-    fn step_b(self, right: Cell, down: Cell) -> Self {
+    pub fn mat(&self) -> &M {
+        &self.mat
+    }
+    pub fn mat_mut(&mut self) -> &mut M {
+        &mut self.mat
+    }
+}
+
+impl<M: Material + Clone> Cell<M> {
+    fn step_b(self, right: Self, down: Self, _delta_t: f32) -> Self {
         // Maxwell's equation: del x E = -dB/dt
         // Expand: dE_y/dx - dE_x/dy = -dB_z/dt
         // Rearrange: dB_z/dt = dE_x/dy - dE_y/dx
@@ -136,29 +151,62 @@ impl Cell {
             ex: self.ex,
             ey: self.ey,
             bz: self.bz + delta_bz,
+            mat: self.mat,
         }
     }
 
-    fn step_e(self, left: Cell, up: Cell) -> Self {
+    /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
-        // Maxwell's equation: del x B = mu_0 eps_0 dE/dt
+    /// since the simulation spends half a timestep in step_b and the other half in step_e.
-        // N.B: c = 1/sqrt(mu_0 eps_0) so:
+    /// delta_x and delta_y are derived from delta_t (so, make sure delta_t is constant across all calls if the grid spacing is also constant!)
-        // Rearrange: dE/dt = c^2 del x B
+    fn step_e(self, left: Self, up: Self, delta_t: f32) -> Self {
-        // Expand: dE_x/dt = c^2 dB_z/dy   (1);  dE_y/dt = -c^2 dB_z/dx   (2)
+        // Maxwell's equation: \del x B = \mu_0 (J + \eps_0 dE/dt)    where J = current density = \sigma E, \sigma being a material parameter
+        // Expand: \del x B = \mu_0 \sigma E + \mu_0 \eps_0 dE/dt
+        // Substitute: \del x B = S + 1/c^2 dE/dt  where c = 1/\sqrt{\mu_0 \eps_0} is the speed of light, and S = \mu_0 \sigma E for convenience
+        // Rearrange: dE/dt = c^2 (\del x B - S)
+        // Expand: dE_x/dt = c^2 (dB_z/dy - S_x)   (1);   dE_y/dt = c^2 (-dB_z/dx - S_y)   (2)
         // 
-        // Discretize (1): (delta E_x)/(delta t) = c^2 (delta B_z)/(delta y)
+        // Discretize (1): (\delta E_x)/(\delta t) = c^2 (\delta B_z / \delta y - S_x)
-        // Recall: (delta y)/(delta t) = 2c, as from step_b
+        //   Information is propagated across 1/2 \delta x where \delta x = grid spacing of cells.
-        // Substitute: (delta E_x) = c/2 (delta B_z,y)
+        //   Therefore 1/2 \delta x = c \delta t  or  \delta t / \delta x = 1/(2c)
+        // Rearrange: \delta E_x = c^2 (\delta B_z \delta t / \delta y  -  \delta t S_x)
+        // Rearrange: \delta E_x = c (\delta B_z/2 - c \delta t S_x)
         // 
-        // Discretize (2): (delta E_y)/(delta t) = -c^2 (delta B_z)/(delta x)
+        // Discretize (2): (\delta E_y)/(\delta t) = c^2 (-\delta B_z / \delta x - S_y)
-        // Substitute c: (delta E_y) = -c/2 (delta B_z,x)
+        // Rearrange: \delta E_y = c (-\delta B_z / 2 - c \delta_t S_y)
+        
         let delta_bz_y = self.bz - up.bz;
-        let delta_ex = (0.5f32*consts::C) * delta_bz_y;
+        let static_ex = consts::MU0 * self.mat.conductivity() * self.ex;
+        let delta_ex = consts::C * (0.5 * delta_bz_y - consts::C * delta_t * static_ex);
+
         let delta_bz_x = self.bz - left.bz;
-        let delta_ey = -(0.5f32*consts::C) * delta_bz_x;
+        let static_ey = consts::MU0 * self.mat.conductivity() * self.ey;
+        let delta_ey = consts::C * (-0.5 * delta_bz_x - consts::C * delta_t * static_ey);
+        
         Cell {
             ex: self.ex + delta_ex,
             ey: self.ey + delta_ey,
             bz: self.bz,
+            mat: self.mat,
         }
     }
 }
+
+
+pub trait Material {
+    /// Return \sigma, the electrical conductivity.
+    /// For a vacuum, this is zero. For a perfect conductor, \inf.
+    fn conductivity(&self) -> f32 {
+        0.0
+    }
+}
+
+#[derive(Clone, Default)]
+pub struct GenericMaterial {
+    pub conductivity: f32,
+}
+
+impl Material for GenericMaterial {
+    fn conductivity(&self) -> f32 {
+        self.conductivity
+    }
+}

src/render.rs

@@ -55,7 +55,7 @@ impl ColorTermRenderer {
                 //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);
+                let g = norm_color((cell.bz() * 3.0e8) as _);
                 write!(&mut buf, "{}", RGB(r, g, b).paint(square));
             }
             write!(&mut buf, "\n");