Skip to content

DOC_FD-TD.md

Latest commit

 

History

History
659 lines (515 loc) · 47.8 KB

DOC_FD-TD.md

File metadata and controls

659 lines (515 loc) · 47.8 KB
Raw

FD-TD — GPU-Accelerated 2D Finite-Difference Time-Domain Electromagnetic Simulator

Location: python/NOTEBOKS/FD-TD/ · Language: Python 3 + WGSL (WebGPU Shading Language) · Documented: 2026-07-11

Abstract

This program solves the two-dimensional transverse-magnetic (TM$_z$) Maxwell equations with the Finite-Difference Time-Domain (FDTD) method of Yee [1], entirely on the GPU through the WebGPU compute API (wgpu-py). Absorbing boundaries are implemented with Bérenger's split-field Perfectly Matched Layer (PML) [3]. The user selects one of several pre-built scenes (lenses, optical fibers, refraction interfaces, a Shack–Hartmann-style wavefront sensor, etc.) defined in config.py; the simulation runs in a Jupyter notebook (main_gpu.ipynb) and displays a live false-color animation of the electric field $E_z$ superimposed on a rendering of the material distribution. All quantities are in normalized (dimensionless) units with $c = \varepsilon_0 = \mu_0 = 1$.

1. Purpose and Scope

The program is a visual, interactive 2D electromagnetic wave simulator. It is intended for demonstrating and exploring wave-optics phenomena — refraction, focusing by lenses, guiding in optical fibers, anti-reflection structures, atmospheric distortion and wavefront sensing — rather than for quantitative engineering analysis.

What it does:

  • Integrates the 2D TM$_z$ Maxwell equations on a uniform grid (default $800 \times 416$ cells) with heterogeneous relative permittivity $\varepsilon_r(x,y)$, relative permeability $\mu_r(x,y)$, and electric conductivity $\sigma(x,y)$.
  • Runs the update kernels on the GPU (WGSL compute shaders), keeping all field state in GPU buffers; only $E_z$ is read back each frame for display.
  • Absorbs outgoing waves with a split-field PML [3] (an "anechoic chamber").
  • Injects point sources: continuous sine, Gaussian pulse, or Gaussian-modulated sine burst.

What it does NOT do:

  • No 3D, no dispersive or nonlinear media, no far-field transforms, no physical-unit input/output (everything is normalized).
  • No TE$_z$ polarization (only $E_z$, $H_x$, $H_y$ are evolved).
  • No file output of field data; the only outputs are the on-screen animation and timing printouts (the output/ folder holds manually saved screen recordings).

2. Theory and Mathematical Model

2.1 Governing equations

In a 2D problem invariant along $z$, Maxwell's curl equations [5] decouple into two polarizations. The program solves the TM$_z$ set, involving the scalar field components $E_z(x,y,t)$, $H_x(x,y,t)$, $H_y(x,y,t)$:

$$ \varepsilon \frac{\partial E_z}{\partial t} + \sigma E_z = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} $$

$$ \mu \frac{\partial H_x}{\partial t} = -\frac{\partial E_z}{\partial y}, \qquad \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} $$

where $\varepsilon = \varepsilon_0 \varepsilon_r$ is the permittivity, $\mu = \mu_0 \mu_r$ the permeability and $\sigma$ the electric conductivity. Unit system: the code uses normalized units $\varepsilon_0 = \mu_0 = 1$, hence the speed of light $c = 1/\sqrt{\varepsilon_0\mu_0} = 1$; the shorter side of the domain has length $1$ (the grid step is dy = 1/max(Nx,Ny), config.py:97). All symbols below are therefore dimensionless; where a "unit" is quoted it is the normalized one (lengths in units of the domain height, times in units of a light-crossing time of the domain height, frequencies in inverse crossing times, $\sigma$ in inverse times since $\varepsilon_0=1$).

2.2 Bérenger split-field PML

To terminate the grid without reflections, the outermost cells use Bérenger's perfectly matched layer [3]. $E_z$ is split into two non-physical sub-components,

$$ E_z = E_{zx} + E_{zy}, $$

each damped by an independent, direction-selective conductivity $\sigma_x(x)$ (ramped at the left/right boundaries) or $\sigma_y(y)$ (ramped at the top/bottom boundaries):

$$ \varepsilon \frac{\partial E_{zx}}{\partial t} + \sigma_x E_{zx} = \frac{\partial H_y}{\partial x}, \qquad \varepsilon \frac{\partial E_{zy}}{\partial t} + \sigma_y E_{zy} = -\frac{\partial H_x}{\partial y}. $$

A wave impinging on a vertical boundary varies along $x$ and thus feeds mainly $E_{zx}$, which $\sigma_x$ absorbs; a match of impedance at the PML interface makes the layer reflectionless in the continuum limit [3]. The split fields are stored in Ezix/Eziy (see the comment at wgpu_helper.py:99-101). The magnetic updates carry no magnetic loss in this implementation ($\sigma^{(*)} = 0$; coefficients c3 = c4 = dt/ur, util.py:163-164), so the PML is not perfectly matched in the discrete sense — it works because the conductivity profile is graded smoothly (Section 3.3) and the incident energy is absorbed over many cells.

The conductivity profile inside a layer of $d$ cells follows the standard polynomial grading [2, 3]

$$ \sigma(u) = \sigma_{\max}\left(\frac{u}{d}\right)^{m}, $$

with $u$ the depth into the layer in cells. The defaults are $d = 60$ cells, $\sigma_{\max} = 1000$, $m = 4$ (config.py:115, util.py:13-28).

2.3 Material model

Four $N_y \times N_x$ arrays describe the medium:

  • $\varepsilon_r[i,j]$ (er) — relative permittivity; lenses, fibers and interfaces are built by painting geometric shapes or PNG images into this array (generateEpsilon(), util.py:254).
  • $\mu_r[i,j]$ (ur) — relative permeability; always $1$ in the provided scenes.
  • $\sigma_x[i,j], \sigma_y[i,j]$ (sx, sy) — split conductivities. Besides the PML they are also used to build opaque/absorbing objects (apertures, screens, detector masks) by adding large values ($\sim 150$–$550$) inside shapes.
  • $w[i,j]$ (wire) — a multiplicative mask applied to the neighbor fields inside the curl terms; $w = 1$ is free space, $w = 0$ suppresses the cell's contribution to its neighbors' updates, approximating a perfect electric conductor. All provided scenes use $w \equiv 1$.

2.4 Sources

Sources are soft (additive) point sources: a time function $v(t)$ is added to $E_z$ at a fixed grid point each step (kernel.wgsl:103-123, half to each split component). Three kinds are implemented in generate_source() (kernel.wgsl:43-60), selected by the integer kind with parameter vector $(p_1,p_2,p_3,p_4)$:

kind Model Formula Parameters
1 continuous sine $v(t) = p_1 \sin\left(2\pi p_2 t\right)$ $p_1$ amplitude, $p_2$ frequency
2 Gaussian pulse $v(t) = p_1 \exp\left(-\dfrac{(t-p_2)^2}{2 p_3}\right)$ $p_1$ amplitude, $p_2$ center time $t_0$, $p_3$ variance $s^2$
3 Gaussian-modulated sine (burst) $v(t) = p_1 \sin\left(2\pi p_2 t\right)\exp\left(-\dfrac{(t-p_3)^2}{2 p_4}\right)$ $p_1$ amplitude, $p_2$ frequency, $p_3$ center time, $p_4$ variance

Note that $p_3$ (kind 2) and $p_4$ (kind 3) are the Gaussian variance $s^2$, not the standard deviation. A source of frequency $f$ in these units radiates a wavelength $\lambda = c/f = 1/f$ domain-height units; e.g. $f = 40$ gives $\lambda = 1/40$ of the domain height.

2.5 Symbol ↔ code variable table

Symbol Code variable Meaning Units (normalized)
$E_{zx}, E_{zy}$ Ezix, Eziy (Ezx,Ezy in the notebook) split components of $E_z$ [3] field units
$E_z$ Ezi (= Ezix + Eziy) out-of-plane electric field field units
$H_x, H_y$ Hyi, Hxi (see note) in-plane magnetic field field units
$\varepsilon_r$ er relative permittivity – (range 1–10 in scenes)
$\mu_r$ ur relative permeability – (always 1)
$\sigma_x, \sigma_y$ sx, sy split electric conductivities 1/time
$w$ wire PEC mask – (0 or 1)
$\Delta x, \Delta y$ dx, dy grid steps (always equal) length; default $1/800$
$\Delta t$ dt time step time; default $\approx 4.42\times10^{-4}$
$N_x, N_y, N_t$ Nx, Ny, Nt grid columns, rows, time steps –
$c_{1x}, c_{1y}$ c1x, c1y lossy-update decay coefficients –
$c_{2x}, c_{2y}$ c2x, c2y lossy-update curl coefficients time
$c_3, c_4$ c3, c4 magnetic update coefficients $\Delta t/\mu_r$ time
$t$ t simulation time (GPU-side accumulator) time
$d, \sigma_{\max}, m$ depth, sigma_max, m in PML() PML thickness (cells), peak conductivity, grading order –, 1/time, –

Naming note. The kernel's arrays Hxi, Hyi correspond to $-H_y$ and $-H_x$ of the conventional orientation (rows = $y$, columns = $x$): Hxi is updated from the $E_z$ difference along columns and Hyi from the difference along rows (kernel.wgsl:92-98). Likewise the labels dx/dy inside the kernel are swapped with respect to row/column differences. Because the code always sets dx = dy and the media are isotropic, this relabeling has no numerical effect, but it must be kept in mind when modifying the kernel.

3. Numerical Methods

3.1 Yee leapfrog scheme

The equations of Section 2 are discretized with Yee's staggered leapfrog scheme [1, 2]: second-order centered differences in space and time, with $E$ and $H$ evaluated at alternating half steps. Writing $\Delta = \Delta x = \Delta y$, row index $i$, column index $j$, and time level superscript $(n)$, the electric updates (leap_frog1, kernel.wgsl:63-80) are

$$ E_{zx}^{(n+1)}[i,j] = c_{1x}[i,j], E_{zx}^{(n)}[i,j] - c_{2x}[i,j], \frac{\tilde H^{(n+\frac{1}{2})}[i,j] - \tilde H^{(n+\frac{1}{2})}[i,j-1]}{\Delta} $$

$$ E_{zy}^{(n+1)}[i,j] = c_{1y}[i,j], E_{zy}^{(n)}[i,j] + c_{2y}[i,j], \frac{\bar H^{(n+\frac{1}{2})}[i,j] - \bar H^{(n+\frac{1}{2})}[i-1,j]}{\Delta} $$

for $1 \le i \le N_y-1$, $1 \le j \le N_x-1$, where $\tilde H$ = Hxi $= -H_y$ and $\bar H$ = Hyi $= -H_x$. The magnetic updates (leap_frog2, kernel.wgsl:83-100) then use the freshly updated $E_z = E_{zx}+E_{zy}$:

$$ \tilde H^{(n+\frac{3}{2})}[i,j] = \tilde H^{(n+\frac{1}{2})}[i,j] + c_3[i,j], \frac{E_z^{(n+1)}[i,j] - E_z^{(n+1)}[i,j+1]}{\Delta} $$

$$ \bar H^{(n+\frac{3}{2})}[i,j] = \bar H^{(n+\frac{1}{2})}[i,j] + c_4[i,j], \frac{E_z^{(n+1)}[i+1,j] - E_z^{(n+1)}[i,j]}{\Delta} $$

for $0 \le i \le N_y-2$, $0 \le j \le N_x-2$. In the shader, one thread per cell computes exactly these expressions; the wire mask multiplies each neighbor read:

// kernel.wgsl:71-79  (leap_frog1 — E update)
let G1x =  c1x[idx_row+1][idx_col+1]*Ezix_1[idx_row+1][idx_col+1];
let G4 = wire[idx_row+1][idx_col+1]*Hxi_1[idx_row+1][idx_col+1];
let G5 = wire[idx_row+1][idx_col+0]*Hxi_1[idx_row+1][idx_col+0];
Ezix_2[idx_row+1][idx_col+1] = G1x + c2x[idx_row+1][idx_col+1]*( -(G4 - G5)/dy );

3.2 Semi-implicit treatment of loss and the update coefficients

The conduction term $\sigma E$ is time-averaged over the two levels (semi-implicit, unconditionally stable for the loss term [2]):

$$ \varepsilon_r \frac{E^{(n+1)} - E^{(n)}}{\Delta t} + \sigma \frac{E^{(n+1)} + E^{(n)}}{2} = (\nabla \times \mathbf{h})_z $$

Solving for $E^{(n+1)}$ yields the coefficients computed once on the CPU in leapFrog2D_constants() (util.py:151-166):

$$ c_{1x} = \frac{1 - \dfrac{\sigma_x \Delta t}{2\varepsilon_r}}{1 + \dfrac{\sigma_x \Delta t}{2\varepsilon_r}}, \qquad c_{2x} = \frac{\Delta t / \varepsilon_r}{1 + \dfrac{\sigma_x \Delta t}{2\varepsilon_r}}, \qquad c_3 = c_4 = \frac{\Delta t}{\mu_r}, $$

and identically for $c_{1y}, c_{2y}$ with $\sigma_y$:

# util.py:153-164
ax = (1.0 - dt*sx/(2.0*er));  bx = (1.0 + dt*sx/(2.0*er))
c1x = ax/bx
c2x = (dt/er)/bx
c3 = dt/ur

3.3 Stability (CFL condition)

The 2D Yee scheme is stable when the Courant–Friedrichs–Lewy condition [4, 6] holds:

$$ c,\Delta t \sqrt{\frac{1}{\Delta x^{2}} + \frac{1}{\Delta y^{2}}} \le 1 . $$

The default configuration uses half the limit — Courant number $S = 0.5$ — to keep the high-$\sigma$ PML cells stable (config.py:100-101):

CLF = np.sqrt(1/S['dx']**2 + 1/S['dy']**2)*2.0  # *2.0 to prevent instability in the PML borders
S['dt'] = 1.0 / CLF

The optical ckt scene (config.py:519-521) runs at the exact limit $S = 1$ instead. With $\Delta = 1/800$ the default step is $\Delta t = 1/(2\sqrt{2}\cdot 800) \approx 4.42\times10^{-4}$. The scheme is second-order accurate in space and time [2]; numerical dispersion follows the standard 2D FDTD dispersion relation [2, ch. 4].

3.4 PML discretization

PML() (util.py:13-28) builds the two graded conductivity maps by mirroring a 1D ramp $\left(u/d\right)$, raised to the power $m$ and scaled by $\sigma_{\max}$, into the four boundary strips:

# util.py:23-24
sigma_x = sigma_max*sx**m     # ramps in the first/last `depth` rows
sigma_y = sigma_max*sy**m     # ramps in the first/last `depth` columns

Callers assign the returned pair as S['sy'], S['sx'] = PML(...) (config.py:115) — i.e. the first returned array (row-ramped) becomes $\sigma_y$ and the second (column-ramped) becomes $\sigma_x$, consistent with the kernel convention of Section 2.5.

3.5 Time-marching, double buffering and GPU dispatch

One time step (WGPUHelper.compute(), wgpu_helper.py:66-71) dispatches three compute passes in order:

  1. update_source — single thread; advances $t \mathrel{+}= \Delta t$ and adds $v(t)/2$ to both split components of $E_z$ in both buffer copies (kernel.wgsl:102-123).
  2. leap_frog1 — one thread per cell; reads $E^{(n)}$/$H^{(n+1/2)}$ from buffer set A, writes $E^{(n+1)}$ to buffer set B.
  3. leap_frog2 — reads the new $E^{(n+1)}$ (set B) and old $H$ (set A), writes $H^{(n+3/2)}$ to set B.

To avoid read/write data races inside a pass, every field is stored twice (wgpu_Ezix_1/_2, etc.) and two pre-built bind groups swap the read/write roles on even/odd step number n (wgpu_helper.py:330-341) — a classic ping-pong (double-buffer) scheme. Threads are launched in 2D workgroups of workgroup_size × workgroup_size (16 by default, 4 in the notebook), with 1 + Ny//wg × 1 + Nx//wg workgroups; out-of-range threads return immediately (kernel.wgsl:67-69). Grid dimensions [[N]], [[M]] and [[workgoup_size]] are placeholders substituted into the WGSL source at pipeline-build time (wgpu_helper.py:122-125), because WGSL array sizes must be compile-time constants.

3.6 Scene construction (rasterization of materials)

Material maps are built on the CPU by painting primitives into $N_y \times N_x$ arrays:

  • Lens (generateEpsilon, type 'lens', util.py:259-279): a biconvex lens of aperture radius $R$ and edge sagittae $r_\mathrm{left}, r_\mathrm{right}$ is the intersection of two disks; the surface-circle radii are recovered from the sagitta relation $x_c = \dfrac{R^2 - r^2}{2r}$, $C_r = x_c + r$, and the region $\left{(x-x_\text{center})^2 + (y-c_y)^2 < C_r^2\right}$ is filled with $\varepsilon_r$.
  • Lens2 (util.py:327-346): a lens defined by thickness and height, as the intersection of two circles of radius $R = \dfrac{\tau^2 + h^2}{2\tau}$ ($\tau$ = lensTickness, $h$ = lensHeight), optionally rotated by angle.
  • Box, interface, image: axis-aligned rectangles; oriented half-planes $-\left(p_x\cos a - p_y\sin a\right) < 0$; and grayscale PNG images (red channel, flipped vertically, scaled) placed by placeSubImage().
  • Smoothing: edges are anti-aliased by $n$ passes of a normalized $6\times6$ disk-shaped convolution kernel (smooth_image(), util.py:31-39), which reduces staircase scattering artifacts.

4. Code Architecture

File Purpose
main_gpu.ipynb Entry point (Jupyter). Selects the scene, builds initial arrays and coefficients, sets up the GPU, runs the time loop, animates $E_z$.
config.py Scene library: configSimulation(name) dispatches to one config function per scene; each fills the dictionary S with grid, medium, PML and source definitions.
util.py CPU-side numerics and helpers: PML builder, update coefficients, material rasterizers (generateEpsilon, generateConductivity), legacy CPU source generator, image smoothing, colormaps and frame composition. Defines the GPUSource NumPy dtype matching the WGSL Source struct.
wgpu_helper.py WGPUHelper class: WebGPU device setup, buffer creation, shader compilation (with placeholder substitution), pipelines, double-buffer bind groups, dispatch, and host↔GPU transfer.
kernel.wgsl GPU compute shaders: leap_frog1 ($E$ update), leap_frog2 ($H$ update), update_source (time accumulator and source injection), generate_source (waveform library).
images/*.png Grayscale masks rasterized into $\varepsilon_r$ or $\sigma$ (waveguides, atmosphere phase screens, detector, anti-reflex crystal, …).
output/ Manually saved animation recordings (.mp4), not written by the code.

Data flow:

flowchart LR
    A["configSimulation(name)<br/>config.py"] -->|"S: Nx,Ny,Nt,dx,dy,dt,<br/>er,ur,sx,sy,wire,sources"| B["leapFrog2D_constants()<br/>util.py"]
    A -->|images/*.png| A
    B -->|"c1x,c1y,c2x,c2y,c3,c4"| C["WGPUHelper.setup_compute()<br/>wgpu_helper.py"]
    C -->|"buffers + pipelines"| D["GPU: kernel.wgsl<br/>update_source → leap_frog1 → leap_frog2"]
    E["notebook time loop<br/>n = 1..Nt"] -->|"compute(n)"| D
    D -->|"sync_gpu_host_data(): Ez"| F["addEzToComposite()<br/>composeImage(), colormaps<br/>util.py"]
    F -->|RGB frame| G["matplotlib window<br/>(every 10 steps)"]
Loading

5. All Functions Details

5.1 config.py

def configSimulation(name, anecroic=True)
  • description: Top-level scene dispatcher. Builds the default simulation dictionary via config_default() and then applies the named scene configurator (switch-case on name). Prints a warning for unknown names and returns the default scene.
  • parameters:
    • name : str — scene identifier; one of 'single lens', 'laser pulse', 'antireflex', 'opticalfiber', 'optical ckt', 'dual optcal fiber', 'refraction', 'wavefront sensor' (the branches 'single lens plane' and 'yagi_antenna' reference functions that are currently commented out and would raise NameError).
    • anecroic : bool — passed to config_default(); kept for backward compatibility with the legacy cosine absorber (the current PML ignores it).
  • return: S : dict — complete simulation description with keys Nx, Ny, Nt (int), dx, dy, dt (float), er, ur, sx, sy, wire (np.ndarray (Ny,Nx) float), sources (np.ndarray of dtype GPUSource).
  • example: S = configSimulation('refraction') builds the double-slit + prism-interface demo used in the notebook.
def config_default(anecroic=True)
  • description: Creates the baseline scene: $800\times416$ grid, $N_t = 3500$ steps, $\Delta x = \Delta y = 1/800$, $\Delta t$ at Courant number 0.5 (config.py:100-101), PML of 60 cells with $\sigma_{\max}=1000$, $m=4$, vacuum medium (er = ur = wire = 1), no sources.
  • parameters:
    • anecroic : bool — legacy switch for the disabled cosine boundary absorber; has no effect on the active PML path.
  • return: S : dict — default scene (same key layout as above; sources must be added by a scene function).
def config_single_lens(S)
  • description: Adds one biconvex lens ($\varepsilon_r = 3$, heavily smoothed with 100 passes) centered at normalized $(c_x, c_y) = (1.0, 0.5)$ and one continuous sine point source (kind 1, amplitude 1.2, frequency 40) at grid point $(0.5 N_y,\ 0.2 N_x)$.
  • parameters:
    • S : dict — scene dictionary from config_default(); modified in place.
  • return: S : dict — updated scene.
def config_laser_pulse(S)
  • description: Single Gaussian-modulated sine source (kind 3): amplitude 2.0, frequency 40, pulse center $t_0 = 0.1$, variance $6\times10^{-4}$; empty medium. Demonstrates a propagating pulse absorbed by the PML.
  • parameters / return: as config_single_lens.
def config_antireflex(S)
  • description: Anti-reflection demo: an image-defined absorbing divider (images/antireflex_s.png, $\sigma$ scale 150), four conductive boxes forming two apertures, two weak lenses ($\varepsilon_r = 1.1$), an image-defined crystal slab (images/antireflex_er2.png, $\varepsilon_r$ scale 1.47), and two identical sine sources (kind 1, amplitude 1.2, frequency 40) illuminating the two halves.
  • parameters / return: as config_single_lens.
def config_optical_fiber(S)
  • description: Step-index fiber demo, $N_t = 2500$: two conductive boxes ($\sigma = 550$) form the entrance aperture, an image lens (images/lens.png, $\varepsilon_r$ scale 4) focuses into a rectangular core ($\varepsilon_r = 3$) extending to the right edge; one sine source (kind 1, amplitude 1.2, frequency 40).
  • parameters / return: as config_single_lens.
def config_optical_ckt(S)
  • description: Largest scene: rebuilds the grid at $1802\times802$, $N_t = 25000$, $\Delta t$ at the exact CFL limit, PML depth 80. Waveguide layout from images/waveguide.png ($\varepsilon_r$ scale 17) and a curved fiber from images/fiber_curve.png (scale 6); aperture boxes with $\sigma = 350$; one sine source (kind 1, amplitude 3.0, frequency 20).
  • parameters / return: as config_single_lens (note: it replaces Nx, Ny, dt and all medium arrays).
def config_dual_optical_fiber(S)
  • description: Two crossing step-index fibers ($\varepsilon_r = 2.2$) fed through a small lens ($\varepsilon_r = 4$) behind a three-box aperture ($\sigma = 550$); one sine source (kind 1, amplitude 1.2, frequency 40). Demonstrates evanescent coupling between cores.
  • parameters / return: as config_single_lens.
def config_refraction(S)
  • description: The notebook's "mega demo", $N_t = 6000$: a double-aperture screen ($\sigma = 550$), two lenses ($\varepsilon_r = 3$) collimating the beam, then a $45°$ tilted half-plane interface raising $\varepsilon_r$ by a factor 10 and a second tilted interface restoring it (Section 3.6, type 'interface'). Two co-located sine sources of different frequencies (kind 1: amplitudes 2.2 and 1.2, frequencies 40 and 40 — the comments suggest 10 and 80 for a "rainbow" dispersion variant). Demonstrates Snell refraction of the beam.
  • parameters / return: as config_single_lens.
def wavefront_sensor(S)
  • description: Shack–Hartmann-style wavefront sensor: a column of five plano-convex lenslets ('lens2', $\varepsilon_r = 2$) held by conductive support boxes ($\sigma = 550$) at $x = 1.6$, an image-defined detector screen (images/detector.png, $\sigma$ scale 550) and an atmospheric phase screen (images/atmosfere2.png, $\varepsilon_r$ scale 2). Grid rescaled by factor sc = 0.95 ($760\times395$). One sine source (kind 1, amplitude 0.8, frequency 80). The dictionary is mutated in place; configSimulation ignores its return value, which works because dict is mutable.
  • parameters:
    • S : dict — scene dictionary; rebuilt in place (grid, PML, medium and sources are all replaced).
  • return: S : dict — the same dictionary.

5.2 util.py

GPUSource = np.dtype([('position', np.int32, (2,)), ('kind', np.int32),
                      ('_', np.int32), ('parameters', np.float32, (4,))])
  • description: NumPy structured dtype mirroring the WGSL struct Source byte layout (kernel.wgsl:13-18): vec2<i32> position (row, column), i32 kind, one i32 of explicit padding so parameters lands on a 16-byte boundary, vec4<f32> parameters. Total 24 bytes → padded struct compatible with the storage buffer.
def PML(Nx, Ny, depth, sigma_max, m)
  • description: Builds the two polynomial-graded PML conductivity maps of Section 3.4: $\sigma(u) = \sigma_{\max}(u/d)^m$ ramping from 0 at the inner PML edge to $\sigma_{\max}$ at the outer grid edge, in all four boundary strips.
  • parameters:
    • Nx : int — grid columns.
    • Ny : int — grid rows.
    • depth : int — PML thickness $d$ in cells (60–80 in the scenes).
    • sigma_max : float — peak conductivity $\sigma_{\max}$ (1000 in the scenes), normalized units.
    • m : int — grading polynomial order $m$ (4 in the scenes; 3–4 is the standard recommendation [2]).
  • return: (sigma_x, sigma_y) : tuple of np.ndarray (Ny,Nx) — first array ramps in the first/last depth rows (absorbs at top/bottom), second in the first/last depth columns (left/right). Callers assign them as S['sy'], S['sx'] = PML(...).
  • example: sy, sx = PML(800, 416, 60, 1000, 4).
def smooth_image(I, n)
  • description: Applies n passes of a normalized $6\times6$ disk-shaped moving-average convolution (via scipy.signal.convolve2d, 'same' mode) to soften material edges and reduce staircase diffraction artifacts.
  • parameters:
    • I : np.ndarray (Ny,Nx) — image/material map.
    • n : int — number of smoothing passes (0 = no-op).
  • return: np.ndarray (Ny,Nx) float — smoothed map.
def fastInd2RGB(Iidx, M)
  • description: Vectorized indexed-color lookup: clips an index image to $[0,255]$ and maps each pixel through a $255\times3$ colormap by flat fancy-indexing. Note the cast is to np.int8 (signed); values above 127 wrap negative and index the LUT from its end, which the arctan-saturated colormaps tolerate visually.
  • parameters:
    • Iidx : np.ndarray (Ny,Nx) — scalar index image.
    • M : np.ndarray (255,3) float — RGB colormap in $[0,1]$.
  • return: np.ndarray (Ny,Nx,3) float — RGB image.
def addEzToComposite(Ii, Ez, M)
  • description: Produces one animation frame: rescales the field as $E_s = 40 E_z + 130$ (bias 130 ≈ mid-LUT for $E_z=0$), maps it through colormap M, and adds it to the static material image Ii with clipping to $[0,255]$.
  • parameters:
    • Ii : np.ndarray (Ny,Nx,3) uint8 — precomputed material composite (from composeImage).
    • Ez : np.ndarray (Ny,Nx) float — field snapshot (the notebook passes 100.0*Ezi).
    • M : np.ndarray (255,3) float — diverging colormap (M1: red for negative, blue for positive).
  • return: np.ndarray (Ny,Nx,3) uint8 — displayable RGB frame.
  • example: I = addEzToComposite(Imaterial, 100.0*Ezi, M1) (notebook display cell).
def composeImage(s, wire, er, M2, M3, M4)
  • description: Renders the static material overlay shown under the field: conductivity as $s/2 + 128$ through M2 (yellow/blue), wire mask as $255(1-w)$ through M3 (grayscale), permittivity as $255(\varepsilon_r - 1)/2$ through M4 (cyan); averages and clips the three layers.
  • parameters:
    • s : np.ndarray (Ny,Nx) — conductivity map to display (the notebook passes (sx+sy)/2).
    • wire : np.ndarray (Ny,Nx) — PEC mask.
    • er : np.ndarray (Ny,Nx) — relative permittivity.
    • M2, M3, M4 : np.ndarray (255,3) — colormaps from customColormaps().
  • return: np.ndarray (Ny,Nx,3) uint8 — material image.
def customColormaps()
  • description: Builds four 255-entry RGB lookup tables from arctangent sigmoids $\dfrac{\pi/2 + \arctan\left(\beta, (i - i_0)/255\right)}{\pi}$ with sharpness $\beta = 8$ and pivots $i_0 \in {60, 128, 196}$, then normalizes each by its global maximum (MATLAB-style). M1: field colormap (red↔blue diverging); M2: conductivity (blue→yellow); M3: grayscale; M4: permittivity (black→cyan).
  • parameters: none.
  • return: (M1, M2, M3, M4) : tuple of np.ndarray (255,3) float — colormaps in $[0,1]$.
def placeSubImage(Nx, Ny, cx, cy, Ii)
  • description: Pastes a (vertically flipped) sub-image centered at pixel $(c_x, c_y)$ into a zero background of size $N_y\times N_x$. If the sub-image would overflow any border, prints 'Image place forbidden' and returns an empty array.
  • parameters:
    • Nx, Ny : int — target grid size.
    • cx, cy : float — center in pixels (callers pass cx*Nx, cy*Ny).
    • Ii : np.ndarray (Hi,Wi) — sub-image to place.
  • return: np.ndarray (Ny,Nx) float — background with the sub-image inserted (or empty array on overflow).
def leapFrog2D_constants(er, ur, sx, sy, dt)
  • description: Precomputes the update coefficients of Section 3.2 from the material maps: $c_{1x,y}$ (decay), $c_{2x,y}$ (curl gain), $c_3 = c_4 = \Delta t/\mu_r$. This is where $\varepsilon_r$, $\sigma$ and $\Delta t$ enter the GPU scheme; the kernel itself only sees these six arrays.
  • parameters:
    • er, ur : np.ndarray (Ny,Nx) — relative permittivity/permeability.
    • sx, sy : np.ndarray (Ny,Nx) — split conductivities (normalized).
    • dt : float or np.ndarray shape (1,) — time step.
  • return: (c1x, c1y, c2x, c2y, c3, c4) : tuple of np.ndarray (Ny,Nx) — see the symbol table (Section 2.5).
  • example: c1x, c1y, c2x, c2y, c3, c4 = leapFrog2D_constants(er, ur, sx, sy, np.array([dt])).
def generateConductivity(Nx, Ny, cond_type, parameters)
  • description: Rasterizes a conductivity map $\sigma[i,j]$ of one of three types. 'boundary': legacy cosine-taper frame absorber (superseded by PML()); 'image': grayscale PNG (red channel, flipped, resized by imScale, smoothed) scaled to scale and placed at $(c_x N_x, c_y N_y)$; 'box': axis-aligned rectangle of value value with optional smoothing. Box coordinates are normalized to the domain height ($\times N_y$) for both axes, so $x \in [0, N_x/N_y]$ (≈ $[0,1.92]$ on the default grid).
  • parameters:
    • Nx, Ny : int — grid size.
    • cond_type : str — 'boundary' | 'image' | 'box'.
    • parameters : dict — per type: boundary: L, maxS, value; image: name (path), cx, cy (normalized center), scale (σ of a white pixel), imScale (image resize factor), optional smooth; box: upper_left_x/y, lower_right_x/y (normalized by $N_y$), value, optional smooth (passes).
  • return: np.ndarray (Ny,Nx) float — conductivity contribution (callers add it to S['sx']/S['sy']).
  • example: S['sx'] += generateConductivity(800, 416, 'box', {'upper_left_x':0.5, 'upper_left_y':0.01, 'lower_right_x':0.57, 'lower_right_y':0.42, 'value':550, 'smooth':10}) builds the upper aperture blade of the refraction scene.
def generateEpsilon(Nx, Ny, eps_type, parameters)
  • description: Rasterizes a relative-permittivity map. Types: 'lens' (biconvex lens from two-disk intersection, Section 3.6; background 1, lens body er); 'image' (PNG placed as $1 + \text{scale}\cdot I/255$); 'box' (rectangle of value on a background of 1); 'interface' (tilted half-plane step of height amplitude plus 1); 'lens2' (thickness/height-defined lens, background 0 — callers add it). Several other types ('plane_x', 'percolation') exist only as commented-out MATLAB legacy.
  • parameters:
    • Nx, Ny : int — grid size.
    • eps_type : str — 'lens' | 'image' | 'box' | 'interface' | 'lens2'.
    • parameters : dict — per type: lens: cx, cy, r_left, r_right, R, er, optional smooth; image: name, cx, cy, scale, imScale, optional smooth; box: corner keys as in generateConductivity plus value; interface: point_x, point_y, angle (rad), amplitude, optional smooth; lens2: cx, cy, lensTickness, lensHeight, angle, er, optional smooth. All coordinates normalized by $N_y$.
  • return: np.ndarray (Ny,Nx) float — permittivity map (multiplied, added or assigned by the caller depending on the scene).
  • example: er = generateEpsilon(800, 416, 'lens', {'cx':1.0,'cy':0.5,'r_left':0.15,'r_right':0.15,'R':0.5,'er':3,'smooth':100}).
def generateSource(Nx, Ny, source_type, parameters)
  • description: Legacy CPU source generator (unused in the GPU pipeline, kept from the pre-GPU version): returns a full grid with the source value written at one point (or one column for the plane types), to be added to the field each step. Types: 'pulse', 'burst', 'sine', 'plain_square_x', 'plain_sine_x'; several more remain as commented MATLAB code. In the current program sources are instead injected on the GPU (update_source, kernel.wgsl:102).
  • parameters:
    • Nx, Ny : int — grid size.
    • source_type : str — waveform/geometry type.
    • parameters : dict — cx, cy (normalized by $N_y$), amplitude, time, and per type frequency, t0, width (std-dev), offset (row margin for the plane types).
  • return: np.ndarray (Ny,Nx) float — field increment for the current time.

5.3 wgpu_helper.py — class WGPUHelper

def __init__(self, workgroup_size=16, shader_source='kernel.wgsl')
  • description: Stores configuration and initializes all host- and GPU-side field references to None. No GPU work happens here.
  • parameters:
    • workgroup_size : int — side of the square GPU workgroup (threads per group = workgroup_size**2; the notebook uses 4). Must divide the hardware limit (typically ≤ 256 threads/group).
    • shader_source : str — path to the WGSL kernel file.
  • return: WGPUHelper instance (constructor).
def setup_compute(self, c1x, c1y, c2x, c2y, c3, c4, dx, dy, dt, t, wire,
                  Ezix, Eziy, Hxi, Hyi, sources)
  • description: One-time GPU initialization: (1) casts and stores all arrays as float32 copies; (2) requests a high-performance WebGPU adapter and device, printing the adapter summary; (3) loads kernel.wgsl, substitutes the [[N]], [[M]], [[workgoup_size]] placeholders with the actual grid/workgroup sizes and compiles the shader module; (4) creates one storage buffer per array (coefficients, fields ×2 for double buffering, t, sources) and uniform buffers for dx, dy, dt; (5) builds three compute pipelines (leap_frog1, leap_frog2, update_source) over a shared 20-binding layout; (6) builds the two ping-pong bind groups that swap the _1/_2 field buffers (Section 3.5).
  • parameters:
    • c1x … c4 : np.ndarray (Ny,Nx) — update coefficients from leapFrog2D_constants().
    • dx, dy, dt, t : np.ndarray shape (1,) float — grid steps, time step, initial time (scalars wrapped in arrays so they map to GPU buffers).
    • wire : np.ndarray (Ny,Nx) — PEC mask.
    • Ezix, Eziy, Hxi, Hyi : np.ndarray (Ny,Nx) — initial fields (normally zeros).
    • sources : np.ndarray dtype=GPUSource — source table (uploaded via .tobytes()).
  • return: None (state stored on the instance).
  • example: wgpu_helper.setup_compute(c1x, c1y, c2x, c2y, c3, c4, dx, dy, dt, t, wire, Ezx, Ezy, Hx, Hy, S['sources']) (notebook cell).
def compute(self, n)
  • description: Advances the simulation one full time step by dispatching, in order: update_source (1 thread), leap_frog1 (E update, one thread/cell), leap_frog2 (H update). The step parity n selects which bind group (buffer set) is read vs. written.
  • parameters:
    • n : int — current step number (controls ping-pong parity; the source pass is always dispatched with parity 0 because it writes both buffer sets anyway).
  • return: None (GPU state advanced; work is submitted asynchronously to the queue).
def call_kernel(self, N_workers, M_workers, n, entry_point_idx)
  • description: Low-level dispatch: creates a command encoder and compute pass, selects pipeline entry_point_idx (0 = leap_frog1, 1 = leap_frog2, 2 = update_source), binds bind group 1 for even n / 2 for odd n, and dispatches (1 + N_workers//wg, 1 + M_workers//wg) workgroups.
  • parameters:
    • N_workers : int — number of rows to cover ($N_y$).
    • M_workers : int — number of columns to cover ($N_x$).
    • n : int — step parity selector.
    • entry_point_idx : int — pipeline index (see above).
  • return: None.
def sync_gpu_host_data(self)
  • description: Blocking read-back of the two split-field buffers (wgpu_Ezix_1, wgpu_Eziy_1), reshaped to the grid shape and summed to reconstruct $E_z = E_{zx} + E_{zy}$. Only the display path uses this; $H$ is never read back. Note it always reads buffer set 1, so on even steps the returned field is one step behind the freshest set (irrelevant for visualization).
  • parameters: none.
  • return: np.ndarray (Ny,Nx) float32 — current $E_z$ field.
def sync_host_gpu_data(self)
  • description: Uploads every host-side array (coefficients, dx, dy, dt, t, wire, and both copies of all four fields) to its GPU buffer via queue.write_buffer. Called once after setup_compute(); can be reused to hot-reset the simulation state.
  • parameters: none.
  • return: None.

5.4 kernel.wgsl (GPU compute shaders)

fn generate_source(kind : i32, parameters : vec4<f32>, t : f32) -> f32
  • description: Waveform library evaluating $v(t)$ for the three source kinds of Section 2.4; returns 0.0 for unknown kinds.
  • parameters:
    • kind : i32 — 1 sine, 2 Gaussian, 3 Gaussian-modulated sine.
    • parameters : vec4<f32> — $(p_1, p_2, p_3, p_4)$ as in the table of Section 2.4 (unused slots are NaN in the configs and never read).
    • t : f32 — current simulation time.
  • return: f32 — source amplitude $v(t)$.
@compute fn leap_frog1(@builtin(global_invocation_id) index: vec3<u32>)
  • description: Electric-field update (first leapfrog half, Section 3.1). Thread $(\texttt{index.x}, \texttt{index.y}) = (i,j)$ updates cell $[i{+}1, j{+}1]$ of Ezix_2/Eziy_2 from Ezix_1/Eziy_1 and the wire-masked magnetic neighbors in Hxi_1/Hyi_1. Threads with $i \ge N-1$ or $j \ge M-1$ return immediately (boundary guard).
  • parameters: index : vec3<u32> — global thread id (.x = row, .y = column).
  • return: none (writes storage buffers at bindings 15–16).
@compute fn leap_frog2(@builtin(global_invocation_id) index: vec3<u32>)
  • description: Magnetic-field update (second leapfrog half). Reconstructs $E_z = E_{zx}+E_{zy}$ from the freshly written set-2 buffers (masked by wire) and increments Hxi_2/Hyi_2 from Hxi_1/Hyi_1 with the centered differences of Section 3.1. (K4 recomputes the same value as K1; redundant but harmless.)
  • parameters / return: as leap_frog1 (writes bindings 17–18).
@compute fn update_source(@builtin(global_invocation_id) index: vec3<u32>)
  • description: Single-thread pass run once per step: advances the GPU-resident clock t += dt, then loops over the source table and adds $v(t)/2$ to Ezix and Eziy at each source's grid point in both buffer sets (so the injection survives the ping-pong swap). The loop bound is hard-coded to 2 sources (for (var i=0; i<2; i++), kernel.wgsl:110); see Section 9.
  • parameters: index : vec3<u32> — thread id (all but thread (0,0) return).
  • return: none (writes bindings 10–12, 15–16).

6. Inputs and Outputs

Inputs

  • Scene name — string passed to configSimulation() in the notebook's third cell; the active scene is chosen by (un)commenting lines. Available names and what they demonstrate:

    Name Demo Notes
    'laser pulse' Gaussian burst in vacuum shortest run
    'single lens' point source focused by one lens
    'antireflex' anti-reflection coating comparison uses 2 PNG masks
    'opticalfiber' step-index fiber coupling $N_t = 2500$
    'optical ckt' integrated optical circuit $1802\times802$ grid, $N_t = 25000$, longest
    'dual optcal fiber' coupled fibers (sic — note the typo in the key)
    'refraction' beam refraction at tilted interfaces notebook default, $N_t = 6000$
    'wavefront sensor' Shack–Hartmann lenslet array + atmosphere uses 2 PNG masks
    'single lens plane', 'yagi_antenna' — broken: their config functions are commented out (config.py:213, config.py:310)

  • PNG masks (images/*.png) — grayscale images whose red channel is rasterized into $\varepsilon_r$ or $\sigma$ (Section 3.6). White (255) maps to the full scale value, black to 0. They are flipped vertically on load.

  • Scene dictionary S — all physical/numerical parameters (see Sections 2.5 and 5.1). No CLI arguments and no input data files other than the PNGs.

Outputs

  • Live animation (matplotlib Qt window): every 10 steps the frame $I = \mathrm{clip}\left(255,\mathrm{LUT}!\left(40\cdot 100,E_z + 130\right) + I_\text{material}\right)$ is drawn — red/blue = negative/positive $E_z$, yellow = conductive objects, cyan = dielectric objects (Section 5.2). Nothing is written to disk.
  • Console: GPU adapter info at setup; total elapsed wall time at the end; optional per-stage timings when show_time = True in the main loop cell.
  • output/*.mp4 — pre-recorded screen captures of past runs, not produced by the code.

7. Dependencies and Environment

Dependency Used for Import site
Python 3 (≥3.9 recommended) host language –
numpy arrays, dtypes, coefficients all .py files
scipy signal.convolve2d (smoothing) util.py:3
matplotlib animation window notebook, util.py:2
Pillow (PIL) PNG loading/resizing util.py:4
wgpu (wgpu-py [8]) WebGPU device, buffers, pipelines wgpu_helper.py:1
rendercanvas imported by wgpu_helper.py:3 (unused in the compute path but required for the import to succeed) wgpu_helper.py:3
PyQt5 (or another Qt binding) %matplotlib qt5 backend notebook cell 2
Jupyter running main_gpu.ipynb –

No versions are pinned anywhere in the repository (no requirements.txt). The program is OS-agnostic wherever wgpu-py finds a WebGPU backend (Metal on macOS — the development platform, Vulkan on Linux, D3D12/Vulkan on Windows). No compilation step: the WGSL shader is compiled at runtime by the wgpu driver.

Install:

pip install numpy scipy matplotlib pillow wgpu rendercanvas pyqt5 jupyter

8. How to Run

  1. Open the notebook:

    cd "python/NOTEBOKS/FD-TD"
jupyter notebook main_gpu.ipynb

  2. Run the cells in order:

    • Cell 1 — imports.
    • Cell 2 — %matplotlib qt5 (pops the animation into a native window; replace with qt or another interactive backend if PyQt5 is unavailable).
    • Cell 3 — scene selection: uncomment exactly one S = configSimulation(...) line (default: 'refraction'); the same cell zeroes the fields and computes the update coefficients.
    • Cell 4 (WGPUHelper cell) — wgpu_helper = WGPUHelper(workgroup_size=4), setup_compute(...), sync_host_gpu_data(). Prints the selected GPU adapter.
    • Last cell — builds the colormaps and material composite, then runs the while n < Nt loop, drawing every 10th frame.

  3. Minimal working example (script form, no notebook):

    from wgpu_helper import *
from config import *
from util import *

S = configSimulation('laser pulse')
Z = lambda: np.zeros((S['Ny'], S['Nx']))
c = leapFrog2D_constants(S['er'], S['ur'], S['sx'], S['sy'], np.array([S['dt']]))

h = WGPUHelper(workgroup_size=4)
h.setup_compute(*c, np.array([S['dx']]), np.array([S['dy']]),
                np.array([S['dt']]), np.array([0.0]),
                S['wire'], Z(), Z(), Z(), Z(), S['sources'])
h.sync_host_gpu_data()
for n in range(1, S['Nt'] + 1):
    h.compute(n)
Ez = h.sync_gpu_host_data()   # final field, shape (Ny, Nx)

Expected runtime: on the development machine (Apple-silicon Mac, Metal backend) the default scenes ($800\times416$, 2500–6000 steps, read-back every step, draw every 10th) run in the order of tens of seconds to a few minutes; 'optical ckt' ($1802\times802$, 25000 steps) is roughly an order of magnitude longer. The per-step read-back of $E_z$ (sync_gpu_host_data() is called every iteration, not only on drawing steps) is a major part of the cost.

9. Validation and Limitations

Qualitative validation cases built into the scenes: PML absorption ('laser pulse' — the pulse should vanish at the borders without visible reflection), focusing consistent with the lensmaker geometry ('single lens'), Snell-law beam deviation at the $45°$ interfaces ('refraction'), total internal reflection guiding ('opticalfiber'). No quantitative regression tests exist.

Known pitfalls and limitations found in the code:

  1. Hard-coded source count. update_source loops over exactly 2 sources (kernel.wgsl:110). Scenes with one source (e.g. 'single lens', 'laser pulse', 'wavefront sensor') read sources[1] out of bounds; WGSL robust-buffer-access clamps the index, so the single source is effectively injected twice (doubled amplitude). Scenes with more than 2 sources would silently ignore the extras. The loop bound should be the source-array length (e.g. via arrayLength).
  2. Broken scene names. 'single lens plane' and 'yagi_antenna' dispatch to functions that are commented out and raise NameError (config.py:50-57).
  3. CFL safety factor. The default $\Delta t$ uses Courant number 0.5 explicitly to keep the strong PML stable (config.py:100); 'optical ckt' runs at Courant number 1.0 (config.py:519), the exact 2D stability limit — any added lossy structure there operates with no margin.
  4. Smoothing vs. background. As warned in config.py:35, smoothing $\varepsilon_r$ maps is dangerous because the background is 1, not 0; the disk kernel bleeds the background into objects and vice versa. Scene functions therefore combine maps carefully with + ... - 1 and * patterns.
  5. PEC mask is approximate. wire masks neighbor contributions inside the curls but does not zero the field storage itself; it is adequate for visualization, not a rigorous PEC boundary.
  6. No magnetic PML conductivity. $\sigma^{(*)} = 0$ (Section 2.2), so the layer is not discretely matched; residual reflection is suppressed by the deep ($60$–$80$ cell), high-order ($m=4$) grading.
  7. fastInd2RGB signed cast. The np.int8 cast (util.py:47) wraps indices above 127; visually masked by the saturated colormaps but incorrect as written (should be np.uint8).
  8. Display latency. sync_gpu_host_data() always reads buffer set 1 (wgpu_helper.py:346-347), which on even steps is one time level stale, and it is called every step even when no frame is drawn — a straightforward performance optimization target.
  9. Normalized units only. Mapping to physical problems requires rescaling lengths by the chosen domain size $L_0$ and times by $L_0/c$.

10. References

[1] K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302–307, 1966. DOI: 10.1109/TAP.1966.1138693.

[2] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed., Artech House, Norwood, MA, 2005.

[3] J.-P. Bérenger, "A perfectly matched layer for the absorption of electromagnetic waves," Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994. DOI: 10.1006/jcph.1994.1159.

[4] R. Courant, K. Friedrichs, and H. Lewy, "Über die partiellen Differenzengleichungen der mathematischen Physik," Mathematische Annalen, vol. 100, pp. 32–74, 1928. DOI: 10.1007/BF01448839.

[5] J. C. Maxwell, "A dynamical theory of the electromagnetic field," Philosophical Transactions of the Royal Society of London, vol. 155, pp. 459–512, 1865. DOI: 10.1098/rstl.1865.0008.

[6] A. Taflove and M. E. Brodwin, "Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell's equations," IEEE Transactions on Microwave Theory and Techniques, vol. 23, no. 8, pp. 623–630, 1975. DOI: 10.1109/TMTT.1975.1128640.

[7] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007. (General background on finite-difference schemes and stability.)

[8] wgpu-py developers, wgpu-py: WebGPU for Python, https://wgpu-py.readthedocs.io (accessed 2026-07-11); W3C, WebGPU Shading Language (WGSL) Specification, https://www.w3.org/TR/WGSL/.

No authoritative published reference was found for the arctangent-sigmoid colormap construction (customColormaps(), util.py:87) or for the two-disk lens rasterization; both appear to be original constructions of the author.