Built a real-time gravitational lensing shader pair for a personal project: single Schwarzschild black hole, and a binary inspiral/merger set-piece. Sharing the code and, more importantly, a full technical writeup of a numerical instability I hit and didn’t fully solve — figured this community would find the failure mode as interesting as the working part.
Repo: GitHub - guruprasadregar1-afk/blackhole-physics · GitHub (MIT)
Writeup: [repo]/docs/WRITEUP.md
Physics model
Single-body: standard Schwarzschild null geodesic, integrated with fixed-step RK4 per pixel in GLSL:
d²u/dφ² + u = 3Mu² where u = 1/r
Validated against known results in a CPU-side reference (weak-field deflection δφ ≈ 4M/b, photon sphere at r = 3M, event horizon capture at r = 2M) — 50 unit tests in the package, no GPU required to check the math.
Binary extension: two-body lensing via summed Schwarzschild deflections (documented as a heuristic superposition, not true nonlinear two-body GR — no closed-form solution exists), inspiral timing from the Peters (1964) GW-decay formula, and a quadrupole strain visualization for the merger.
The interesting failure
The binary merger vista showed hard radial spoke artifacts — looked exactly like dome-mesh faceting at first (48×32 sphere, vWorldPos interpolation). Rewrote it to a true per-pixel fullscreen ray-cast path (uInverseProjection/uInverseView, same pattern as reconstructing rays from NDC). Confirmed via GPU shader-source capture that the fullscreen path was actually compiled and running.
The spokes didn’t go away.
Built a float64 CPU mirror of the exact GLSL integrator to test it properly rather than guessing from screenshots. Findings, on a 128×72 grid:
- ~9.8% of adjacent screen-space rays flip to a different discrete outcome class (disk-hit vs. capture vs. escape) despite near-identical input directions
- Two-body shows ~2x the adjacent-outcome-flip rate of the equivalent single-body trace (8.8% vs 4.3%) under the same fixed-step RK4 pattern
- Flip rate spikes to ~29% near the horizon vs ~10% in the far field — classic strong-field sensitivity
Root cause: hard threshold predicates (u >= uHorizon, disk-plane-crossing sign flip, escape condition) compound with fixed-step RK4’s lack of adaptation near strong-field regions — tiny differences in adjacent ray paths cross those thresholds differently, producing sharp discrete boundaries in the render.
Tried softening the thresholds into smooth blended weights instead of hard returns. Measured result: adjacent-pixel discontinuity metrics got slightly worse, not better (median luminance diff 1.17 → 1.59). Reverted to the dome path, which incidentally reduces the visible artifact by quantizing ray angles to mesh vertices — a mask, not a fix. Documented all of this plainly in the writeup rather than shipping the fullscreen version with the bug hidden by camera angle luck.
If anyone’s solved something like this for multi-body geodesic tracing — adaptive stepping strategies, a better classification scheme than hard thresholds, or a different way to structure the two-body integration — genuinely curious how you approached it. My best current guess is the fix needs either adaptive step size near strong-field regions or a restructured 3D path integration (currently locked to a 2D orbital plane relative to body 1, which seems to be part of what’s amplifying sensitivity near the secondary body).
Not claiming this is research-grade or competitive with offline VFX — it’s explicitly scoped as an educational/outreach reference, and the writeup says so plainly. Screenshots below: both show the faceted-edge artifact discussed above — present on the single-hole dome as well as the binary case, though the underlying cause differs (single-hole is dome-mesh tessellation; binary layers a genuine numerical instability in the two-body integrator on top of that). Left in deliberately rather than cropped to hide it.
