Simulating black holes computationally presents one of the most challenging intersections of theoretical physics and computer graphics. So when I first stumbled across Kavan's black hole simulation video on YouTube, I was immediately struck by how he managed to bridge the gap between Einstein's abstract equations and tangible, visual reality. As someone fascinated by both physics and computer graphics, I found myself diving deep into his codebase, trying to understand how theoretical concepts translate into real-time rendering.
The project progresses elegantly from basic 2D ray tracing to a full 3D GPU-accelerated simulation. Along the way, it demonstrates gravitational lensing effects, accretion disks, and curved spacetime visualization through deformed grid rendering. Each component builds upon the previous one, making complex concepts accessible through careful implementation. Let me walk you through each component and explain why his approach works so well.
Before we can appreciate the programming implementation, I need to establish the physics foundation that drives everything. At the heart of any realistic black hole simulation lies the Schwarzschild metric, which describes how massive objects curve spacetime around them. Kavan's code correctly implements the Schwarzschild radius formula:
r_s = 2.0 * G * mass / (c*c);
This elegant equation represents the event horizon, the crucial boundary beyond which nothing, not even light, can escape. What I find particularly thoughtful about Kavan's approach is how he treats this as the visual radius of the black hole, which turns out to be physically accurate for a non-rotating black hole.
The key insight that makes his simulation work so effectively involves understanding that light rays follow geodesics: the straightest possible paths in curved spacetime. Rather than applying simple Newtonian gravitational forces (which would be incorrect for light), the simulation calculates how spacetime curvature itself affects light propagation. This distinction proves crucial and separates his work from simpler, less accurate approaches.
Moving from theory to practice, the 2D simulation (2D_lensing.cpp) implements the geodesic equation for null rays through what Kavan calls the geodesicRHS function.
void geodesicRHS(const Ray& ray, double rhs[4], double rs) {
double r = ray.r;
double dr = ray.dr;
double dphi = ray.dphi;
double E = ray.E;
double f = 1.0 - rs/r;
double dt_dλ = E / f;
rhs[2] = - (rs/(2*r*r)) * f * (dt_dλ*dt_dλ)
+ (rs/(2*r*r*f)) * (dr*dr)
+ (r - rs) * (dphi*dphi);
rhs[3] = -2.0 * dr * dphi / r;
}
Breaking down these equations reveals their physical significance. The first term captures gravitational redshift effects, while the second handles coordinate singularities inherent in radial coordinates. Meanwhile, the third term accounts for centrifugal effects in curved spacetime, and the final equation ensures angular momentum conservation. This attention to physical detail makes Kavan's simulation convincing and accurate.
What really impressed me about the implementation was how carefully Kavan handles conserved quantities. In relativistic mechanics, certain quantities remain constant along geodesic paths, and his code properly implements these:
L = r*r * dphi; *// Angular momentum per unit mass*
double f = 1.0 - SagA.r_s/r;
double dt_dλ = sqrt((dr*dr)/(f*f) + (r*r*dphi*dphi)/f);
E = f * dt_dλ; *// Energy per unit mass*
These conservation laws serve dual purposes: they ensure physical accuracy while providing numerical stability. Without them, the simulation would gradually drift away from realistic behavior as small computational errors accumulate over time.
Kavan's approach structures the project as a clear evolution from proof-of-concept to production-ready simulation. This progression makes complex concepts much more approachable than jumping directly into 3D implementation.