Fix the conductivity error (was calculating the wrong timestep)

Also, use f64 and force them to be real.

Cargo.toml

@@ -8,4 +8,5 @@ edition = "2018"
 
 [dependencies]
 ansi_term = "0.12"
+decorum = "0.3"
 ndarray = "0.13"

examples/coremem.rs

@@ -6,8 +6,10 @@ use std::{thread, time};
 fn main() {
     let mut state = SimState::new(101, 101);
 
-    for x in 0..100 {
+    for y in 70..100 {
-        state.get_mut(x, 70).mat_mut().conductivity = 0.00000001;
+        for x in 0..100 {
+            state.get_mut(x, y).mat_mut().conductivity = 1.0.into();
+        }
     }
 
     let mut step = 0u64;

src/lib.rs

@@ -5,6 +5,7 @@
 //!
 //! [1] https://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf
 
+use decorum::R64;
 use ndarray::Array2;
 
 pub mod render;
@@ -19,12 +20,21 @@ pub mod render;
 pub mod consts {
     /// Speed of light in a vacuum; m/s.
     /// Also equal to 1/sqrt(epsilon_0 mu_0)
-    pub const C: f32 = 299792458f32;
+    pub const C: f64 = 299792458.0;
-    // pub const Z0: f32 = 376.73031366857f32;
+    // pub const Z0: f64 = 376.73031366857f32;
     /// Vacuum Permeability
-    pub const MU0: f32 = 1.2566370621219e-6; // H/m
+    pub const MU0: f64 = 1.2566370621219e-6; // H/m
     // Vacuum Permittivity
-    // pub const Eps0: f32 = 8.854187812813e-12 // F⋅m−1
+    // pub const Eps0: f64 = 8.854187812813e-12 // F⋅m−1
+    pub mod real {
+        use crate::R64;
+        pub fn C() -> R64 {
+            super::C.into()
+        }
+        pub fn MU0() -> R64 {
+            super::MU0.into()
+        }
+    }
 }
 
 #[derive(Default)]
@@ -41,7 +51,7 @@ impl SimState {
 
     pub fn step(&mut self) {
         // feature size: 1mm.
-        let half_time_step = 0.0005 * consts::C;
+        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() {
@@ -49,7 +59,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);
+                working_cells[[down_y-1, right_x-1]] = cell.step_b(right_cell, down_cell, half_time_step.into());
             }
         }
         std::mem::swap(&mut working_cells, &mut self.cells);
@@ -60,19 +70,19 @@ 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, half_time_step);
+                working_cells[[y, x]] = cell.step_e(left_cell, up_cell, half_time_step.into());
             }
         }
         std::mem::swap(&mut working_cells, &mut self.cells);
     }
 
-    pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f32) {
+    pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f64) {
         self.cells[[y, x]].ex += ex;
     }
-    pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f32) {
+    pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f64) {
         self.cells[[y, x]].ey += ey;
     }
-    pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f32) {
+    pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f64) {
         self.cells[[y, x]].bz += bz;
     }
 
@@ -110,21 +120,21 @@ impl SimState {
 /// the pluses.
 #[derive(Copy, Clone, Default)]
 pub struct Cell<M> {
-    ex: f32,
+    ex: R64,
-    ey: f32,
+    ey: R64,
-    bz: f32,
+    bz: R64,
     mat: M,
 }
 
 impl<M> Cell<M> {
-    pub fn ex(&self) -> f32 {
+    pub fn ex(&self) -> f64 {
-        self.ex
+        self.ex.into()
     }
-    pub fn ey(&self) -> f32 {
+    pub fn ey(&self) -> f64 {
-        self.ey
+        self.ey.into()
     }
-    pub fn bz(&self) -> f32 {
+    pub fn bz(&self) -> f64 {
-        self.bz
+        self.bz.into()
     }
     pub fn mat(&self) -> &M {
         &self.mat
@@ -135,7 +145,7 @@ impl<M> Cell<M> {
 }
 
 impl<M: Material + Clone> Cell<M> {
-    fn step_b(self, right: Self, down: Self, _delta_t: f32) -> Self {
+    fn step_b(self, right: Self, down: Self, _delta_t: R64) -> 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
@@ -146,7 +156,7 @@ impl<M: Material + Clone> Cell<M> {
         // 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) / (2f32*consts::C);
+        let delta_bz = (delta_ex - delta_ey) / (2.0*consts::C);
         Cell {
             ex: self.ex,
             ey: self.ey,
@@ -158,7 +168,7 @@ impl<M: Material + Clone> Cell<M> {
     /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
     /// since the simulation spends half a timestep in step_b and the other half in step_e.
     /// 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!)
-    fn step_e(self, left: Self, up: Self, delta_t: f32) -> Self {
+    fn step_e(self, left: Self, up: Self, delta_t: R64) -> Self {
         // 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
@@ -174,13 +184,17 @@ impl<M: Material + Clone> Cell<M> {
         // Discretize (2): (\delta E_y)/(\delta t) = c^2 (-\delta B_z / \delta x - S_y)
         // Rearrange: \delta E_y = c (-\delta B_z / 2 - c \delta_t S_y)
         
+        use consts::real::{C, MU0};
+
+        let half = R64::from_inner(0.5);
+        
         let delta_bz_y = self.bz - up.bz;
-        let static_ex = consts::MU0 * self.mat.conductivity() * self.ex;
+        let static_ex: R64 = MU0() * self.mat.conductivity() * self.ex;
-        let delta_ex = consts::C * (0.5 * delta_bz_y - consts::C * delta_t * static_ex);
+        let delta_ex: R64 = C() * (half * delta_bz_y - C() * delta_t * static_ex);
 
         let delta_bz_x = self.bz - left.bz;
-        let static_ey = consts::MU0 * self.mat.conductivity() * self.ey;
+        let static_ey: R64 = MU0() * self.mat.conductivity() * self.ey;
-        let delta_ey = consts::C * (-0.5 * delta_bz_x - consts::C * delta_t * static_ey);
+        let delta_ey: R64 = C() * (-half * delta_bz_x - C() * delta_t * static_ey);
         
         Cell {
             ex: self.ex + delta_ex,
@@ -195,18 +209,18 @@ impl<M: Material + Clone> Cell<M> {
 pub trait Material {
     /// Return \sigma, the electrical conductivity.
     /// For a vacuum, this is zero. For a perfect conductor, \inf.
-    fn conductivity(&self) -> f32 {
+    fn conductivity(&self) -> f64 {
-        0.0
+        0.0.into()
     }
 }
 
 #[derive(Clone, Default)]
 pub struct GenericMaterial {
-    pub conductivity: f32,
+    pub conductivity: f64,
 }
 
 impl Material for GenericMaterial {
-    fn conductivity(&self) -> f32 {
+    fn conductivity(&self) -> f64 {
         self.conductivity
     }
 }

src/render.rs

@@ -29,8 +29,8 @@ fn clamp(v: f32, range: f32) -> f32 {
     v.min(range).max(-range)
 }
 
-fn norm_color(v: f32) -> u8 {
+fn norm_color(v: f64) -> u8 {
-    (v * 64.0 + 128.0).max(0f32).min(255f32) as u8
+    (v * 64.0 + 128.0).max(0.0).min(255.0) as u8
 }
 
 fn curl(x: f32, y: f32) -> f32 {
@@ -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) as _);
+                let g = norm_color((cell.bz() * 3.0e8).into());
                 write!(&mut buf, "{}", RGB(r, g, b).paint(square));
             }
             write!(&mut buf, "\n");