path: root/docs/physics.md

# 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 tabs
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<br/>through curved spacetime}
    C -->|falls in| D[Event horizon<br/>black shadow]
    C -->|hits disk| E[Accretion disk<br/>redshifted blackbody]
    C -->|hits object| F[Placed sphere<br/>shaded]
    C -->|escapes| G[Background sky<br/>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<br/>r = rₛ = 2M"] --- P["Photon sphere<br/>r = 1.5 rₛ = 3M"] --- I["ISCO<br/>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 tabs and the export pipeline.