diff options
Diffstat (limited to 'docs/physics.md')
| -rw-r--r-- | docs/physics.md | 373 |
1 files changed, 373 insertions, 0 deletions
diff --git a/docs/physics.md b/docs/physics.md new file mode 100644 index 0000000..4877f64 --- /dev/null +++ b/docs/physics.md @@ -0,0 +1,373 @@ +# 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 $\text{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 = \operatorname{clamp}\!\left(0.02\,\max(r - r_\text{photon},\,0),\; \Delta_\text{min},\; \Delta_\text{max}\right), +\qquad +\begin{aligned} +\Delta_\text{min} &= 10^6\\ +\Delta_\text{max} &= 2\times10^{10} +\end{aligned} +$$ + +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. |