# the physics of donut everything donut draws comes from tracing light backward through the curved spacetime around a black hole. this document works through the physics and maths of that trace, in roughly the order the shader applies it, with references to the code in [`assets/shaders/geodesic.slang`](../assets/shaders/geodesic.slang) — symbol names below (`InitRay`, `GeodesicRHS`, `DiskEmission`) all live in that file. for the software side — how the shader gets fed, the render backends, the docking ui and the export path — see [`architecture.md`](architecture.md). ## contents - [overview](#overview) - [units and scale](#units-and-scale) - [the schwarzschild metric](#the-schwarzschild-metric) - [null geodesics and conserved quantities](#null-geodesics-and-conserved-quantities) - [the equations of motion](#the-equations-of-motion) - [numerical integration](#numerical-integration) - [the three critical radii](#the-three-critical-radii) - [the accretion disk](#the-accretion-disk) - [redshift, doppler beaming and colour](#redshift-doppler-beaming-and-colour) - [the impact parameter](#the-impact-parameter) - [observable channels](#observable-channels) ## overview a black hole isn't drawn like ordinary geometry. for each pixel donut casts a ray from the camera and follows it *backward* until one of four things happens: it crosses the event horizon, it strikes the accretion disk, it hits a placed object, or it escapes to the background sky. mass bends the path of light, so the rays curve, and that one effect produces the whole picture: the dark shadow, the bright ring wrapped around it, the far side of the disk folded up over the top of the hole, and the doppler-brightened leading edge. ```mermaid flowchart LR A[Camera pixel] --> B[Build ray direction] B --> C{March the geodesic
through curved spacetime} C -->|falls in| D[Event horizon
black shadow] C -->|hits disk| E[Accretion disk
redshifted blackbody] C -->|hits object| F[Placed sphere
shaded] C -->|escapes| G[Background sky
HDRI environment] D --> H[Pixel colour] E --> H F --> H G --> H ``` ## units and scale donut works in geometric units, $G = c = 1$. mass then carries units of length, and the schwarzschild radius reduces to $$ r_s = \frac{2GM}{c^2} = 2M, \qquad\text{so}\qquad M = \frac{r_s}{2}. $$ one number describes the hole. for sagittarius a* the shader fixes it as ``` static const float SagA_rs = 1.269e10; // metres (M ≈ 4.3×10⁶ M☉) ``` every distance the integrator handles is a physical length in metres, written as a multiple of `SagA_rs`, so the critical radii come out as constants: ``` R_PHOTON = 1.5 * SagA_rs // photon sphere (3M) R_ISCO = 3.0 * SagA_rs // ISCO (6M) ``` the scene editor uses a friendlier grid. the constant `SCENE_UNITS_PER_RS = 3.0` (in [`src/scene/scene_types.h`](../src/scene/scene_types.h)) sets three grid units to one schwarzschild radius. when the renderer hands a placed object to the shader it scales the position by $r_s/3$ (the code's `SagA_rs / 3`) to get metres, so the editor and the simulation always agree on where things sit. ## the schwarzschild metric sgr a* is treated as a non-rotating, uncharged black hole, whose spacetime is the exact schwarzschild solution of einstein's equations. in spherical coordinates $(t, r, \theta, \phi)$ the line element is $$ ds^2 = -\left(1-\frac{r_s}{r}\right)dt^2 + \left(1-\frac{r_s}{r}\right)^{-1}dr^2 + r^2\left(d\theta^2 + \sin^2\theta\, d\phi^2\right). $$ the factor that keeps recurring is abbreviated $$f(r) = 1 - \frac{r_s}{r},$$ which is `float f = 1.0 - SagA_rs / r;` in the code. as $r \to r_s$, $f \to 0$ and the metric coefficients diverge. that divergence is a coordinate artifact rather than a real singularity, but it is why the integrator stops a ray once it reaches $r \le r_s$ instead of pushing through. ## null geodesics and conserved quantities light follows null geodesics, the curves with $ds^2 = 0$. the metric has no explicit dependence on $t$ or $\phi$ (a time-translation symmetry and an axial rotation symmetry), so two quantities stay constant along every ray: $$ E = f(r)\,\frac{dt}{d\lambda} \qquad\text{(energy)}, \qquad\qquad L_z = r^2\sin^2\theta\,\frac{d\phi}{d\lambda} \qquad\text{(axial angular momentum)}, $$ with $\lambda$ an affine parameter along the ray. `InitRay` sets both at the camera: it turns the camera-space ray direction into the spherical components $(\dot r, \dot\theta, \dot\phi)$, then computes ``` ray.L = r*r * sin(theta) * dphi; // angular momentum dt_dL = sqrt(dr*dr/f + r*r*(dtheta² + sin²θ·dphi²)); ray.E = f * dt_dL; // energy ``` `E` is put to work during integration: `GeodesicRHS` reads the time component $\dot t = E/f$ from it, so the $t$ coordinate never has to be integrated on its own — one fewer equation per step. `L` is computed at initialisation as the ray's angular momentum but isn't fed back into the equations of motion; the azimuthal motion is carried directly by $\dot\phi$. for a null geodesic the affine parameter has an arbitrary overall scale, and the ray's shape — which is all the image depends on — doesn't change with it, so the exact normalisation of `E` is only a convention. ## the equations of motion marching a ray means solving the geodesic equation $\ddot x^\mu + \Gamma^\mu_{\alpha\beta}\dot x^\alpha \dot x^\beta = 0$ for the schwarzschild metric. `GeodesicRHS` writes it as a first-order system in the six ray variables $(r,\theta,\phi,\dot r,\dot\theta,\dot\phi)$. the three positions advance by their velocities, $$\dot r,\qquad \dot\theta,\qquad \dot\phi,$$ and the three velocities accelerate with the curvature: $$ \ddot r = -\frac{r_s}{2r^2}\,f\,\dot t^2 + \frac{r_s}{2r^2 f}\,\dot r^2 + r\left(\dot\theta^2 + \sin^2\theta\,\dot\phi^2\right), \qquad \dot t = \frac{E}{f}, $$ $$ \ddot\theta = -\frac{2}{r}\,\dot r\,\dot\theta + \sin\theta\cos\theta\,\dot\phi^2, $$ $$ \ddot\phi = -\frac{2}{r}\,\dot r\,\dot\phi - 2\cot\theta\,\dot\theta\,\dot\phi. $$ the terms are the christoffel symbols of the metric. in $\ddot r$ the first term is the inward pull of gravity (it carries $\dot t^2$, hence the energy); the rest are the centrifugal contributions from angular motion. the $\theta$ and $\phi$ equations are the angular-momentum couplings that hold the ray to its orbital plane and sweep it around the hole. the code is a direct transcription: ``` d2.x = -(SagA_rs/(2r²))·f·dt_dL² + (SagA_rs/(2r²f))·dr² + r·(dtheta² + sin²θ·dphi²); d2.y = -2·dr·dtheta/r + sin(theta)·cos(theta)·dphi²; d2.z = -2·dr·dphi/r - 2·(cos/sin)(theta)·dtheta·dphi; ``` ## numerical integration there is no closed form for a general ray, so the integrator advances it in steps. `RK4Step` takes one step. it evaluates `GeodesicRHS` once, advances the six variables by `dL` times their rates, and recomputes the cartesian position from the new spherical coordinates. that is a single forward-euler stage, despite the name: only the first slope `k1` is evaluated, where a genuine fourth-order step would also compute `k2`, `k3` and `k4` at intermediate points. moving to real rk4 is the obvious accuracy upgrade; as it stands, almost all of the accuracy comes from the step-size control instead. `CalculateAdaptiveStepSize` chooses the step length. a fixed step would waste time far from the hole and lose the trajectory near it, so the step scales with distance from the photon sphere: $$ \Delta\lambda = \mathrm{clamp}\!\left(0.02\,\max(r - r_\text{photon},\,0),\; \Delta_\text{min},\; \Delta_\text{max}\right), $$ with $\Delta_\text{min} = 10^6$ and $\Delta_\text{max} = 2\times10^{10}$ metres. far out, the ray is in near-flat space and crosses it in a handful of long strides. near the photon sphere, where the path bends hardest and mistakes show the most, the step shrinks to follow the curve. a second clamp forces the step down to the disk's half-thickness whenever the ray is near the disk plane, so a thin, nearly edge-on disk is never stepped straight over. a ray's march ends on the first of these: | condition | meaning | | --- | --- | | $r \le r_s$ (`Intercept`) | fell through the horizon → shadow (black) | | crossed / entered the disk slab | hit the opaque disk → emit its colour | | `InterceptObject` (every 5 steps) | hit a placed sphere → shade it | | $r > $ `earlyExitDistance` ($2\times10^{12}$) | left the rendered region → sample the sky | | $\dot r > 0$ and $r > 50\,r_s$ | outbound in flat space, direction frozen → sample the sky early | | step count exceeds the budget | up to `quality_steps` (default 15000, clamped 1000–15000) | the step budget is the same whether the camera is moving or settled. at steep, strongly-lensed poses the disk only resolves with a high step count, so cutting it during motion would make the disk flicker. responsiveness during a drag comes from the rendering resolution and sample count instead — a smaller target and one sample per pixel while moving, sharpening to full resolution and 4× supersampling once the camera settles (see [`architecture.md`](architecture.md#progressive-resolution-and-supersampling)). ## the three critical radii three radii set up everything you see: ```mermaid flowchart LR subgraph one[" "] direction LR H["Event horizon
r = rₛ = 2M"] --- P["Photon sphere
r = 1.5 rₛ = 3M"] --- I["ISCO
r = 3 rₛ = 6M"] end ``` the **event horizon** at $r_s = 2M$ is the point of no return; the set of directions whose rays end there is the black shadow. the **photon sphere** at $\tfrac{3}{2}r_s = 3M$ is where light can circle the hole on unstable orbits, so rays passing near it loop around once or more before escaping — this makes the thin photon ring against the shadow and the folded multiple images of the disk. the **isco** at $3r_s = 6M$ is the innermost stable circular orbit, inside which matter can't hold a steady orbit; it is the disk's inner edge, and `DiskEmission` clamps the inner radius with `max(disk.disk_r1, R_ISCO)`. ## the accretion disk donut models the disk as a thin, opaque, self-luminous slab in the equatorial plane ($y = 0$), not a volumetric cloud. a ray hits it the first time it crosses the midplane (or grazes into the slab of half-thickness `disk.thickness`) inside the radial band $[r_\text{in}, r_\text{out}]$, and that surface's emission is the pixel colour. the default band runs from $3\,r_s$ to $12\,r_s$. a steady thin accretion disk radiates with a flux that rises from zero at the inner edge, peaks just outside it, and tails off with radius: $$ F(r) \;\propto\; \frac{1}{r^3}\left(1 - \sqrt{\frac{r_\text{in}}{r}}\right). $$ this is the novikov–thorne / shakura–sunyaev thin-disk profile. its peak sits at $r/r_\text{in} \approx 1.36$ with value `FLUX_PEAK = 0.0569`, which normalises it. a blackbody's flux goes as $T^4$ (stefan–boltzmann), so the local temperature is $$ T(r) = T_\text{peak}\left(\frac{F(r)}{F_\text{peak}}\right)^{1/4}, $$ where $T_\text{peak}$ is the tunable `disk.temperature`, 4800 k by default: ``` flux = max((1 - sqrt(1/xr)) / (xr*xr*xr), 0); xr = rc / r_in Tn = pow(flux / FLUX_PEAK, 0.25); // normalised temperature, peak ≈ 1 Temit = disk.temperature * Tn; ``` an optional turbulence overlay (the `turbulence` parameter, `disk.disk_num`) modulates the brightness with animated fractal noise to suggest churning gas. it never changes the fact that the disk is an opaque surface. ## redshift, doppler beaming and colour the disk is hot gas on relativistic orbits, deep in the gravity well. two effects shift its light on the way to the camera, and both collapse into a single redshift factor $g$ (observed frequency over emitted). the gas moves on prograde circular geodesics. for schwarzschild, the locally measured orbital speed is $$ v = \sqrt{\frac{M}{r - 2M}} = \sqrt{\frac{r_s/2}{r - r_s}}, $$ which is exactly $0.5\,c$ at the isco. the velocity vector is $\boldsymbol\beta = v\,\hat\phi$, tangent to the orbit. combining the gravitational and time-dilation shift of a circular orbit with the relativistic doppler shift from that motion gives $$ g = \frac{\sqrt{\,1 - \tfrac{3}{2}\,\dfrac{r_s}{r_c}\,}}{1 - \boldsymbol\beta\cdot\hat n}, $$ where $\hat n$ points along the photon toward the observer and $r_c$ is the cylindrical radius of the emission point. the numerator is the gravitational part (it vanishes at the photon sphere $r_c = \tfrac{3}{2}r_s$, where even orbiting light is infinitely redshifted); the denominator is the doppler part, which brightens and blueshifts the side turning toward the camera and dims and redshifts the receding side. as a check, $g \to \sqrt{1/2}$ at the isco, matching the code. two things follow from $g$, both physical: $$ T_\text{obs} = g\,T_\text{emit} \qquad\text{(colour: a redshifted blackbody)}, \qquad\qquad I_\text{obs} = g^4\,I_\text{emit} \qquad\text{(relativistic beaming)}. $$ the colour is the planckian blackbody colour at the observed temperature, `Blackbody(g · Temit)`, using a tanner-helland fit to the planckian locus. the brightness keeps the physical $g^4$ beaming — the real approaching/receding asymmetry — while the enormous $T^4$ radial range is compressed to $T_n^2$ for display, so the colour gradient across the disk stays visible instead of collapsing to a single saturated ring: ``` bright = pow(Tn, 2.0) * pow(g, 4.0) * edge; // edge = soft inner/outer falloff colour = Blackbody(g * Temit) * bright; ``` that $T_n^2$ in place of the physical $T_n^4 = F$ is the one intentional concession to legibility; the rest of the disk model is the genuine relativistic result. ## the impact parameter a ray's impact parameter $b$ is the perpendicular distance from the hole's centre to the straight line the ray would have followed with no gravity — the quantity that sets how strongly it deflects. donut reads it straight off the camera geometry (in units of $r_s$): $$ b = \frac{\lVert \mathbf{r}_\text{cam} \times \hat d\,\rVert}{r_s}, $$ with $\mathbf{r}_\text{cam}$ the camera position relative to the hole and $\hat d$ the pixel's ray direction. rays whose $b$ is near the critical value (about $\tfrac{3\sqrt3}{2}r_s$) are the ones that skim the photon sphere and build the ring. ## observable channels the renderer already computes these physical quantities while tracing, so it can output them directly instead of only the final colour. the shader's `outputChannel` picks which quantity each pixel reports, and `rawOutput` picks whether to write the raw floating-point value (for analysis) or a false-coloured / tone-mapped version (for viewing): | channel | quantity | notes | | --- | --- | --- | | 0 | colour | the final tone-mapped hdr radiance — the normal image | | 1 | redshift $g$ | disk pixels only; validity flagged in alpha | | 2 | emission temperature $T_\text{emit}$ (k) | disk pixels only; validity in alpha | | 3 | impact parameter $b$ ($r_s$) | a per-ray geometric quantity, defined everywhere | the raw channels are what make the export usable as data rather than just imagery; [`architecture.md`](architecture.md#the-export-pipeline) covers how they are rendered off-screen and written to pfm or csv. --- see also [`architecture.md`](architecture.md) for how the renderer is built, from the portable gpu layer up through the panels and the export pipeline.