jan/meanas
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

Publish docs for local build

Browse source
This commit is contained in:
Forgejo Actions 2026-04-19 16:40:30 -07:00
commit 693e3af8fa
9 changed files with 2042 additions and 817 deletions
Download patch file Download diff file Expand all files Collapse all files

95
api/fdfd/index.html
View file

@ -1345,58 +1345,63 @@ matrix equation (Ax=b) or eigenvalue problem.</p>
<p>From the "Frequency domain" section of <code>meanas.fdmath</code>, we have</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{E}_{l, \vec{r}} &amp;= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\
\tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &amp;= \tilde{H}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\
\tilde{J}_{l, \vec{r}} &amp;= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\
\tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &amp;= \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}} &amp;= -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\
\Omega &amp;= 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}
\]</div>
<p>resulting in</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{\partial}_t &amp;\Rightarrow -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\
\hat{\partial}_t &amp;\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}
\]</div>
<p>Maxwell's equations are then</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{\nabla} \times \tilde{E}_{\vec{r}} &amp;=
\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}} &amp;=
\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}} &amp;= 0 \\
\hat{\nabla} \cdot \tilde{D}_{\vec{r}} &amp;= \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}
\]</div>
<p>With <span class="arithmatex">\(\Delta_t \to 0\)</span>, this simplifies to</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{E}_{l, \vec{r}} &amp;\to \tilde{E}_{\vec{r}} \\
\tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &amp;\to \tilde{H}_{\vec{r} + \frac{1}{2}} \\
\tilde{J}_{l, \vec{r}} &amp;\to \tilde{J}_{\vec{r}} \\
\tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &amp;\to \tilde{M}_{\vec{r} + \frac{1}{2}} \\
\Omega &amp;\to \omega \\
\tilde{\partial}_t &amp;\to -\imath \omega \\
\hat{\partial}_t &amp;\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}
\]</div>
<p>and then</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{\nabla} \times \tilde{E}_{\vec{r}} &amp;=
\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}} &amp;=
\hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &=
-\imath \omega \tilde{D}_{\vec{r}}
+ \tilde{J}_{\vec{r}} \\
\end{aligned}
\]</div>
<div class="arithmatex">\[
\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</p
<div class="arithmatex">\[
\frac{A Q - Q A}{-i \omega} E.
\]</div>
<p>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 <code>A</code> and <code>Q</code> fail to commute. The returned current is therefore the
@ -2416,11 +2422,15 @@ the <code>meanas.fdmath.types</code> submodule for details.</p>
<div class="doc doc-contents ">
<p>Wave operator
$$ \nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon $$</p>
<div class="arithmatex">\[ \nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon \]</div>
</p>
<div class="highlight"><pre><span></span><code>del x (1/mu * del x) - omega**2 * epsilon
</code></pre></div>
<p>for use with the E-field, with wave equation
$$ (\nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon) E = -\imath \omega J $$</p>
<div class="arithmatex">\[ (\nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon) E = -\imath \omega J \]</div>
</p>
<div class="highlight"><pre><span></span><code>(del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
</code></pre></div>
<p>To make this matrix symmetric, use the preconditioners from <code>e_full_preconditioners()</code>.</p>
@ -2678,11 +2688,15 @@ The PMC is applied per-field-component (i.e. <code>pmc.size == epsilon.size</cod
<div class="doc doc-contents ">
<p>Wave operator
$$ \nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu $$</p>
<div class="arithmatex">\[ \nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu \]</div>
</p>
<div class="highlight"><pre><span></span><code>del x (1/epsilon * del x) - omega**2 * mu
</code></pre></div>
<p>for use with the H-field, with wave equation
$$ (\nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu) E = \imath \omega M $$</p>
<div class="arithmatex">\[ (\nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu) E = \imath \omega M \]</div>
</p>
<div class="highlight"><pre><span></span><code>(del x (1/epsilon * del x) - omega**2 * mu) E = i * omega * M
</code></pre></div>
@ -2855,25 +2869,30 @@ The PMC is applied per-field-component (i.e. <code>pmc.size == epsilon.size</cod
<div class="doc doc-contents ">
<p>Wave operator for <code>[E, H]</code> field representation. This operator implements Maxwell's
equations without cancelling out either E or H. The operator is
$$ \begin{bmatrix}
-\imath \omega \epsilon &amp; \nabla \times \
\nabla \times &amp; \imath \omega \mu
\end{bmatrix} $$</p>
equations without cancelling out either E or H. The operator is</p>
<div class="arithmatex">\[
\begin{bmatrix}
-\imath \omega \epsilon & \nabla \times \\
\nabla \times & \imath \omega \mu
\end{bmatrix}
\]</div>
<div class="highlight"><pre><span></span><code>[[-i * omega * epsilon, del x ],
[del x, i * omega * mu]]
</code></pre></div>
<p>for use with a field vector of the form <code>cat(vec(E), vec(H))</code>:
$$ \begin{bmatrix}
-\imath \omega \epsilon &amp; \nabla \times \
\nabla \times &amp; \imath \omega \mu
<p>for use with a field vector of the form <code>cat(vec(E), vec(H))</code>:</p>
<div class="arithmatex">\[
\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} $$</p>
\end{bmatrix}
\]</div>
<p><span class="doc-section-title">Parameters:</span></p>
@ -3408,6 +3427,7 @@ usual antisymmetry of the cross product,</p>
<div class="arithmatex">\[
H \times E = -(E \times H),
\]</div>
<p>once the same staggered field placement is used on both sides.</p>
@ -3527,6 +3547,7 @@ Then the TFSF current operator is the commutator</p>
<div class="arithmatex">\[
\frac{A Q - Q A}{-i \omega}.
\]</div>
<p>Inside regions where <code>Q</code> is locally constant, <code>A</code> and <code>Q</code> 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

203
api/fdmath/index.html
View file

@ -1903,8 +1903,9 @@ series of other matrices).</p>
which covers a superset of this material with similar notation and more detail.</p>
<h4 id="meanas.fdmath--scalar-derivatives-and-cell-shifts">Scalar derivatives and cell shifts<a class="headerlink" href="#meanas.fdmath--scalar-derivatives-and-cell-shifts" title="Permanent link">&para;</a></h4>
<p>Define the discrete forward derivative as
$$ [\tilde{\partial}<em _="+" _frac_1="\frac{1" m="m">x f]</em> - f_m) $$
where }{2}} = \frac{1}{\Delta_{x, m}} (f_{m + 1<span class="arithmatex">\(f\)</span> is a function defined at discrete locations on the x-axis (labeled using <span class="arithmatex">\(m\)</span>).
<div class="arithmatex">\[ [\tilde{\partial}_x f]_{m + \frac{1}{2}} = \frac{1}{\Delta_{x, m}} (f_{m + 1} - f_m) \]</div></p>
<p>where <span class="arithmatex">\(f\)</span> is a function defined at discrete locations on the x-axis (labeled using <span class="arithmatex">\(m\)</span>).
The value at <span class="arithmatex">\(m\)</span> occupies a length <span class="arithmatex">\(\Delta_{x, m}\)</span> along the x-axis. Note that <span class="arithmatex">\(m\)</span>
is an index along the x-axis, <em>not</em> necessarily an x-coordinate, since each length
<span class="arithmatex">\(\Delta_{x, m}, \Delta_{x, m+1}, ...\)</span> is independently chosen.</p>
@ -1913,8 +1914,9 @@ along the x-axis, the forward derivative is</p>
<div class="highlight"><pre><span></span><code>deriv_forward(f)[i] = (f[i + 1] - f[i]) / dx[i]
</code></pre></div>
<p>Likewise, discrete reverse derivative is
$$ [\hat{\partial}<em -="-" _frac_1="\frac{1" m="m">x f ]</em>) $$
or}{2}} = \frac{1}{\Delta_{x, m}} (f_{m} - f_{m - 1</p>
<div class="arithmatex">\[ [\hat{\partial}_x f ]_{m - \frac{1}{2}} = \frac{1}{\Delta_{x, m}} (f_{m} - f_{m - 1}) \]</div></p>
<p>or</p>
<div class="highlight"><pre><span></span><code>deriv_back(f)[i] = (f[i] - f[i - 1]) / dx[i]
</code></pre></div>
<p>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.</p>
<p>The fractional subscript <span class="arithmatex">\(m + \frac{1}{2}\)</span> is used to indicate values defined
at shifted locations relative to the original <span class="arithmatex">\(m\)</span>, with corresponding lengths
$$ \Delta_{x, m + \frac{1}{2}} = \frac{1}{2} * (\Delta_{x, m} + \Delta_{x, m + 1}) $$</p>
<div class="arithmatex">\[ \Delta_{x, m + \frac{1}{2}} = \frac{1}{2} * (\Delta_{x, m} + \Delta_{x, m + 1}) \]</div>
</p>
<p>Just as <span class="arithmatex">\(m\)</span> is not itself an x-coordinate, neither is <span class="arithmatex">\(m + \frac{1}{2}\)</span>;
carefully note the positions of the various cells in the above figure vs their labels.
If the positions labeled with <span class="arithmatex">\(m\)</span> are considered the "base" or "original" grid,
@ -1961,12 +1965,18 @@ See the <code>Grid description</code> section below for additional information o
and generalization to three dimensions.</p>
<h4 id="meanas.fdmath--gradients-and-fore-vectors">Gradients and fore-vectors<a class="headerlink" href="#meanas.fdmath--gradients-and-fore-vectors" title="Permanent link">&para;</a></h4>
<p>Expanding to three dimensions, we can define two gradients
$$ [\tilde{\nabla} f]<em _="+" _frac_1="\frac{1" m="m">{m,n,p} = \vec{x} [\tilde{\partial}_x f]</em> +
\vec{y} [\tilde{\partial}}{2},n,p<em _="+" _frac_1="\frac{1" m_n="m,n">y f]</em> +
\vec{z} [\tilde{\partial}}{2},p<em _="+" _frac_1="\frac{1" m_n_p="m,n,p">z f]</em> $$
$$ [\hat{\nabla} f]}{2}<em _="+" _frac_1="\frac{1" m="m">{m,n,p} = \vec{x} [\hat{\partial}_x f]</em> +
\vec{y} [\hat{\partial}}{2},n,p<em _="+" _frac_1="\frac{1" m_n="m,n">y f]</em> +
\vec{z} [\hat{\partial}}{2},p<em _="+" _frac_1="\frac{1" m_n_p="m,n,p">z f]</em> $$}{2}</p>
<br />
<div class="arithmatex">\[
[\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}}
\]</div></p>
<div class="arithmatex">\[
[\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}}
\]</div>
<p>or</p>
<div class="highlight"><pre><span></span><code>[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.</p>
<p>We call the resulting object a "fore-vector" or "back-vector", depending
on the direction of the shift. We write it as
$$ \tilde{g}<em _="+" _frac_1="\frac{1" m="m">{m,n,p} = \vec{x} g^x</em> +
<br />
<div class="arithmatex">\[
\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<em -="-" _frac_1="\frac{1" m="m">{m,n,p} = \vec{x} g^x</em> +
\vec{z} g^z_{m,n,p + \frac{1}{2}}
\]</div></p>
<div class="arithmatex">\[
\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</p>
\vec{z} g^z_{m,n,p - \frac{1}{2}}
\]</div>
<div class="highlight"><pre><span></span><code>[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
</code></pre></div>
<h4 id="meanas.fdmath--divergences">Divergences<a class="headerlink" href="#meanas.fdmath--divergences" title="Permanent link">&para;</a></h4>
<p>There are also two divergences,</p>
<p>$$ d_{n,m,p} = [\tilde{\nabla} \cdot \hat{g}]<em m_n_p="m,n,p">{n,m,p}
= [\tilde{\partial}_x g^x]</em> +
[\tilde{\partial}<em m_n_p="m,n,p">y g^y]</em> +
[\tilde{\partial}<em m_n_p="m,n,p">z g^z]</em> $$</p>
<p>$$ d_{n,m,p} = [\hat{\nabla} \cdot \tilde{g}]<em m_n_p="m,n,p">{n,m,p}
= [\hat{\partial}_x g^x]</em> +
[\hat{\partial}<em m_n_p="m,n,p">y g^y]</em> +
[\hat{\partial}<em m_n_p="m,n,p">z g^z]</em> $$</p>
<div class="arithmatex">\[
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}
\]</div>
<div class="arithmatex">\[
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}
\]</div>
<p>or</p>
<div class="highlight"><pre><span></span><code>[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 <span class="arithmatex
</code></pre></div>
<h4 id="meanas.fdmath--curls">Curls<a class="headerlink" href="#meanas.fdmath--curls" title="Permanent link">&para;</a></h4>
<p>The two curls are then</p>
<p>$$ \begin{aligned}
\hat{h}<em _="+" _frac_1="\frac{1" m="m">{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &amp;= \
[\tilde{\nabla} \times \tilde{g}]</em> &amp;=
\vec{x} (\tilde{\partial}}{2}, n + \frac{1}{2}, p + \frac{1}{2}<em _="+" _frac_1="\frac{1" m_n_p="m,n,p">y g^z</em>}{2}} - \tilde{\partial<em _="+" _frac_1="\frac{1" m_n="m,n">z g^y</em>) \
&amp;+ \vec{y} (\tilde{\partial}}{2},p<em _="+" _frac_1="\frac{1" m="m">z g^x</em>}{2},n,p} - \tilde{\partial<em _="+" _frac_1="\frac{1" m_n_p="m,n,p">x g^z</em>) \
&amp;+ \vec{z} (\tilde{\partial}}{2}<em _="+" _frac_1="\frac{1" m_n="m,n">x g^y</em>}{2},p} - \tilde{\partial<em _="+" _frac_1="\frac{1" m="m">y g^z</em>)
\end{aligned} $$}{2},n,p</p>
<div class="arithmatex">\[
\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}
\]</div>
<p>and</p>
<p>$$ \tilde{h}<em -="-" _frac_1="\frac{1" m="m">{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} =
[\hat{\nabla} \times \hat{g}]</em> $$}{2}, n - \frac{1}{2}, p - \frac{1}{2}</p>
<div class="arithmatex">\[
\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}}
\]</div>
<p>where <span class="arithmatex">\(\hat{g}\)</span> and <span class="arithmatex">\(\tilde{g}\)</span> are located at <span class="arithmatex">\((m,n,p)\)</span>
with components at <span class="arithmatex">\((m \pm \frac{1}{2},n,p)\)</span> etc.,
while <span class="arithmatex">\(\hat{h}\)</span> and <span class="arithmatex">\(\tilde{h}\)</span> are located at <span class="arithmatex">\((m \pm \frac{1}{2}, n \pm \frac{1}{2}, p \pm \frac{1}{2})\)</span>
@ -2103,19 +2131,25 @@ curl_back(g)[i,j,k] = [Dy_back(gz)[i, j, k] - Dz_back(gy)[i, j, k],
</code></pre></div>
<h3 id="meanas.fdmath--maxwells-equations">Maxwell's Equations<a class="headerlink" href="#meanas.fdmath--maxwells-equations" title="Permanent link">&para;</a></h3>
<p>If we discretize both space (m,n,p) and time (l), Maxwell's equations become</p>
<p>$$ \begin{aligned}
\tilde{\nabla} \times \tilde{E}<em l-_frac_1="l-\frac{1">{l,\vec{r}} &amp;= -\tilde{\partial}_t \hat{B}</em>
- \hat{M}}{2}, \vec{r} + \frac{1}{2}<em l-_frac_1="l-\frac{1">{l, \vec{r} + \frac{1}{2}} \
\hat{\nabla} \times \hat{H}</em>}{2},\vec{r} + \frac{1}{2}} &amp;= \hat{\partial<em _vec_r="\vec{r" l_="l,">t \tilde{D}</em>
+ \tilde{J}}<em l-_frac_1="l-\frac{1">{l-\frac{1}{2},\vec{r}} \
\tilde{\nabla} \cdot \hat{B}</em> &amp;= 0 \
\hat{\nabla} \cdot \tilde{D}}{2}, \vec{r} + \frac{1}{2}<em l_vec_r="l,\vec{r">{l,\vec{r}} &amp;= \rho</em>
\end{aligned} $$}</p>
<div class="arithmatex">\[
\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}
\]</div>
<p>with</p>
<p>$$ \begin{aligned}
\hat{B}<em _vec_r="\vec{r">{\vec{r}} &amp;= \mu</em>} + \frac{1}{2}} \cdot \hat{H<em _vec_r="\vec{r">{\vec{r} + \frac{1}{2}} \
\tilde{D}</em>
\end{aligned} $$}} &amp;= \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}</p>
<div class="arithmatex">\[
\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}
\]</div>
<p>where the spatial subscripts are abbreviated as <span class="arithmatex">\(\vec{r} = (m, n, p)\)</span> and
<span class="arithmatex">\(\vec{r} + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2})\)</span>,
<span class="arithmatex">\(\tilde{E}\)</span> and <span class="arithmatex">\(\hat{H}\)</span> are the electric and magnetic fields,
@ -2177,94 +2211,108 @@ r - 1/2 = | / | /
</code></pre></div>
<p>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.</p>
<div class="arithmatex">\[ \hat{\nabla} \cdot \tilde{J} + \hat{\partial}_t \rho = 0 \]</div></p>
<p>implying that the discrete Maxwell's equations do not produce spurious charges.</p>
<h4 id="meanas.fdmath--wave-equation">Wave equation<a class="headerlink" href="#meanas.fdmath--wave-equation" title="Permanent link">&para;</a></h4>
<p>Taking the backward curl of the <span class="arithmatex">\(\tilde{\nabla} \times \tilde{E}\)</span> equation and
replacing the resulting <span class="arithmatex">\(\hat{\nabla} \times \hat{H}\)</span> term using its respective equation,
and setting <span class="arithmatex">\(\hat{M}\)</span> to zero, we can form the discrete wave equation:</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &amp;=
\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}} &amp;=
\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}}) &amp;=
\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}}) &amp;=
\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}}) &amp;=
\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}}
&amp;= \tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}}
&= \tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}}
\end{aligned}
\]</div>
<h4 id="meanas.fdmath--frequency-domain">Frequency domain<a class="headerlink" href="#meanas.fdmath--frequency-domain" title="Permanent link">&para;</a></h4>
<p>We can substitute in a time-harmonic fields</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{E}_{l, \vec{r}} &amp;= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\
\tilde{J}_{l, \vec{r}} &amp;= \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}
\]</div>
<p>resulting in</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{\partial}_t &amp;\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 &amp;\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 &amp;= 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}
\]</div>
<p>This gives the frequency-domain wave equation,</p>
<div class="arithmatex">\[
\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} \\
\]</div>
<h4 id="meanas.fdmath--plane-waves-and-dispersion-relation">Plane waves and Dispersion relation<a class="headerlink" href="#meanas.fdmath--plane-waves-and-dispersion-relation" title="Permanent link">&para;</a></h4>
<p>With uniform material distribution and no sources</p>
<div class="arithmatex">\[
\begin{aligned}
\mu_{\vec{r} + \frac{1}{2}} &amp;= \mu \\
\epsilon_{\vec{r}} &amp;= \epsilon \\
\tilde{J}_{\vec{r}} &amp;= 0 \\
\mu_{\vec{r} + \frac{1}{2}} &= \mu \\
\epsilon_{\vec{r}} &= \epsilon \\
\tilde{J}_{\vec{r}} &= 0 \\
\end{aligned}
\]</div>
<p>the frequency domain wave equation simplifies to</p>
<div class="arithmatex">\[ \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}} - \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 \]</div>
<p>Since <span class="arithmatex">\(\hat{\nabla} \cdot \tilde{E}_{\vec{r}} = 0\)</span>, we can simplify</p>
<div class="arithmatex">\[
\begin{aligned}
\hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}}
&amp;= \tilde{\nabla}(\hat{\nabla} \cdot \tilde{E}_{\vec{r}}) - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\
&amp;= - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\
&amp;= - \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}
\]</div>
<p>and we get</p>
<div class="arithmatex">\[ \tilde{\nabla}^2 \tilde{E}_{\vec{r}} + \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 \]</div>
<div class="arithmatex">\[ \tilde{\nabla}^2 \tilde{E}_{\vec{r}} + \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 \]</div>
<p>We can convert this to three scalar-wave equations of the form</p>
<div class="arithmatex">\[ (\tilde{\nabla}^2 + K^2) \phi_{\vec{r}} = 0 \]</div>
<p>with <span class="arithmatex">\(K^2 = \Omega^2 \mu \epsilon\)</span>. Now we let</p>
<div class="arithmatex">\[ \phi_{\vec{r}} = A e^{\imath (k_x m \Delta_x + k_y n \Delta_y + k_z p \Delta_z)} \]</div>
<div class="arithmatex">\[ \phi_{\vec{r}} = A e^{\imath (k_x m \Delta_x + k_y n \Delta_y + k_z p \Delta_z)} \]</div>
<p>resulting in</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{\partial}_x &amp;\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 &amp;\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 &amp;= 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}
\]</div>
<p>with similar expressions for the y and z dimnsions (and <span class="arithmatex">\(K_y, K_z\)</span>).</p>
<p>This implies</p>
<div class="arithmatex">\[
\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
\]</div>
<p>where <span class="arithmatex">\(c = \sqrt{\mu \epsilon}\)</span>.</p>
<p>Assuming real <span class="arithmatex">\((k_x, k_y, k_z), \omega\)</span> will be real only if</p>
<div class="arithmatex">\[ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} &lt; 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) \]</div>
<div class="arithmatex">\[ 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}) \]</div>
<p>If <span class="arithmatex">\(\Delta_x = \Delta_y = \Delta_z\)</span>, this simplifies to <span class="arithmatex">\(c \Delta_t &lt; \Delta_x / \sqrt{3}\)</span>.
This last form can be interpreted as enforcing causality; the distance that light
travels in one timestep (i.e., <span class="arithmatex">\(c \Delta_t\)</span>) must be less than the diagonal
@ -2461,15 +2509,16 @@ mu = [mu_xx, mu_yy, mu_zz]
</code></pre></div>
<p>or</p>
<div class="arithmatex">\[
\epsilon = \begin{bmatrix} \epsilon_{xx} &amp; 0 &amp; 0 \\
0 &amp; \epsilon_{yy} &amp; 0 \\
0 &amp; 0 &amp; \epsilon_{zz} \end{bmatrix}
$$
$$
\mu = \begin{bmatrix} \mu_{xx} &amp; 0 &amp; 0 \\
0 &amp; \mu_{yy} &amp; 0 \\
0 &amp; 0 &amp; \mu_{zz} \end{bmatrix}
\epsilon = \begin{bmatrix} \epsilon_{xx} & 0 & 0 \\
0 & \epsilon_{yy} & 0 \\
0 & 0 & \epsilon_{zz} \end{bmatrix}
\]</div>
<div class="arithmatex">\[
\mu = \begin{bmatrix} \mu_{xx} & 0 & 0 \\
0 & \mu_{yy} & 0 \\
0 & 0 & \mu_{zz} \end{bmatrix}
\]</div>
<p>where the off-diagonal terms (e.g. <code>epsilon_xy</code>) are assumed to be zero.</p>
<p>High-accuracy volumetric integration of shapes on multiple grids can be performed
by the <a href="https://mpxd.net/code/jan/gridlock">gridlock</a> module.</p>

503
api/fdtd/index.html
View file

@ -809,6 +809,50 @@
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.temporal_phasor" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;temporal_phasor
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.temporal_phasor_scale" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;temporal_phasor_scale
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.real_injection_scale" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;real_injection_scale
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real
</span>
</a>
</li>
<li class="md-nav__item">
@ -842,6 +886,39 @@
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real_e" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real_e
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real_h" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real_h
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real_j" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real_j
</span>
</a>
</li>
</ul>
@ -1223,6 +1300,50 @@
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.temporal_phasor" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;temporal_phasor
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.temporal_phasor_scale" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;temporal_phasor_scale
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.real_injection_scale" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;real_injection_scale
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real
</span>
</a>
</li>
<li class="md-nav__item">
@ -1256,6 +1377,39 @@
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real_e" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real_e
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real_h" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real_h
</span>
</a>
</li>
<li class="md-nav__item">
<a href="#meanas.fdtd.phasor.reconstruct_real_j" class="md-nav__link">
<span class="md-ellipsis">
<code class="doc-symbol doc-symbol-toc doc-symbol-function"></code>&nbsp;reconstruct_real_j
</span>
</a>
</li>
</ul>
@ -1301,7 +1455,8 @@ mathematical background.</p>
<h3 id="meanas.fdtd--timestep">Timestep<a class="headerlink" href="#meanas.fdtd--timestep" title="Permanent link">&para;</a></h3>
<p>From the discussion of "Plane waves and the Dispersion relation" in <code>meanas.fdmath</code>,
we have</p>
<div class="arithmatex">\[ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} &lt; 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) \]</div>
<div class="arithmatex">\[ 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}) \]</div>
<p>or, if <span class="arithmatex">\(\Delta_x = \Delta_y = \Delta_z\)</span>, then <span class="arithmatex">\(c \Delta_t &lt; \frac{\Delta_x}{\sqrt{3}}\)</span>.</p>
<p>Based on this, we can set</p>
<div class="highlight"><pre><span></span><code>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</p>
<p>The <code>dx_min</code>, <code>dy_min</code>, <code>dz_min</code> should be the minimum value across both the base and derived grids.</p>
<h3 id="meanas.fdtd--poynting-vector-and-energy-conservation">Poynting Vector and Energy Conservation<a class="headerlink" href="#meanas.fdtd--poynting-vector-and-energy-conservation" title="Permanent link">&para;</a></h3>
<p>Let</p>
<div class="arithmatex">\[ \begin{aligned}
\tilde{S}_{l, l', \vec{r}} &amp;=&amp; &amp;\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}} \\
&amp;=&amp; &amp;\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}}) \\
&amp; &amp;+ &amp;\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}}) \\
&amp; &amp;+ &amp;\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}})
<div class="arithmatex">\[
\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}
\]</div>
<p>where <span class="arithmatex">\(\vec{r} = (m, n, p)\)</span> and <span class="arithmatex">\(\otimes\)</span> is a modified cross product
in which the <span class="arithmatex">\(\tilde{E}\)</span> terms are shifted as indicated.</p>
<p>By taking the divergence and rearranging terms, we can show that</p>
<div class="arithmatex">\[
\begin{aligned}
\hat{\nabla} \cdot \tilde{S}_{l, l', \vec{r}}
&amp;= \hat{\nabla} \cdot (\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}}) \\
&amp;= \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}} \\
&amp;= \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}}) \\
&amp;= \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}
\]</div>
<p>where in the last line the spatial subscripts have been dropped to emphasize
the time subscripts <span class="arithmatex">\(l, l'\)</span>, i.e.</p>
<div class="arithmatex">\[
\begin{aligned}
\tilde{E}_l &amp;= \tilde{E}_{l, \vec{r}} \\
\hat{H}_l &amp;= \tilde{H}_{l, \vec{r} + \frac{1}{2}} \\
\tilde{\epsilon} &amp;= \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}
\]</div>
<p>etc.
For <span class="arithmatex">\(l' = l + \frac{1}{2}\)</span> we get</p>
<div class="arithmatex">\[
\begin{aligned}
\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}}
&amp;= \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}} \\
&amp;= (-\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}} \\
&amp;= -(\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 <span class="arithmatex">\(l' = l + \frac{1}{2}\)</span> we get</p>
- \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
\end{aligned}
\]</div>
<p>and for <span class="arithmatex">\(l' = l - \frac{1}{2}\)</span>,</p>
<div class="arithmatex">\[
\begin{aligned}
\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}}
&amp;= (\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 <span class="arithmatex">\(l' = l + \frac{1}{2}\)</span> we get</p>
- \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
\end{aligned}
\]</div>
<p>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:</p>
<div class="arithmatex">\[
\begin{aligned}
U_l &amp;= \epsilon \tilde{E}^2_l + \mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} \\
U_{l + \frac{1}{2}} &amp;= \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}
\]</div>
<p>Rewriting the Poynting theorem in terms of the energy expressions,</p>
<div class="arithmatex">\[
\begin{aligned}
(U_{l+\frac{1}{2}} - U_l) / \Delta_t
&amp;= -\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
&amp;= -\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}
\]</div>
<p>This result is exact and should practically hold to within numerical precision. No time-
or spatial-averaging is necessary.</p>
<p>Note that each value of <span class="arithmatex">\(J\)</span> 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.</p>
<p><code>accumulate_phasor</code> in <code>meanas.fdtd.phasor</code> performs the phase accumulation for one
or more target frequencies, while leaving source normalization and simulation-loop
policy to the caller. Convenience wrappers <code>accumulate_phasor_e</code>,
<code>accumulate_phasor_h</code>, and <code>accumulate_phasor_j</code> apply the usual Yee time offsets.
policy to the caller. <code>temporal_phasor(...)</code> and <code>temporal_phasor_scale(...)</code>
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. <code>real_injection_scale(...)</code> is the corresponding helper for the
common real-valued injection pattern <code>numpy.real(scale * waveform)</code>. Convenience
wrappers <code>accumulate_phasor_e</code>, <code>accumulate_phasor_h</code>, and <code>accumulate_phasor_j</code>
apply the usual Yee time offsets. <code>reconstruct_real(...)</code> and the corresponding
<code>reconstruct_real_e/h/j(...)</code> wrappers perform the inverse operation, converting
phasors back into real-valued field snapshots at explicit Yee-aligned times.
For scalar <code>omega</code>, the reconstruction helpers accept either a plain field
phasor or the batched <code>(1, *sample_shape)</code> form used by the accumulator helpers.
The helpers accumulate</p>
<div class="arithmatex">\[ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l \]</div>
<p>with caller-provided sample times and weights. In this codebase the matching
electric-current convention is typically <code>E -= dt * J / epsilon</code> in FDTD and
<code>-i \omega J</code> on the right-hand side of the FDFD wave equation.</p>
<p>For FDTD/FDFD equivalence, this phasor path is the primary comparison workflow.
It isolates the guided <code>+\omega</code> response that the frequency-domain solver
targets directly, regardless of whether the underlying FDTD run used real- or
complex-valued fields.</p>
<p>For exact pulsed FDTD/FDFD equivalence it is often simplest to keep the
injected source, fields, and CPML auxiliary state complex-valued. The
<code>real_injection_scale(...)</code> helper is instead for the more ordinary one-run
real-valued source path, where the intended positive-frequency waveform is
injected through <code>numpy.real(scale * waveform)</code> and any remaining negative-
frequency leakage is controlled by the pulse bandwidth and run length.</p>
<p><code>reconstruct_real(...)</code> 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,</p>
<div class="arithmatex">\[ E_{\text{monitor}} = E_{\text{guided}} + E_{\text{residual}} . \]</div>
<p>Phasor/modal comparisons mostly validate the guided <code>+\omega</code> 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.</p>
<p>The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse
shape. It can be written</p>
<div class="arithmatex">\[ f_r(t) = (1 - \frac{1}{2} (\omega (t - \tau))^2) e^{-(\frac{\omega (t - \tau)}{2})^2} \]</div>
<p>with <span class="arithmatex">\(\tau &gt; \frac{2 * \pi}{\omega}\)</span> as a minimum delay to avoid a discontinuity at
t=0 (assuming the source is off for t&lt;0 this gives <span class="arithmatex">\(\sim 10^{-3}\)</span> error at t=0).</p>
<h3 id="meanas.fdtd--boundary-conditions">Boundary conditions<a class="headerlink" href="#meanas.fdtd--boundary-conditions" title="Permanent link">&para;</a></h3>
@ -1447,6 +1640,7 @@ coordinate system in both places:</p>
<div class="arithmatex">\[
E \leftarrow E - \Delta_t J / \epsilon
\]</div>
<p>which matches the FDFD right-hand side</p>
<div class="arithmatex">\[
-i \omega J.
@ -2390,19 +2584,20 @@ boundaries of the corresponding <code>U</code> cell. For example,</p>
<a id="__codelineno-0-15" name="__codelineno-0-15" href="#__codelineno-0-15"></a> assert_close(sum(planes), du_half_h2e[mask])
</code></pre></div>
<p>(see <code>meanas.tests.test_fdtd.test_poynting_planes</code>)</p>
<p>The full relationship is
$$
<p>The full relationship is</p>
<div class="arithmatex">\[
\begin{aligned}
(U_{l+\frac{1}{2}} - U_l) / \Delta_t
&amp;= -\hat{\nabla} \cdot \tilde{S}<em l_frac_1="l+\frac{1">{l, l + \frac{1}{2}} \
- \hat{H}</em>}{2}} \cdot \hat{M<em l_frac_1="l+\frac{1">l \
- \tilde{E}_l \cdot \tilde{J}</em> \
&= -\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
&amp;= -\hat{\nabla} \cdot \tilde{S}}{2}<em l-_frac_1="l-\frac{1">{l, l - \frac{1}{2}} \
- \hat{H}</em>}{2}} \cdot \hat{M<em l-_frac_1="l-\frac{1">l \
- \tilde{E}_l \cdot \tilde{J}</em> \
&= -\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}</p>
\]</div>
<p>These equalities are exact and should practically hold to within numerical precision.
No time- or spatial-averaging is necessary. (See <code>meanas.fdtd</code> docs for derivation.)</p>
@ -3579,6 +3774,18 @@ the sampling policy. The accumulated quantity is</p>
<li><code>accumulate_phasor_h(..., step=l)</code> for <code>H_{l + 1/2}</code></li>
<li><code>accumulate_phasor_j(..., step=l)</code> for <code>J_{l + 1/2}</code></li>
</ul>
<p><code>temporal_phasor(...)</code> and <code>temporal_phasor_scale(...)</code> 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.
<code>real_injection_scale(...)</code> is a companion helper for the common real-valued
injection pattern <code>numpy.real(scale * waveform)</code>, where <code>waveform</code> is the
analytic positive-frequency signal and the injected real current should still
produce a desired phasor response.
<code>reconstruct_real(...)</code> and its <code>E/H/J</code> 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 <code>(1, *sample_shape)</code> form used elsewhere in
this module.</p>
<p>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 <code>E -= dt * J / epsilon</code>,
@ -3650,6 +3857,154 @@ factor was derived from a discrete sum that already includes <code>dt</code>, pa
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.temporal_phasor" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">temporal_phasor</span>
<a href="#meanas.fdtd.phasor.temporal_phasor" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">temporal_phasor</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">samples</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="o">*</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a> <span class="n">start_step</span><span class="p">:</span> <span class="nb">int</span> <span class="o">=</span> <span class="mi">0</span><span class="p">,</span>
<a id="__codelineno-0-10" name="__codelineno-0-10" href="#__codelineno-0-10"></a> <span class="n">offset_steps</span><span class="p">:</span> <span class="nb">float</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">,</span>
<a id="__codelineno-0-11" name="__codelineno-0-11" href="#__codelineno-0-11"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">complexfloating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Fourier-project a 1D temporal waveform onto one or more angular frequencies.</p>
<p>The returned quantity is</p>
<div class="highlight"><pre><span></span><code>dt * sum(exp(-1j * omega * t_step) * samples[step_index])
</code></pre></div>
<p>where <code>t_step = (start_step + step_index + offset_steps) * dt</code>.</p>
</div>
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.temporal_phasor_scale" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">temporal_phasor_scale</span>
<a href="#meanas.fdtd.phasor.temporal_phasor_scale" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">temporal_phasor_scale</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">samples</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="o">*</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a> <span class="n">start_step</span><span class="p">:</span> <span class="nb">int</span> <span class="o">=</span> <span class="mi">0</span><span class="p">,</span>
<a id="__codelineno-0-10" name="__codelineno-0-10" href="#__codelineno-0-10"></a> <span class="n">offset_steps</span><span class="p">:</span> <span class="nb">float</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">,</span>
<a id="__codelineno-0-11" name="__codelineno-0-11" href="#__codelineno-0-11"></a> <span class="n">target</span><span class="p">:</span> <span class="n">ArrayLike</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">,</span>
<a id="__codelineno-0-12" name="__codelineno-0-12" href="#__codelineno-0-12"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">complexfloating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Return the scalar multiplier that gives a desired temporal phasor response.</p>
<p>The returned scale satisfies</p>
<div class="highlight"><pre><span></span><code>temporal_phasor(scale * samples, omegas, dt, ...) == target
</code></pre></div>
<p>for each target frequency. The result keeps a leading frequency axis even
when <code>omegas</code> is scalar.</p>
</div>
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.real_injection_scale" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">real_injection_scale</span>
<a href="#meanas.fdtd.phasor.real_injection_scale" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">real_injection_scale</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">samples</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="o">*</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a> <span class="n">start_step</span><span class="p">:</span> <span class="nb">int</span> <span class="o">=</span> <span class="mi">0</span><span class="p">,</span>
<a id="__codelineno-0-10" name="__codelineno-0-10" href="#__codelineno-0-10"></a> <span class="n">offset_steps</span><span class="p">:</span> <span class="nb">float</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">,</span>
<a id="__codelineno-0-11" name="__codelineno-0-11" href="#__codelineno-0-11"></a> <span class="n">target</span><span class="p">:</span> <span class="n">ArrayLike</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">,</span>
<a id="__codelineno-0-12" name="__codelineno-0-12" href="#__codelineno-0-12"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">complexfloating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Return the scale for a real-valued injection built from an analytic waveform.</p>
<p>If the time-domain source is applied as</p>
<div class="highlight"><pre><span></span><code>numpy.real(scale * samples[step])
</code></pre></div>
<p>then the desired positive-frequency phasor is obtained by compensating for
the 1/2 factor between the real-valued source and its analytic component:</p>
<div class="highlight"><pre><span></span><code>scale = 2 * target / temporal_phasor(samples, ...)
</code></pre></div>
<p>This helper normalizes only the intended positive-frequency component. Any
residual negative-frequency leakage is controlled by the waveform design and
the accumulation window.</p>
</div>
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.reconstruct_real" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">reconstruct_real</span>
<a href="#meanas.fdtd.phasor.reconstruct_real" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">reconstruct_real</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">phasors</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="n">step</span><span class="p">:</span> <span class="nb">int</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a> <span class="o">*</span><span class="p">,</span>
<a id="__codelineno-0-10" name="__codelineno-0-10" href="#__codelineno-0-10"></a> <span class="n">offset_steps</span><span class="p">:</span> <span class="nb">float</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">,</span>
<a id="__codelineno-0-11" name="__codelineno-0-11" href="#__codelineno-0-11"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">floating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Reconstruct a real-valued field snapshot from one or more phasors.</p>
<p>The returned quantity is</p>
<div class="highlight"><pre><span></span><code>real(phasor * exp(1j * omega * t_step))
</code></pre></div>
<p>where <code>t_step = (step + offset_steps) * dt</code>.</p>
<p>For multi-frequency inputs, the leading frequency axis is preserved. For a
scalar <code>omega</code>, callers may pass either <code>(1, *sample_shape)</code> or
<code>sample_shape</code>; the return shape matches that choice.</p>
</div>
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.accumulate_phasor_e" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">accumulate_phasor_e</span>
@ -3740,6 +4095,90 @@ factor was derived from a discrete sum that already includes <code>dt</code>, pa
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.reconstruct_real_e" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">reconstruct_real_e</span>
<a href="#meanas.fdtd.phasor.reconstruct_real_e" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">reconstruct_real_e</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">phasors</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="n">step</span><span class="p">:</span> <span class="nb">int</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">floating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Reconstruct a real E-field snapshot taken at integer timestep <code>step</code>.</p>
</div>
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.reconstruct_real_h" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">reconstruct_real_h</span>
<a href="#meanas.fdtd.phasor.reconstruct_real_h" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">reconstruct_real_h</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">phasors</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="n">step</span><span class="p">:</span> <span class="nb">int</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">floating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Reconstruct a real H-field snapshot corresponding to <code>H_{step + 1/2}</code>.</p>
</div>
</div>
<div class="doc doc-object doc-function">
<h3 id="meanas.fdtd.phasor.reconstruct_real_j" class="doc doc-heading">
<code class="doc-symbol doc-symbol-heading doc-symbol-function"></code> <span class="doc doc-object-name doc-function-name">reconstruct_real_j</span>
<a href="#meanas.fdtd.phasor.reconstruct_real_j" class="headerlink" title="Permanent link">&para;</a></h3>
<div class="doc-signature highlight"><pre><span></span><code><a id="__codelineno-0-1" name="__codelineno-0-1" href="#__codelineno-0-1"></a><span class="nf">reconstruct_real_j</span><span class="p">(</span>
<a id="__codelineno-0-2" name="__codelineno-0-2" href="#__codelineno-0-2"></a> <span class="n">phasors</span><span class="p">:</span> <span class="n">ArrayLike</span><span class="p">,</span>
<a id="__codelineno-0-3" name="__codelineno-0-3" href="#__codelineno-0-3"></a> <span class="n">omegas</span><span class="p">:</span> <span class="nb">float</span>
<a id="__codelineno-0-4" name="__codelineno-0-4" href="#__codelineno-0-4"></a> <span class="o">|</span> <span class="nb">complex</span>
<a id="__codelineno-0-5" name="__codelineno-0-5" href="#__codelineno-0-5"></a> <span class="o">|</span> <span class="n">Sequence</span><span class="p">[</span><span class="nb">float</span> <span class="o">|</span> <span class="nb">complex</span><span class="p">]</span>
<a id="__codelineno-0-6" name="__codelineno-0-6" href="#__codelineno-0-6"></a> <span class="o">|</span> <span class="n">NDArray</span><span class="p">,</span>
<a id="__codelineno-0-7" name="__codelineno-0-7" href="#__codelineno-0-7"></a> <span class="n">dt</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<a id="__codelineno-0-8" name="__codelineno-0-8" href="#__codelineno-0-8"></a> <span class="n">step</span><span class="p">:</span> <span class="nb">int</span><span class="p">,</span>
<a id="__codelineno-0-9" name="__codelineno-0-9" href="#__codelineno-0-9"></a><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">NDArray</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">floating</span><span class="p">]</span>
</code></pre></div>
<div class="doc doc-contents ">
<p>Reconstruct a real current snapshot corresponding to <code>J_{step + 1/2}</code>.</p>
</div>
</div>
</div>

192
api/waveguides/index.html
View file

@ -1402,132 +1402,144 @@ a structure with some (x, y) cross-section extending uniformly into the z dimens
with a diagonal <span class="arithmatex">\(\epsilon\)</span> tensor, we have</p>
<div class="arithmatex">\[
\begin{aligned}
\nabla \times \vec{E}(x, y, z) &amp;= -\imath \omega \mu \vec{H} \\
\nabla \times \vec{H}(x, y, z) &amp;= \imath \omega \epsilon \vec{E} \\
\vec{E}(x,y,z) &amp;= (\vec{E}_t(x, y) + E_z(x, y)\vec{z}) e^{-\imath \beta z} \\
\vec{H}(x,y,z) &amp;= (\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}
\]</div>
<p>Expanding the first two equations into vector components, we get</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{xx} H_x &amp;= \partial_y E_z - \partial_z E_y \\
-\imath \omega \mu_{yy} H_y &amp;= \partial_z E_x - \partial_x E_z \\
-\imath \omega \mu_{zz} H_z &amp;= \partial_x E_y - \partial_y E_x \\
\imath \omega \epsilon_{xx} E_x &amp;= \partial_y H_z - \partial_z H_y \\
\imath \omega \epsilon_{yy} E_y &amp;= \partial_z H_x - \partial_x H_z \\
\imath \omega \epsilon_{zz} E_z &amp;= \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}
\]</div>
<p>Substituting in our expressions for <span class="arithmatex">\(\vec{E}\)</span>, <span class="arithmatex">\(\vec{H}\)</span> and discretizing:</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{xx} H_x &amp;= \tilde{\partial}_y E_z + \imath \beta E_y \\
-\imath \omega \mu_{yy} H_y &amp;= -\imath \beta E_x - \tilde{\partial}_x E_z \\
-\imath \omega \mu_{zz} H_z &amp;= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\
\imath \omega \epsilon_{xx} E_x &amp;= \hat{\partial}_y H_z + \imath \beta H_y \\
\imath \omega \epsilon_{yy} E_y &amp;= -\imath \beta H_x - \hat{\partial}_x H_z \\
\imath \omega \epsilon_{zz} E_z &amp;= \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}
\]</div>
<p>Rewrite the last three equations as</p>
<div class="arithmatex">\[
\begin{aligned}
\imath \beta H_y &amp;= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\
\imath \beta H_x &amp;= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\
\imath \omega E_z &amp;= \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}
\]</div>
<p>Now apply <span class="arithmatex">\(\imath \beta \tilde{\partial}_x\)</span> to the last equation,
then substitute in for <span class="arithmatex">\(\imath \beta H_x\)</span> and <span class="arithmatex">\(\imath \beta H_y\)</span>:</p>
<div class="arithmatex">\[
\begin{aligned}
\imath \beta \tilde{\partial}_x \imath \omega E_z &amp;= \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 \\
&amp;= \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) \\
&amp;= \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 &amp;= \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}
\]</div>
<p>With a similar approach (but using <span class="arithmatex">\(\imath \beta \tilde{\partial}_y\)</span> instead), we can get</p>
<div class="arithmatex">\[
\begin{aligned}
\imath \beta \tilde{\partial}_y E_z &amp;= \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}
\]</div>
<p>We can combine this equation for <span class="arithmatex">\(\imath \beta \tilde{\partial}_y E_z\)</span> with
the unused <span class="arithmatex">\(\imath \omega \mu_{xx} H_x\)</span> and <span class="arithmatex">\(\imath \omega \mu_{yy} H_y\)</span> equations to get</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{xx} \imath \beta H_x &amp;= -\beta^2 E_y + \imath \beta \tilde{\partial}_y E_z \\
-\imath \omega \mu_{xx} \imath \beta H_x &amp;= -\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}
\]</div>
<p>and</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{yy} \imath \beta H_y &amp;= \beta^2 E_x - \imath \beta \tilde{\partial}_x E_z \\
-\imath \omega \mu_{yy} \imath \beta H_y &amp;= \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}
\]</div>
<p>However, based on our rewritten equation for <span class="arithmatex">\(\imath \beta H_x\)</span> and the so-far unused
equation for <span class="arithmatex">\(\imath \omega \mu_{zz} H_z\)</span> we can also write</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{xx} (\imath \beta H_x) &amp;= -\imath \omega \mu_{xx} (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\
&amp;= -\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)) \\
&amp;= -\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}
\]</div>
<p>and, similarly,</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{yy} (\imath \beta H_y) &amp;= \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}
\]</div>
<p>By combining both pairs of expressions, we get</p>
<div class="arithmatex">\[
\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)
) &amp;= \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)
) &amp;= -\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}
\]</div>
<p>Using these, we can construct the eigenvalue problem</p>
<div class="arithmatex">\[
\beta^2 \begin{bmatrix} E_x \\
E_y \end{bmatrix} =
(\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} &amp; 0 \\
0 &amp; \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 &amp; \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} &amp; \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}
\]</div>
<p>In the literature, <span class="arithmatex">\(\beta\)</span> is usually used to denote the lossless/real part of the propagation constant,
but in <code>meanas</code> it is allowed to be complex.</p>
<p>An equivalent eigenvalue problem can be formed using the <span class="arithmatex">\(H_x\)</span> and <span class="arithmatex">\(H_y\)</span> fields, if those are more convenient.</p>
@ -1580,15 +1592,16 @@ mu * [[-Dy], [Dx]] / mu * [-Dy, Dx] +
<p>for use with a field vector of the form <code>cat([E_x, E_y])</code>.</p>
<p>More precisely, the operator is</p>
<div class="arithmatex">\[
\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} &amp; 0 \\
0 &amp; \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 &amp; \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} &amp; \hat{\partial}_y \epsilon_{yy} \end{bmatrix}
\begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix}
\]</div>
<p><span class="arithmatex">\(\tilde{\partial}_x\)</span> and <span class="arithmatex">\(\hat{\partial}_x\)</span> are the forward and backward derivatives along x,
and each <span class="arithmatex">\(\epsilon_{xx}\)</span>, <span class="arithmatex">\(\mu_{yy}\)</span>, etc. is a diagonal matrix containing the vectorized material
property distribution.</p>
@ -1735,15 +1748,16 @@ epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +
<p>for use with a field vector of the form <code>cat([H_x, H_y])</code>.</p>
<p>More precisely, the operator is</p>
<div class="arithmatex">\[
\omega^2 \begin{bmatrix} \epsilon_{yy} \mu_{xx} &amp; 0 \\
0 &amp; \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 &amp; \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} &amp; \tilde{\partial}_y \mu_{yy} \end{bmatrix}
\begin{bmatrix} \tilde{\partial}_x \mu_{xx} & \tilde{\partial}_y \mu_{yy} \end{bmatrix}
\]</div>
<p><span class="arithmatex">\(\tilde{\partial}_x\)</span> and <span class="arithmatex">\(\hat{\partial}_x\)</span> are the forward and backward derivatives along x,
and each <span class="arithmatex">\(\epsilon_{xx}\)</span>, <span class="arithmatex">\(\mu_{yy}\)</span>, etc. is a diagonal matrix containing the vectorized material
property distribution.</p>
@ -2066,6 +2080,7 @@ and <code>exy2h(...)</code>, then normalizes them to unit forward power using
<div class="arithmatex">\[
\Re\left[\mathrm{inner\_product}(e, h, \mathrm{conj\_h}=True)\right] = 1,
\]</div>
<p>so the returned fields represent a unit-power forward mode under the
discrete Yee-grid Poynting inner product.</p>
</details>
@ -2721,21 +2736,23 @@ identical representatives for nearly symmetric modes.</p>
<div class="arithmatex">\[
\imath \omega \epsilon_{zz} E_z = \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
\]</div>
<p>as well as the intermediate equations</p>
<div class="arithmatex">\[
\begin{aligned}
\imath \beta H_y &amp;= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\
\imath \beta H_x &amp;= -\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}
\]</div>
<p>Combining these, we get</p>
<div class="arithmatex">\[
\begin{aligned}
E_z &amp;= \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))
&amp;= \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}
\]</div>
@ -3654,10 +3671,12 @@ E_z &amp;= \frac{1}{- \omega \beta \epsilon_{zz}} ((
<span class="arithmatex">\(\beta\)</span> to changes in the dielectric structure <span class="arithmatex">\(\epsilon\)</span>.</p>
<p>The output is a vector of the same size as <code>vec(epsilon)</code>, with each element specifying the
sensitivity of <code>wavenumber</code> to changes in the corresponding element in <code>vec(epsilon)</code>, i.e.</p>
<div class="arithmatex">\[sens_{i} = \frac{\partial\beta}{\partial\epsilon_i}\]</div>
<div class="arithmatex">\[ sens_{i} = \frac{\partial\beta}{\partial\epsilon_i} \]</div>
<p>An adjoint approach is used to calculate the sensitivity; the derivation is provided here:</p>
<p>Starting with the eigenvalue equation</p>
<div class="arithmatex">\[\beta^2 E_{xy} = A_E E_{xy}\]</div>
<div class="arithmatex">\[ \beta^2 E_{xy} = A_E E_{xy} \]</div>
<p>where <span class="arithmatex">\(A_E\)</span> is the waveguide operator from <code>operator_e()</code>, and <span class="arithmatex">\(E_{xy} = \begin{bmatrix} E_x \\
E_y \end{bmatrix}\)</span>,
we can differentiate with respect to one of the <span class="arithmatex">\(\epsilon\)</span> elements (i.e. at one Yee grid point), <span class="arithmatex">\(\epsilon_i\)</span>:</p>
@ -3665,11 +3684,13 @@ we can differentiate with respect to one of the <span class="arithmatex">\(\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}
\]</div>
<p>We then multiply by <span class="arithmatex">\(H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}\)</span> from the left:</p>
<div class="arithmatex">\[
(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}
\]</div>
<p>However, <span class="arithmatex">\(H_{yx}^\star\)</span> is actually a left-eigenvector of <span class="arithmatex">\(A_E\)</span>. This can be verified by inspecting
the form of <code>operator_h</code> (<span class="arithmatex">\(A_H\)</span>) and comparing its conjugate transpose to <code>operator_e</code> (<span class="arithmatex">\(A_E\)</span>). Also, note
<span class="arithmatex">\(H_{yx}^\star \cdot E_{xy} = H^\star \times E\)</span> recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201
@ -3677,11 +3698,13 @@ for a similar approach. Therefore,</p>
<div class="arithmatex">\[
H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
\]</div>
<p>and we can simplify to</p>
<div class="arithmatex">\[
\partial_{\epsilon_i}(\beta)
= \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}}
\]</div>
<p>This expression can be quickly calculated for all <span class="arithmatex">\(i\)</span> by writing out the various terms of
<span class="arithmatex">\(\partial_{\epsilon_i} A_E\)</span> and recognizing that the vector-matrix-vector products (i.e. scalars)
<span class="arithmatex">\(sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}\)</span>, indexed by <span class="arithmatex">\(i\)</span>, can be expressed as
@ -4161,6 +4184,7 @@ longitudinal Poynting flux,</p>
<div class="arithmatex">\[
\frac{1}{2}\int (E_x H_y - E_y H_x) \, dx \, dy
\]</div>
<p>with the Yee-grid staggering and optional propagation-phase adjustment used
by the waveguide helpers in this module.</p>
@ -4866,6 +4890,7 @@ same sign convention used elsewhere in the package:</p>
<div class="arithmatex">\[
\sum \mathrm{overlap\_e} \; E_\mathrm{other}^*
\]</div>
<p>where the sum is over the full Yee-grid field storage.</p>
<p>The construction uses a two-cell window immediately upstream of the selected
slice:</p>
@ -4881,6 +4906,7 @@ overlap relation involving</p>
<div class="arithmatex">\[
\int ((E \times H_\mathrm{mode}) + (E_\mathrm{mode} \times H)) \cdot dn.
\]</div>
<p>E x H_mode + E_mode x H
-&gt; 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</p>
<div class="arithmatex">\[
e^{-i \, \mathrm{polarity} \, wavenumber \, \Delta z}
\]</div>
<p>to each copied slice.</p>
</details>
@ -5254,11 +5281,13 @@ and the fields are assumed to have the form</p>
<div class="arithmatex">\[
\vec{E}(r, y, \theta), \vec{H}(r, y, \theta) \propto e^{-\imath m \theta},
\]</div>
<p>where <code>m</code> is the angular wavenumber returned by <code>solve_mode(s)</code>. It is often
convenient to introduce the corresponding linear wavenumber</p>
<div class="arithmatex">\[
\beta = \frac{m}{r_{\min}},
\]</div>
<p>so that the cylindrical problem resembles the straight-waveguide problem with
additional metric factors.</p>
<p>Those metric factors live on the staggered radial Yee grids. If the left edge of
@ -5266,33 +5295,36 @@ the computational window is at <code>r = r_{\min}</code>, define the electric-gr
magnetic-grid radial sample locations by</p>
<div class="arithmatex">\[
\begin{aligned}
r_a(n) &amp;= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\
r_b\!\left(n + \tfrac{1}{2}\right) &amp;= r_{\min} + \tfrac{1}{2}\Delta r_{e, n}
+ \sum_{j &lt; 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}
\]</div>
<p>and from them the diagonal metric matrices</p>
<div class="arithmatex">\[
\begin{aligned}
T_a &amp;= \operatorname{diag}(r_a / r_{\min}), \\
T_b &amp;= \operatorname{diag}(r_b / r_{\min}).
T_a &= \operatorname{diag}(r_a / r_{\min}), \\
T_b &= \operatorname{diag}(r_b / r_{\min}).
\end{aligned}
\]</div>
<p>With the same forward/backward derivative notation used in <code>waveguide_2d</code>, the
coordinate-transformed discrete curl equations used here are</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{rr} H_r &amp;= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
-\imath \omega \mu_{yy} H_y &amp;= -\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 &amp;= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\
\imath \beta H_y &amp;= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\
\imath \beta H_r &amp;= \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 &amp;= 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}
\]</div>
<p>The first three equations are the cylindrical analogue of the straight-guide
relations for <code>H_r</code>, <code>H_y</code>, and <code>H_z</code>. The next two are the metric-weighted
versions of the straight-guide identities for <code>\imath \beta H_y</code> and
@ -5306,6 +5338,7 @@ and then eliminate <code>H_z</code> with</p>
H_z = \frac{1}{-\imath \omega \mu_{zz}}
\left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right).
\]</div>
<p>This yields the transverse electric eigenproblem implemented by
<code>cylindrical_operator(...)</code>:</p>
<div class="arithmatex">\[
@ -5315,8 +5348,8 @@ H_z = \frac{1}{-\imath \omega \mu_{zz}}
\left(
\omega^2
\begin{bmatrix}
T_b^2 \mu_{yy} \epsilon_{xx} &amp; 0 \\
0 &amp; 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} &amp; 0 \\
\end{bmatrix}
T_b \mu_{zz}^{-1}
\begin{bmatrix}
-\tilde{\partial}_y &amp; \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} &amp;
\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}.
\]</div>
<p>Since <code>\beta = m / r_{\min}</code>, the solver implemented in this file returns the
angular wavenumber <code>m</code>, while the operator itself is most naturally written in
terms of the linear quantity <code>\beta</code>. The helpers below reconstruct the full
@ -5389,17 +5423,18 @@ product used by <code>waveguide_2d</code>.</p>
<p>Cylindrical coordinate waveguide operator of the form</p>
<div class="arithmatex">\[
(\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} &amp; 0 \\
0 &amp; 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 &amp; \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} &amp; \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}
\]</div>
<p>for use with a field vector of the form <code>[E_r, E_y]</code>.</p>
<p>This operator can be used to form an eigenvalue problem of the form
A @ [E_r, E_y] = beta**2 * [E_r, E_y]</p>
@ -6288,15 +6323,15 @@ from E_r and E_y in order to then calculate E_z.</p>
<p>Returns an operator which, when applied to a vectorized E eigenfield, produces
the vectorized H eigenfield.</p>
<p>This operator is created directly from the initial coordinate-transformed equations:
$$
<p>This operator is created directly from the initial coordinate-transformed equations:</p>
<div class="arithmatex">\[
\begin{aligned}
-\imath \omega \mu_{rr} H_r &amp;= \tilde{\partial}<em yy="yy">y E_z + \imath \beta T_a^{-1} E_y, \
-\imath \omega \mu</em> E_r
- T_b^{-1} \tilde{\partial}} H_y &amp;= -\imath \beta T_b^{-1<em zz="zz">r (T_a E_z), \
-\imath \omega \mu</em>_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 &amp;= \tilde{\partial}_r E_y - \tilde{\partial</p>
\]</div>
<p><span class="doc-section-title">Parameters:</span></p>
@ -6745,6 +6780,7 @@ enforces unit forward power under the discrete inner product</p>
<div class="arithmatex">\[
\frac{1}{2}\int (E_r H_y^* - E_y H_r^*) \, dr \, dy.
\]</div>
<p>The angular wavenumber <code>m</code> 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.</p>