Black Hole Simulation

A dependency-free C99 renderer that traces null geodesics around a black hole and returns the shadow, the lensed accretion disk, the photon ring and the relativistic beaming asymmetry at interactive frame rates. The same physics was later compressed into a single WebGL2 fragment shader, which is what renders behind the hero of this portfolio.

Brief

The problem

Every public black-hole image I could find was either a real observation (EHT), an offline render from a research pipeline, or a GLSL screensaver that had the shape roughly right and the physics entirely wrong. I wanted the second thing to look like the first thing, and I wanted to be able to explain every pixel.

The interesting constraint was that I refused to use a library. No Eigen, no OpenGL, no stb, no ray-tracing framework. That decision sounds like machismo, but it was pedagogical: the moment you import a vector library you stop knowing which quantities are conserved and which are just convenient. The whole point of the exercise was to be forced to know.

If the physics is a lie, the picture is just a screensaver. So the rule was: nothing in the image that does not come out of the metric.

That rules out the obvious shortcut, which is to draw a black disk, paste a bright ellipse on it, and call it lensing. The shadow, the ring, and the way the disk appears to stand up and wrap over the top are all consequences of one thing: light rays are not straight.

How I got there

Attempt one: bend the ray once. The textbook weak-field deflection for a light ray passing a mass at impact parameter b is about 4M/b. My first renderer marched the ray in straight segments and applied a single angular kick toward the hole at each step. It produced a dark circle with a smeared halo. There was no photon ring, no secondary image of the far side of the disk, and the shadow was the wrong size. The failure was structural, not numerical: a single-kick approximation is a first-order statement about a nearly-flat spacetime, and everything worth looking at happens in the strong field, where the deflection diverges logarithmically as the impact parameter approaches the critical value. Rays that should loop the hole once and come back out were instead either captured or flung away, because the approximation had no mechanism for the loop at all.

Attempt two: integrate the geodesic equation with the metric connection. More correct, worse in practice. Writing out the Christoffel symbols for Schwarzschild in standard coordinates gives an integrator that behaves fine far away and catastrophically badly near the horizon, where the coordinate speed of light goes to zero. Fixed step sizes either wasted most of the budget in the vacuum outside or overshot the photon sphere entirely. I was fighting the coordinates, not the physics.

Attempt three: use the conserved quantities. This is the version that shipped, in both implementations. For a static, spherically symmetric metric, a null geodesic has two constants of motion, energy E and angular momentum L, and the orbit equation collapses to a single vector acceleration on the spatial position,

a = -(3/2) · rs · L² · r^-5 · p

in units where G = c = 1, with rs the Schwarzschild radius. Everything else - the bending, the capture cone, the photon sphere at 1.5 rs - falls out of that one term. The WebGL version in this repository uses exactly this expression. What remains is an integration problem, which is where the actual engineering time went.

The integrator

Each pixel is one ray. Each ray is marched for a fixed budget of 125 steps along the acceleration above, using a semi-implicit update (velocity integrated first, then position, which keeps the orbit from inflating the way explicit Euler does).

The step size is adaptive and banded on radius, because the information content of the trajectory is not uniform:

| Region | Step length | Why | | --- | ---: | | Front vacuum, r > 2.4 | 0.040 - 0.280 | nearly straight, do not waste steps | | Disk envelope, 1.0 < r < 2.4 | 0.008 - 0.022 | this is where the image is made | | Near-horizon, r < 1.0 | 0.010 - 0.028 | capture decision is binary but must not tunnel | | Default / far side | 0.036 - 0.080 | cheap |

Two early-outs keep the budget honest. If the ray crosses the horizon, the pixel is black and the march stops. If it escapes past the outer bound, the march stops and the remaining direction is used to sample the background sky. Termination by radius, not by step count, is what makes 125 steps enough.

The payoff is the shadow size. The critical impact parameter for a Schwarzschild hole is b_crit = (sqrt(27)/2) · rs = 2.598 · rs, so with rs = 0.52 the silhouette radius is 1.352 in scene units, against an inner disk edge at 1.00. The renderer reports the ratio as 1.35 : 1, and it is not a tunable art parameter - it is a prediction of the integrator. Getting it right was the moment I trusted the march.

The disk

The accretion disk is a thin Keplerian structure between the innermost stable circular orbit and an outer edge. The inner edge is placed by the Bardeen-Press-Teukolsky ISCO condition rather than eyeballed; the angular velocity is Keplerian, omega proportional to r^-3/2, which is what makes the inner regions rotate visibly faster than the outer ones.

Temperature follows a Shakura-Sunyaev style radial profile, T(r) proportional to r^-3/4 times the standard (1 - sqrt(r_in/r))^1/4 stress-free inner boundary term, multiplied by a Gaussian in the vertical direction so the disk has a soft thickness rather than a hard plane. Colour comes from a five-stage Planck blackbody ramp, hottest at the inner edge.

The trouble with a smooth analytic disk is that it looks like soup. Real disks are turbulent. My first fix was to evaluate 3D value noise inside the fragment loop at several octaves, which halved the frame rate on a laptop GPU. The second fix was the one that stuck: bake the turbulence once into a 512 x 256 framebuffer texture at startup, in cylindrical coordinates (angle, radius), with shear, high-frequency filaments and a dust channel packed into the RGBA channels. The march then reads the disk with two texture fetches instead of forty noise evaluations. The texture is sampled with REPEAT wrapping in azimuth and CLAMP_TO_EDGE in radius, so the seam at angle zero is invisible.

Animation is a two-phase crossfade between the same texture read at two time offsets over a 36 second cycle, which gives long-lived structure that slowly reorganises instead of a boiling surface.

Relativistic optics

Three effects are applied per sample, and each one is cheap enough to survive inside the march.

Doppler beaming. The shift factor is delta = sqrt(1 - beta^2) / (1 - beta · cos theta), where theta is the angle between the emitting material's velocity and the outgoing ray. Observed intensity goes as delta^3. Orbital speed is clamped at beta = 0.72 rather than the true ISCO value, and the beaming factor is clamped to [0.42, 1.45] - unclamped, the receding limb goes to black and the image loses all structure on one side. Even clamped, the approaching limb is roughly 5x the receding limb. Beaming is applied before tone mapping, so the ACES curve compresses the highlight instead of clipping it.

Gravitational redshift. The factor g = sqrt(1 - 1.5 · rs/r) is applied to the colour temperature only, not to the intensity. This is a deliberate compromise and I want to be honest about it: shifting the temperature is what produces the visible reddening of the inner disk, and folding the full g^4 intensity penalty in as well makes the inner half of the image disappear into the shadow. The physically complete version is on the list of things to do properly.

Photon ring. The march tracks the closest approach radius of each ray. Rays whose minimum radius sits just outside the capture threshold are the ones that looped the hole, and they get a narrow bright contribution, an exp(-dMin^2 · 750) core over an exp(-36 · dMin) glow. That is a caustic approximation, not a proper ray-splitting or backward-tracing scheme, but it produces the thin ring inside the ejecta that the real images have and that naive renderers do not.

The background is a lensed procedural starfield - hashed cells with core plus halo, a slow twinkle, and a drift of 0.024 radians per second so the sky moves under the camera. It is sampled with the ray's final direction after the march, which is the entire point: the stars bend around the hole for free.

Results

Runtime behaviour in the WebGL version is guarded by a watchdog rather than a fixed quality setting: any frame slower than 20 ms decrements the resolution scale by 0.05 down to a floor of 0.65, and 60 consecutive frames faster than 16.8 ms raise it back by 0.05. Device pixel ratio is capped at 1.25 on desktop and 1.0 on touch devices. The canvas pauses when off-screen or when the tab is hidden, and a MutationObserver re-initialises the whole pipeline after a theme swap, because the drawing buffer is cleared.

What I'd change

The WebGL hero is not a Kerr renderer. It uses the Schwarzschild orbit equation and borrows Kerr-flavoured disk kinematics; frame dragging, the oblate horizon, and the asymmetric ISCO of a spinning hole are approximated, not solved. The C version is closer to honest but still treats the disk as a test fluid.

I would also separate the physics step from the shading step. Right now the beaming clamp, the redshift compromise, and the tone curve are entangled, which means changing one silently changes the look of the others. A proper spectral integral over the blackbody against the sensor response would replace both the five-stop palette and the redshift cheat, at a cost I have not measured.

And the caustic is a hack. The correct way to render a photon ring is to split rays near the critical impact parameter and accumulate the higher-order images, which is what the research pipelines do. Mine gets the visual gist in one extra exponential.

Credits

• Andrew Hamilton, University of Colorado - the Riemann-Solver/Cormac black hole explainer, and the beaming derivations I worked from
• Bardeen, Press and Teukolsky (1972), The kerr metric - ISCO and the rotating-hole reference values
• Shakura and Sunyaev (1973) - the accretion disk temperature profile
• Schneider, Ehler & Falcke - ray tracing in curved spacetime
• John Nash's A Beautiful Mind for the Monty Python bar scene, which is the intuition behind the beaming asymmetry
• The EHT collaboration's M87* and Sgr A* images, which set the target