From a82c2f60e73d654afa60e15c0a5ab33c54ebfb11 Mon Sep 17 00:00:00 2001 From: Colin Date: Tue, 14 Jul 2020 23:26:27 -0700 Subject: [PATCH] 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 | 1 + examples/coremem.rs | 18 ++--- src/lib.rs | 156 +++++++++++++++++++++++++++----------------- src/render.rs | 37 ++++++----- 4 files changed, 126 insertions(+), 86 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 3d7ea70..9003afc 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,3 +8,4 @@ edition = "2018" [dependencies] ansi_term = "0.12" +ndarray = "0.13" diff --git a/examples/coremem.rs b/examples/coremem.rs index c9920a6..ab9aed3 100644 --- a/examples/coremem.rs +++ b/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)); } diff --git a/src/lib.rs b/src/lib.rs index e258955..d2ae3d4 100644 --- a/src/lib.rs +++ b/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, + cells: Array2, } 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, } } } diff --git a/src/render.rs b/src/render.rs index 48c949f..4b623fa 100644 --- a/src/render.rs +++ b/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"); } }