From 693e3af8fa767294f79f822885df2f7a580ae4c8 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Sun, 19 Apr 2026 16:40:30 -0700 Subject: [PATCH] Publish docs for local build --- api/fdfd/index.html | 95 ++-- api/fdmath/index.html | 203 ++++--- api/fdtd/index.html | 503 ++++++++++++++-- api/waveguides/index.html | 192 ++++--- index.html | 20 +- objects.inv | Bin 1378 -> 1475 bytes print_page/index.html | 1135 +++++++++++++++++++++++++------------ search/search_index.json | 2 +- standalone.html | 709 +++++++++++++++-------- 9 files changed, 2042 insertions(+), 817 deletions(-) diff --git a/api/fdfd/index.html b/api/fdfd/index.html index 07b78ef..43fe48f 100644 --- a/api/fdfd/index.html +++ b/api/fdfd/index.html @@ -1345,58 +1345,63 @@ matrix equation (Ax=b) or eigenvalue problem.

From the "Frequency domain" section of meanas.fdmath, we have

\[ \begin{aligned} - \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{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{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{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} \\ \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 + -\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

\[ \begin{aligned} - \tilde{\partial}_t &\Rightarrow -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\ - \hat{\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}\\ \end{aligned} \]
+

Maxwell's equations are then

\[ \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}} - \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}} + \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}} + \tilde{\nabla} \cdot \hat{B}_{\vec{r} + \frac{1}{2}} &= 0 \\ + \hat{\nabla} \cdot \tilde{D}_{\vec{r}} &= \rho_{\vec{r}} \end{aligned} \]
+

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 \\ - \hat{\partial}_t &\to -\imath \omega \\ + \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 \\ + \hat{\partial}_t &\to -\imath \omega \\ \end{aligned} \]
+

and then

\[ \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}} - \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}} + \tilde{J}_{\vec{r}} \\ \end{aligned} \]
+
\[ \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}} \\ @@ -2064,6 +2069,7 @@ mask selecting the total-field region, then the TFSF source is the commutator

\[ \frac{A Q - Q A}{-i \omega} E. \]
+

This vanishes in the interior of the total-field and scattered-field regions and is supported only at their shared boundary, where the mask discontinuity makes A and Q fail to commute. The returned current is therefore the @@ -2416,11 +2422,15 @@ the meanas.fdmath.types submodule for details.

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

To make this matrix symmetric, use the preconditioners from e_full_preconditioners().

@@ -2678,11 +2688,15 @@ The PMC is applied per-field-component (i.e. pmc.size == epsilon.size

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
@@ -2855,25 +2869,30 @@ The PMC is applied per-field-component (i.e. pmc.size == epsilon.size

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} - -\imath \omega \epsilon & \nabla \times \ - \nabla \times & \imath \omega \mu - \end{bmatrix} $$

+ equations without cancelling out either E or H. The operator is

+
\[ +\begin{bmatrix} + -\imath \omega \epsilon & \nabla \times \\ + \nabla \times & \imath \omega \mu +\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} - -\imath \omega \epsilon & \nabla \times \ - \nabla \times & \imath \omega \mu +

for use with a field vector of the form cat(vec(E), vec(H)):

+
\[ +\begin{bmatrix} + -\imath \omega \epsilon & \nabla \times \\ + \nabla \times & \imath \omega \mu \end{bmatrix} - \begin{bmatrix} E \ + \begin{bmatrix} E \\ H \end{bmatrix} - = \begin{bmatrix} J \ + = \begin{bmatrix} J \\ -M - \end{bmatrix} $$

+\end{bmatrix} +\]

Parameters:

@@ -3408,6 +3427,7 @@ usual antisymmetry of the cross product,

\[ H \times E = -(E \times H), \]
+

once the same staggered field placement is used on both sides.

@@ -3527,6 +3547,7 @@ Then the TFSF current operator is the commutator

\[ \frac{A Q - Q A}{-i \omega}. \]
+

Inside regions where Q is locally constant, A and Q commute and the source vanishes. Only cells at the TF/SF boundary contribute nonzero current, which is exactly the desired distributed source for injecting the chosen diff --git a/api/fdmath/index.html b/api/fdmath/index.html index 61228f3..d687628 100644 --- a/api/fdmath/index.html +++ b/api/fdmath/index.html @@ -1903,8 +1903,9 @@ series of other matrices).

which covers a superset of this material with similar notation and more detail.

Scalar derivatives and cell shifts

Define the discrete forward derivative as - $$ [\tilde{\partial}x f] - f_m) $$ - where }{2}} = \frac{1}{\Delta_{x, m}} (f_{m + 1\(f\) is a function defined at discrete locations on the x-axis (labeled using \(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\) 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.

@@ -1913,8 +1914,9 @@ along the x-axis, the forward derivative is

deriv_forward(f)[i] = (f[i + 1] - f[i]) / dx[i]

Likewise, discrete reverse derivative is - $$ [\hat{\partial}x f ]) $$ - or}{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 @@ -1949,7 +1951,9 @@ Periodic boundaries are used here and elsewhere unless otherwise noted. etc.

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}\); 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, @@ -1961,12 +1965,18 @@ See the Grid description section below for additional information o and generalization to three dimensions.

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] + - \vec{y} [\tilde{\partial}}{2},n,py f] + - \vec{z} [\tilde{\partial}}{2},pz f] $$ - $$ [\hat{\nabla} f]}{2}{m,n,p} = \vec{x} [\hat{\partial}_x f] + - \vec{y} [\hat{\partial}}{2},n,py f] + - \vec{z} [\hat{\partial}}{2},pz f] $$}{2}

+
+
\[ + [\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} + + \vec{y} [\hat{\partial}_y f]_{m,n + \frac{1}{2},p} + + \vec{z} [\hat{\partial}_z f]_{m,n,p + \frac{1}{2}} + \]
+

or

[code: gradients]
 grad_forward(f)[i,j,k] = [Dx_forward(f)[i, j, k],
@@ -1989,12 +1999,18 @@ at different points: the x-component is shifted in the x-direction,
 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 +
+

+
\[
+  \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{z} g^z_{m,n,p + \frac{1}{2}} $$
-  $$ \hat{g}}{2},n,p{m,n,p} = \vec{x} g^x +
+  \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} +
                        \vec{y} g^y_{m,n - \frac{1}{2},p} +
-                       \vec{z} g^z_{m,n,p - \frac{1}{2}} $$}{2},n,p

+  \vec{z} g^z_{m,n,p - \frac{1}{2}}
+  \]

+
 
[figure: gradient / fore-vector]
    (m, n+1, p+1) ______________ (m+1, n+1, p+1)
                 /:            /|
@@ -2010,14 +2026,20 @@ on the direction of the shift. We write it as
 

 
Divergences

 
There are also two divergences,

-
$$ d_{n,m,p} = [\tilde{\nabla} \cdot \hat{g}]{n,m,p}
-               = [\tilde{\partial}_x g^x] +
-                 [\tilde{\partial}y g^y] +
-                 [\tilde{\partial}z g^z]   $$

-
$$ d_{n,m,p} = [\hat{\nabla} \cdot \tilde{g}]{n,m,p}
-               = [\hat{\partial}_x g^x] +
-                 [\hat{\partial}y g^y] +
-                 [\hat{\partial}z g^z]  $$

+
\[
+  d_{n,m,p} = [\tilde{\nabla} \cdot \hat{g}]_{n,m,p}
+               = [\tilde{\partial}_x g^x]_{m,n,p} +
+                 [\tilde{\partial}_y g^y]_{m,n,p} +
+  [\tilde{\partial}_z g^z]_{m,n,p}
+  \]

+
+
\[
+  d_{n,m,p} = [\hat{\nabla} \cdot \tilde{g}]_{n,m,p}
+               = [\hat{\partial}_x g^x]_{m,n,p} +
+                 [\hat{\partial}_y g^y]_{m,n,p} +
+  [\hat{\partial}_z g^z]_{m,n,p}
+  \]

+
 
or

 
[code: divergences]
 div_forward(g)[i,j,k] = Dx_forward(gx)[i, j, k] +
@@ -2056,16 +2078,22 @@ is defined at the back-vector's (fore-vector's) location Curls
 
The two curls are then

-
$$ \begin{aligned}
-     \hat{h}{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &= \
-     [\tilde{\nabla} \times \tilde{g}] &=
-        \vec{x} (\tilde{\partial}}{2}, n + \frac{1}{2}, p + \frac{1}{2}y g^z}{2}} - \tilde{\partialz g^y) \
-     &+ \vec{y} (\tilde{\partial}}{2},pz g^x}{2},n,p} - \tilde{\partialx g^z) \
-     &+ \vec{z} (\tilde{\partial}}{2}x g^y}{2},p} - \tilde{\partialy g^z)
-     \end{aligned} $$}{2},n,p

+
\[
+  \begin{aligned}
+     \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}} &=
+        \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})
+  \end{aligned}
+  \]

+
 
and

-
$$ \tilde{h}{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} =
-     [\hat{\nabla} \times \hat{g}] $$}{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}}
+  \]

+
 
where \(\hat{g}\) and \(\tilde{g}\) are located at \((m,n,p)\)
   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})\)
@@ -2103,19 +2131,25 @@ curl_back(g)[i,j,k] = [Dy_back(gz)[i, j, k] - Dz_back(gy)[i, j, k],
 

 
Maxwell's Equations

 
If we discretize both space (m,n,p) and time (l), Maxwell's equations become

-
$$ \begin{aligned}
-  \tilde{\nabla} \times \tilde{E}{l,\vec{r}} &= -\tilde{\partial}_t \hat{B}
-                                                                   - \hat{M}}{2}, \vec{r} + \frac{1}{2}{l, \vec{r} + \frac{1}{2}}  \
-  \hat{\nabla} \times \hat{H}}{2},\vec{r} + \frac{1}{2}} &= \hat{\partialt \tilde{D}
-                                                                   + \tilde{J}}{l-\frac{1}{2},\vec{r}} \
-  \tilde{\nabla} \cdot \hat{B} &= 0 \
-  \hat{\nabla} \cdot \tilde{D}}{2}, \vec{r} + \frac{1}{2}{l,\vec{r}} &= \rho
- \end{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}}
+                                                                   - \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}}
+                                                                   + \tilde{J}_{l-\frac{1}{2},\vec{r}} \\
+  \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}}
+ \end{aligned}
+ \]

+
 
with

-
$$ \begin{aligned}
-  \hat{B}{\vec{r}} &= \mu} + \frac{1}{2}} \cdot \hat{H{\vec{r} + \frac{1}{2}} \
-  \tilde{D}
- \end{aligned} $$}} &= \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}

+
\[
+ \begin{aligned}
+  \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}}
+ \end{aligned}
+ \]

+
 
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})\),
 \(\tilde{E}\) and \(\hat{H}\) are the electric and magnetic fields,
@@ -2177,94 +2211,108 @@ r - 1/2 =    | /                      | /

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 $$ - implying that the discrete Maxwell's equations do not produce spurious charges.

+ +
\[ \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 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:

\[ \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}} - \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}} \\ - \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 (\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}} \\ - \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}} \\ \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} \]
+

Frequency domain

We can substitute in a time-harmonic fields

\[ \begin{aligned} - \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{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} \end{aligned} \]
+

resulting in

\[ \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}\\ - \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 + \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}\\ + \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t \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}}) -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\ \]
+

Plane waves and Dispersion relation

With uniform material distribution and no sources

\[ \begin{aligned} - \mu_{\vec{r} + \frac{1}{2}} &= \mu \\ - \epsilon_{\vec{r}} &= \epsilon \\ - \tilde{J}_{\vec{r}} &= 0 \\ + \mu_{\vec{r} + \frac{1}{2}} &= \mu \\ + \epsilon_{\vec{r}} &= \epsilon \\ + \tilde{J}_{\vec{r}} &= 0 \\ \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 \]
+

Since \(\hat{\nabla} \cdot \tilde{E}_{\vec{r}} = 0\), we can simplify

\[ \begin{aligned} \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}} \\ - &= - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\ - &= - \tilde{\nabla}^2 \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}} \\ + &= - \tilde{\nabla}^2 \tilde{E}_{\vec{r}} \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 \]
+

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} - \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}\\ - K_x &= 2 \sin(k_x \Delta_x / 2) / \Delta_x \\ + \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}\\ + 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\)).

This implies

\[ \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 \]
+

where \(c = \sqrt{\mu \epsilon}\).

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}\). 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 @@ -2461,15 +2509,16 @@ mu = [mu_xx, mu_yy, mu_zz]

or

\[ - \epsilon = \begin{bmatrix} \epsilon_{xx} & 0 & 0 \\ - 0 & \epsilon_{yy} & 0 \\ - 0 & 0 & \epsilon_{zz} \end{bmatrix} -$$ -$$ - \mu = \begin{bmatrix} \mu_{xx} & 0 & 0 \\ - 0 & \mu_{yy} & 0 \\ - 0 & 0 & \mu_{zz} \end{bmatrix} + \epsilon = \begin{bmatrix} \epsilon_{xx} & 0 & 0 \\ + 0 & \epsilon_{yy} & 0 \\ + 0 & 0 & \epsilon_{zz} \end{bmatrix} \]
+
\[ + \mu = \begin{bmatrix} \mu_{xx} & 0 & 0 \\ + 0 & \mu_{yy} & 0 \\ + 0 & 0 & \mu_{zz} \end{bmatrix} +\]
+

where the off-diagonal terms (e.g. epsilon_xy) are assumed to be zero.

High-accuracy volumetric integration of shapes on multiple grids can be performed by the gridlock module.

diff --git a/api/fdtd/index.html b/api/fdtd/index.html index 6511297..2b544da 100644 --- a/api/fdtd/index.html +++ b/api/fdtd/index.html @@ -809,6 +809,50 @@ + + +
  • + + + +  temporal_phasor + + + + +
    • + +
  • + + + +  temporal_phasor_scale + + + + +
    • + +
  • + + + +  real_injection_scale + + + + +
    • + +
  • + + + +  reconstruct_real + + + +
  • @@ -842,6 +886,39 @@ +
    • + +
  • + + + +  reconstruct_real_e + + + + +
    • + +
  • + + + +  reconstruct_real_h + + + + +
    • + +
  • + + + +  reconstruct_real_j + + + +
    • @@ -1223,6 +1300,50 @@ + + +
  • + + + +  temporal_phasor + + + + +
    • + +
  • + + + +  temporal_phasor_scale + + + + +
    • + +
  • + + + +  real_injection_scale + + + + +
    • + +
  • + + + +  reconstruct_real + + + +
  • @@ -1256,6 +1377,39 @@ +
    • + +
  • + + + +  reconstruct_real_e + + + + +
    • + +
  • + + + +  reconstruct_real_h + + + + +
    • + +
  • + + + +  reconstruct_real_j + + + +
    • @@ -1301,7 +1455,8 @@ mathematical background.

    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}}\).

    Based on this, we can set

    dt = sqrt(mu.min() * epsilon.min()) / sqrt(1/dx_min**2 + 1/dy_min**2 + 1/dz_min**2)
@@ -1309,54 +1464,58 @@ we have
    
 
    The dx_min, dy_min, dz_min should be the minimum value across both the base and derived grids.
    
 
    Poynting Vector and Energy Conservation
    
 
    Let
    
-
    \[ \begin{aligned}
-  \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{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}})
+
    \[
+\begin{aligned}
+  \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{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}})
    \end{aligned}
 \]
    
+
 
    where \(\vec{r} = (m, n, p)\) and \(\otimes\) is a modified cross product
 in which the \(\tilde{E}\) terms are shifted as indicated.
    
 
    By taking the divergence and rearranging terms, we can show that
    
 
    \[
   \begin{aligned}
   \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{H}_{l', \vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l, \vec{r}} -
+   &= \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}} -
       \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}} -
               \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{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}})
       - \hat{H}_{l'} \cdot \hat{M}_{l} - \tilde{E}_l \cdot \tilde{J}_{l'} \\
   \end{aligned}
 \]
    
+
 
    where in the last line the spatial subscripts have been dropped to emphasize
 the time subscripts \(l, l'\), i.e.
    
 
    \[
   \begin{aligned}
-  \tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\
-  \hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}}  \\
-  \tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}}  \\
+  \tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\
+  \hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}}  \\
+  \tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}}  \\
   \end{aligned}
 \]
    
+
 
    etc.
 For \(l' = l + \frac{1}{2}\) we get
    
 
    \[
   \begin{aligned}
   \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}}) -
       \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}} \\
-   &= (-\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)
       - \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 \\
       +(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
        +\epsilon \tilde{E}^2_l) / \Delta_t \\
@@ -1364,11 +1523,12 @@ For \(l' = l + \frac{1}{2}\) we get
    
       - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
   \end{aligned}
 \]
    
+
 
    and for \(l' = l - \frac{1}{2}\),
    
 
    \[
   \begin{aligned}
   \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 \\
       -(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
        +\epsilon \tilde{E}^2_l) / \Delta_t \\
@@ -1376,28 +1536,31 @@ For \(l' = l + \frac{1}{2}\) we get
    
       - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
   \end{aligned}
 \]
    
+
 
    These two results form the discrete time-domain analogue to Poynting's theorem.
 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}
- 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 &= \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}} \\
  \end{aligned}
 \]
    
+
 
    Rewriting the Poynting theorem in terms of the energy expressions,
    
 
    \[
   \begin{aligned}
   (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 \\
       - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
   (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 \\
       - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
  \end{aligned}
 \]
    
+
 
    This result is exact and should practically hold to within numerical precision. No time-
 or spatial-averaging is necessary.
    
 
    Note that each value of \(J\) contributes to the energy twice (i.e. once per field update)
@@ -1408,16 +1571,46 @@ extract the frequency-domain response by performing an on-the-fly Fourier transf
 of the time-domain fields.
    
 
    accumulate_phasor in meanas.fdtd.phasor performs the phase accumulation for one
 or more target frequencies, while leaving source normalization and simulation-loop
-policy to the caller. Convenience wrappers accumulate_phasor_e,
-accumulate_phasor_h, and accumulate_phasor_j apply the usual Yee time offsets.
+policy to the caller. temporal_phasor(...) and temporal_phasor_scale(...)
+apply the same Fourier sum to a scalar waveform, which is useful when a pulsed
+source envelope must be normalized before being applied to a point source or
+mode source. real_injection_scale(...) is the corresponding helper for the
+common real-valued injection pattern numpy.real(scale * waveform). Convenience
+wrappers accumulate_phasor_e, accumulate_phasor_h, and accumulate_phasor_j
+apply the usual Yee time offsets. reconstruct_real(...) and the corresponding
+reconstruct_real_e/h/j(...) wrappers perform the inverse operation, converting
+phasors back into real-valued field snapshots at explicit Yee-aligned times.
+For scalar omega, the reconstruction helpers accept either a plain field
+phasor or the batched (1, *sample_shape) form used by the accumulator helpers.
 The helpers accumulate
    
 
    \[ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l \]
    
+
 
    with caller-provided sample times and weights. In this codebase the matching
 electric-current convention is typically E -= dt * J / epsilon in FDTD and
 -i \omega J on the right-hand side of the FDFD wave equation.
    
+
    For FDTD/FDFD equivalence, this phasor path is the primary comparison workflow.
+It isolates the guided +\omega response that the frequency-domain solver
+targets directly, regardless of whether the underlying FDTD run used real- or
+complex-valued fields.
    
+
    For exact pulsed FDTD/FDFD equivalence it is often simplest to keep the
+injected source, fields, and CPML auxiliary state complex-valued. The
+real_injection_scale(...) helper is instead for the more ordinary one-run
+real-valued source path, where the intended positive-frequency waveform is
+injected through numpy.real(scale * waveform) and any remaining negative-
+frequency leakage is controlled by the pulse bandwidth and run length.
    
+
    reconstruct_real(...) is for a different question: given a phasor, what late
+real-valued field snapshot should it produce? That raw-snapshot comparison is
+stricter and noisier because a monitor plane generally contains both the guided
+field and the remaining orthogonal content,
    
+
    \[ E_{\text{monitor}} = E_{\text{guided}} + E_{\text{residual}} . \]
    
+
+
    Phasor/modal comparisons mostly validate the guided +\omega term. Raw
+real-field comparisons expose both terms at once, so they should be treated as
+secondary diagnostics rather than the main solver-equivalence benchmark.
    
 
    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} \]
    
+
 
    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).
    
 
    Boundary conditions
    
@@ -1447,6 +1640,7 @@ coordinate system in both places:
    
 
    \[
 E \leftarrow E - \Delta_t J / \epsilon
 \]
    
+
 
    which matches the FDFD right-hand side
    
 
    \[
 -i \omega J.
@@ -2390,19 +2584,20 @@ boundaries of the corresponding U cell. For example,
    
     assert_close(sum(planes), du_half_h2e[mask])

    (see meanas.tests.test_fdtd.test_poynting_planes)

    -

    The full relationship is -$$ +

    The full relationship is

    +
    \[ \begin{aligned} (U_{l+\frac{1}{2}} - U_l) / \Delta_t - &= -\hat{\nabla} \cdot \tilde{S}{l, l + \frac{1}{2}} \ - - \hat{H}}{2}} \cdot \hat{Ml \ - - \tilde{E}_l \cdot \tilde{J} \ + &= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\ + - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\ + - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\ (U_l - U_{l-\frac{1}{2}}) / \Delta_t - &= -\hat{\nabla} \cdot \tilde{S}}{2}{l, l - \frac{1}{2}} \ - - \hat{H}}{2}} \cdot \hat{Ml \ - - \tilde{E}_l \cdot \tilde{J} \ + &= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\ + - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\ + - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\ \end{aligned} -$$}{2}

    +\]
    +

    These equalities are exact and should practically hold to within numerical precision. No time- or spatial-averaging is necessary. (See meanas.fdtd docs for derivation.)

    @@ -3579,6 +3774,18 @@ the sampling policy. The accumulated quantity is

  • accumulate_phasor_h(..., step=l) for H_{l + 1/2}
  • accumulate_phasor_j(..., step=l) for J_{l + 1/2}
    • +

    temporal_phasor(...) and temporal_phasor_scale(...) apply the same Fourier +sum to a 1D scalar waveform. They are useful for normalizing a pulsed source +before that scalar waveform is applied to a point source or spatial mode source. +real_injection_scale(...) is a companion helper for the common real-valued +injection pattern numpy.real(scale * waveform), where waveform is the +analytic positive-frequency signal and the injected real current should still +produce a desired phasor response. +reconstruct_real(...) and its E/H/J wrappers perform the inverse operation: +they turn one or more phasors back into real-valued field snapshots at explicit +Yee-aligned sample times. For a scalar target frequency they accept either a +plain field phasor or the batched (1, *sample_shape) form used elsewhere in +this module.

    These helpers do not choose warmup/accumulation windows or pulse-envelope normalization. They also do not impose a current sign convention. In this codebase, electric-current injection normally appears as E -= dt * J / epsilon, @@ -3650,6 +3857,154 @@ factor was derived from a discrete sum that already includes dt, pa

    +

    + temporal_phasor + + +

    +
    temporal_phasor(
+    samples: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    *,
+    start_step: int = 0,
+    offset_steps: float = 0.0,
+) -> NDArray[numpy.complexfloating]
+
    + +
    + +

    Fourier-project a 1D temporal waveform onto one or more angular frequencies.

    +

    The returned quantity is

    +
    dt * sum(exp(-1j * omega * t_step) * samples[step_index])
+
    +

    where t_step = (start_step + step_index + offset_steps) * dt.

    + + +
    + +
    + +
    + + +

    + temporal_phasor_scale + + +

    +
    temporal_phasor_scale(
+    samples: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    *,
+    start_step: int = 0,
+    offset_steps: float = 0.0,
+    target: ArrayLike = 1.0,
+) -> NDArray[numpy.complexfloating]
+
    + +
    + +

    Return the scalar multiplier that gives a desired temporal phasor response.

    +

    The returned scale satisfies

    +
    temporal_phasor(scale * samples, omegas, dt, ...) == target
+
    +

    for each target frequency. The result keeps a leading frequency axis even +when omegas is scalar.

    + + +
    + +
    + +
    + + +

    + real_injection_scale + + +

    +
    real_injection_scale(
+    samples: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    *,
+    start_step: int = 0,
+    offset_steps: float = 0.0,
+    target: ArrayLike = 1.0,
+) -> NDArray[numpy.complexfloating]
+
    + +
    + +

    Return the scale for a real-valued injection built from an analytic waveform.

    +

    If the time-domain source is applied as

    +
    numpy.real(scale * samples[step])
+
    +

    then the desired positive-frequency phasor is obtained by compensating for +the 1/2 factor between the real-valued source and its analytic component:

    +
    scale = 2 * target / temporal_phasor(samples, ...)
+
    +

    This helper normalizes only the intended positive-frequency component. Any +residual negative-frequency leakage is controlled by the waveform design and +the accumulation window.

    + + +
    + +
    + +
    + + +

    + reconstruct_real + + +

    +
    reconstruct_real(
+    phasors: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    step: int,
+    *,
+    offset_steps: float = 0.0,
+) -> NDArray[numpy.floating]
+
    + +
    + +

    Reconstruct a real-valued field snapshot from one or more phasors.

    +

    The returned quantity is

    +
    real(phasor * exp(1j * omega * t_step))
+
    +

    where t_step = (step + offset_steps) * dt.

    +

    For multi-frequency inputs, the leading frequency axis is preserved. For a +scalar omega, callers may pass either (1, *sample_shape) or +sample_shape; the return shape matches that choice.

    + + +
    + +
    + +
    + +

    accumulate_phasor_e @@ -3740,6 +4095,90 @@ factor was derived from a discrete sum that already includes dt, pa

    +
    + + +

    + reconstruct_real_e + + +

    +
    reconstruct_real_e(
+    phasors: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    step: int,
+) -> NDArray[numpy.floating]
+
    + +
    + +

    Reconstruct a real E-field snapshot taken at integer timestep step.

    + + +
    + +
    + +
    + + +

    + reconstruct_real_h + + +

    +
    reconstruct_real_h(
+    phasors: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    step: int,
+) -> NDArray[numpy.floating]
+
    + +
    + +

    Reconstruct a real H-field snapshot corresponding to H_{step + 1/2}.

    + + +
    + +
    + +
    + + +

    + reconstruct_real_j + + +

    +
    reconstruct_real_j(
+    phasors: ArrayLike,
+    omegas: float
+    | complex
+    | Sequence[float | complex]
+    | NDArray,
+    dt: float,
+    step: int,
+) -> NDArray[numpy.floating]
+
    + +
    + +

    Reconstruct a real current snapshot corresponding to J_{step + 1/2}.

    + + +
    + +
    + diff --git a/api/waveguides/index.html b/api/waveguides/index.html index 64260c7..468141f 100644 --- a/api/waveguides/index.html +++ b/api/waveguides/index.html @@ -1402,132 +1402,144 @@ a structure with some (x, y) cross-section extending uniformly into the z dimens with a diagonal \(\epsilon\) tensor, we have

    \[ \begin{aligned} -\nabla \times \vec{E}(x, y, z) &= -\imath \omega \mu \vec{H} \\ -\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^{-\imath \beta z} \\ -\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\imath \beta z} \\ +\nabla \times \vec{E}(x, y, z) &= -\imath \omega \mu \vec{H} \\ +\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^{-\imath \beta z} \\ +\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\imath \beta z} \\ \end{aligned} \]
    +

    Expanding the first two equations into vector components, we get

    \[ \begin{aligned} --\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_{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_{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 \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_{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_{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 \\ \end{aligned} \]
    +

    Substituting in our expressions for \(\vec{E}\), \(\vec{H}\) and discretizing:

    \[ \begin{aligned} --\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta E_y \\ --\imath \omega \mu_{yy} H_y &= -\imath \beta 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 \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta H_y \\ -\imath \omega \epsilon_{yy} E_y &= -\imath \beta 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 \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta E_y \\ +-\imath \omega \mu_{yy} H_y &= -\imath \beta 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 \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta H_y \\ +\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \hat{\partial}_x H_z \\ +\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ \end{aligned} \]
    +

    Rewrite the last three equations as

    \[ \begin{aligned} -\imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ -\imath \beta 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 \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ +\imath \beta 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 \\ \end{aligned} \]
    +

    Now apply \(\imath \beta \tilde{\partial}_x\) to the last equation, then substitute in for \(\imath \beta H_x\) and \(\imath \beta H_y\):

    \[ \begin{aligned} -\imath \beta \tilde{\partial}_x \imath \omega E_z &= \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y +\imath \beta \tilde{\partial}_x \imath \omega E_z &= \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \imath \beta \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}_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) \\ -\imath \beta \tilde{\partial}_x E_z &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) +\imath \beta \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) \\ \end{aligned} \]
    +

    With a similar approach (but using \(\imath \beta \tilde{\partial}_y\) instead), we can get

    \[ \begin{aligned} -\imath \beta \tilde{\partial}_y E_z &= \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) +\imath \beta \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) \\ \end{aligned} \]
    +

    We can combine this equation for \(\imath \beta \tilde{\partial}_y E_z\) with the unused \(\imath \omega \mu_{xx} H_x\) and \(\imath \omega \mu_{yy} H_y\) equations to get

    \[ \begin{aligned} --\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \imath \beta \tilde{\partial}_y E_z \\ --\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \tilde{\partial}_y ( +-\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \imath \beta \tilde{\partial}_y E_z \\ +-\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \tilde{\partial}_y ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ \end{aligned} \]
    +

    and

    \[ \begin{aligned} --\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \imath \beta \tilde{\partial}_x E_z \\ --\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \tilde{\partial}_x ( +-\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \imath \beta \tilde{\partial}_x E_z \\ +-\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \tilde{\partial}_x ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ \end{aligned} \]
    +

    However, based on our rewritten equation for \(\imath \beta H_x\) and the so-far unused equation for \(\imath \omega \mu_{zz} H_z\) we can also write

    \[ \begin{aligned} --\imath \omega \mu_{xx} (\imath \beta 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 + \imath \omega \mu_{xx} \hat{\partial}_x ( +-\imath \omega \mu_{xx} (\imath \beta 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 + \imath \omega \mu_{xx} \hat{\partial}_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) \\ \end{aligned} \]
    +

    and, similarly,

    \[ \begin{aligned} --\imath \omega \mu_{yy} (\imath \beta H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +-\imath \omega \mu_{yy} (\imath \beta 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) \\ \end{aligned} \]
    +

    By combining both pairs of expressions, we get

    \[ \begin{aligned} \beta^2 E_x - \tilde{\partial}_x ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \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) \\ -\beta^2 E_y + \tilde{\partial}_y ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \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) \\ \end{aligned} \]
    +

    Using these, we can construct the eigenvalue problem

    \[ \beta^2 \begin{bmatrix} E_x \\ E_y \end{bmatrix} = - (\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ - 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + + (\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ + 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\ \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 \\ \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 \\ E_y \end{bmatrix} \]
    +

    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.

    @@ -1580,15 +1592,16 @@ mu * [[-Dy], [Dx]] / mu * [-Dy, Dx] +

    for use with a field vector of the form cat([E_x, E_y]).

    More precisely, the operator is

    \[ -\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ - 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + +\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ + 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\ \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 \\ \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, and each \(\epsilon_{xx}\), \(\mu_{yy}\), etc. is a diagonal matrix containing the vectorized material property distribution.

    @@ -1735,15 +1748,16 @@ epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +

    for use with a field vector of the form cat([H_x, H_y]).

    More precisely, the operator is

    \[ -\omega^2 \begin{bmatrix} \epsilon_{yy} \mu_{xx} & 0 \\ - 0 & \epsilon_{xx} \mu_{yy} \end{bmatrix} + +\omega^2 \begin{bmatrix} \epsilon_{yy} \mu_{xx} & 0 \\ + 0 & \epsilon_{xx} \mu_{yy} \end{bmatrix} + \begin{bmatrix} -\epsilon_{yy} \tilde{\partial}_y \\ \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 \\ \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, and each \(\epsilon_{xx}\), \(\mu_{yy}\), etc. is a diagonal matrix containing the vectorized material property distribution.

    @@ -2066,6 +2080,7 @@ and exy2h(...), then normalizes them to unit forward power using
    \[ \Re\left[\mathrm{inner\_product}(e, h, \mathrm{conj\_h}=True)\right] = 1, \]
    +

    so the returned fields represent a unit-power forward mode under the discrete Yee-grid Poynting inner product.

    @@ -2721,21 +2736,23 @@ identical representatives for nearly symmetric modes.

    \[ \imath \omega \epsilon_{zz} E_z = \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ \]
    +

    as well as the intermediate equations

    \[ \begin{aligned} -\imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ -\imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ +\imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ +\imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ \end{aligned} \]
    +

    Combining these, we get

    \[ \begin{aligned} -E_z &= \frac{1}{- \omega \beta \epsilon_{zz}} (( +E_z &= \frac{1}{- \omega \beta \epsilon_{zz}} (( \hat{\partial}_y \hat{\partial}_x H_z -\hat{\partial}_x \hat{\partial}_y H_z) + \imath \omega (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y)) - &= \frac{1}{\imath \beta \epsilon_{zz}} (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y) + &= \frac{1}{\imath \beta \epsilon_{zz}} (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y) \end{aligned} \]
    @@ -3654,10 +3671,12 @@ E_z &= \frac{1}{- \omega \beta \epsilon_{zz}} (( \(\beta\) to changes in the dielectric structure \(\epsilon\).

    The output is a vector of the same size as vec(epsilon), with each element specifying the sensitivity of wavenumber to changes in the corresponding element in vec(epsilon), i.e.

    -
    \[sens_{i} = \frac{\partial\beta}{\partial\epsilon_i}\]
    +
    \[ sens_{i} = \frac{\partial\beta}{\partial\epsilon_i} \]
    +

    An adjoint approach is used to calculate the sensitivity; the derivation is provided here:

    Starting with the eigenvalue equation

    -
    \[\beta^2 E_{xy} = A_E E_{xy}\]
    +
    \[ \beta^2 E_{xy} = A_E E_{xy} \]
    +

    where \(A_E\) is the waveguide operator from operator_e(), and \(E_{xy} = \begin{bmatrix} E_x \\ E_y \end{bmatrix}\), we can differentiate with respect to one of the \(\epsilon\) elements (i.e. at one Yee grid point), \(\epsilon_i\):

    @@ -3665,11 +3684,13 @@ we can differentiate with respect to one of the \(\epsi (2 \beta) \partial_{\epsilon_i}(\beta) E_{xy} + \beta^2 \partial_{\epsilon_i} E_{xy} = \partial_{\epsilon_i}(A_E) E_{xy} + A_E \partial_{\epsilon_i} E_{xy} \] +

    We then multiply by \(H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}\) from the left:

    \[ (2 \beta) \partial_{\epsilon_i}(\beta) H_{yx}^\star E_{xy} + \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy} = H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} \]
    +

    However, \(H_{yx}^\star\) is actually a left-eigenvector of \(A_E\). This can be verified by inspecting the form of operator_h (\(A_H\)) and comparing its conjugate transpose to operator_e (\(A_E\)). Also, note \(H_{yx}^\star \cdot E_{xy} = H^\star \times E\) recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201 @@ -3677,11 +3698,13 @@ for a similar approach. Therefore,

    \[ H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy} \]
    +

    and we can simplify to

    \[ \partial_{\epsilon_i}(\beta) = \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}} \]
    +

    This expression can be quickly calculated for all \(i\) by writing out the various terms of \(\partial_{\epsilon_i} A_E\) and recognizing that the vector-matrix-vector products (i.e. scalars) \(sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}\), indexed by \(i\), can be expressed as @@ -4161,6 +4184,7 @@ longitudinal Poynting flux,

    \[ \frac{1}{2}\int (E_x H_y - E_y H_x) \, dx \, dy \]
    +

    with the Yee-grid staggering and optional propagation-phase adjustment used by the waveguide helpers in this module.

    @@ -4866,6 +4890,7 @@ same sign convention used elsewhere in the package:

    \[ \sum \mathrm{overlap\_e} \; E_\mathrm{other}^* \]
    +

    where the sum is over the full Yee-grid field storage.

    The construction uses a two-cell window immediately upstream of the selected slice:

    @@ -4881,6 +4906,7 @@ overlap relation involving

    \[ \int ((E \times H_\mathrm{mode}) + (E_\mathrm{mode} \times H)) \cdot dn. \]
    +

    E x H_mode + E_mode x H -> Ex Hmy - EyHmx + Emx Hy - Emy Hx (Z-prop) Ex we/B Emx + Ex i/B dy Hmz - Ey (-we/B Emy) - Ey i/B dx Hmz @@ -5219,6 +5245,7 @@ along the propagation axis and applies the phase factor

    \[ e^{-i \, \mathrm{polarity} \, wavenumber \, \Delta z} \]
    +

    to each copied slice.

    @@ -5254,11 +5281,13 @@ and the fields are assumed to have the form

    \[ \vec{E}(r, y, \theta), \vec{H}(r, y, \theta) \propto e^{-\imath m \theta}, \]
    +

    where m is the angular wavenumber returned by solve_mode(s). It is often convenient to introduce the corresponding linear wavenumber

    \[ \beta = \frac{m}{r_{\min}}, \]
    +

    so that the cylindrical problem resembles the straight-waveguide problem with additional metric factors.

    Those metric factors live on the staggered radial Yee grids. If the left edge of @@ -5266,33 +5295,36 @@ the computational window is at r = r_{\min}, define the electric-gr magnetic-grid radial sample locations by

    \[ \begin{aligned} -r_a(n) &= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\ -r_b\!\left(n + \tfrac{1}{2}\right) &= r_{\min} + \tfrac{1}{2}\Delta r_{e, n} - + \sum_{j < n} \Delta r_{h, j}, +r_a(n) &= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\ +r_b\!\left(n + \tfrac{1}{2}\right) &= r_{\min} + \tfrac{1}{2}\Delta r_{e, n} + + \sum_{j < n} \Delta r_{h, j}, \end{aligned} \]
    +

    and from them the diagonal metric matrices

    \[ \begin{aligned} -T_a &= \operatorname{diag}(r_a / r_{\min}), \\ -T_b &= \operatorname{diag}(r_b / r_{\min}). +T_a &= \operatorname{diag}(r_a / r_{\min}), \\ +T_b &= \operatorname{diag}(r_b / r_{\min}). \end{aligned} \]
    +

    With the same forward/backward derivative notation used in waveguide_2d, the coordinate-transformed discrete curl equations used here are

    \[ \begin{aligned} --\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\ --\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r +-\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\ +-\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\ --\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\ -\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\ -\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y +-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\ +\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\ +\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\ -\imath \omega E_z &= T_a \epsilon_{zz}^{-1} +\imath \omega E_z &= T_a \epsilon_{zz}^{-1} \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right). \end{aligned} \]
    +

    The first three equations are the cylindrical analogue of the straight-guide relations for H_r, H_y, and H_z. The next two are the metric-weighted versions of the straight-guide identities for \imath \beta H_y and @@ -5306,6 +5338,7 @@ and then eliminate H_z with

    H_z = \frac{1}{-\imath \omega \mu_{zz}} \left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right). \] +

    This yields the transverse electric eigenproblem implemented by cylindrical_operator(...):

    \[ @@ -5315,8 +5348,8 @@ H_z = \frac{1}{-\imath \omega \mu_{zz}} \left( \omega^2 \begin{bmatrix} -T_b^2 \mu_{yy} \epsilon_{xx} & 0 \\ -0 & T_a^2 \mu_{xx} \epsilon_{yy} +T_b^2 \mu_{yy} \epsilon_{xx} & 0 \\ +0 & T_a^2 \mu_{xx} \epsilon_{yy} \end{bmatrix} + \begin{bmatrix} @@ -5325,7 +5358,7 @@ T_b^2 \mu_{yy} \epsilon_{xx} & 0 \\ \end{bmatrix} T_b \mu_{zz}^{-1} \begin{bmatrix} --\tilde{\partial}_y & \tilde{\partial}_x +-\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + \begin{bmatrix} @@ -5334,12 +5367,13 @@ T_b \mu_{zz}^{-1} \end{bmatrix} T_a \epsilon_{zz}^{-1} \begin{bmatrix} -\hat{\partial}_x T_b \epsilon_{xx} & +\hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix} \right) \begin{bmatrix} E_r \\ E_y \end{bmatrix}. \]
    +

    Since \beta = m / r_{\min}, the solver implemented in this file returns the angular wavenumber m, while the operator itself is most naturally written in terms of the linear quantity \beta. The helpers below reconstruct the full @@ -5389,17 +5423,18 @@ product used by waveguide_2d.

    Cylindrical coordinate waveguide operator of the form

    \[ - (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\ - 0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} + + (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\ + 0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} + \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\ T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \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 \\ \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1} - \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix}) + \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix}) \begin{bmatrix} E_r \\ E_y \end{bmatrix} \]
    +

    for use with a field vector of the form [E_r, E_y].

    This operator can be used to form an eigenvalue problem of the form A @ [E_r, E_y] = beta**2 * [E_r, E_y]

    @@ -6288,15 +6323,15 @@ from E_r and E_y in order to then calculate E_z.

    Returns an operator which, when applied to a vectorized E eigenfield, produces the vectorized H eigenfield.

    -

    This operator is created directly from the initial coordinate-transformed equations: -$$ +

    This operator is created directly from the initial coordinate-transformed equations:

    +
    \[ \begin{aligned} --\imath \omega \mu_{rr} H_r &= \tilde{\partial}y E_z + \imath \beta T_a^{-1} E_y, \ --\imath \omega \mu E_r - - T_b^{-1} \tilde{\partial}} H_y &= -\imath \beta T_b^{-1r (T_a E_z), \ --\imath \omega \mu_y E_r, +-\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\ +-\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r + - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\ +-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \end{aligned} -$$} H_z &= \tilde{\partial}_r E_y - \tilde{\partial

    +\]

    Parameters:

    @@ -6745,6 +6780,7 @@ enforces unit forward power under the discrete inner product

    \[ \frac{1}{2}\int (E_r H_y^* - E_y H_r^*) \, dr \, dy. \]
    +

    The angular wavenumber m is first converted into the full three-component fields, then the overall complex phase and sign are fixed so the result is reproducible for symmetric modes.

    diff --git a/index.html b/index.html index 9b144f6..fbce33a 100644 --- a/index.html +++ b/index.html @@ -647,10 +647,26 @@ conventions.

    For most users, the tracked examples under examples/ are the right entry point. They show the intended combinations of tools for solving complete problems.

    +

    Relevant starting examples:

    + +

    For solver equivalence, prefer the phasor-based examples first. They compare +the extracted +\omega content of the FDTD run directly against the FDFD +solution and are the main accuracy benchmarks in the test suite.

    +

    examples/waveguide_real.py answers a different, stricter question: how well a +late raw real snapshot matches Re(E_\omega e^{i\omega t}) on a monitor plane. +That diagnostic is useful, but it also includes orthogonal residual structure +that the phasor comparison intentionally filters out.

    The API pages are better read as a toolbox map and derivation reference: