From a4008dcc1dfae42ce9f63d1994b44eced1fed320 Mon Sep 17 00:00:00 2001
From: colin <colin@uninsane.org>
Date: Fri, 26 Aug 2022 00:41:41 -0700
Subject: [PATCH] cross: step: document the step_h math

---
 crates/cross/src/step/mod.rs | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/crates/cross/src/step/mod.rs b/crates/cross/src/step/mod.rs
index c6fd315..95e900d 100644
--- a/crates/cross/src/step/mod.rs
+++ b/crates/cross/src/step/mod.rs
@@ -63,7 +63,7 @@ impl<'a, R: Real, M: Material<R>> StepEContext<'a, R, M> {
         e_matrix[idx] = new_e;
     }
 
-    /// return the e field for this cell at `t=t_{prev} + 2 \Delta t`
+    /// return the e field for this cell at one time step from where it was previously evaluated.
     pub fn step_e(self) -> Vec3<R> {
         // ```tex
         // Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
@@ -175,15 +175,21 @@ impl<'a, R: Real, M: Material<R>> StepHContext<'a, R, M> {
         m_matrix[idx] = new_m;
     }
 
-    /// return the `(h, m)` fields for this cell at `t=t_{now} + \delta t`
+    /// return the `(h, m)` fields for this cell at one time-step from where it was previously
+    /// evaluated.
     pub fn step_h(self) -> (Vec3<R>, Vec3<R>) {
+        // ```tex
+        // Maxwell-Faraday equation:
+        //     $\nabla x E = -dB/dt$  (1)
+        //
+        // where
+        //     $B = \mu_0 (H + M)$  (2)
+        // ```
         let mu0 = R::mu0();
         let mu0_inv = R::mu0_inv();
         let deltas = self.in_e.delta_e();
-        // println!("spirv-step_h delta_e_struct: {:?}", deltas);
         // \nabla x E
         let nabla_e = deltas.nabla() * self.inv_feature_size;
-        // println!("spirv-step_h nabla_e: {:?}", nabla_e);
         let delta_b = nabla_e * (-self.time_step);
 
         // Relation between these is: B = mu0*(H + M)
@@ -191,11 +197,15 @@ impl<'a, R: Real, M: Material<R>> StepHContext<'a, R, M> {
         let old_m = self.in_m;
         let old_b = (old_h + old_m) * mu0;
 
+        // evaluate the next B value.
+        // eq (1) enforces this relation. H and M can vary freely within that context,
+        // so long as they sum to B as per (2).
+        // we let the Material dictate M, and from that compute H.
         let new_b = old_b + delta_b + self.stim_h * mu0;
-        let mat = self.mat;
+        let mat = self.mat;  // XXX: copy off of self to avoid a rust-gpu bug around ZSTs.
         let new_m = mat.move_b_vec(old_m, new_b);
         let new_h = new_b * mu0_inv - new_m;
-        // println!("spirv-step_h delta_h: {:?}", delta_h);
+
         (new_h, new_m)
     }
 }