make_docs.sh

@@ -12,7 +12,7 @@ cd ~/projects/meanas
 rm -rf _doc_mathimg
 pdoc --pdf --force --template-dir pdoc_templates -o doc . > doc.md
 pandoc --metadata=title:"meanas" --from=markdown+abbreviations --to=html --output=doc.htex --gladtex -s --css pdoc_templates/pdoc.css doc.md
-gladtex -a -n -d _doc_mathimg -c white -b black doc.htex
+gladtex -a -n -d _doc_mathimg -c white doc.htex
 
 # Approach 3: html with gladtex
 #pdoc3 --html --force --template-dir pdoc_templates -o doc .

meanas/fdfd/__init__.py

@@ -1,4 +1,4 @@
-r"""
+"""
 Tools for finite difference frequency-domain (FDFD) simulations and calculations.
 
 These mostly involve picking a single frequency, then setting up and solving a
@@ -19,71 +19,71 @@ Submodules:
 From the "Frequency domain" section of `meanas.fdmath`, we have
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\
+ \\tilde{E}_{l, \\vec{r}} &= \\tilde{E}_{\\vec{r}} e^{-\\imath \\omega l \\Delta_t} \\\\
- \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{H}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\
+ \\tilde{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= \\tilde{H}_{\\vec{r} + \\frac{1}{2}} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t} \\\\
- \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\
+ \\tilde{J}_{l, \\vec{r}} &= \\tilde{J}_{\\vec{r}} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t} \\\\
- \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{M}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega l \Delta_t} \\
+ \\tilde{M}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= \\tilde{M}_{\\vec{r} + \\frac{1}{2}} e^{-\\imath \\omega l \\Delta_t} \\\\
- \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}})
+ \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}})
-    -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} &= -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\
+    -\\Omega^2 \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}} &= -\\imath \\Omega \\tilde{J}_{\\vec{r}} e^{\\imath \\omega \\Delta_t / 2} \\\\
- \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t
+ \\Omega &= 2 \\sin(\\omega \\Delta_t / 2) / \\Delta_t
- \end{aligned}
+ \\end{aligned}
 $$
 
 resulting in
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \tilde{\partial}_t &\Rightarrow -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\
+ \\tilde{\\partial}_t &\\Rightarrow -\\imath \\Omega e^{-\\imath \\omega \\Delta_t / 2}\\\\
-   \hat{\partial}_t &\Rightarrow -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\
+   \\hat{\\partial}_t &\\Rightarrow -\\imath \\Omega e^{ \\imath \\omega \\Delta_t / 2}\\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 Maxwell's equations are then
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \tilde{\nabla} \times \tilde{E}_{\vec{r}} &=
+  \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}} &=
-         \imath \Omega e^{-\imath \omega \Delta_t / 2} \hat{B}_{\vec{r} + \frac{1}{2}}
+         \\imath \\Omega e^{-\\imath \\omega \\Delta_t / 2} \\hat{B}_{\\vec{r} + \\frac{1}{2}}
-                                                     - \hat{M}_{\vec{r} + \frac{1}{2}}  \\
+                                                          - \\hat{M}_{\\vec{r} + \\frac{1}{2}}  \\\\
-  \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &=
+  \\hat{\\nabla} \\times \\hat{H}_{\\vec{r} + \\frac{1}{2}} &=
-        -\imath \Omega e^{ \imath \omega \Delta_t / 2} \tilde{D}_{\vec{r}}
+        -\\imath \\Omega e^{ \\imath \\omega \\Delta_t / 2} \\tilde{D}_{\\vec{r}}
-                                                     + \tilde{J}_{\vec{r}} \\
+                                                          + \\tilde{J}_{\\vec{r}} \\\\
-  \tilde{\nabla} \cdot \hat{B}_{\vec{r} + \frac{1}{2}} &= 0 \\
+  \\tilde{\\nabla} \\cdot \\hat{B}_{\\vec{r} + \\frac{1}{2}} &= 0 \\\\
-  \hat{\nabla} \cdot \tilde{D}_{\vec{r}} &= \rho_{\vec{r}}
+  \\hat{\\nabla} \\cdot \\tilde{D}_{\\vec{r}} &= \\rho_{\\vec{r}}
- \end{aligned}
+ \\end{aligned}
 $$
 
-With $\Delta_t \to 0$, this simplifies to
+With $\\Delta_t \\to 0$, this simplifies to
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \tilde{E}_{l, \vec{r}} &\to \tilde{E}_{\vec{r}} \\
+ \\tilde{E}_{l, \\vec{r}} &\\to \\tilde{E}_{\\vec{r}} \\\\
- \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{H}_{\vec{r} + \frac{1}{2}} \\
+ \\tilde{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &\\to \\tilde{H}_{\\vec{r} + \\frac{1}{2}} \\\\
- \tilde{J}_{l, \vec{r}} &\to \tilde{J}_{\vec{r}} \\
+ \\tilde{J}_{l, \\vec{r}} &\\to \\tilde{J}_{\\vec{r}} \\\\
- \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{M}_{\vec{r} + \frac{1}{2}} \\
+ \\tilde{M}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &\\to \\tilde{M}_{\\vec{r} + \\frac{1}{2}} \\\\
- \Omega &\to \omega \\
+ \\Omega &\\to \\omega \\\\
- \tilde{\partial}_t &\to -\imath \omega \\
+ \\tilde{\\partial}_t &\\to -\\imath \\omega \\\\
-   \hat{\partial}_t &\to -\imath \omega \\
+   \\hat{\\partial}_t &\\to -\\imath \\omega \\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 and then
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \tilde{\nabla} \times \tilde{E}_{\vec{r}} &=
+  \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}} &=
-         \imath \omega \hat{B}_{\vec{r} + \frac{1}{2}}
+         \\imath \\omega \\hat{B}_{\\vec{r} + \\frac{1}{2}}
-                     - \hat{M}_{\vec{r} + \frac{1}{2}}  \\
+                       - \\hat{M}_{\\vec{r} + \\frac{1}{2}}  \\\\
-  \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &=
+  \\hat{\\nabla} \\times \\hat{H}_{\\vec{r} + \\frac{1}{2}} &=
-        -\imath \omega \tilde{D}_{\vec{r}}
+        -\\imath \\omega \\tilde{D}_{\\vec{r}}
-                     + \tilde{J}_{\vec{r}} \\
+                       + \\tilde{J}_{\\vec{r}} \\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 $$
- \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}})
+ \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}})
-    -\omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \omega \tilde{J}_{\vec{r}} \\
+    -\\omega^2 \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}} = -\\imath \\omega \\tilde{J}_{\\vec{r}} \\\\
 $$
 
 # TODO FDFD?

meanas/fdfd/functional.py

@@ -189,10 +189,10 @@ def e_tfsf_source(
 
 
 def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield_t, cfdfield_t], cfdfield_t]:
-    r"""
+    """
     Generates a function that takes the single-frequency `E` and `H` fields
-    and calculates the cross product `E` x `H` = $E \times H$ as required
+    and calculates the cross product `E` x `H` = $E \\times H$ as required
-    for the Poynting vector, $S = E \times H$
+    for the Poynting vector, $S = E \\times H$
 
     Note:
         This function also shifts the input `E` field by one cell as required

meanas/fdfd/operators.py

@@ -45,14 +45,14 @@ def e_full(
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    r"""
+    """
     Wave operator
-     $$ \nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon $$
+     $$ \\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\Omega^2 \\epsilon $$
 
         del x (1/mu * del x) - omega**2 * epsilon
 
      for use with the E-field, with wave equation
-     $$ (\nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon) E = -\imath \omega J $$
+     $$ (\\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\Omega^2 \\epsilon) E = -\\imath \\omega J $$
 
         (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
 
@@ -131,14 +131,14 @@ def h_full(
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    r"""
+    """
     Wave operator
-     $$ \nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu $$
+     $$ \\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu $$
 
         del x (1/epsilon * del x) - omega**2 * mu
 
      for use with the H-field, with wave equation
-     $$ (\nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu) E = \imath \omega M $$
+     $$ (\\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu) E = \\imath \\omega M $$
 
         (del x (1/epsilon * del x) - omega**2 * mu) E = i * omega * M
 
@@ -188,28 +188,28 @@ def eh_full(
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    r"""
+    """
     Wave operator for `[E, H]` field representation. This operator implements Maxwell's
      equations without cancelling out either E or H. The operator is
-    $$  \begin{bmatrix}
+    $$  \\begin{bmatrix}
-        -\imath \omega \epsilon  &  \nabla \times      \\
+        -\\imath \\omega \\epsilon  &  \\nabla \\times      \\\\
-        \nabla \times            &  \imath \omega \mu
+        \\nabla \\times             &  \\imath \\omega \\mu
-        \end{bmatrix} $$
+        \\end{bmatrix} $$
 
         [[-i * omega * epsilon,  del x         ],
          [del x,                 i * omega * mu]]
 
     for use with a field vector of the form `cat(vec(E), vec(H))`:
-    $$  \begin{bmatrix}
+    $$  \\begin{bmatrix}
-        -\imath \omega \epsilon  &  \nabla \times      \\
+        -\\imath \\omega \\epsilon  &  \\nabla \\times      \\\\
-        \nabla \times            &  \imath \omega \mu
+        \\nabla \\times             &  \\imath \\omega \\mu
-        \end{bmatrix}
+        \\end{bmatrix}
-        \begin{bmatrix} E \\
+        \\begin{bmatrix} E \\\\
-                        H
+                         H
-        \end{bmatrix}
+        \\end{bmatrix}
-        = \begin{bmatrix} J \\
+        = \\begin{bmatrix} J \\\\
-                         -M
+                          -M
-          \end{bmatrix} $$
+          \\end{bmatrix} $$
 
     Args:
         omega: Angular frequency of the simulation

meanas/fdfd/waveguide_2d.py

@@ -1,4 +1,4 @@
-r"""
+"""
 Operators and helper functions for waveguides with unchanging cross-section.
 
 The propagation direction is chosen to be along the z axis, and all fields
@@ -12,166 +12,166 @@ As the z-dependence is known, all the functions in this file assume a 2D grid
 
 Consider Maxwell's equations in continuous space, in the frequency domain. Assuming
 a structure with some (x, y) cross-section extending uniformly into the z dimension,
-with a diagonal $\epsilon$ tensor, we have
+with a diagonal $\\epsilon$ tensor, we have
 
 $$
-\begin{aligned}
+\\begin{aligned}
-\nabla \times \vec{E}(x, y, z) &= -\imath \omega \mu \vec{H} \\
+\\nabla \\times \\vec{E}(x, y, z) &= -\\imath \\omega \\mu \\vec{H} \\\\
-\nabla \times \vec{H}(x, y, z) &=  \imath \omega \epsilon \vec{E} \\
+\\nabla \\times \\vec{H}(x, y, z) &=  \\imath \\omega \\epsilon \\vec{E} \\\\
-\vec{E}(x,y,z) &= (\vec{E}_t(x, y) + E_z(x, y)\vec{z}) e^{-\gamma z} \\
+\\vec{E}(x,y,z) = (\\vec{E}_t(x, y) + E_z(x, y)\\vec{z}) e^{-\\gamma z} \\\\
-\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\gamma z} \\
+\\vec{H}(x,y,z) = (\\vec{H}_t(x, y) + H_z(x, y)\\vec{z}) e^{-\\gamma z} \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
 Expanding the first two equations into vector components, we get
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\imath \omega \mu_{xx} H_x &= \partial_y E_z - \partial_z E_y \\
+-\\imath \\omega \\mu_{xx} H_x &= \\partial_y E_z - \\partial_z E_y \\\\
--\imath \omega \mu_{yy} H_y &= \partial_z E_x - \partial_x E_z \\
+-\\imath \\omega \\mu_{yy} H_y &= \\partial_z E_x - \\partial_x E_z \\\\
--\imath \omega \mu_{zz} H_z &= \partial_x E_y - \partial_y E_x \\
+-\\imath \\omega \\mu_{zz} H_z &= \\partial_x E_y - \\partial_y E_x \\\\
-\imath \omega \epsilon_{xx} E_x &= \partial_y H_z - \partial_z H_y \\
+\\imath \\omega \\epsilon_{xx} E_x &= \\partial_y H_z - \\partial_z H_y \\\\
-\imath \omega \epsilon_{yy} E_y &= \partial_z H_x - \partial_x H_z \\
+\\imath \\omega \\epsilon_{yy} E_y &= \\partial_z H_x - \\partial_x H_z \\\\
-\imath \omega \epsilon_{zz} E_z &= \partial_x H_y - \partial_y H_x \\
+\\imath \\omega \\epsilon_{zz} E_z &= \\partial_x H_y - \\partial_y H_x \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
-Substituting in our expressions for $\vec{E}$, $\vec{H}$ and discretizing:
+Substituting in our expressions for $\\vec{E}$, $\\vec{H}$ and discretizing:
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \gamma E_y \\
+-\\imath \\omega \\mu_{xx} H_x &= \\tilde{\\partial}_y E_z + \\gamma E_y \\\\
--\imath \omega \mu_{yy} H_y &= -\gamma E_x - \tilde{\partial}_x E_z \\
+-\\imath \\omega \\mu_{yy} H_y &= -\\gamma E_x - \\tilde{\\partial}_x E_z \\\\
--\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\
+-\\imath \\omega \\mu_{zz} H_z &= \\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x \\\\
-\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \gamma H_y \\
+\\imath \\omega \\epsilon_{xx} E_x &= \\hat{\\partial}_y H_z + \\gamma H_y \\\\
-\imath \omega \epsilon_{yy} E_y &= -\gamma H_x - \hat{\partial}_x H_z \\
+\\imath \\omega \\epsilon_{yy} E_y &= -\\gamma H_x - \\hat{\\partial}_x H_z \\\\
-\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
+\\imath \\omega \\epsilon_{zz} E_z &= \\hat{\\partial}_x H_y - \\hat{\\partial}_y H_x \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
 Rewrite the last three equations as
 $$
-\begin{aligned}
+\\begin{aligned}
-\gamma H_y &=  \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\
+\\gamma H_y &=  \\imath \\omega \\epsilon_{xx} E_x - \\hat{\\partial}_y H_z \\\\
-\gamma H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\
+\\gamma H_x &= -\\imath \\omega \\epsilon_{yy} E_y - \\hat{\\partial}_x H_z \\\\
-\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
+\\imath \\omega E_z &= \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x H_y - \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y H_x \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
-Now apply $\gamma \tilde{\partial}_x$ to the last equation,
+Now apply $\\gamma \\tilde{\\partial}_x$ to the last equation,
-then substitute in for $\gamma H_x$ and $\gamma H_y$:
+then substitute in for $\\gamma H_x$ and $\\gamma H_y$:
 
 $$
-\begin{aligned}
+\\begin{aligned}
-\gamma \tilde{\partial}_x \imath \omega E_z &= \gamma \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y
+\\gamma \\tilde{\\partial}_x \\imath \\omega E_z &= \\gamma \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x H_y
-                                             - \gamma \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
+                                                - \\gamma \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y H_x \\\\
-        &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z)
+        &= \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x ( \\imath \\omega \\epsilon_{xx} E_x - \\hat{\\partial}_y H_z)
-         - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z)  \\
+         - \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (-\\imath \\omega \\epsilon_{yy} E_y - \\hat{\\partial}_x H_z)  \\\\
-        &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x)
+        &= \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x ( \\imath \\omega \\epsilon_{xx} E_x)
-         - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y)  \\
+         - \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (-\\imath \\omega \\epsilon_{yy} E_y)  \\\\
-\gamma \tilde{\partial}_x E_z &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+\\gamma \\tilde{\\partial}_x E_z &= \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                               + \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\
+                                  + \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y) \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
-With a similar approach (but using $\gamma \tilde{\partial}_y$ instead), we can get
+With a similar approach (but using $\\gamma \\tilde{\\partial}_y$ instead), we can get
 
 $$
-\begin{aligned}
+\\begin{aligned}
-\gamma \tilde{\partial}_y E_z &= \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+\\gamma \\tilde{\\partial}_y E_z &= \\tilde{\\partial}_y \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                               + \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\
+                                  + \\tilde{\\partial}_y \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y) \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
-We can combine this equation for $\gamma \tilde{\partial}_y E_z$ with
+We can combine this equation for $\\gamma \\tilde{\\partial}_y E_z$ with
-the unused $\imath \omega \mu_{xx} H_x$ and $\imath \omega \mu_{yy} H_y$ equations to get
+the unused $\\imath \\omega \\mu_{xx} H_x$ and $\\imath \\omega \\mu_{yy} H_y$ equations to get
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\imath \omega \mu_{xx} \gamma H_x &=  \gamma^2 E_y + \gamma \tilde{\partial}_y E_z \\
+-\\imath \\omega \\mu_{xx} \\gamma H_x &=  \\gamma^2 E_y + \\gamma \\tilde{\\partial}_y E_z \\\\
--\imath \omega \mu_{xx} \gamma H_x &=  \gamma^2 E_y + \tilde{\partial}_y (
+-\\imath \\omega \\mu_{xx} \\gamma H_x &=  \\gamma^2 E_y + \\tilde{\\partial}_y (
-                                      \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+                                      \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                                    + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+                                    + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-                                    )\\
+                                    )\\\\
-\end{aligned}
+\\end{aligned}
 $$
 
 and
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\imath \omega \mu_{yy} \gamma H_y &= -\gamma^2 E_x - \gamma \tilde{\partial}_x E_z \\
+-\\imath \\omega \\mu_{yy} \\gamma H_y &= -\\gamma^2 E_x - \\gamma \\tilde{\\partial}_x E_z \\\\
--\imath \omega \mu_{yy} \gamma H_y &= -\gamma^2 E_x - \tilde{\partial}_x (
+-\\imath \\omega \\mu_{yy} \\gamma H_y &= -\\gamma^2 E_x - \\tilde{\\partial}_x (
-                                      \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+                                      \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                                    + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+                                    + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-                                    )\\
+                                    )\\\\
-\end{aligned}
+\\end{aligned}
 $$
 
-However, based on our rewritten equation for $\gamma H_x$ and the so-far unused
+However, based on our rewritten equation for $\\gamma H_x$ and the so-far unused
-equation for $\imath \omega \mu_{zz} H_z$ we can also write
+equation for $\\imath \\omega \\mu_{zz} H_z$ we can also write
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\imath \omega \mu_{xx} (\gamma H_x) &= -\imath \omega \mu_{xx} (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\
+-\\imath \\omega \\mu_{xx} (\\gamma H_x) &= -\\imath \\omega \\mu_{xx} (-\\imath \\omega \\epsilon_{yy} E_y - \\hat{\\partial}_x H_z) \\\\
-                                     &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y
+                                         &= -\\omega^2 \\mu_{xx} \\epsilon_{yy} E_y
-                                        +\imath \omega \mu_{xx} \hat{\partial}_x (
+                                            +\\imath \\omega \\mu_{xx} \\hat{\\partial}_x (
-                                           \frac{1}{-\imath \omega \mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x)) \\
+                                                \\frac{1}{-\\imath \\omega \\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x)) \\\\
-                                     &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y
+                                         &= -\\omega^2 \\mu_{xx} \\epsilon_{yy} E_y
-                                        -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+                                            -\\mu_{xx} \\hat{\\partial}_x \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
 and, similarly,
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\imath \omega \mu_{yy} (\gamma H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x
+-\\imath \\omega \\mu_{yy} (\\gamma H_y) &= \\omega^2 \\mu_{yy} \\epsilon_{xx} E_x
-                                           +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+                                           +\\mu_{yy} \\hat{\\partial}_y \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
 By combining both pairs of expressions, we get
 
 $$
-\begin{aligned}
+\\begin{aligned}
--\gamma^2 E_x - \tilde{\partial}_x (
+-\\gamma^2 E_x - \\tilde{\\partial}_x (
-    \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+    \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-  + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+  + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-    ) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x
+    ) &= \\omega^2 \\mu_{yy} \\epsilon_{xx} E_x
-        +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+        +\\mu_{yy} \\hat{\\partial}_y \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\gamma^2 E_y + \tilde{\partial}_y (
+\\gamma^2 E_y + \\tilde{\\partial}_y (
-    \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+    \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-  + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+  + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-    ) &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y
+    ) &= -\\omega^2 \\mu_{xx} \\epsilon_{yy} E_y
-         -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+         -\\mu_{xx} \\hat{\\partial}_x \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\end{aligned}
+\\end{aligned}
 $$
 
 Using these, we can construct the eigenvalue problem
 
 $$
-\beta^2 \begin{bmatrix} E_x \\
+\\beta^2 \\begin{bmatrix} E_x \\\\
-                        E_y \end{bmatrix} =
+                          E_y \\end{bmatrix} =
-    (\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\
+    (\\omega^2 \\begin{bmatrix} \\mu_{yy} \\epsilon_{xx} & 0 \\\\
-                                                   0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+                                                       0 & \\mu_{xx} \\epsilon_{yy} \\end{bmatrix} +
-              \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\
+                 \\begin{bmatrix} -\\mu_{yy} \\hat{\\partial}_y \\\\
-                               \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1}
+                                   \\mu_{xx} \\hat{\\partial}_x \\end{bmatrix} \\mu_{zz}^{-1}
-              \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+                 \\begin{bmatrix} -\\tilde{\\partial}_y & \\tilde{\\partial}_x \\end{bmatrix} +
-      \begin{bmatrix} \tilde{\partial}_x \\
+      \\begin{bmatrix} \\tilde{\\partial}_x \\\\
-                      \tilde{\partial}_y \end{bmatrix} \epsilon_{zz}^{-1}
+                       \\tilde{\\partial}_y \\end{bmatrix} \\epsilon_{zz}^{-1}
-                 \begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix})
+                 \\begin{bmatrix} \\hat{\\partial}_x \\epsilon_{xx} & \\hat{\\partial}_y \\epsilon_{yy} \\end{bmatrix})
-    \begin{bmatrix} E_x \\
+    \\begin{bmatrix} E_x \\\\
-                    E_y \end{bmatrix}
+                     E_y \\end{bmatrix}
 $$
 
-where $\gamma = \imath\beta$. In the literature, $\beta$ is usually used to denote
+where $\\gamma = \\imath\\beta$. In the literature, $\\beta$ is usually used to denote
 the lossless/real part of the propagation constant, but in `meanas` it is allowed to
 be complex.
 
 An equivalent eigenvalue problem can be formed using the $H_x$ and $H_y$ fields, if those are more convenient.
 
-Note that $E_z$ was never discretized, so $\gamma$ and $\beta$ will need adjustment
+Note that $E_z$ was never discretized, so $\\gamma$ and $\\beta$ will need adjustment
 to account for numerical dispersion if the result is introduced into a space with a discretized z-axis.
 
 
@@ -198,7 +198,7 @@ def operator_e(
         epsilon: vfdfield_t,
         mu: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    r"""
+    """
     Waveguide operator of the form
 
         omega**2 * mu * epsilon +
@@ -210,18 +210,18 @@ def operator_e(
     More precisely, the operator is
 
     $$
-    \omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\
+    \\omega^2 \\begin{bmatrix} \\mu_{yy} \\epsilon_{xx} & 0 \\\\
-                                                      0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+                                                       0 & \\mu_{xx} \\epsilon_{yy} \\end{bmatrix} +
-                 \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\
+                 \\begin{bmatrix} -\\mu_{yy} \\hat{\\partial}_y \\\\
-                                   \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1}
+                                   \\mu_{xx} \\hat{\\partial}_x \\end{bmatrix} \\mu_{zz}^{-1}
-                 \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+                 \\begin{bmatrix} -\\tilde{\\partial}_y & \\tilde{\\partial}_x \\end{bmatrix} +
-      \begin{bmatrix} \tilde{\partial}_x \\
+      \\begin{bmatrix} \\tilde{\\partial}_x \\\\
-                       \tilde{\partial}_y \end{bmatrix} \epsilon_{zz}^{-1}
+                       \\tilde{\\partial}_y \\end{bmatrix} \\epsilon_{zz}^{-1}
-                 \begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix}
+                 \\begin{bmatrix} \\hat{\\partial}_x \\epsilon_{xx} & \\hat{\\partial}_y \\epsilon_{yy} \\end{bmatrix}
     $$
 
-    $\tilde{\partial}_x$ and $\hat{\partial}_x$ are the forward and backward derivatives along x,
+    $\\tilde{\\partial}_x$ and $\\hat{\\partial}_x$ are the forward and backward derivatives along x,
-    and each $\epsilon_{xx}$, $\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
+    and each $\\epsilon_{xx}$, $\\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
     property distribution.
 
     This operator can be used to form an eigenvalue problem of the form
@@ -253,10 +253,9 @@ def operator_e(
     mu_yx = sparse.diags(numpy.hstack((mu_parts[1], mu_parts[0])))
     mu_z_inv = sparse.diags(1 / mu_parts[2])
 
-    op = (omega * omega * mu_yx @ eps_xy
+    op = omega * omega * mu_yx @ eps_xy + \
-        + mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx))
+        mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) + \
-        + sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy
+        sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy
-        )
     return op
 
 
@@ -266,7 +265,7 @@ def operator_h(
         epsilon: vfdfield_t,
         mu: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    r"""
+    """
     Waveguide operator of the form
 
         omega**2 * epsilon * mu +
@@ -278,18 +277,18 @@ def operator_h(
     More precisely, the operator is
 
     $$
-    \omega^2 \begin{bmatrix} \epsilon_{yy} \mu_{xx} & 0 \\
+    \\omega^2 \\begin{bmatrix} \\epsilon_{yy} \\mu_{xx} & 0 \\\\
-                                                  0 & \epsilon_{xx} \mu_{yy} \end{bmatrix} +
+                                                      0 & \\epsilon_{xx} \\mu_{yy} \\end{bmatrix} +
-                 \begin{bmatrix} -\epsilon_{yy} \tilde{\partial}_y \\
+                 \\begin{bmatrix} -\\epsilon_{yy} \\tilde{\\partial}_y \\\\
-                                  \epsilon_{xx} \tilde{\partial}_x \end{bmatrix} \epsilon_{zz}^{-1}
+                                   \\epsilon_{xx} \\tilde{\\partial}_x \\end{bmatrix} \\epsilon_{zz}^{-1}
-                 \begin{bmatrix} -\hat{\partial}_y & \hat{\partial}_x \end{bmatrix} +
+                 \\begin{bmatrix} -\\hat{\\partial}_y & \\hat{\\partial}_x \\end{bmatrix} +
-      \begin{bmatrix} \hat{\partial}_x \\
+      \\begin{bmatrix} \\hat{\\partial}_x \\\\
-                      \hat{\partial}_y \end{bmatrix} \mu_{zz}^{-1}
+                       \\hat{\\partial}_y \\end{bmatrix} \\mu_{zz}^{-1}
-                 \begin{bmatrix} \tilde{\partial}_x \mu_{xx} & \tilde{\partial}_y \mu_{yy} \end{bmatrix}
+                 \\begin{bmatrix} \\tilde{\\partial}_x \\mu_{xx} & \\tilde{\\partial}_y \\mu_{yy} \\end{bmatrix}
     $$
 
-    $\tilde{\partial}_x$ and $\hat{\partial}_x$ are the forward and backward derivatives along x,
+    $\\tilde{\\partial}_x$ and $\\hat{\\partial}_x$ are the forward and backward derivatives along x,
-    and each $\epsilon_{xx}$, $\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
+    and each $\\epsilon_{xx}$, $\\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
     property distribution.
 
     This operator can be used to form an eigenvalue problem of the form
@@ -321,10 +320,10 @@ def operator_h(
     mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
     mu_z_inv = sparse.diags(1 / mu_parts[2])
 
-    op = (omega * omega * eps_yx @ mu_xy
+    op = omega * omega * eps_yx @ mu_xy + \
-        + eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx))
+        eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) + \
-        + sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy
+        sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy
-        )
+
     return op
 
 

meanas/fdfd/waveguide_cyl.py

@@ -73,10 +73,8 @@ def cylindrical_operator(
 
     omega2 = omega * omega
 
-    op = (
+    op = (omega2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
-        (omega2 * diag((Tx, Ty)) + pa) @ diag((a0, a1))
         - (sparse.bmat(((None, Ty), (Tx, None))) + pb / omega2) @ diag((b0, b1))
-        )
     return op
 
 

meanas/fdmath/__init__.py

@@ -1,4 +1,4 @@
-r"""
+"""
 
 Basic discrete calculus for finite difference (fd) simulations.
 
@@ -43,11 +43,11 @@ Scalar derivatives and cell shifts
 ----------------------------------
 
 Define the discrete forward derivative as
- $$ [\tilde{\partial}_x f]_{m + \frac{1}{2}} = \frac{1}{\Delta_{x, m}} (f_{m + 1} - f_m) $$
+ $$ [\\tilde{\\partial}_x f]_{m + \\frac{1}{2}} = \\frac{1}{\\Delta_{x, m}} (f_{m + 1} - f_m) $$
  where $f$ is a function defined at discrete locations on the x-axis (labeled using $m$).
- The value at $m$ occupies a length $\Delta_{x, m}$ along the x-axis. Note that $m$
+ The value at $m$ occupies a length $\\Delta_{x, m}$ along the x-axis. Note that $m$
  is an index along the x-axis, _not_ necessarily an x-coordinate, since each length
- $\Delta_{x, m}, \Delta_{x, m+1}, ...$ is independently chosen.
+ $\\Delta_{x, m}, \\Delta_{x, m+1}, ...$ is independently chosen.
 
 If we treat `f` as a 1D array of values, with the `i`-th value `f[i]` taking up a length `dx[i]`
 along the x-axis, the forward derivative is
@@ -56,13 +56,13 @@ along the x-axis, the forward derivative is
 
 
 Likewise, discrete reverse derivative is
- $$ [\hat{\partial}_x f ]_{m - \frac{1}{2}} = \frac{1}{\Delta_{x, m}} (f_{m} - f_{m - 1}) $$
+ $$ [\\hat{\\partial}_x f ]_{m - \\frac{1}{2}} = \\frac{1}{\\Delta_{x, m}} (f_{m} - f_{m - 1}) $$
  or
 
     deriv_back(f)[i] = (f[i] - f[i - 1]) / dx[i]
 
 The derivatives' values are shifted by a half-cell relative to the original function, and
-will have different cell widths if all the `dx[i]` ( $\Delta_{x, m}$ ) are not
+will have different cell widths if all the `dx[i]` ( $\\Delta_{x, m}$ ) are not
 identical:
 
     [figure: derivatives and cell sizes]
@@ -88,19 +88,19 @@ identical:
 
 In the above figure,
  `f0 =` $f_0$, `f1 =` $f_1$
- `Df0 =` $[\tilde{\partial}f]_{0 + \frac{1}{2}}$
+ `Df0 =` $[\\tilde{\\partial}f]_{0 + \\frac{1}{2}}$
- `Df1 =` $[\tilde{\partial}f]_{1 + \frac{1}{2}}$
+ `Df1 =` $[\\tilde{\\partial}f]_{1 + \\frac{1}{2}}$
- `df0 =` $[\hat{\partial}f]_{0 - \frac{1}{2}}$
+ `df0 =` $[\\hat{\\partial}f]_{0 - \\frac{1}{2}}$
  etc.
 
-The fractional subscript $m + \frac{1}{2}$ is used to indicate values defined
+The fractional subscript $m + \\frac{1}{2}$ is used to indicate values defined
  at shifted locations relative to the original $m$, with corresponding lengths
- $$ \Delta_{x, m + \frac{1}{2}} = \frac{1}{2} * (\Delta_{x, m} + \Delta_{x, m + 1}) $$
+ $$ \\Delta_{x, m + \\frac{1}{2}} = \\frac{1}{2} * (\\Delta_{x, m} + \\Delta_{x, m + 1}) $$
 
-Just as $m$ is not itself an x-coordinate, neither is $m + \frac{1}{2}$;
+Just as $m$ is not itself an x-coordinate, neither is $m + \\frac{1}{2}$;
 carefully note the positions of the various cells in the above figure vs their labels.
 If the positions labeled with $m$ are considered the "base" or "original" grid,
-the positions labeled with $m + \frac{1}{2}$ are said to lie on a "dual" or
+the positions labeled with $m + \\frac{1}{2}$ are said to lie on a "dual" or
 "derived" grid.
 
 For the remainder of the `Discrete calculus` section, all figures will show
@@ -113,12 +113,12 @@ Gradients and fore-vectors
 --------------------------
 
 Expanding to three dimensions, we can define two gradients
-  $$ [\tilde{\nabla} f]_{m,n,p} = \vec{x} [\tilde{\partial}_x f]_{m + \frac{1}{2},n,p} +
+  $$ [\\tilde{\\nabla} f]_{m,n,p} = \\vec{x} [\\tilde{\\partial}_x f]_{m + \\frac{1}{2},n,p} +
-                                  \vec{y} [\tilde{\partial}_y f]_{m,n + \frac{1}{2},p} +
+                                    \\vec{y} [\\tilde{\\partial}_y f]_{m,n + \\frac{1}{2},p} +
-                                  \vec{z} [\tilde{\partial}_z f]_{m,n,p + \frac{1}{2}}  $$
+                                    \\vec{z} [\\tilde{\\partial}_z f]_{m,n,p + \\frac{1}{2}}  $$
-  $$ [\hat{\nabla} f]_{m,n,p} = \vec{x} [\hat{\partial}_x f]_{m + \frac{1}{2},n,p} +
+  $$ [\\hat{\\nabla} f]_{m,n,p} = \\vec{x} [\\hat{\\partial}_x f]_{m + \\frac{1}{2},n,p} +
-                                \vec{y} [\hat{\partial}_y f]_{m,n + \frac{1}{2},p} +
+                                  \\vec{y} [\\hat{\\partial}_y f]_{m,n + \\frac{1}{2},p} +
-                                \vec{z} [\hat{\partial}_z f]_{m,n,p + \frac{1}{2}}  $$
+                                  \\vec{z} [\\hat{\\partial}_z f]_{m,n,p + \\frac{1}{2}}  $$
 
  or
 
@@ -144,12 +144,12 @@ y in y, and z in z.
 
 We call the resulting object a "fore-vector" or "back-vector", depending
 on the direction of the shift. We write it as
-  $$ \tilde{g}_{m,n,p} = \vec{x} g^x_{m + \frac{1}{2},n,p} +
+  $$ \\tilde{g}_{m,n,p} = \\vec{x} g^x_{m + \\frac{1}{2},n,p} +
-                         \vec{y} g^y_{m,n + \frac{1}{2},p} +
+                          \\vec{y} g^y_{m,n + \\frac{1}{2},p} +
-                         \vec{z} g^z_{m,n,p + \frac{1}{2}} $$
+                          \\vec{z} g^z_{m,n,p + \\frac{1}{2}} $$
-  $$ \hat{g}_{m,n,p} = \vec{x} g^x_{m - \frac{1}{2},n,p} +
+  $$ \\hat{g}_{m,n,p} = \\vec{x} g^x_{m - \\frac{1}{2},n,p} +
-                       \vec{y} g^y_{m,n - \frac{1}{2},p} +
+                        \\vec{y} g^y_{m,n - \\frac{1}{2},p} +
-                       \vec{z} g^z_{m,n,p - \frac{1}{2}} $$
+                        \\vec{z} g^z_{m,n,p - \\frac{1}{2}} $$
 
 
     [figure: gradient / fore-vector]
@@ -172,15 +172,15 @@ Divergences
 
 There are also two divergences,
 
-  $$ d_{n,m,p} = [\tilde{\nabla} \cdot \hat{g}]_{n,m,p}
+  $$ d_{n,m,p} = [\\tilde{\\nabla} \\cdot \\hat{g}]_{n,m,p}
-               = [\tilde{\partial}_x g^x]_{m,n,p} +
+               = [\\tilde{\\partial}_x g^x]_{m,n,p} +
-                 [\tilde{\partial}_y g^y]_{m,n,p} +
+                 [\\tilde{\\partial}_y g^y]_{m,n,p} +
-                 [\tilde{\partial}_z g^z]_{m,n,p}   $$
+                 [\\tilde{\\partial}_z g^z]_{m,n,p}   $$
 
-  $$ d_{n,m,p} = [\hat{\nabla} \cdot \tilde{g}]_{n,m,p}
+  $$ d_{n,m,p} = [\\hat{\\nabla} \\cdot \\tilde{g}]_{n,m,p}
-               = [\hat{\partial}_x g^x]_{m,n,p} +
+               = [\\hat{\\partial}_x g^x]_{m,n,p} +
-                 [\hat{\partial}_y g^y]_{m,n,p} +
+                 [\\hat{\\partial}_y g^y]_{m,n,p} +
-                 [\hat{\partial}_z g^z]_{m,n,p}  $$
+                 [\\hat{\\partial}_z g^z]_{m,n,p}  $$
 
  or
 
@@ -203,7 +203,7 @@ where `g = [gx, gy, gz]` is a fore- or back-vector field.
 
 Since we applied the forward divergence to the back-vector (and vice-versa), the resulting scalar value
 is defined at the back-vector's (fore-vector's) location $(m,n,p)$ and not at the locations of its components
-$(m \pm \frac{1}{2},n,p)$ etc.
+$(m \\pm \\frac{1}{2},n,p)$ etc.
 
     [figure: divergence]
                                     ^^
@@ -227,23 +227,23 @@ Curls
 
 The two curls are then
 
-  $$ \begin{aligned}
+  $$ \\begin{aligned}
-     \hat{h}_{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &= \\
+     \\hat{h}_{m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2}} &= \\\\
-     [\tilde{\nabla} \times \tilde{g}]_{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &=
+     [\\tilde{\\nabla} \\times \\tilde{g}]_{m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2}} &=
-        \vec{x} (\tilde{\partial}_y g^z_{m,n,p + \frac{1}{2}} - \tilde{\partial}_z g^y_{m,n + \frac{1}{2},p}) \\
+        \\vec{x} (\\tilde{\\partial}_y g^z_{m,n,p + \\frac{1}{2}} - \\tilde{\\partial}_z g^y_{m,n + \\frac{1}{2},p}) \\\\
-     &+ \vec{y} (\tilde{\partial}_z g^x_{m + \frac{1}{2},n,p} - \tilde{\partial}_x g^z_{m,n,p + \frac{1}{2}}) \\
+     &+ \\vec{y} (\\tilde{\\partial}_z g^x_{m + \\frac{1}{2},n,p} - \\tilde{\\partial}_x g^z_{m,n,p + \\frac{1}{2}}) \\\\
-     &+ \vec{z} (\tilde{\partial}_x g^y_{m,n + \frac{1}{2},p} - \tilde{\partial}_y g^z_{m + \frac{1}{2},n,p})
+     &+ \\vec{z} (\\tilde{\\partial}_x g^y_{m,n + \\frac{1}{2},p} - \\tilde{\\partial}_y g^z_{m + \\frac{1}{2},n,p})
-     \end{aligned} $$
+     \\end{aligned} $$
 
  and
 
-  $$ \tilde{h}_{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} =
+  $$ \\tilde{h}_{m - \\frac{1}{2}, n - \\frac{1}{2}, p - \\frac{1}{2}} =
-     [\hat{\nabla} \times \hat{g}]_{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} $$
+     [\\hat{\\nabla} \\times \\hat{g}]_{m - \\frac{1}{2}, n - \\frac{1}{2}, p - \\frac{1}{2}} $$
 
-  where $\hat{g}$ and $\tilde{g}$ are located at $(m,n,p)$
+  where $\\hat{g}$ and $\\tilde{g}$ are located at $(m,n,p)$
-  with components at $(m \pm \frac{1}{2},n,p)$ etc.,
+  with components at $(m \\pm \\frac{1}{2},n,p)$ etc.,
-  while $\hat{h}$ and $\tilde{h}$ are located at $(m \pm \frac{1}{2}, n \pm \frac{1}{2}, p \pm \frac{1}{2})$
+  while $\\hat{h}$ and $\\tilde{h}$ are located at $(m \\pm \\frac{1}{2}, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})$
-  with components at $(m, n \pm \frac{1}{2}, p \pm \frac{1}{2})$ etc.
+  with components at $(m, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})$ etc.
 
 
     [code: curls]
@@ -287,27 +287,27 @@ Maxwell's Equations
 
 If we discretize both space (m,n,p) and time (l), Maxwell's equations become
 
- $$ \begin{aligned}
+ $$ \\begin{aligned}
-  \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &= -\tilde{\partial}_t \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}
+  \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &= -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
-                                                                   - \hat{M}_{l, \vec{r} + \frac{1}{2}}  \\
+                                                                         - \\hat{M}_{l, \\vec{r} + \\frac{1}{2}}  \\\\
-  \hat{\nabla} \times \hat{H}_{l-\frac{1}{2},\vec{r} + \frac{1}{2}} &= \hat{\partial}_t \tilde{D}_{l, \vec{r}}
+  \\hat{\\nabla} \\times \\hat{H}_{l-\\frac{1}{2},\\vec{r} + \\frac{1}{2}} &= \\hat{\\partial}_t \\tilde{D}_{l, \\vec{r}}
-                                                                   + \tilde{J}_{l-\frac{1}{2},\vec{r}} \\
+                                                                             + \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\
-  \tilde{\nabla} \cdot \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} &= 0 \\
+  \\tilde{\\nabla} \\cdot \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= 0 \\\\
-  \hat{\nabla} \cdot \tilde{D}_{l,\vec{r}} &= \rho_{l,\vec{r}}
+  \\hat{\\nabla} \\cdot \\tilde{D}_{l,\\vec{r}} &= \\rho_{l,\\vec{r}}
- \end{aligned} $$
+ \\end{aligned} $$
 
  with
 
- $$ \begin{aligned}
+ $$ \\begin{aligned}
-  \hat{B}_{\vec{r}} &= \mu_{\vec{r} + \frac{1}{2}} \cdot \hat{H}_{\vec{r} + \frac{1}{2}} \\
+  \\hat{B}_{\\vec{r}} &= \\mu_{\\vec{r} + \\frac{1}{2}} \\cdot \\hat{H}_{\\vec{r} + \\frac{1}{2}} \\\\
-  \tilde{D}_{\vec{r}} &= \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}}
+  \\tilde{D}_{\\vec{r}} &= \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}}
- \end{aligned} $$
+ \\end{aligned} $$
 
-where the spatial subscripts are abbreviated as $\vec{r} = (m, n, p)$ and
+where the spatial subscripts are abbreviated as $\\vec{r} = (m, n, p)$ and
-$\vec{r} + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2})$,
+$\\vec{r} + \\frac{1}{2} = (m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2})$,
-$\tilde{E}$ and $\hat{H}$ are the electric and magnetic fields,
+$\\tilde{E}$ and $\\hat{H}$ are the electric and magnetic fields,
-$\tilde{J}$ and $\hat{M}$ are the electric and magnetic current distributions,
+$\\tilde{J}$ and $\\hat{M}$ are the electric and magnetic current distributions,
-and $\epsilon$ and $\mu$ are the dielectric permittivity and magnetic permeability.
+and $\\epsilon$ and $\\mu$ are the dielectric permittivity and magnetic permeability.
 
 The above is Yee's algorithm, written in a form analogous to Maxwell's equations.
 The time derivatives can be expanded to form the update equations:
@@ -369,34 +369,34 @@ Each component forms its own grid, offset from the others:
 
 The divergence equations can be derived by taking the divergence of the curl equations
 and combining them with charge continuity,
- $$ \hat{\nabla} \cdot \tilde{J} + \hat{\partial}_t \rho = 0 $$
+ $$ \\hat{\\nabla} \\cdot \\tilde{J} + \\hat{\\partial}_t \\rho = 0 $$
  implying that the discrete Maxwell's equations do not produce spurious charges.
 
 
 Wave equation
 -------------
 
-Taking the backward curl of the $\tilde{\nabla} \times \tilde{E}$ equation and
+Taking the backward curl of the $\\tilde{\\nabla} \\times \\tilde{E}$ equation and
-replacing the resulting $\hat{\nabla} \times \hat{H}$ term using its respective equation,
+replacing the resulting $\\hat{\\nabla} \\times \\hat{H}$ term using its respective equation,
-and setting $\hat{M}$ to zero, we can form the discrete wave equation:
+and setting $\\hat{M}$ to zero, we can form the discrete wave equation:
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &=
+  \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
-     -\tilde{\partial}_t \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}
+      -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
-                       - \hat{M}_{l-1, \vec{r} + \frac{1}{2}}  \\
+                          - \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}}  \\\\
-  \mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &=
+  \\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
-   -\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}  \\
+    -\\tilde{\\partial}_t \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}  \\\\
-  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
+  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
-   \hat{\nabla} \times (-\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}})  \\
+    \\hat{\\nabla} \\times (-\\tilde{\\partial}_t \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}})  \\\\
-  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
+  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
-   -\tilde{\partial}_t \hat{\nabla} \times \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}  \\
+    -\\tilde{\\partial}_t \\hat{\\nabla} \\times \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}  \\\\
-  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
+  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
-   -\tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \tilde{E}_{l, \vec{r}} + \hat{\partial}_t \tilde{J}_{l-\frac{1}{2},\vec{r}} \\
+    -\\tilde{\\partial}_t \\hat{\\partial}_t \\epsilon_{\\vec{r}} \\tilde{E}_{l, \\vec{r}} + \\hat{\\partial}_t \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\
-  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}})
+  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}})
-           + \tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \cdot \tilde{E}_{l, \vec{r}}
+            + \\tilde{\\partial}_t \\hat{\\partial}_t \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{l, \\vec{r}}
-           &= \tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}}
+            &= \\tilde{\\partial}_t \\tilde{J}_{l - \\frac{1}{2}, \\vec{r}}
-  \end{aligned}
+  \\end{aligned}
 $$
 
 
@@ -406,27 +406,27 @@ Frequency domain
 We can substitute in a time-harmonic fields
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\
+ \\tilde{E}_{l, \\vec{r}} &= \\tilde{E}_{\\vec{r}} e^{-\\imath \\omega l \\Delta_t} \\\\
- \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t}
+ \\tilde{J}_{l, \\vec{r}} &= \\tilde{J}_{\\vec{r}} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t}
- \end{aligned}
+ \\end{aligned}
 $$
 
 resulting in
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \tilde{\partial}_t &\Rightarrow (e^{ \imath \omega \Delta_t} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{-\imath \omega \Delta_t / 2} = -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\
+ \\tilde{\\partial}_t &\\Rightarrow (e^{ \\imath \\omega \\Delta_t} - 1) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_t} \\sin(\\omega \\Delta_t / 2) e^{-\\imath \\omega \\Delta_t / 2} = -\\imath \\Omega e^{-\\imath \\omega \\Delta_t / 2}\\\\
-   \hat{\partial}_t &\Rightarrow (1 - e^{-\imath \omega \Delta_t}) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{ \imath \omega \Delta_t / 2} = -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\
+   \\hat{\\partial}_t &\\Rightarrow (1 - e^{-\\imath \\omega \\Delta_t}) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_t} \\sin(\\omega \\Delta_t / 2) e^{ \\imath \\omega \\Delta_t / 2} = -\\imath \\Omega e^{ \\imath \\omega \\Delta_t / 2}\\\\
- \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t
+ \\Omega &= 2 \\sin(\\omega \\Delta_t / 2) / \\Delta_t
- \end{aligned}
+ \\end{aligned}
 $$
 
 This gives the frequency-domain wave equation,
 
 $$
- \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}})
+ \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}})
-    -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\
+    -\\Omega^2 \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}} = -\\imath \\Omega \\tilde{J}_{\\vec{r}} e^{\\imath \\omega \\Delta_t / 2} \\\\
 $$
 
 
@@ -436,48 +436,48 @@ Plane waves and Dispersion relation
 With uniform material distribution and no sources
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \mu_{\vec{r} + \frac{1}{2}} &= \mu \\
+ \\mu_{\\vec{r} + \\frac{1}{2}} &= \\mu \\\\
- \epsilon_{\vec{r}} &= \epsilon \\
+ \\epsilon_{\\vec{r}} &= \\epsilon \\\\
- \tilde{J}_{\vec{r}} &= 0 \\
+ \\tilde{J}_{\\vec{r}} &= 0 \\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 the frequency domain wave equation simplifies to
 
-$$ \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}} - \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 $$
+$$ \\hat{\\nabla} \\times \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}} - \\Omega^2 \\epsilon \\mu \\tilde{E}_{\\vec{r}} = 0 $$
 
-Since $\hat{\nabla} \cdot \tilde{E}_{\vec{r}} = 0$, we can simplify
+Since $\\hat{\\nabla} \\cdot \\tilde{E}_{\\vec{r}} = 0$, we can simplify
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}}
+ \\hat{\\nabla} \\times \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}}
-  &= \tilde{\nabla}(\hat{\nabla} \cdot \tilde{E}_{\vec{r}}) - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\
+   &= \\tilde{\\nabla}(\\hat{\\nabla} \\cdot \\tilde{E}_{\\vec{r}}) - \\hat{\\nabla} \\cdot \\tilde{\\nabla} \\tilde{E}_{\\vec{r}} \\\\
-  &= - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\
+   &= - \\hat{\\nabla} \\cdot \\tilde{\\nabla} \\tilde{E}_{\\vec{r}} \\\\
-  &= - \tilde{\nabla}^2 \tilde{E}_{\vec{r}}
+   &= - \\tilde{\\nabla}^2 \\tilde{E}_{\\vec{r}}
- \end{aligned}
+ \\end{aligned}
 $$
 
 and we get
 
-$$  \tilde{\nabla}^2 \tilde{E}_{\vec{r}} + \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 $$
+$$  \\tilde{\\nabla}^2 \\tilde{E}_{\\vec{r}} + \\Omega^2 \\epsilon \\mu \\tilde{E}_{\\vec{r}} = 0 $$
 
 We can convert this to three scalar-wave equations of the form
 
-$$ (\tilde{\nabla}^2 + K^2) \phi_{\vec{r}} = 0 $$
+$$ (\\tilde{\\nabla}^2 + K^2) \\phi_{\\vec{r}} = 0 $$
 
-with $K^2 = \Omega^2 \mu \epsilon$. Now we let
+with $K^2 = \\Omega^2 \\mu \\epsilon$. Now we let
 
-$$  \phi_{\vec{r}} = A e^{\imath (k_x m \Delta_x + k_y n \Delta_y + k_z p \Delta_z)}  $$
+$$  \\phi_{\\vec{r}} = A e^{\\imath (k_x m \\Delta_x + k_y n \\Delta_y + k_z p \\Delta_z)}  $$
 
 resulting in
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- \tilde{\partial}_x &\Rightarrow (e^{ \imath k_x \Delta_x} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{ \imath k_x \Delta_x / 2} = \imath K_x e^{ \imath k_x \Delta_x / 2}\\
+ \\tilde{\\partial}_x &\\Rightarrow (e^{ \\imath k_x \\Delta_x} - 1) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_x} \\sin(k_x \\Delta_x / 2) e^{ \\imath k_x \\Delta_x / 2} = \\imath K_x e^{ \\imath k_x \\Delta_x / 2}\\\\
-   \hat{\partial}_x &\Rightarrow (1 - e^{-\imath k_x \Delta_x}) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{-\imath k_x \Delta_x / 2} = \imath K_x e^{-\imath k_x \Delta_x / 2}\\
+   \\hat{\\partial}_x &\\Rightarrow (1 - e^{-\\imath k_x \\Delta_x}) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_x} \\sin(k_x \\Delta_x / 2) e^{-\\imath k_x \\Delta_x / 2} = \\imath K_x e^{-\\imath k_x \\Delta_x / 2}\\\\
- K_x &= 2 \sin(k_x \Delta_x / 2) / \Delta_x \\
+ K_x &= 2 \\sin(k_x \\Delta_x / 2) / \\Delta_x \\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 with similar expressions for the y and z dimnsions (and $K_y, K_z$).
@@ -485,20 +485,20 @@ with similar expressions for the y and z dimnsions (and $K_y, K_z$).
 This implies
 
 $$
-  \tilde{\nabla}^2 = -(K_x^2 + K_y^2 + K_z^2) \phi_{\vec{r}} \\
+  \\tilde{\\nabla}^2 = -(K_x^2 + K_y^2 + K_z^2) \\phi_{\\vec{r}} \\\\
-  K_x^2 + K_y^2 + K_z^2 = \Omega^2 \mu \epsilon = \Omega^2 / c^2
+  K_x^2 + K_y^2 + K_z^2 = \\Omega^2 \\mu \\epsilon = \\Omega^2 / c^2
 $$
 
-where $c = \sqrt{\mu \epsilon}$.
+where $c = \\sqrt{\\mu \\epsilon}$.
 
-Assuming real $(k_x, k_y, k_z), \omega$ will be real only if
+Assuming real $(k_x, k_y, k_z), \\omega$ will be real only if
 
-$$ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) $$
+$$ c^2 \\Delta_t^2 = \\frac{\\Delta_t^2}{\\mu \\epsilon} < 1/(\\frac{1}{\\Delta_x^2} + \\frac{1}{\\Delta_y^2} + \\frac{1}{\\Delta_z^2}) $$
 
-If $\Delta_x = \Delta_y = \Delta_z$, this simplifies to $c \Delta_t < \Delta_x / \sqrt{3}$.
+If $\\Delta_x = \\Delta_y = \\Delta_z$, this simplifies to $c \\Delta_t < \\Delta_x / \\sqrt{3}$.
 This last form can be interpreted as enforcing causality; the distance that light
-travels in one timestep (i.e., $c \Delta_t$) must be less than the diagonal
+travels in one timestep (i.e., $c \\Delta_t$) must be less than the diagonal
-of the smallest cell ( $\Delta_x / \sqrt{3}$ when on a uniform cubic grid).
+of the smallest cell ( $\\Delta_x / \\sqrt{3}$ when on a uniform cubic grid).
 
 
 Grid description
@@ -515,8 +515,8 @@ to make the illustration simpler; we need at least two cells in the x dimension
 demonstrate how nonuniform `dx` affects the various components.
 
 Place the E fore-vectors at integer indices $r = (m, n, p)$ and the H back-vectors
-at fractional indices $r + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2},
+at fractional indices $r + \\frac{1}{2} = (m + \\frac{1}{2}, n + \\frac{1}{2},
-p + \frac{1}{2})$. Remember that these are indices and not coordinates; they can
+p + \\frac{1}{2})$. Remember that these are indices and not coordinates; they can
 correspond to arbitrary (monotonically increasing) coordinates depending on the cell widths.
 
 Draw lines to denote the planes on which the H components and back-vectors are defined.
@@ -718,14 +718,14 @@ composed of the three diagonal tensor components:
 or
 
 $$
- \epsilon = \begin{bmatrix} \epsilon_{xx} & 0 & 0 \\
+ \\epsilon = \\begin{bmatrix} \\epsilon_{xx} & 0 & 0 \\\\
-                            0 & \epsilon_{yy} & 0 \\
+                              0 & \\epsilon_{yy} & 0 \\\\
-                            0 & 0 & \epsilon_{zz} \end{bmatrix}
+                              0 & 0 & \\epsilon_{zz} \\end{bmatrix}
 $$
 $$
- \mu = \begin{bmatrix} \mu_{xx} & 0 & 0 \\
+ \\mu = \\begin{bmatrix} \\mu_{xx} & 0 & 0 \\\\
-                         0 & \mu_{yy} & 0 \\
+                         0 & \\mu_{yy} & 0 \\\\
-                         0 & 0 & \mu_{zz} \end{bmatrix}
+                         0 & 0 & \\mu_{zz} \\end{bmatrix}
 $$
 
 where the off-diagonal terms (e.g. `epsilon_xy`) are assumed to be zero.

meanas/fdmath/functional.py

@@ -62,7 +62,7 @@ def deriv_back(
 def curl_forward(
         dx_e: Sequence[NDArray[numpy.float_]] | None = None,
         ) -> fdfield_updater_t:
-    r"""
+    """
     Curl operator for use with the E field.
 
     Args:
@@ -71,7 +71,7 @@ def curl_forward(
 
     Returns:
         Function `f` for taking the discrete forward curl of a field,
-        `f(E)` -> curlE $= \nabla_f \times E$
+        `f(E)` -> curlE $= \\nabla_f \\times E$
     """
     Dx, Dy, Dz = deriv_forward(dx_e)
 
@@ -91,7 +91,7 @@ def curl_forward(
 def curl_back(
         dx_h: Sequence[NDArray[numpy.float_]] | None = None,
         ) -> fdfield_updater_t:
-    r"""
+    """
     Create a function which takes the backward curl of a field.
 
     Args:
@@ -100,7 +100,7 @@ def curl_back(
 
     Returns:
         Function `f` for taking the discrete backward curl of a field,
-        `f(H)` -> curlH $= \nabla_b \times H$
+        `f(H)` -> curlH $= \\nabla_b \\times H$
     """
     Dx, Dy, Dz = deriv_back(dx_h)
 

meanas/fdtd/__init__.py

@@ -1,4 +1,4 @@
-r"""
+"""
 Utilities for running finite-difference time-domain (FDTD) simulations
 
 See the discussion of `Maxwell's Equations` in `meanas.fdmath` for basic
@@ -11,9 +11,9 @@ Timestep
 From the discussion of "Plane waves and the Dispersion relation" in `meanas.fdmath`,
 we have
 
-$$ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) $$
+$$ c^2 \\Delta_t^2 = \\frac{\\Delta_t^2}{\\mu \\epsilon} < 1/(\\frac{1}{\\Delta_x^2} + \\frac{1}{\\Delta_y^2} + \\frac{1}{\\Delta_z^2}) $$
 
-or, if $\Delta_x = \Delta_y = \Delta_z$, then $c \Delta_t < \frac{\Delta_x}{\sqrt{3}}$.
+or, if $\\Delta_x = \\Delta_y = \\Delta_z$, then $c \\Delta_t < \\frac{\\Delta_x}{\\sqrt{3}}$.
 
 Based on this, we can set
 
@@ -27,81 +27,81 @@ Poynting Vector and Energy Conservation
 
 Let
 
-$$ \begin{aligned}
+$$ \\begin{aligned}
-  \tilde{S}_{l, l', \vec{r}} &=& &\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}}  \\
+  \\tilde{S}_{l, l', \\vec{r}} &=& &\\tilde{E}_{l, \\vec{r}} \\otimes \\hat{H}_{l', \\vec{r} + \\frac{1}{2}}  \\\\
-  &=&  &\vec{x} (\tilde{E}^y_{l,m+1,n,p} \hat{H}^z_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^z_{l,m+1,n,p} \hat{H}^y_{l', \vec{r} + \frac{1}{2}}) \\
+  &=&  &\\vec{x} (\\tilde{E}^y_{l,m+1,n,p} \\hat{H}^z_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^z_{l,m+1,n,p} \\hat{H}^y_{l', \\vec{r} + \\frac{1}{2}}) \\\\
-  & &+ &\vec{y} (\tilde{E}^z_{l,m,n+1,p} \hat{H}^x_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^x_{l,m,n+1,p} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \\
+  & &+ &\\vec{y} (\\tilde{E}^z_{l,m,n+1,p} \\hat{H}^x_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^x_{l,m,n+1,p} \\hat{H}^z_{l', \\vec{r} + \\frac{1}{2}}) \\\\
-  & &+ &\vec{z} (\tilde{E}^x_{l,m,n,p+1} \hat{H}^y_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^y_{l,m,n,p+1} \hat{H}^z_{l', \vec{r} + \frac{1}{2}})
+  & &+ &\\vec{z} (\\tilde{E}^x_{l,m,n,p+1} \\hat{H}^y_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^y_{l,m,n,p+1} \\hat{H}^z_{l', \\vec{r} + \\frac{1}{2}})
-   \end{aligned}
+   \\end{aligned}
 $$
 
-where $\vec{r} = (m, n, p)$ and $\otimes$ is a modified cross product
+where $\\vec{r} = (m, n, p)$ and $\\otimes$ is a modified cross product
-in which the $\tilde{E}$ terms are shifted as indicated.
+in which the $\\tilde{E}$ terms are shifted as indicated.
 
 By taking the divergence and rearranging terms, we can show that
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \hat{\nabla} \cdot \tilde{S}_{l, l', \vec{r}}
+  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l', \\vec{r}}
-   &= \hat{\nabla} \cdot (\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}})  \\
+   &= \\hat{\\nabla} \\cdot (\\tilde{E}_{l, \\vec{r}} \\otimes \\hat{H}_{l', \\vec{r} + \\frac{1}{2}})  \\\\
-   &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l, \vec{r}} -
+   &= \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l, \\vec{r}} -
-      \tilde{E}_{l, \vec{r}} \cdot \hat{\nabla} \times \hat{H}_{l', \vec{r} + \frac{1}{2}} \\
+      \\tilde{E}_{l, \\vec{r}} \\cdot \\hat{\\nabla} \\times \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\\\
-   &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot
+   &= \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\cdot
-          (-\tilde{\partial}_t \mu_{\vec{r} + \frac{1}{2}} \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} -
+            (-\\tilde{\\partial}_t \\mu_{\\vec{r} + \\frac{1}{2}} \\hat{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} -
-              \hat{M}_{l, \vec{r} + \frac{1}{2}}) -
+                \\hat{M}_{l, \\vec{r} + \\frac{1}{2}}) -
-      \tilde{E}_{l, \vec{r}} \cdot (\hat{\partial}_t \tilde{\epsilon}_{\vec{r}} \tilde{E}_{l'+\frac{1}{2}, \vec{r}} +
+      \\tilde{E}_{l, \\vec{r}} \\cdot (\\hat{\\partial}_t \\tilde{\\epsilon}_{\\vec{r}} \\tilde{E}_{l'+\\frac{1}{2}, \\vec{r}} +
-              \tilde{J}_{l', \vec{r}}) \\
+                \\tilde{J}_{l', \\vec{r}}) \\\\
-   &= \hat{H}_{l'} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) -
+   &= \\hat{H}_{l'} \\cdot (-\\mu / \\Delta_t)(\\hat{H}_{l + \\frac{1}{2}} - \\hat{H}_{l - \\frac{1}{2}}) -
-      \tilde{E}_l \cdot (\epsilon / \Delta_t )(\tilde{E}_{l'+\frac{1}{2}} - \tilde{E}_{l'-\frac{1}{2}})
+      \\tilde{E}_l \\cdot (\\epsilon / \\Delta_t )(\\tilde{E}_{l'+\\frac{1}{2}} - \\tilde{E}_{l'-\\frac{1}{2}})
-      - \hat{H}_{l'} \cdot \hat{M}_{l} - \tilde{E}_l \cdot \tilde{J}_{l'} \\
+      - \\hat{H}_{l'} \\cdot \\hat{M}_{l} - \\tilde{E}_l \\cdot \\tilde{J}_{l'} \\\\
-  \end{aligned}
+  \\end{aligned}
 $$
 
 where in the last line the spatial subscripts have been dropped to emphasize
 the time subscripts $l, l'$, i.e.
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\
+  \\tilde{E}_l &= \\tilde{E}_{l, \\vec{r}} \\\\
-  \hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}}  \\
+  \\hat{H}_l &= \\tilde{H}_{l, \\vec{r} + \\frac{1}{2}}  \\\\
-  \tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}}  \\
+  \\tilde{\\epsilon} &= \\tilde{\\epsilon}_{\\vec{r}}  \\\\
-  \end{aligned}
+  \\end{aligned}
 $$
 
 etc.
-For $l' = l + \frac{1}{2}$ we get
+For $l' = l + \\frac{1}{2}$ we get
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}}
+  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}}
-   &= \hat{H}_{l + \frac{1}{2}} \cdot
+   &= \\hat{H}_{l + \\frac{1}{2}} \\cdot
-            (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) -
+            (-\\mu / \\Delta_t)(\\hat{H}_{l + \\frac{1}{2}} - \\hat{H}_{l - \\frac{1}{2}}) -
-      \tilde{E}_l \cdot (\epsilon / \Delta_t)(\tilde{E}_{l+1} - \tilde{E}_l)
+      \\tilde{E}_l \\cdot (\\epsilon / \\Delta_t)(\\tilde{E}_{l+1} - \\tilde{E}_l)
-      - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\
+      - \\hat{H}_{l'} \\cdot \\hat{M}_l - \\tilde{E}_l \\cdot \\tilde{J}_{l + \\frac{1}{2}} \\\\
-   &= (-\mu / \Delta_t)(\hat{H}^2_{l + \frac{1}{2}} - \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}) -
+   &= (-\\mu / \\Delta_t)(\\hat{H}^2_{l + \\frac{1}{2}} - \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}) -
-      (\epsilon / \Delta_t)(\tilde{E}_{l+1} \cdot \tilde{E}_l - \tilde{E}^2_l)
+      (\\epsilon / \\Delta_t)(\\tilde{E}_{l+1} \\cdot \\tilde{E}_l - \\tilde{E}^2_l)
-      - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\
+      - \\hat{H}_{l'} \\cdot \\hat{M}_l - \\tilde{E}_l \\cdot \\tilde{J}_{l + \\frac{1}{2}} \\\\
-   &= -(\mu \hat{H}^2_{l + \frac{1}{2}}
+   &= -(\\mu \\hat{H}^2_{l + \\frac{1}{2}}
-       +\epsilon \tilde{E}_{l+1} \cdot \tilde{E}_l) / \Delta_t \\
+       +\\epsilon \\tilde{E}_{l+1} \\cdot \\tilde{E}_l) / \\Delta_t \\ \\
-      +(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
+      +(\\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}
-       +\epsilon \tilde{E}^2_l) / \Delta_t \\
+       +\\epsilon \\tilde{E}^2_l) / \\Delta_t \\ \\
-      - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
+      - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
-  \end{aligned}
+  \\end{aligned}
 $$
 
-and for $l' = l - \frac{1}{2}$,
+and for $l' = l - \\frac{1}{2}$,
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  \hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}}
+  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}}
-   &=  (\mu \hat{H}^2_{l - \frac{1}{2}}
+   &=  (\\mu \\hat{H}^2_{l - \\frac{1}{2}}
-       +\epsilon \tilde{E}_{l-1} \cdot \tilde{E}_l) / \Delta_t \\
+       +\\epsilon \\tilde{E}_{l-1} \\cdot \\tilde{E}_l) / \\Delta_t \\ \\
-      -(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
+      -(\\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}
-       +\epsilon \tilde{E}^2_l) / \Delta_t \\
+       +\\epsilon \\tilde{E}^2_l) / \\Delta_t \\ \\
-      - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
+      - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
-  \end{aligned}
+  \\end{aligned}
 $$
 
 These two results form the discrete time-domain analogue to Poynting's theorem.
@@ -109,25 +109,25 @@ They hint at the expressions for the energy, which can be calculated at the same
 time-index as either the E or H field:
 
 $$
- \begin{aligned}
+ \\begin{aligned}
- U_l &= \epsilon \tilde{E}^2_l + \mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} \\
+ U_l &= \\epsilon \\tilde{E}^2_l + \\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}} \\\\
- U_{l + \frac{1}{2}} &= \epsilon \tilde{E}_l \cdot \tilde{E}_{l + 1} + \mu \hat{H}^2_{l + \frac{1}{2}} \\
+ U_{l + \\frac{1}{2}} &= \\epsilon \\tilde{E}_l \\cdot \\tilde{E}_{l + 1} + \\mu \\hat{H}^2_{l + \\frac{1}{2}} \\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 Rewriting the Poynting theorem in terms of the energy expressions,
 
 $$
-  \begin{aligned}
+  \\begin{aligned}
-  (U_{l+\frac{1}{2}} - U_l) / \Delta_t
+  (U_{l+\\frac{1}{2}} - U_l) / \\Delta_t
-   &= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\
+   &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}} \\ \\
-      - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
+      - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
-  (U_l - U_{l-\frac{1}{2}}) / \Delta_t
+  (U_l - U_{l-\\frac{1}{2}}) / \\Delta_t
-   &= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\
+   &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}} \\ \\
-      - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
+      - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
- \end{aligned}
+ \\end{aligned}
 $$
 
 This result is exact and should practically hold to within numerical precision. No time-
@@ -147,10 +147,10 @@ of the time-domain fields.
 The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse
 shape. It can be written
 
-$$ f_r(t) = (1 - \frac{1}{2} (\omega (t - \tau))^2) e^{-(\frac{\omega (t - \tau)}{2})^2} $$
+$$ f_r(t) = (1 - \\frac{1}{2} (\\omega (t - \\tau))^2) e^{-(\\frac{\\omega (t - \\tau)}{2})^2} $$
 
-with $\tau > \frac{2 * \pi}{\omega}$ as a minimum delay to avoid a discontinuity at
+with $\\tau > \\frac{2 * \\pi}{\\omega}$ as a minimum delay to avoid a discontinuity at
-t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0).
+t=0 (assuming the source is off for t<0 this gives $\\sim 10^{-3}$ error at t=0).
 
 
 

meanas/fdtd/energy.py

@@ -12,7 +12,7 @@ def poynting(
         h: fdfield_t,
         dxes: dx_lists_t | None = None,
         ) -> fdfield_t:
-    r"""
+    """
     Calculate the poynting vector `S` ($S$).
 
     This is the energy transfer rate (amount of energy `U` per `dt` transferred
@@ -44,16 +44,16 @@ def poynting(
 
     The full relationship is
     $$
-      \begin{aligned}
+      \\begin{aligned}
-      (U_{l+\frac{1}{2}} - U_l) / \Delta_t
+      (U_{l+\\frac{1}{2}} - U_l) / \\Delta_t
-       &= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\
+       &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}} \\ \\
-          - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
+          - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-          - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
+          - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
-      (U_l - U_{l-\frac{1}{2}}) / \Delta_t
+      (U_l - U_{l-\\frac{1}{2}}) / \\Delta_t
-       &= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\
+       &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}} \\ \\
-          - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
+          - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-          - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
+          - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
-     \end{aligned}
+     \\end{aligned}
     $$
 
     These equalities are exact and should practically hold to within numerical precision.

pdoc_templates/html_helpers.py

@@ -321,6 +321,7 @@ class _ToMarkdown:
         """Wrap URLs in Python-Markdown-compatible <angle brackets>."""
         return re.sub(r'(?<![<"\'])(\s*)((?:http|ftp)s?://[^>)\s]+)(\s*)', r'\1<\2>\3', text)
 
+import subprocess
 
 class _MathPattern(InlineProcessor):
     NAME = 'pdoc-math'