diff --git a/src/lib.rs b/src/lib.rs index 3a2ab92..4a6a9a9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -79,7 +79,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.into()); + working_cells[[down_y-1, right_x-1]] = cell.step_h(right_cell, down_cell, half_time_step.into()); } } std::mem::swap(&mut working_cells, &mut self.cells); @@ -98,13 +98,13 @@ impl SimState { } pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f64) { - self.cells[[y, x]].ex += ex; + self.cells[[y, x]].state.ex += ex; } pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f64) { - self.cells[[y, x]].ey += ey; + self.cells[[y, x]].state.ey += ey; } pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f64) { - self.cells[[y, x]].bz += bz; + self.cells[[y, x]].impulse_bz(bz); } pub fn width(&self) -> usize { @@ -145,21 +145,19 @@ impl SimState { /// the pluses. #[derive(Copy, Clone, Default)] pub struct Cell { - ex: R64, - ey: R64, - bz: R64, + state: CellState, mat: M, } impl Cell { pub fn ex(&self) -> f64 { - self.ex.into() + self.state.ex() } pub fn ey(&self) -> f64 { - self.ey.into() + self.state.ey() } - pub fn bz(&self) -> f64 { - self.bz.into() + pub fn hz(&self) -> f64 { + self.state.hz() } pub fn mat(&self) -> &M { &self.mat @@ -170,7 +168,15 @@ impl Cell { } impl Cell { - fn step_b(self, right: Self, down: Self, _delta_t: R64) -> Self { + pub fn bz(&self) -> f64 { + consts::MU0 * (self.hz() + self.mat.mz()) + } + + fn impulse_bz(&mut self, delta_bz: f64) { + self.state.hz = self.mat.step_h(&self.state, delta_bz).into(); + } + + fn step_h(mut self, right: Self, down: Self, _delta_t: R64) -> Self { // ```tex // Maxwell's equation: $\nabla x E = -dB/dt$ // Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$ @@ -181,13 +187,15 @@ impl Cell { // therefore $(\Delta x)/(2 \Delta t) = (\Delta y)/(2 \Delta t) = c$ if we use SI units. // 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_ex = down.ex() - self.ex(); + let delta_ey = right.ey() - self.ey(); let delta_bz = (delta_ex - delta_ey) / (2.0*consts::C); + let hz = self.mat.step_h(&self.state, delta_bz); Cell { - ex: self.ex, - ey: self.ey, - bz: self.bz + delta_bz, + state: CellState { + hz: hz.into(), + ..self.state + }, mat: self.mat, } } @@ -222,29 +230,58 @@ impl Cell { let sigma = self.mat.conductivity(); - let delta_bz_y = self.bz - up.bz; - let ex_rhs = self.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size; + let delta_bz_y = self.bz() - up.bz(); + let ex_rhs = self.state.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size; let ex_next = ex_rhs / (ONE() + C2()*MU0()*sigma*delta_t); - let delta_bz_x = self.bz - left.bz; - let ey_rhs = self.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size; + let delta_bz_x = self.bz() - left.bz(); + let ey_rhs = self.state.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size; let ey_next = ey_rhs / (ONE() + C2()*MU0()*sigma*delta_t); Cell { - ex: ex_next, - ey: ey_next, - bz: self.bz, + state: CellState { + ex: ex_next, + ey: ey_next, + hz: self.state.hz, + }, mat: self.mat, } } } +#[derive(Copy, Clone, Default)] +pub struct CellState { + ex: R64, + ey: R64, + hz: R64, +} + +impl CellState { + pub fn ex(&self) -> f64 { + self.ex.into() + } + pub fn ey(&self) -> f64 { + self.ey.into() + } + pub fn hz(&self) -> f64 { + self.hz.into() + } +} pub trait Material { /// Return \sigma, the electrical conductivity. /// For a vacuum, this is zero. For a perfect conductor, \inf. fn conductivity(&self) -> f64 { - 0.0.into() + 0.0 + } + /// Return the magnetization in the z direction (M_z). + fn mz(&self) -> f64 { + 0.0 + } + /// Return the new h_z, and optionally change any internal state (e.g. magnetization). + fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + // B = mu0*(H+M), and in a default material M=0. + context.hz() + delta_bz / consts::MU0 } }