diff options
Diffstat (limited to 'docs/physics.md')
| -rw-r--r-- | docs/physics.md | 369 |
1 files changed, 0 insertions, 369 deletions
diff --git a/docs/physics.md b/docs/physics.md deleted file mode 100644 index d4599a8..0000000 --- a/docs/physics.md +++ /dev/null @@ -1,369 +0,0 @@ -# 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<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 panels and the export pipeline. |