Format these math comments so that they render properly for a tex-aware editor

Also document an error in the current logic

examples/coremem.rs

@@ -4,23 +4,47 @@ use coremem::consts;
 use std::{thread, time};
 
 fn main() {
-    let mut state = SimState::new(101, 101, 1e-3 /* feature size */);
+    let width = 201;
+    let mut state = SimState::new(width, 101, 1e-3 /* feature size */);
 
-    for y in 70..100 {
+    for inset in 0..20 {
-        for x in 0..100 {
+        for x in 0..width {
-            state.get_mut(x, y).mat_mut().conductivity = 1.0.into();
+            state.get_mut(x, inset).mat_mut().conductivity += (0.1*(20.0 - inset as f64));
+        }
+        for y in 0..100 {
+            state.get_mut(inset, y).mat_mut().conductivity += (0.1*(20.0 - inset as f64));
+            state.get_mut(width-1 - inset, y).mat_mut().conductivity += (0.1*(20.0 - inset as f64));
+        }
+    }
+    for y in 75..100 {
+        for x in 0..width {
+            // from https://www.thoughtco.com/table-of-electrical-resistivity-conductivity-608499
+            // NB: different sources give pretty different values for this
+            // NB: Simulation misbehaves for values > 10... Proably this model isn't so great.
+            // Maybe use \eps or \xi instead of conductivity.
+            // No, I don't think that gives us what's needed. Rather, the calculation of static_ex
+            // is wrong (it's computed for the wrong value of time)
+            state.get_mut(x, y).mat_mut().conductivity = (2.17).into();
         }
     }
 
     let mut step = 0u64;
     loop {
         step += 1;
-        let imp = 50.0 * ((step as f64)*0.05).sin();
+        let imp = if step < 50 { 
+            2.5 * ((step as f64)*0.04*std::f64::consts::PI).sin()
+        } else {
+            0.0
+        };
         // state.impulse_ex(50, 50, imp);
         // state.impulse_ey(50, 50, imp);
-        state.impulse_bz(50, 50, (imp / 3.0e8) as _);
+        // state.impulse_bz(20, 20, (imp / 3.0e8) as _);
+        // state.impulse_bz(80, 20, (imp / 3.0e8) as _);
+        for x in 0..width {
+            state.impulse_bz(x, 20, (imp / 3.0e8) as _);
+        }
         Renderer.render(&state);
         state.step();
-        thread::sleep(time::Duration::from_millis(33));
+        thread::sleep(time::Duration::from_millis(67));
     }
 }

src/lib.rs

@@ -164,14 +164,16 @@ impl<M> Cell<M> {
 
 impl<M: Material + Clone> Cell<M> {
     fn step_b(self, right: Self, down: Self, _delta_t: R64) -> Self {
-        // Maxwell's equation: del x E = -dB/dt
+        // ```tex
-        // Expand: dE_y/dx - dE_x/dy = -dB_z/dt
+        // Maxwell's equation: $\nabla x E = -dB/dt$
-        // Rearrange: dB_z/dt = dE_x/dy - dE_y/dx
+        // Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$
-        // Discretize: (delta B_z)/(delta t) = (delta E_x)/(delta y) - (delta E_y)/(delta x)
+        // 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, and it we're only stepping half a timestep here (the other half when we advance E),
-        //   therefore (delta x)/(2 delta t) = (delta y)/(2 delta t) = c if we use SI units.
+        //   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)
+        // 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) / (2.0*consts::C);
@@ -187,24 +189,33 @@ impl<M: Material + Clone> Cell<M> {
     /// 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: 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
+        // ```tex
-        // Expand: \del x B = \mu_0 \sigma E + \mu_0 \eps_0 dE/dt
+        // Ampere's circuital law with Maxwell's addition, in SI units:
-        // 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
+        //   $\nabla x B = \mu_0 (J + \epsilon_0 dE/dt)$    where J = current density = $\sigma E$, $\sigma$ being a material parameter
-        // Rearrange: dE/dt = c^2 (\del x B - S)
+        // Expand: $\nabla x B = \mu_0 \sigma E + \mu_0 \epsilon_0 dE/dt$
-        // Expand: dE_x/dt = c^2 (dB_z/dy - S_x)   (1);   dE_y/dt = c^2 (-dB_z/dx - S_y)   (2)
+        // Substitute: $\nabla x B = S + 1/c^2 dE/dt$  where $c = 1/\sqrt{\mu_0 \epsilon_0}$ is the speed of light, and $S = \mu_0 \sigma E$ for convenience
+        // Rearrange: $dE/dt = c^2 (\nabla x B - S)$
+        // Expand: $dE_x/dt = c^2 (dB_z/dy - S_x)$   (1);   $dE_y/dt = c^2 (-dB_z/dx - S_y)$   (2)
         // 
-        // Discretize (1): (\delta E_x)/(\delta t) = c^2 (\delta B_z / \delta y - S_x)
+        // Discretize (1): $(\Delta E_x)/(\Delta t) = c^2 (\Delta B_z / \Delta y - S_x)$
-        //   Information is propagated across 1/2 \delta x where \delta x = grid spacing of cells.
+        //   Information is propagated across $1/2 \Delta x$ where $\Delta x$ = grid spacing of cells.
-        //   Therefore 1/2 \delta x = c \delta t  or  \delta t / \delta x = 1/(2c)
+        //   Therefore $1/2 \Delta x = c \Delta t$  or  $\Delta t / \Delta x = 1/(2c)$
-        // Rearrange: \delta E_x = c^2 (\delta B_z \delta t / \delta y  -  \delta t S_x)
+        // Rearrange: $\Delta E_x = c^2 (\Delta B_z \Delta t / \Delta y  -  \Delta t S_x)$
-        // Rearrange: \delta E_x = c (\delta B_z/2 - c \delta t S_x)
+        // Rearrange: $\Delta E_x = c (\Delta B_z/2 - c \Delta t S_x)$
         // 
-        // Discretize (2): (\delta E_y)/(\delta t) = c^2 (-\delta B_z / \delta x - S_y)
+        // 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)
+        // Rearrange: $\Delta E_y = c (-\Delta B_z / 2 - c \Delta_t S_y)$
+        // ```
         
         use consts::real::{C, HALF, MU0};
 
         let delta_bz_y = self.bz - up.bz;
+        // TYPO: self.ex here should actually be replaced with 0.5(self.ex + self.ex+delta_ex).
+        // i.e. it should be the Ex halfway through the timestep.
+        // BUT: meep defines metals as a medium with epsilon=-\inf. No mention of
+        // conductivity/sigma. Maybe that route is simpler (if equivalent?)
+        // Note that metals don't really seem to have 'bound' current: just 'free' current:
+        //   https://physics.stackexchange.com/questions/227014/are-the-conducting-electrons-in-a-metal-counted-as-free-or-bound-charges
         let static_ex: R64 = MU0() * self.mat.conductivity() * self.ex;
         let delta_ex: R64 = C() * (HALF() * delta_bz_y - C() * delta_t * static_ex);
 

src/render.rs

@@ -51,7 +51,8 @@ impl ColorTermRenderer {
                 let cell = state.get(x, y);
                 //let r = norm_color(cell.bz() * consts::C);
                 let r = 0;
-                let b = 0;
+                let b = (55.0*cell.mat().conductivity).min(255.0) as u8;
+                //let b = 0;
                 //let g = norm_color(cell.ex());
                 //let b = norm_color(cell.ey());
                 //let g = norm_color(curl(cell.ex(), cell.ey()));