diff --git a/Cargo.toml b/Cargo.toml index 9003afc..445c6da 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,4 +8,5 @@ edition = "2018" [dependencies] ansi_term = "0.12" +decorum = "0.3" ndarray = "0.13" diff --git a/examples/coremem.rs b/examples/coremem.rs index 58ec234..63f706d 100644 --- a/examples/coremem.rs +++ b/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 { - state.get_mut(x, 70).mat_mut().conductivity = 0.00000001; + for y in 70..100 { + for x in 0..100 { + state.get_mut(x, y).mat_mut().conductivity = 1.0.into(); + } } let mut step = 0u64; diff --git a/src/lib.rs b/src/lib.rs index 540b509..bad5781 100644 --- a/src/lib.rs +++ b/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 Z0: f32 = 376.73031366857f32; + pub const C: f64 = 299792458.0; + // 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 { - ex: f32, - ey: f32, - bz: f32, + ex: R64, + ey: R64, + bz: R64, mat: M, } impl Cell { - pub fn ex(&self) -> f32 { - self.ex + pub fn ex(&self) -> f64 { + self.ex.into() } - pub fn ey(&self) -> f32 { - self.ey + pub fn ey(&self) -> f64 { + self.ey.into() } - pub fn bz(&self) -> f32 { - self.bz + pub fn bz(&self) -> f64 { + self.bz.into() } pub fn mat(&self) -> &M { &self.mat @@ -135,7 +145,7 @@ impl Cell { } impl Cell { - 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 Cell { // 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 Cell { /// 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 Cell { // 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 delta_ex = consts::C * (0.5 * delta_bz_y - consts::C * delta_t * static_ex); + let static_ex: R64 = MU0() * self.mat.conductivity() * self.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 delta_ey = consts::C * (-0.5 * delta_bz_x - consts::C * delta_t * static_ey); + let static_ey: R64 = MU0() * self.mat.conductivity() * self.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 Cell { 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 + fn conductivity(&self) -> f64 { + 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 } } diff --git a/src/render.rs b/src/render.rs index 52b743d..81c074b 100644 --- a/src/render.rs +++ b/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 { - (v * 64.0 + 128.0).max(0f32).min(255f32) as u8 +fn norm_color(v: f64) -> 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");