fdtd¶

meanas.fdtd ¶

Utilities for running finite-difference time-domain (FDTD) simulations

See the discussion of Maxwell's Equations in meanas.fdmath for basic mathematical background.

Timestep¶

From the discussion of "Plane waves and the Dispersion relation" in meanas.fdmath, we have

\[ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) \]

or, if \(\Delta_x = \Delta_y = \Delta_z\), then \(c \Delta_t < \frac{\Delta_x}{\sqrt{3}}\).

Based on this, we can set

dt = sqrt(mu.min() * epsilon.min()) / sqrt(1/dx_min**2 + 1/dy_min**2 + 1/dz_min**2)

The dx_min, dy_min, dz_min should be the minimum value across both the base and derived grids.

Poynting Vector and Energy Conservation¶

Let

\[ \begin{aligned} \tilde{S}_{l, l', \vec{r}} &=& &\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}} \\ &=& &\vec{x} (\tilde{E}^y_{l,m+1,n,p} \hat{H}^z_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^z_{l,m+1,n,p} \hat{H}^y_{l', \vec{r} + \frac{1}{2}}) \\ & &+ &\vec{y} (\tilde{E}^z_{l,m,n+1,p} \hat{H}^x_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^x_{l,m,n+1,p} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \\ & &+ &\vec{z} (\tilde{E}^x_{l,m,n,p+1} \hat{H}^y_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^y_{l,m,n,p+1} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \end{aligned} \]

where \(\vec{r} = (m, n, p)\) and \(\otimes\) is a modified cross product in which the \(\tilde{E}\) terms are shifted as indicated.

By taking the divergence and rearranging terms, we can show that

\[ \begin{aligned} \hat{\nabla} \cdot \tilde{S}_{l, l', \vec{r}} &= \hat{\nabla} \cdot (\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}}) \\ &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l, \vec{r}} - \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 (-\tilde{\partial}_t \mu_{\vec{r} + \frac{1}{2}} \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} - \hat{M}_{l, \vec{r} + \frac{1}{2}}) - \tilde{E}_{l, \vec{r}} \cdot (\hat{\partial}_t \tilde{\epsilon}_{\vec{r}} \tilde{E}_{l'+\frac{1}{2}, \vec{r}} + \tilde{J}_{l', \vec{r}}) \\ &= \hat{H}_{l'} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) - \tilde{E}_l \cdot (\epsilon / \Delta_t )(\tilde{E}_{l'+\frac{1}{2}} - \tilde{E}_{l'-\frac{1}{2}}) - \hat{H}_{l'} \cdot \hat{M}_{l} - \tilde{E}_l \cdot \tilde{J}_{l'} \\ \end{aligned} \]

where in the last line the spatial subscripts have been dropped to emphasize the time subscripts \(l, l'\), i.e.

\[ \begin{aligned} \tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\ \hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}} \\ \tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}} \\ \end{aligned} \]

etc. For \(l' = l + \frac{1}{2}\) we get

\[ \begin{aligned} \hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} &= \hat{H}_{l + \frac{1}{2}} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) - \tilde{E}_l \cdot (\epsilon / \Delta_t)(\tilde{E}_{l+1} - \tilde{E}_l) - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\ &= (-\mu / \Delta_t)(\hat{H}^2_{l + \frac{1}{2}} - \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}) - (\epsilon / \Delta_t)(\tilde{E}_{l+1} \cdot \tilde{E}_l - \tilde{E}^2_l) - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\ &= -(\mu \hat{H}^2_{l + \frac{1}{2}} +\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 \\ - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\ - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\ \end{aligned} \]

and for \(l' = l - \frac{1}{2}\),

\[ \begin{aligned} \hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} &= (\mu \hat{H}^2_{l - \frac{1}{2}} +\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 \\ - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\ - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\ \end{aligned} \]

These two results form the discrete time-domain analogue to Poynting's theorem. They hint at the expressions for the energy, which can be calculated at the same time-index as either the E or H field:

\[ \begin{aligned} U_l &= \epsilon \tilde{E}^2_l + \mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} \\ U_{l + \frac{1}{2}} &= \epsilon \tilde{E}_l \cdot \tilde{E}_{l + 1} + \mu \hat{H}^2_{l + \frac{1}{2}} \\ \end{aligned} \]

Rewriting the Poynting theorem in terms of the energy expressions,

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

This result is exact and should practically hold to within numerical precision. No time- or spatial-averaging is necessary.

Note that each value of \(J\) contributes to the energy twice (i.e. once per field update) despite only causing the value of \(E\) to change once (same for \(M\) and \(H\)).

Sources¶

It is often useful to excite the simulation with an arbitrary broadband pulse and then extract the frequency-domain response by performing an on-the-fly Fourier transform of the time-domain fields.

accumulate_phasor in meanas.fdtd.phasor performs the phase accumulation for one or more target frequencies, while leaving source normalization and simulation-loop policy to the caller. temporal_phasor(...) and temporal_phasor_scale(...) apply the same Fourier sum to a scalar waveform, which is useful when a pulsed source envelope must be normalized before being applied to a point source or mode source. real_injection_scale(...) is the corresponding helper for the common real-valued injection pattern numpy.real(scale * waveform). Convenience wrappers accumulate_phasor_e, accumulate_phasor_h, and accumulate_phasor_j apply the usual Yee time offsets. reconstruct_real(...) and the corresponding reconstruct_real_e/h/j(...) wrappers perform the inverse operation, converting phasors back into real-valued field snapshots at explicit Yee-aligned times. For scalar omega, the reconstruction helpers accept either a plain field phasor or the batched (1, *sample_shape) form used by the accumulator helpers. The helpers accumulate

\[ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l \]

with caller-provided sample times and weights. In this codebase the matching electric-current convention is typically E -= dt * J / epsilon in FDTD and -i \omega J on the right-hand side of the FDFD wave equation.

For FDTD/FDFD equivalence, this phasor path is the primary comparison workflow. It isolates the guided +\omega response that the frequency-domain solver targets directly, regardless of whether the underlying FDTD run used real- or complex-valued fields.

For exact pulsed FDTD/FDFD equivalence it is often simplest to keep the injected source, fields, and CPML auxiliary state complex-valued. The real_injection_scale(...) helper is instead for the more ordinary one-run real-valued source path, where the intended positive-frequency waveform is injected through numpy.real(scale * waveform) and any remaining negative- frequency leakage is controlled by the pulse bandwidth and run length.

reconstruct_real(...) is for a different question: given a phasor, what late real-valued field snapshot should it produce? That raw-snapshot comparison is stricter and noisier because a monitor plane generally contains both the guided field and the remaining orthogonal content,

\[ E_{\text{monitor}} = E_{\text{guided}} + E_{\text{residual}} . \]

Phasor/modal comparisons mostly validate the guided +\omega term. Raw real-field comparisons expose both terms at once, so they should be treated as secondary diagnostics rather than the main solver-equivalence benchmark.

The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse shape. It can be written

\[ f_r(t) = (1 - \frac{1}{2} (\omega (t - \tau))^2) e^{-(\frac{\omega (t - \tau)}{2})^2} \]

with \(\tau > \frac{2 * \pi}{\omega}\) as a minimum delay to avoid a discontinuity at t=0 (assuming the source is off for t<0 this gives \(\sim 10^{-3}\) error at t=0).

Boundary conditions¶

meanas.fdtd exposes two boundary-related building blocks:

conducting_boundary(...) for simple perfect-electric-conductor style field clamping at one face of the domain.
cpml_params(...) and updates_with_cpml(...) for convolutional perfectly matched layers (CPMLs) attached to one or more faces of the Yee grid.

updates_with_cpml(...) accepts a three-by-two table of CPML parameter blocks:

cpml_params[axis][polarity_index]

where axis is 0, 1, or 2 and polarity_index corresponds to (-1, +1). Passing None for one entry disables CPML on that face while leaving the other directions unchanged. This is how mixed boundary setups such as "absorbing in x, periodic in y/z" are expressed.

When comparing an FDTD run against an FDFD solve, use the same stretched coordinate system in both places:

Build the FDTD update with the desired CPML parameters.
Stretch the FDFD dxes with the matching SCPML transform.
Compare the extracted phasor against the FDFD field or residual on those stretched dxes.

The electric-current sign convention used throughout the examples and tests is

\[ E \leftarrow E - \Delta_t J / \epsilon \]

which matches the FDFD right-hand side

\[ -i \omega J. \]

Core update and analysis helpers¶

meanas.fdtd.base ¶

Basic FDTD field updates

maxwell_e ¶

maxwell_e(
    dt: float, dxes: dx_lists_t | None = None
) -> fdfield_updater_t

Build a function which performs a portion the time-domain E-field update,

E += curl_back(H[t]) / epsilon

The full update should be

E += (curl_back(H[t]) + J) / epsilon

which requires an additional step of E += J / epsilon which is not performed by the generated function.

See meanas.fdmath for descriptions of

This update step: "Maxwell's equations" section
dxes: "Datastructure: dx_lists_t" section
epsilon: "Permittivity and Permeability" section

Also see the "Timestep" section of meanas.fdtd for a discussion of the dt parameter.

Parameters:

Name	Type	Description	Default
`dt`	`float`	Timestep. See `meanas.fdtd` for details.	required
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Type	Description
`fdfield_updater_t`	Function `f(E_old, H_old, epsilon) -> E_new`.

maxwell_h ¶

maxwell_h(
    dt: float, dxes: dx_lists_t | None = None
) -> fdfield_updater_t

Build a function which performs part of the time-domain H-field update,

H -= curl_forward(E[t]) / mu

The full update should be

H -= (curl_forward(E[t]) + M) / mu

which requires an additional step of H -= M / mu which is not performed by the generated function; this step can be omitted if there is no magnetic current M.

See meanas.fdmath for descriptions of

This update step: "Maxwell's equations" section
dxes: "Datastructure: dx_lists_t" section
mu: "Permittivity and Permeability" section

Also see the "Timestep" section of meanas.fdtd for a discussion of the dt parameter.

Parameters:

Name	Type	Description	Default
`dt`	`float`	Timestep. See `meanas.fdtd` for details.	required
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Type	Description
`fdfield_updater_t`	Function `f(E_old, H_old, epsilon) -> E_new`.

meanas.fdtd.pml ¶

Convolutional perfectly matched layer (CPML) support for FDTD updates.

The helpers in this module construct per-face CPML parameters and then wrap the standard Yee updates with the additional auxiliary psi fields needed by the CPML recurrence.

The intended call pattern is:

Build a cpml_params[axis][polarity_index] table with cpml_params(...).
Pass that table into updates_with_cpml(...) together with dt, dxes, and epsilon.
Advance the returned update_E / update_H closures in the simulation loop.

Each face can be enabled or disabled independently by replacing one table entry with None.

cpml_params ¶

cpml_params(
    axis: int,
    polarity: int,
    dt: float,
    thickness: int = 8,
    ln_R_per_layer: float = -1.6,
    epsilon_eff: float = 1,
    mu_eff: float = 1,
    m: float = 3.5,
    ma: float = 1,
    cfs_alpha: float = 0,
) -> dict[str, Any]

Construct the parameter block for one CPML face.

Parameters:

Name	Type	Description	Default
`axis`	`int`	Which Cartesian axis the CPML is normal to (`0`, `1`, or `2`).	required
`polarity`	`int`	Which face along that axis (`-1` for the low-index face, `+1` for the high-index face).	required
`dt`	`float`	Timestep used by the Yee update.	required
`thickness`	`int`	Number of Yee cells occupied by the CPML region.	`8`
`ln_R_per_layer`	`float`	Logarithmic attenuation target per layer.	`-1.6`
`epsilon_eff`	`float`	Effective permittivity used when choosing the CPML scaling.	`1`
`mu_eff`	`float`	Effective permeability used when choosing the CPML scaling.	`1`
`m`	`float`	Polynomial grading exponent for `sigma` and `kappa`.	`3.5`
`ma`	`float`	Polynomial grading exponent for the complex-frequency shift `alpha`.	`1`
`cfs_alpha`	`float`	Maximum complex-frequency shift parameter.	`0`

Returns:

Type	Description
`dict[str, Any]`	Dictionary with:
`dict[str, Any]`	`param_e`: `(p0, p1, p2)` arrays for the E update
`dict[str, Any]`	`param_h`: `(p0, p1, p2)` arrays for the H update
`dict[str, Any]`	`region`: slice tuple selecting the CPML cells on that face

updates_with_cpml ¶

updates_with_cpml(
    cpml_params: Sequence[Sequence[dict[str, Any] | None]],
    dt: float,
    dxes: dx_lists_t,
    epsilon: fdfield,
    *,
    dtype: DTypeLike = numpy.float32,
) -> tuple[
    Callable[[fdfield_t, fdfield_t, fdfield_t], None],
    Callable[[fdfield_t, fdfield_t, fdfield_t], None],
]

Build Yee-step update closures augmented with CPML terms.

Parameters:

Name	Type	Description	Default
`cpml_params`	`Sequence[Sequence[dict[str, Any] \| None]]`	Three-by-two sequence indexed as `[axis][polarity_index]`. Entries are the dictionaries returned by `cpml_params(...)`; use `None` to disable CPML on one face.	required
`dt`	`float`	Timestep.	required
`dxes`	`dx_lists_t`	Yee-grid spacing lists `[dx_e, dx_h]`.	required
`epsilon`	`fdfield`	Electric material distribution used by the E update.	required
`dtype`	`DTypeLike`	Storage dtype for the auxiliary CPML state arrays.	`float32`

Returns:

Type	Description
`Callable[[fdfield_t, fdfield_t, fdfield_t], None]`	`(update_E, update_H)` closures with the same call shape as the basic
`Callable[[fdfield_t, fdfield_t, fdfield_t], None]`	Yee updates:
`tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]`	`update_E(e, h, epsilon)`
`tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]`	`update_H(e, h, mu)`
`tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]`	The closures retain the CPML auxiliary state internally.

meanas.fdtd.boundaries ¶

Boundary conditions

TODO conducting boundary documentation¶

meanas.fdtd.energy ¶

poynting ¶

poynting(
    e: fdfield, h: fdfield, dxes: dx_lists_t | None = None
) -> fdfield_t

Calculate the poynting vector S (\(S\)).

This is the energy transfer rate (amount of energy U per dt transferred between adjacent cells) in each direction that happens during the half-step bounded by the two provided fields.

The returned vector field S is the energy flow across +x, +y, and +z boundaries of the corresponding U cell. For example,

    mx = numpy.roll(mask, -1, axis=0)
    my = numpy.roll(mask, -1, axis=1)
    mz = numpy.roll(mask, -1, axis=2)

    u_hstep = fdtd.energy_hstep(e0=es[ii - 1], h1=hs[ii], e2=es[ii],     **args)
    u_estep = fdtd.energy_estep(h0=hs[ii],     e1=es[ii], h2=hs[ii + 1], **args)
    delta_j_B = fdtd.delta_energy_j(j0=js[ii], e1=es[ii], dxes=dxes)
    du_half_h2e = u_estep - u_hstep - delta_j_B

    s_h2e = -fdtd.poynting(e=es[ii], h=hs[ii], dxes=dxes) * dt
    planes = [s_h2e[0, mask].sum(), -s_h2e[0, mx].sum(),
              s_h2e[1, mask].sum(), -s_h2e[1, my].sum(),
              s_h2e[2, mask].sum(), -s_h2e[2, mz].sum()]

    assert_close(sum(planes), du_half_h2e[mask])

(see meanas.tests.test_fdtd.test_poynting_planes)

The full relationship is

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

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

Parameters:

Name	Type	Description	Default
`e`	`fdfield`	E-field	required
`h`	`fdfield`	H-field (one half-timestep before or after `e`)	required
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Name	Type	Description
`s`	`fdfield_t`	Vector field. Components indicate the energy transfer rate from the corresponding energy cell into its +x, +y, and +z neighbors during the half-step from the time of the earlier input field until the time of later input field.

poynting_divergence ¶

poynting_divergence(
    s: fdfield | None = None,
    *,
    e: fdfield | None = None,
    h: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Calculate the divergence of the poynting vector.

This is the net energy flow for each cell, i.e. the change in energy U per dt caused by transfer of energy to nearby cells (rather than absorption/emission by currents J or M).

See poynting and meanas.fdtd for more details. Args: s: Poynting vector, as calculated with poynting. Optional; caller can provide e and h instead. e: E-field (optional; need either s or both e and h) h: H-field (optional; need either s or both e and h) dxes: Grid description; see meanas.fdmath.

Returns:

Name	Type	Description
`ds`	`fdfield_t`	Divergence of the poynting vector. Entries indicate the net energy flow out of the corresponding energy cell.

energy_hstep ¶

energy_hstep(
    e0: fdfield,
    h1: fdfield,
    e2: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Calculate energy U at the time of the provided H-field h1.

TODO: Figure out what this means spatially.

Parameters:

Name	Type	Description	Default
`e0`	`fdfield`	E-field one half-timestep before the energy.	required
`h1`	`fdfield`	H-field (at the same timestep as the energy).	required
`e2`	`fdfield`	E-field one half-timestep after the energy.	required
`epsilon`	`fdfield \| None`	Dielectric constant distribution.	`None`
`mu`	`fdfield \| None`	Magnetic permeability distribution.	`None`
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Type	Description
`fdfield_t`	Energy, at the time of the H-field `h1`.

energy_estep ¶

energy_estep(
    h0: fdfield,
    e1: fdfield,
    h2: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Calculate energy U at the time of the provided E-field e1.

TODO: Figure out what this means spatially.

Parameters:

Name	Type	Description	Default
`h0`	`fdfield`	H-field one half-timestep before the energy.	required
`e1`	`fdfield`	E-field (at the same timestep as the energy).	required
`h2`	`fdfield`	H-field one half-timestep after the energy.	required
`epsilon`	`fdfield \| None`	Dielectric constant distribution.	`None`
`mu`	`fdfield \| None`	Magnetic permeability distribution.	`None`
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Type	Description
`fdfield_t`	Energy, at the time of the E-field `e1`.

delta_energy_h2e ¶

delta_energy_h2e(
    dt: float,
    e0: fdfield,
    h1: fdfield,
    e2: fdfield,
    h3: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Change in energy during the half-step from h1 to e2.

This is just from (e2 * e2 + h3 * h1) - (h1 * h1 + e0 * e2)

Parameters:

Name	Type	Description	Default
`e0`	`fdfield`	E-field one half-timestep before the start of the energy delta.	required
`h1`	`fdfield`	H-field at the start of the energy delta.	required
`e2`	`fdfield`	E-field at the end of the energy delta (one half-timestep after `h1`).	required
`h3`	`fdfield`	H-field one half-timestep after the end of the energy delta.	required
`epsilon`	`fdfield \| None`	Dielectric constant distribution.	`None`
`mu`	`fdfield \| None`	Magnetic permeability distribution.	`None`
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Type	Description
`fdfield_t`	Change in energy from the time of `h1` to the time of `e2`.

delta_energy_e2h ¶

delta_energy_e2h(
    dt: float,
    h0: fdfield,
    e1: fdfield,
    h2: fdfield,
    e3: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Change in energy during the half-step from e1 to h2.

This is just from (h2 * h2 + e3 * e1) - (e1 * e1 + h0 * h2)

Parameters:

Name	Type	Description	Default
`h0`	`fdfield`	E-field one half-timestep before the start of the energy delta.	required
`e1`	`fdfield`	H-field at the start of the energy delta.	required
`h2`	`fdfield`	E-field at the end of the energy delta (one half-timestep after `e1`).	required
`e3`	`fdfield`	H-field one half-timestep after the end of the energy delta.	required
`epsilon`	`fdfield \| None`	Dielectric constant distribution.	`None`
`mu`	`fdfield \| None`	Magnetic permeability distribution.	`None`
`dxes`	`dx_lists_t \| None`	Grid description; see `meanas.fdmath`.	`None`

Returns:

Type	Description
`fdfield_t`	Change in energy from the time of `e1` to the time of `h2`.

delta_energy_j ¶

delta_energy_j(
    j0: fdfield, e1: fdfield, dxes: dx_lists_t | None = None
) -> fdfield_t

Calculate the electric-current work term \(J \cdot E\) on the Yee grid.

This is the source contribution that appears beside the flux divergence in the discrete Poynting identities documented in meanas.fdtd.

Note that each value of J contributes twice in a full Yee cycle (once per half-step energy balance) even though it directly changes E only once.

Parameters:

Name	Type	Description	Default
`j0`	`fdfield`	Electric-current density sampled at the same half-step as the current work term.	required
`e1`	`fdfield`	Electric field sampled at the matching integer timestep.	required
`dxes`	`dx_lists_t \| None`	Grid description; defaults to unit spacing.	`None`

Returns:

Type	Description
`fdfield_t`	Per-cell source-work contribution with the scalar field shape.

dxmul ¶

dxmul(
    ee: fdfield,
    hh: fdfield,
    epsilon: fdfield | float | None = None,
    mu: fdfield | float | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Multiply E- and H-like field products by material weights and cell volumes.

Parameters:

Name	Type	Description	Default
`ee`	`fdfield`	Three-component electric-field product, such as `e0 * e2`.	required
`hh`	`fdfield`	Three-component magnetic-field product, such as `h1 * h1`.	required
`epsilon`	`fdfield \| float \| None`	Electric material weight; defaults to `1`.	`None`
`mu`	`fdfield \| float \| None`	Magnetic material weight; defaults to `1`.	`None`
`dxes`	`dx_lists_t \| None`	Grid description; defaults to unit spacing.	`None`

Returns:

Type	Description
`fdfield_t`	Scalar field containing the weighted electric plus magnetic contribution
`fdfield_t`	for each Yee cell.

meanas.fdtd.phasor ¶

Helpers for extracting single- or multi-frequency phasors from FDTD samples.

These helpers are intentionally low-level: callers own the accumulator arrays and the sampling policy. The accumulated quantity is

dt * sum(weight * exp(-1j * omega * t_step) * sample_step)

where t_step = (step + offset_steps) * dt.

The usual Yee offsets are:

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

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

These helpers do not choose warmup/accumulation windows or pulse-envelope normalization. They also do not impose a current sign convention. In this codebase, electric-current injection normally appears as E -= dt * J / epsilon, which matches the FDFD right-hand side -1j * omega * J.

accumulate_phasor ¶

accumulate_phasor(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    offset_steps: float = 0.0,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Add one time-domain sample into a phasor accumulator.

The added quantity is

dt * weight * exp(-1j * omega * t_step) * sample

where t_step = (step + offset_steps) * dt.

Note

This helper already multiplies by dt. If the caller's normalization factor was derived from a discrete sum that already includes dt, pass weight / dt here.

temporal_phasor ¶

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

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

The returned quantity is

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

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

temporal_phasor_scale ¶

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

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

The returned scale satisfies

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

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

real_injection_scale ¶

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

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

If the time-domain source is applied as

numpy.real(scale * samples[step])

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

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

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

reconstruct_real ¶

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

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

The returned quantity is

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

where t_step = (step + offset_steps) * dt.

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

accumulate_phasor_e ¶

accumulate_phasor_e(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Accumulate an E-field sample taken at integer timestep step.

accumulate_phasor_h ¶

accumulate_phasor_h(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Accumulate an H-field sample corresponding to H_{step + 1/2}.

accumulate_phasor_j ¶

accumulate_phasor_j(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Accumulate a current sample corresponding to J_{step + 1/2}.

reconstruct_real_e ¶

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

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

reconstruct_real_h ¶

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

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

reconstruct_real_j ¶

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

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