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

Use raw strings to avoid double backslashes

Browse Source
This commit is contained in:
Jan Petykiewicz 2024-07-15 16:32:31 -07:00
parent 2d48858973
commit 77715da8b4
8 changed files with 291 additions and 291 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
meanas/fdfd/__init__.py
View File

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

6
meanas/fdfd/functional.py
View File

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

42
meanas/fdfd/operators.py
View File

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

2
meanas/fdfd/waveguide_2d.py
View File

@ -211,7 +211,7 @@ def operator_e(
$$ $$
\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ \omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\
0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} +
\begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\ \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\
\mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1} \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1}
\begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +

266
meanas/fdmath/__init__.py
View File

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

8
meanas/fdmath/functional.py
View File

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

148
meanas/fdtd/__init__.py
View File

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

22
meanas/fdtd/energy.py
View File

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