broken: expand to two dimensions.

Broken because there seems to be a scaling issue, and also the previous
1d simulation has no effect when ported to this new 2d world.

Cargo.toml

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

examples/coremem.rs

@@ -1,18 +1,18 @@
 use coremem::SimState;
-use coremem::render::ColorTermRenderer;
+use coremem::render::NumericTermRenderer;
 use coremem::consts;
 use std::{thread, time};
 
 fn main() {
-    let mut state = SimState::new(16);
-    state.impulse_b(0, 1.0/consts::C);
-    state.impulse_e(0, 1.0);
-    state.impulse_b(7, 1.0/consts::C);
-    state.impulse_e(7, 1.0);
-    state.impulse_b(15, 1.0/consts::C);
-    state.impulse_e(15, 1.0);
+    let mut state = SimState::new(16, 4);
+    state.impulse_bz(0, 0, 1.0/consts::C);
+    //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_bz(15, 0, 1.0/consts::C);
+    //state.impulse_ey(15, 0, 1.0);
     loop {
-        ColorTermRenderer.render(&state);
+        NumericTermRenderer.render(&state);
         state.step();
         thread::sleep(time::Duration::from_millis(100));
     }

src/lib.rs

@@ -5,6 +5,8 @@
 //!
 //! [1] https://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf
 
+use ndarray::Array2;
+
 pub mod render;
 
 // Some things to keep in mind:
@@ -23,104 +25,136 @@ pub mod consts {
 
 #[derive(Default)]
 pub struct SimState {
-    cells: Vec<Cell>,
+    cells: Array2<Cell>,
 }
 
 impl SimState {
-    pub fn new(size: usize) -> Self {
+    pub fn new(width: usize, height: usize) -> Self {
         Self {
-            cells: vec![Cell::default(); size],
+            cells: Array2::default((height, width))
         }
     }
 
     pub fn step(&mut self) {
-        let mut working_cells = vec![Cell::default(); self.cells.len()];
+        let mut working_cells = self.cells.clone();
+        //let mut working_cells = vec![Cell::default(); self.cells.len()];
         // 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 down_y in 1..self.height() {
+            for right_x in 1..self.width() {
+                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);
+            }
         }
+        std::mem::swap(&mut working_cells, &mut self.cells);
 
-        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);
+        // now advance electic fields
+        for y in 1..self.height() {
+            for x in 1..self.width() {
+                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);
+            }
         }
+        std::mem::swap(&mut working_cells, &mut self.cells);
     }
 
-    pub fn impulse_e(&mut self, idx: usize, e: f32) {
-        self.cells[idx].ez += e;
+    pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f32) {
+        self.cells[[y, x]].ex += ex;
+    }
+    pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f32) {
+        self.cells[[y, x]].ey += ey;
+    }
+    pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f32) {
+        self.cells[[y, x]].bz += bz;
     }
 
-    pub fn impulse_b(&mut self, idx: usize, y: f32) {
-        self.cells[idx].by += y;
+    pub fn width(&self) -> usize {
+        self.cells.shape()[1]
     }
-
-    pub fn cells(&self) -> &[Cell] {
-        &*self.cells
+    pub fn height(&self) -> usize {
+        self.cells.shape()[0]
+    }
+    pub fn get(&self, x: usize, y: usize) -> Cell {
+        self.cells[[y, x]]
     }
 }
 
 /// Conceptually, one cell looks like this:
 ///
-/// +-------------+
-/// |             |
-/// .Ez    .By    .  next cell
-/// |             |
-/// +-------------+
+/// +-------.-------+
+/// |       Ex      |
+/// |               |
+/// |               |
+/// .Ey     .Bz     .
+/// |               |
+/// |               |
+/// |               |
+/// +-------.-------+
 ///
-/// Where +By points up and +Ez points into the page, and the cell has a unit length of 1.
+/// Where the right hand rule indicates that positive Bz is pointing out of the page, towards 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
+/// (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 {
-    /// electric field
-    ez: f32,
-    /// magnetic field
-    by: f32,
+    ex: f32,
+    ey: f32,
+    bz: f32,
 }
 
 impl Cell {
-    pub fn ez(&self) -> f32 {
-        self.ez
+    pub fn ex(&self) -> f32 {
+        self.ex
     }
-    pub fn by(&self) -> f32 {
-        self.by
+    pub fn ey(&self) -> f32 {
+        self.ey
     }
-    fn step_b(self, right: Cell) -> Self {
+    pub fn bz(&self) -> f32 {
+        self.bz
+    }
+    fn step_b(self, right: Cell, down: 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
+        // Expand: dE_y/dx - dE_x/dy = -dB_z/dt
+        // 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.
+        // Simplify: delta B_z = (delta E_x)/c - (delta E_y)/c
+        let delta_ex = down.ex - self.ex;
+        let delta_ey = right.ey - self.ey;
+        let delta_bz = (delta_ex - delta_ey) / consts::C;
         Cell {
-            ez: self.ez,
-            by: self.by + delta_b,
+            ex: self.ex,
+            ey: self.ey,
+            bz: self.bz + delta_bz,
         }
     }
 
-    fn step_e(self, left: Cell) -> Self {
+    fn step_e(self, left: Cell, up: Cell) -> Self {
         // Maxwell's equation: del x B = mu_0 eps_0 dE/dt
-        // Expands: dB_y/dx = mu_0 eps_0 dE_y/dt
-        // Rearrange: dE_y/dt = 1/(mu_0 eps_0) dB_y/dx
-        // Discretize: (delta E_y)/(delta t) = 1/(mu_0 eps_0) (delta dB_y)/(delta x)
-        // Rearrange: delta E_y = (delta t)/(delta x) 1/(mu_0 eps_0) (delta B_y)
-        // Substitute c as in step_b: delta E_y = (mu_0 eps_0)/c (delta B_y)
-        // Note that c = 1/sqrt(mu_0 eps_0), so this becomes:
-        // delta E_y = c (delta B_y)
-        // XXX once again this differs from [1]
-        let delta_b = self.by - left.by; //< delta B_y
-        let delta_e = delta_b * consts::C; //< delta E_z
+        // N.B: c = 1/sqrt(mu_0 eps_0) so:
+        // Rearrange: dE/dt = c^2 del x B
+        // 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
+        // Substitute: (delta E_x) = c (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)
+        let delta_bz_y = self.bz - up.bz;
+        let delta_ex = consts::C * delta_bz_y;
+        let delta_bz_x = self.bz - left.bz;
+        let delta_ey = consts::C * delta_bz_x;
         Cell {
-            ez: self.ez + delta_e,
-            by: self.by,
+            ex: self.ex + delta_ex,
+            ey: self.ey + delta_ey,
+            bz: self.bz,
         }
     }
 }

src/render.rs

@@ -6,13 +6,17 @@ pub struct NumericTermRenderer;
 
 impl NumericTermRenderer {
     pub fn render(&self, state: &SimState) {
-        print!("B:      ");
-        for cell in state.cells() {
-            print!("{:>10.1e} ", cell.by());
-        }
-        print!("\nE:");
-        for cell in state.cells() {
-            print!("{:>10.1e} ", cell.ez());
+        for y in 0..state.height() {
+            for x in 0..state.width() {
+                let cell = state.get(x, y);
+                print!("           {:>10.1e}", cell.ex());
+            }
+            print!("\n");
+            for x in 0..state.width() {
+                let cell = state.get(x, y);
+                print!("{:>10.1e} {:>10.1e}", cell.ey(), cell.bz());
+            }
+            print!("\n");
         }
         print!("\n");
     }
@@ -26,14 +30,15 @@ fn clamp(v: f32, range: f32) -> f32 {
 
 impl ColorTermRenderer {
     pub fn render(&self, state: &SimState) {
-        let square = "█";
-        for cell in state.cells() {
-            let b_value = clamp(cell.by() * consts::C * 128.0, 255.0);
-            let b_color = RGB(0, b_value.max(0.0) as _, b_value.abs() as _);
-            let e_value = clamp(cell.ez() * 128.0, 255.0);
-            let e_color = RGB(e_value.abs() as _, e_value.max(0.0) as _, 0);
-            print!("{}{}", e_color.paint(square), b_color.paint(square));
-        }
-        print!("\n");
+        unimplemented!()
+        // let square = "█";
+        // for cell in state.cells() {
+        //     let b_value = clamp(cell.by() * consts::C * 128.0, 255.0);
+        //     let b_color = RGB(0, b_value.max(0.0) as _, b_value.abs() as _);
+        //     let e_value = clamp(cell.ez() * 128.0, 255.0);
+        //     let e_color = RGB(e_value.abs() as _, e_value.max(0.0) as _, 0);
+        //     print!("{}{}", e_color.paint(square), b_color.paint(square));
+        // }
+        // print!("\n");
     }
 }