Behind the Sim #1 — EM Waves: How We Put Maxwell's Equations in the Browser

Behind the Sim is a series about the story behind individual simulations: the physics choices, the dead ends, the design decisions, and the moments things suddenly started working. First up: EM Waves — a real-time Finite Difference Time Domain solver for Maxwell's equations, running at interactive frame rates in a plain browser tab.

The Starting Point: What Is an EM Wave?

James Clerk Maxwell unified electricity and magnetism in 1865 with four equations. The key insight for this simulation is the wave equation that falls out of them: a changing electric field induces a magnetic field, and vice versa. The two fields chase each other through space at the speed of light.

In 2D, we work with the transverse magnetic (TM) mode: the electric field has only a z-component Ez, and the magnetic field has x and y components Hx and Hy. This reduces Maxwell's curl equations to three update equations — simple enough to evaluate on every grid cell every frame.

FDTD: The Yee Lattice

The Finite Difference Time Domain (FDTD) method, introduced by Kane Yee in 1966, discretises Maxwell's equations on a staggered space-time grid. Electric and magnetic field components are offset by half a spatial and temporal step, which automatically satisfies Gauss's law (∇·B = 0) without any correction step.

// TM-mode FDTD update equations (2D, uniform grid)
// E and H are interleaved: H updates first, then E

// Update Hx and Hy at time step n+½
Hx[i][j] -= (Cz / dy) * (Ez[i][j+1] - Ez[i][j])
Hy[i][j] += (Cz / dx) * (Ez[i+1][j] - Ez[i][j])

// Update Ez at time step n+1
Ez[i][j] += Cz * (
    (Hy[i][j] - Hy[i-1][j]) / dx -
    (Hx[i][j] - Hx[i][j-1]) / dy
)

// Cz = dt / (ε₀ · μ₀) — Courant number

The Courant stability condition requires dt ≤ dx / (c · √2) for a 2D uniform grid. We set dt = 0.9 · dx / (c · √2) — 90 % of the stability limit for a small safety margin. With a 200×200 grid and dx = 0.01 m, the time step is about 21 ps per update.

The Build Timeline

Day 1 — Bare FDTD on a plain canvas

Started with a 100×100 JavaScript array, pure CPU, single point source. The wave propagated correctly — circular rings spreading outward. But it reflected off the edges perfectly, so the whole domain filled up with standing waves within seconds. Clearly needed absorbing boundaries.

Day 2 — Perfectly Matched Layers (PML)

PML is the standard FDTD absorbing boundary: a lossy sponge region around the domain where the fields are attenuated without reflection. Implementing it doubled the code complexity but eliminated reflections almost perfectly. The domain now felt "open" — waves exit cleanly.

Day 3 — Moving to WebGL

The CPU implementation ran at ~30 FPS on a 100×100 grid. Extending to 200×200 (the minimum grid for useful visualisation) dropped it to 8 FPS. The FDTD update equations are embarrassingly parallel — every cell can be updated independently. Moved the update loop to a WebGL fragment shader writing into a floating-point texture pair (ping-pong scheme). Jumped to 60 FPS at 400×400.

Day 4 — Rendering and sources

Mapped the Ez field to a blue-red colour scale (negative to positive), added a sinusoidal point source modulated by a Gaussian envelope to avoid the initial "click", and implemented a double-slit aperture for the classroom-famous diffraction pattern. Added dipole source, plane wave source, and a movable dielectric region the user can place on the canvas.

Day 5 — Vector field overlay

Overlaying Hx/Hy as arrows using instanced line segments in WebGL. Arrows are subsampled to a 20×20 grid, scaled by field magnitude, and rotated to field direction. The combination of Ez colour map and H-field arrows gives a visceral sense of how the two fields are in quadrature.

The Hardest Part: Floating-Point Textures

WebGL1 does not support reading from a floating-point texture that you've just written to — the "render to float texture" extension (OES_texture_float) exists but isn't universally writable. The solution is a ping-pong buffer: two framebuffers, swapping roles each frame. The shader reads from buffer A and writes to buffer B; next frame it reads from B and writes to A.

On iOS Safari the extension wasn't available at all, requiring a fallback to 16-bit half-float textures with reduced dynamic range. Field values are normalised before storage and rescaled on read — a small precision loss, but invisible at the colour resolution of the display.

The Dielectric Permittivity Texture

To let users place materials (glass, metal, etc.) in the domain, the FDTD update coefficients need to vary spatially. We bake the permittivity εr into a separate R8 texture and sample it during the update shader. Pressing a toolbar button sets a "brush mode" — subsequent mouse events paint permittivity values into the material texture via a small secondary framebuffer.

// In GLSL update shader — reading material texture
float eps_r = texture2D(u_material, vTexCoord).r * 10.0 + 1.0;
// Range: [1.0, 11.0] — air to glass-like dielectric

float Cz_local = dt / (eps_r * eps0 * mu0 * dx * dx);
// Same update formula, but Cz varies with material

Classroom Use: What Works and What Doesn't

We got feedback from several physics teachers after the simulation launched. The features that proved most valuable in the classroom:

What didn't work: the initial free-draw dielectric tool confused students because "drawing glass" looks just like "drawing air" before the wave hits. We added a persistent colour overlay for placed materials (amber tint = dielectric, grey = conductor) which helped enormously.

Try it yourself: Open the EM Waves simulation, hit the "Double Slit" preset, and drag the slit separation slider. You'll see the diffraction pattern change in real time — the same physics as Thomas Young's 1801 experiment, running in your browser.

Next Steps for the Simulation

The 2D TM-mode FDTD is a strong teaching tool, but there are natural extensions I want to build: