Geodesics API

Geodesics

Null-geodesic engine: metric type tree with AD-derived metric derivatives (metrics.jl), plane-parallel camera (camera.jl), the Hamiltonian covariant-momentum integrator (integrator.jl), general analytic and numerical spacetimes, and the exact-Kerr analytic tracer + lensing bands (AART/Krang).

Kerr in Cartesian Kerr-Schild coordinates (blacklight-compatible default).

Flat spacetime.

Any user closure x::SVector{4} -> SMatrix{4,4} (covariant); AD supplies ∂g. Traced through the general AnalyticSpacetime geodesic_rhs/trace_geodesic path, which retains the ∂_t g term — a TIME-DEPENDENT closure integrates correctly. The engine knows no horizon for an arbitrary metric, so trace_geodesic stops only on the escape sphere or λmax.

SuperposedKerrSchild(masses, spins, centers)

Analytic superposed Kerr-Schild binary: two (or one) static Kerr-Schild black holes with spins along ±z,

g_μν = η_μν + Σ_A f_A l^A_μ l^A_ν      (covariant l^A_0 = 1),

exact for one BH, the standard large-separation binary initial-data ansatz for two. The single-KS inverse shortcut is invalid for the sum (cross term fA fB l^A·l^B), so inverse_metric uses the full 4×4 inverse (StaticArrays) and inverse_metric_gradient differentiates THAT with ForwardDiff — reducing exactly to single Kerr when a second mass is zero.

ADMGridSpacetime(xs, ys, zs, data; alpha_min=0.1, offgrid_eta=true)

Frozen-slice numerical spacetime on a REGULAR Cartesian grid. data is (10, nx, ny, nz) holding the ten ADM fields per node in the column order (α, β^x, β^y, β^z, γ_xx, γ_xy, γ_xz, γ_yy, γ_yz, γ_zz); xs,ys,zs are the uniform per-axis node ranges. metric samples the ten fields with a C¹ cubic (Catmull-Rom) interpolant and assembles g_μν via the ADM formula: g_00 = -α² + β_i β^i, g_0i = β_i = γ_ij β^j, g_ij = γ_ij. Off-grid returns Minkowski (offgrid_eta); a ray terminates when the interpolated lapse falls below alpha_min or on escape. Build one from any analytic metric with adm_grid_from_metric.

adm_grid_from_metric(metric_fn, xs, ys, zs; alpha_min=0.1, offgrid_eta=true,
                     excise=nothing)

Sample any analytic covariant metric closure x::SVector{4} -> SMatrix{4,4} onto a regular Cartesian grid, storing the extracted ADM data — the "Kerr-as-grid" factory used to validate the numerical path against the analytic tracer. metric_fn is evaluated on the t=0 slice. excise (optional x::SVector{4} -> Bool) marks nodes to fill with the collapsed-lapse excision data instead of sampling metric_fn — pass the horizon predicate, e.g. x -> ks_radius(st, x) < horizon_radius(st), to give the tracer a grid-alignment-independent inner stopping surface (without it the only excised nodes are those where metric_fn happens to be non-finite — for Kerr the measure-zero z=0 disk nodes, which a shifted/even-count grid misses).

sample_adm(st, x, y, z) -> (α, β::SVector{3}, γ::SMatrix{3,3})

Cubically interpolate the ten ADM fields at Cartesian (x, y, z) (indices may be ForwardDiff.Dual). β is the contravariant shift β^i; γ the spatial metric.

Duck-typed on the grid fields (x0…dz, nx…nz, data) rather than ::ADMGridSpacetime, so the GPU twin param struct (gpu_spacetime.jl DevADMParams, same field names, data::CuDeviceArray) runs this SAME body in-kernel — the CPU interpolant IS the device interpolant.

ADM 4-metric gμν from the cubically interpolated (α, β^i, γij).

Exact Kerr-Schild inverse: g^{ab} = η^{ab} − f l^a l^b with l^a = η^{ab} l_b.

Contravariant metric by direct 4×4 inverse; the generic AD inverse_metric_gradient then supplies ∂g^{μν}.

inverse_metric_gradient(st, x) -> SArray{Tuple{4,4,4}}

dg[a,b,c] = ∂ g^{ab} / ∂ x^c, by forward-mode AD of inverse_metric (compile-time specialized; exact through the KS closed form).

ks_radius(st::KerrSchildCartesian, x::SVector{4}) -> Float64

Kerr–Schild radius $r$ at the Cartesian Kerr–Schild position x.

$r$ is the spheroidal radial coordinate defined implicitly by $\frac{x^2 + y^2}{r^2 + a^2} + \frac{z^2}{r^2} = 1$, not the Euclidean distance $\sqrt{x^2+y^2+z^2}$; the two coincide only as $r \to \infty$ or for $a = 0$. Geometrized units, $M = 1$.

See also horizon_radius, ks_radius_hole.

horizon_radius(st::KerrSchildCartesian) -> Float64

Outer event-horizon radius $r_+ = 1 + \sqrt{1 - a^2}$ in geometrized units ($M = 1$), as a Kerr–Schild radius.

Rays are terminated just outside this radius during integration. Compare against ks_radius, which is the coordinate this is expressed in.

See also horizon_radius_hole for the binary case.

Local KS radius of hole A (for termination / diagnostics).

Per-hole horizon radius MA + sqrt(MA^2 - a_A^2).

christoffel(st, x)           -> SArray{Tuple{4,4,4}}
christoffel(st, x, gcon, dg)  -> SArray{Tuple{4,4,4}}

Γ^μ{αβ} = ½ g^{μλ}(∂α g{λβ} + ∂β g{λα} − ∂λ g{αβ}), derived via AD; index order Γ[μ, α, β]. The four-argument form contracts a PRECOMPUTED inverse metric gcon = g^{μλ} and covariant gradient `dg[a,b,c] = ∂c g{ab}(metricgradient(st, x)) so hot-loop callers that already holdg^{μλ}skip a second closed-form metric evaluation; the two-argument form derives both and delegates. Allocation-free when inlined (@SArray` comprehensions).

metric_gradient(st, x) -> SArray{Tuple{4,4,4}}

dg[a,b,c] = ∂ g_{ab} / ∂ x^c (covariant), by forward-mode AD of metric.

Outer event-horizon radius r₊ = 1 + √(1 − a²) (M = 1, Boyer–Lindquist).

Kerr Δ(r) = r² − 2r + a² (M = 1, Boyer–Lindquist).

kerr_radial_potential(r, a, λ, η) -> R(r)

The Kerr photon radial potential R(r) = (r²+a²−aλ)² − Δ(r)(η+(λ−a)²) (Gralla & Lupsasca PRD 101, P2 Eq. 5); its reciprocal square root ∫dr/√R is the Mino-time integrand. M = 1, Boyer–Lindquist; AD-generic in every argument.

Hamiltonian geodesic RHS for stationary metrics (the exact-zero dp_t/dλ is dropped).

General geodesic RHS for an arbitrary (possibly TIME-DEPENDENT) analytic metric — the metric-agnostic path serving the exported UserMetric. Identical to the stationary form but RETAINS the time term dp_t/dλ = -½ (∂_t g^{μν}) p_μ p_ν. For a stationary closure ∂_t g^{μν} ≡ 0, so this reduces exactly to the stationary RHS (which the more specific StationarySpacetime method drops BY CONSTRUCTION).

trace_geodesic(st, y0; kwargs...) -> ODE solution

Production adaptive solve with event-terminated horizon capture / escape. save_everystep=true retains the trajectory (needed by the radiation stage).

trace_geodesic(st::AnalyticSpacetime, y0; escape, λmax, ...) -> ODE solution

Metric-agnostic driver for an arbitrary analytic metric (the exported UserMetric). Uses the general geodesic_rhs above, so a time-dependent closure integrates correctly. The engine knows no horizon for an arbitrary metric, so integration terminates only on the escape sphere (OUTBOUND crossing) or at λmax — set both so an interior plunge cannot run away.

trace_geodesic(st::ADMGridSpacetime, y0; escape, ...) -> ODE solution

Adaptive geodesic solve terminating when the interpolated lapse α < alpha_min (numerical horizon), on escape, or off-grid. escape should be ≲ the grid's inscribed half-width so a ray leaves through the escape sphere while still on-grid.

trace_rays(st, rays; escape, reltol, abstol) -> (states, dls)

Batch tracer: every camera ray solved with the adaptive trace_geodesic, THREADED across rays, returning per-ray sample trajectories ORDERED SOURCE → CAMERA (reversed from integration order), ready for the radiation stage: states[i] :: Vector{SVector{8}} — samples along ray i dls[i] :: Vector{Float64} — affine step to the NEXT sample (same length; last entry 0) Works for ANY AbstractSpacetime: the per-spacetime trace_geodesic methods own the termination policy (Kerr horizon, binary per-hole capture, ADM lapse).

trace_rays_sampled(st, rays; C=40.0, escape=1.1e3, nmax=6000, R_em=Inf,
                   max_step=Inf, cell_size_fn=nothing, cell_frac=1.0)

Fixed-control RK4 tracer with step h = (r - rh)/C (the blacklight/arXiv:2304.03804 step rule). Unlike the error-adaptive `tracerays, the sample spacing here is also the RADIATION sampling grid — dense where the fluid varies (near the hole), which the error-adaptive tracer under-samples in smooth-geodesic regions. Output contract identical totracerays(source→camera order).Rem(defaultInf, disabled) adds a far-field early-out — once a ray is outbound pastRemit stops sampling the vacuum return leg; only pass it when the emission is confined withinRem`.

OPT-IN AMR-cell step cap (as in RAPTOR — never stride over a fluid cell, which matters for AMR sampling fidelity and slow light): max_step hard-caps h, and cell_size_fn(x::SVector{4})::Real supplies the local fluid cell size, capping h by cell_frac * cell_size_fn(x). Both default to disabled, leaving the sampling BIT-IDENTICAL — the step is only ever reduced, and min(h, Inf) === h.

NOTE (step-density accuracy, measured at a=0.9375; gpu_tracer.jl uses the same rule): C=80 meets image-level budgets (NMSE(I) 1.5e-4 vs C=160, integrated flux 3.8e-4), but the null-constraint drift max|k·k| on plunging rays (r<3) reaches ~1e-7 (median 1.1e-10; the excess is confined to ~2.4% of captured rays' tails), and per-pixel transfer convergence is first-order in 1/C. Use C≥160 for convergence studies (≈2× march cost; max|k·k| ~6e-9).

trace_rays_sampled(st::SuperposedKerrSchild, rays; C=40.0, escape=ESCAPE_RADIUS,
                   nmax=6000, R_em=Inf, max_step=Inf, cell_size_fn=nothing,
                   cell_frac=1.0) -> (states, dls)

Binary sibling of the Kerr fixed-control sampled tracer (integrator.jl): the same step rule and source→camera output contract, generalized to N punctures. Two changes vs the Kerr method — (i) the step uses the distance to the NEAREST horizon, h = min_A(r_A - r_h_A)/C, so sampling densifies toward whichever hole a ray approaches; (ii) capture (at horizon + HORIZON_EPS_SAMPLED) terminates on EITHER hole, matching trace_geodesic(::SuperposedKerrSchild). Escape uses the global coordinate radius. R_em and the opt-in max_step/cell_size_fn AMR caps behave exactly as in the Kerr method, and its step-density accuracy NOTE applies unchanged.

One classical fixed-step RK4 step (AD/GPU mode: static computation graph).

nullify_covariant(st, x, p_spatial) -> SVector{4}

Solve g^{μν} pμ pν = 0 for pt given spatial covariant momentum, taking the past-oriented root (pt < 0 ⇒ E > 0). Outside the ergosphere g^{tt} < 0 and C = g^{ij}pi pj > 0, so the two roots have OPPOSITE signs and exactly one is negative. Inside the ergosphere (f > 1 ⇒ C can be < 0) both roots may be positive: no E > 0 root exists and the contract is unsatisfiable — this errors rather than silently returning a negative-energy (E < 0) covector that a downstream p/(−p_t) normalization would flip.

Conserved energy E = -p_t (exact for stationary metrics in this formulation).

Conserved Lz for CKS axisymmetry (Killing vector ∂φ = x∂y − y∂x): x py − y px.

Null constraint g^{μν} pμ pν (≡ 0 on exact null geodesics; error diagnostic).

Plane-parallel camera at inclination inc_deg [deg], distance d, square fov.

camera_rays(st, cam) -> Vector{SVector{8,Float64}}

Initial states [x^μ; p_μ] for every pixel (row-major over (α, β) pixel centres), null and normalized to p_t = -1. The spatial momentum points INWARD, so forward-λ integration traces camera → strong field (the standard backward- ray-tracing convention). Works for ANY AbstractSpacetime: only metric and nullify_covariant at the camera location are required.

Rays at arbitrary image-plane coordinates [(α, β), ...] (AMR refinement path).

isco_radius(a; prograde=true) -> r_ISCO

Kerr innermost stable circular orbit (Bardeen, Press & Teukolsky 1972), M=1. Returns 6 at a=0, 4.233 at a=0.5, 2.024 at a=0.94 (prograde).

isco_prograde(a) -> r_ISCO

Prograde ISCO radius (M = 1) — AART rms, P1 Eq. B16. Thin alias for isco_radius(a; prograde=true), kept as a public name for the Radiation g-disk source; AD-generic in a.

Prograde photon (marginally bound null) circular orbit radius, M=1.

Marginally bound (E=1) circular-orbit radius (M=1): 2 - a + 2 sqrt(1-a) prograde.

keplerian_omega(a, r; prograde=true) -> Omega_K

Coordinate angular velocity of a circular equatorial (Keplerian) orbit at Boyer–Lindquist radius r in Kerr (M=1): Omega_K = ±1 / (r^{3/2} ± a).

keplerian_specific_angular_momentum(a, r; prograde=true) -> l_K

Specific angular momentum l = -uphi/ut of a circular equatorial geodesic (Bardeen, Press & Teukolsky 1972, sign ±), M=1: lK = ±(r^2 ∓ 2 a sqrt(r) + a^2) / (r^{3/2} - 2 sqrt(r) ± a) — upper sign prograde, lower sign retrograde. At a=0 this is ±r^{3/2}/(r-2) (Schwarzschild). The prograde branch is > 0 with its minimum lms at the ISCO; the retrograde branch carries BPT's overall LEADING minus, so it is NEGATIVE (a counter-rotating orbit has uphi > 0, hence l = -uphi/u_t < 0).

conserved_quantities(α, β, θo, a) -> (λ, η)

Energy-reduced axial angular momentum λ and Carter constant η for a photon landing at Bardeen screen coordinate (α, β) for an observer at inclination θo (rad). P2 Eqs. (58–59). α here is the Bardeen horizontal coordinate; the render seam maps the plane-parallel camera's image-α to this via _cam_to_bardeen.

radial_roots(λ, η, a) -> (r1, r2, r3, r4)::NTuple{4,Complex}

The four roots of the radial potential R(r) = r⁴ + (a²−η−λ²)r² + 2(η+(a−λ)²)r − a²η (P2 Eqs. A1–A8), closed form and complex-safe. Ordered r1 ≤ r2 ≤ r3 ≤ r4 when all real (case II / Type II, four real roots); when only two are real (case III / Type III) r1,r2 are the real pair and r3,r4 a complex-conjugate pair.

angular_turning_points(λ, η, a) -> (u_p, u_m)

Roots u± = cos²θ± of the angular potential (P2 Eq. 11): for equatorial-crossing (non-vortical) rays η ≥ 0 ⇒ up ≥ 0, um ≤ 0. Real for η ≥ 0.

analytic_angular_integrals(n, β, u_p, u_m, θo, a) -> (Gθ, Gφ, Gt)

Angular path integrals for EQUATORIAL emission (θs = π/2 ⇒ source terms vanish), P2 Eqs. (12–22). n is the photon-ring layer (mbar); the effective turning-point count is m = n + H(β) (Heaviside — the disk's far side for β ≥ 0). Gθ is the Mino time τ the radial inversion is evaluated at.

analytic_radial_integrals(rs, ro, roots, a, λ, η, is_case2, νr) -> (I_r, I_φ, I_t)

Definite radial path integrals from source rs to observer ro (P3 Eqs. B1–3). For Type II the ray turns at r4, so the integral is split with the redshift-sign weighting w = (νr == −1); Type III is a plain observer−source difference.

analytic_emission(α, β, θo, a, n; ro=1000.0, tol=1e-12) -> (rs, φs, ts, νr, valid)

Exact-Kerr direct-inversion equatorial ray trace: for a photon landing at the Bardeen screen coordinate (α, β) (observer inclination θo [rad], spin a, M = 1), return the equatorial source point of photon-ring layer n (0 = direct, 1, 2, …) — BL radius rs, source azimuth φs, coordinate time-of-flight ts, the radial-velocity sign νr at emission (∈ {−1, +1}; the redshift sign), and a valid flag (false if the ray is vortical, lands inside the outer horizon, or needs more Mino time than the ray physically has — case II: observer → turn → escape, case III: observer → horizon plunge — so aliased sn/cn ghost landings are rejected and each layer n is confined to its lensing band). ForwardDiff-generic in (α, β, θo, a).

DOMAIN. Exact Kerr symmetries normalise the inputs: θo > π/2 mirrors through the equator (θo → π−θo, β → −β), a < 0 through the azimuth (α → −α; φs, λ mapped back), and |a| < 1e-6 is clamped to a = 1e-6 (the GL angular integrals are 0/0 at a = 0 exactly). α is the Bardeen screen coordinate; the render seam maps the camera image coordinate to it via _cam_to_bardeen. The degenerate β = 0 ray uses the Heaviside convention H(0) = 1 (grazes the disk's far side), shared with lensing_bands.jl.

equatorial_intensity(rs, νr, a; rin=0.0, source=nothing, λ=0, η=0) -> Float64

The equatorial-source SEAM. Given the landing radius rs, the emitted-p^r sign νr and (for source) the ray's conserved quantities (λ, η), return the observed specific intensity, zero off the disk (inside the horizon or rs < rin):

• source === nothing (default) — a placeholder: an optically thin r⁻² surface-brightness falloff (kept so the bare tracer stays self-contained).
• source a callable (rs, νr, a, λ, η) -> I — the injected AART g-disk / g-gas redshift source (Radiation.equatorial_disk_intensity, wired at the render/root scope where Radiation is visible; Geodesics stays dependency-free of Radiation). This is how render_analytic(cfg; source=EquatorialDisk(a)) produces the layered photon-ring images.

The (rs, νr, a; …) positional call shape is preserved for the render seam.

critical_curve(a, θo; n=100) -> (α, β)

The Kerr critical curve (shadow boundary) on the Bardeen screen for spin a, observer inclination θo [rad] (P2 Eqs. 38–39, 55), as a closed parametric loop of ≈2n points (upper branch + mirrored lower branch). For a → 0 it is the Schwarzschild circle of radius 3√3 M, returned analytically (machine precision).

apparent_horizon(a, θo; n=100, ro=1e5) -> (α, β)

The apparent-horizon curve — the n=0 lensing band's INNER boundary (image of the event horizon). For each critical-curve direction the screen point is scaled by s ∈ (0,1) until a ray there originates exactly at the horizon (AART ApparentBH root). For a → 0 the curve is a circle (Schwarzschild); with spin it is D-shaped.

lensing_band(a, θo, n; npts=100, box=25.0) -> (inner, outer)

Inner and outer boundary curves (each an N×2 matrix of (α,β)) of lensing band n for spin a, inclination θo [rad] (θo > π/2 maps through the equatorial mirror β → −β, as in analytic_emission). Band n=0 is bounded by the apparent horizon (inner) and the image box of half-width box (outer). Bands n ≥ 1 are thin annuli straddling the critical curve: for each critical-curve direction the inner edge scales inward (s<1, Type III root of AART nlayers) and the outer edge outward (s>1, Type II root), both hand-rolled bisections.

band_grid(inner, outer; dx=0.02, box=nothing) -> (grid, mask, N)

Adaptive Cartesian (α,β) grid covering the outer boundary's bounding box at spacing dx, with mask[i] true where grid point i lies inside outer but outside inner — i.e. in the lensing band (AART grid_mask, winding-number point-in-polygon). grid is M×2; N is the per-side grid count.

Signed polygon area (shoelace); |area| is the enclosed area for a simple loop.