[docs] correct Yee calculus, source polarity, and phasor conventions

meanas/fdfd/__init__.py

@@ -19,38 +19,39 @@ Submodules:
 
 ================================================================
 
-From the "Frequency domain" section of `meanas.fdmath`, we have
+Following the "Frequency domain" section of `meanas.fdmath`, we use
+$e^{+\imath \omega t}$ time dependence. Phasors are referenced to the
+actual sample times of the staggered fields and sources:
 
 $$
  \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} \\
+ \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \hat{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 - \frac{1}{2}, \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} \\
+ \hat{M}_{l, \vec{r} + \frac{1}{2}} &= \hat{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}})
-    -\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
  \end{aligned}
 $$
 
-resulting in
+The difference multipliers for an integer-indexed sequence are
 
 $$
  \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}
 $$
 
-Maxwell's equations are then
+The half-step factors cancel when the staggered sample times are included.
+Maxwell's equations for the phasors are then
 
 $$
   \begin{aligned}
   \tilde{\nabla} \times \tilde{E}_{\vec{r}} &=
-         \imath \Omega e^{-\imath \omega \Delta_t / 2} \hat{B}_{\vec{r} + \frac{1}{2}}
+        -\imath \Omega \hat{B}_{\vec{r} + \frac{1}{2}}
                                                      - \hat{M}_{\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 \tilde{D}_{\vec{r}}
                                                      + \tilde{J}_{\vec{r}} \\
   \tilde{\nabla} \cdot \hat{B}_{\vec{r} + \frac{1}{2}} &= 0 \\
   \hat{\nabla} \cdot \tilde{D}_{\vec{r}} &= \rho_{\vec{r}}
@@ -61,13 +62,9 @@ With $\Delta_t \to 0$, this simplifies to
 
 $$
  \begin{aligned}
- \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{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}} \\
  \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}
 $$
 
@@ -76,14 +73,17 @@ and then
 $$
   \begin{aligned}
   \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{\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}} \\
  \end{aligned}
 $$
 
+For zero magnetic current, eliminating $H$ gives the electric wave equation
+used by `operators.e_full` and `functional.e_full`:
+
 $$
  \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}} \\

meanas/fdmath/__init__.py

@@ -116,9 +116,9 @@ 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} +
                                   \vec{y} [\tilde{\partial}_y f]_{m,n + \frac{1}{2},p} +
                                   \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
 
@@ -232,7 +232,7 @@ The two curls are then
      [\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{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^x_{m + \frac{1}{2},n,p})
      \end{aligned} $$
 
  and
@@ -314,7 +314,7 @@ The time derivatives can be expanded to form the update equations:
 
     [code: Maxwell's equations updates]
     H[i, j, k] -= dt * (curl_forward(E)[i, j, k] + M[t, i, j, k]) /      mu[i, j, k]
-    E[i, j, k] += dt * (curl_back(   H)[i, j, k] + J[t, i, j, k]) / epsilon[i, j, k]
+    E[i, j, k] += dt * (curl_back(   H)[i, j, k] - J[t, i, j, k]) / epsilon[i, j, k]
 
 Note that the E-field fore-vector and H-field back-vector are offset by a half-cell, resulting
 in distinct locations for all six E- and H-field components:
@@ -384,7 +384,7 @@ $$
   \begin{aligned}
   \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-1, \vec{r} + \frac{1}{2}}  \\
+                       - \hat{M}_{l, \vec{r} + \frac{1}{2}}  \\
   \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}}  \\
   \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
@@ -392,10 +392,10 @@ $$
   \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}}  \\
   \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}} - \tilde{\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}})
            + \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}
 $$
 
@@ -403,37 +403,42 @@ $$
 Frequency domain
 ----------------
 
-We can substitute in a time-harmonic fields
+We use the $e^{+\imath \omega t}$ convention, matching the FDFD operators.
+Substitute time-harmonic fields with phases referenced to their actual sample times:
 
 $$
  \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-\frac{1}{2}, \vec{r}} &= \tilde{J}_{\vec{r}} e^{+\imath \omega (l - \frac{1}{2}) \Delta_t}
  \end{aligned}
 $$
 
-resulting in
+For an integer-indexed sequence proportional to $e^{+\imath \omega l \Delta_t}$,
+the forward and backward difference multipliers are
 
 $$
  \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
  \end{aligned}
 $$
 
-This gives the frequency-domain wave equation,
+The half-step factors account for the shift in the derivative's sample time.
+In particular, $\tilde{\partial}_t J_{l-\frac{1}{2}} = \imath \Omega J e^{\imath \omega l \Delta_t}$.
+Thus, when the staggered fields are expressed using the phasors above, the
+frequency-domain wave equation has no extra half-step phase:
 
 $$
  \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}} \\
 $$
 
 
 Plane waves and Dispersion relation
 ------------------------------------
 
-With uniform material distribution and no sources
+With uniform grid spacings, a uniform material distribution, and no sources
 
 $$
  \begin{aligned}
@@ -474,31 +479,31 @@ resulting in
 
 $$
  \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_x = \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_x = \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 \\
  \end{aligned}
 $$
 
-with similar expressions for the y and z dimnsions (and $K_y, K_z$).
+with similar expressions for the y and z dimensions (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 \phi_{\vec{r}} = -(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
 $$
 
-where $c = \sqrt{\mu \epsilon}$.
+where $c = 1 / \sqrt{\mu \epsilon}$.
 
-Assuming real $(k_x, k_y, k_z), \omega$ will be real only if
+To keep $\omega$ real for every real grid wavevector $(k_x, k_y, k_z)$,
+we use the conservative stability bound
 
 $$ 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}$.
-This last form can be interpreted as enforcing causality; the distance that light
+This limits the distance that light travels in one timestep to less than
-travels in one timestep (i.e., $c \Delta_t$) must be less than the diagonal
+$\Delta_x / \sqrt{3}$ on a uniform cubic grid.
-of the smallest cell ( $\Delta_x / \sqrt{3}$ when on a uniform cubic grid).
 
 
 Grid description

meanas/fdtd/base.py

@@ -17,13 +17,13 @@ def maxwell_e(
     """
     Build a function which performs a portion the time-domain E-field update,
 
-        E += curl_back(H[t]) / epsilon
+        E += dt * curl_back(H[t]) / epsilon
 
     The full update should be
 
-        E += (curl_back(H[t]) + J) / epsilon
+        E += dt * (curl_back(H[t]) - J) / epsilon
 
-    which requires an additional step of `E += J / epsilon` which is not performed
+    which requires an additional step of `E -= dt * J / epsilon` which is not performed
     by the generated function.
 
     See `meanas.fdmath` for descriptions of
@@ -72,13 +72,13 @@ def maxwell_h(
     """
     Build a function which performs part of the time-domain H-field update,
 
-        H -= curl_forward(E[t]) / mu
+        H -= dt * curl_forward(E[t]) / mu
 
     The full update should be
 
-        H -= (curl_forward(E[t]) + M) / mu
+        H -= dt * (curl_forward(E[t]) + M) / mu
 
-    which requires an additional step of `H -= M / mu` which is not performed
+    which requires an additional step of `H -= dt * M / mu` which is not performed
     by the generated function; this step can be omitted if there is no magnetic
     current `M`.