Coding the Cosmos: From Equations to Black Hole Simulations in C++

Introduction

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.

The Physics Foundation

Understanding Schwarzschild Geometry

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.

Implementing the Geodesic Equations

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.

The Role of Conserved Quantities

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.

Programming Architecture Evolution

The Natural Progression from 2D to 3D

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.

The 2D Phase (2D_lensing.cpp) establishes the foundational physics engine. Here, the code works with polar coordinates, implements RK4 integration, and visualizes light ray paths using OpenGL immediate mode rendering. This phase serves as both a learning exercise and validation step, in which debugging physics equations becomes much easier when individual rays can be watched tracing their paths across the screen.

Building on this foundation, the 3D Phase (black_hole.cpp) scales everything up dramatically. This version implements GPU compute shaders for performance, adds intuitive camera controls, includes accretion disk rendering, and provides spacetime grid visualization. The transition between these phases demonstrates how complex simulations can be built incrementally, with each step building confidence and understanding.

Why Numerical Integration Matters

One of the most important technical decisions Kavan makes is transitioning from Euler integration to Runge-Kutta 4th order (RK4). When I first looked at this choice, I wondered if it was really necessary , but then I considered the physics involved.

Euler integration suffers from accumulated errors that become catastrophic when simulating curved paths around black holes. A small directional error can send a light ray on a completely unphysical trajectory. Here's how Kavan implements the more stable RK4 approach:

                  void rk4Step(Ray& ray, double dλ, double rs) {
    double y0[4] = { ray.r, ray.phi, ray.dr, ray.dphi };
    double k1[4], k2[4], k3[4], k4[4], temp[4];

    geodesicRHS(ray, k1, rs);
    // ... four-stage RK4 implementation
    
    ray.r    += (dλ/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
    ray.phi  += (dλ/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);
    // ...
}
                

This RK4's four-stage approach provides fourth-order accuracy, which is essential for maintaining stable orbits and preventing rays from taking unphysical "shortcuts" through the black hole. The difference in visual quality between Euler and RK4 integration is immediately apparent when you see the simulation running.

Tackling the GPU Acceleration Challenge

Perhaps the most impressive aspect of Kavan's implementation lies in addressing the computational bottleneck. The scale becomes clear when calculating the numbers: for an 800x600 resolution with 20,000 integration steps per ray, the system must handle nearly 10 billion calculations per frame. No CPU could manage this workload in real-time.

His solution elegantly leverages GPU computing:

                  void dispatchCompute(const Camera& cam) {
    int cw = cam.moving ? COMPUTE_WIDTH : 200;
    int ch = cam.moving ? COMPUTE_HEIGHT : 150;
    
    glUseProgram(computeProgram);
    uploadCameraUBO(cam);
    uploadDiskUBO();
    uploadObjectsUBO(objects);
    
    glBindImageTexture(0, texture, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA8);
    
    GLuint groupsX = (GLuint)std::ceil(cw / 16.0f);
    GLuint groupsY = (GLuint)std::ceil(ch / 16.0f);
    glDispatchCompute(groupsX, groupsY, 1);
}
                

This implementation uses OpenGL compute shaders with 16x16 thread groups, perfectly leveraging GPU parallelism for the embarrassingly parallel problem of ray tracing. Each pixel can be calculated independently, making this approach ideal for GPU acceleration.

Advanced Graphics Techniques

Efficient Data Management with UBOs

As I dug deeper into the 3D implementation, I discovered that Kavan uses Uniform Buffer Objects (UBOs) for efficient data transfer to the GPU. This might seem like a minor technical detail, but it's actually crucial for performance:

                  struct UBOData {
    vec3 pos; float _pad0;
    vec3 right; float _pad1;
    vec3 up; float _pad2;
    vec3 forward; float _pad3;
    float tanHalfFov;
    float aspect;
    bool moving;
    int _pad4;
} data;
                

Notice the explicit padding throughout the structure. This ensures std140 layout compatibility, preventing subtle but frustrating alignment issues between CPU and GPU memory layouts. These kinds of details separate working code from professional-quality implementations.

Smart Performance Optimizations

One feature that also caught my attention was the adaptive quality system. Rather than always rendering at maximum quality, Kavan implements dynamic quality adjustment based on camera movement:

                  int cw = cam.moving ? COMPUTE_WIDTH : 200;
int ch = cam.moving ? COMPUTE_HEIGHT : 150;
                

This clever optimization provides smooth interaction during camera movement while maintaining high quality for static frames. Such user experience considerations make the difference between a tech demo and a genuinely usable application.

Making Spacetime Visible

Perhaps the most educational feature involves spacetime visualization through grid deformation. This addresses a common challenge in teaching general relativity: making abstract geometric concepts tangible. Kavan's solution proves remarkably effective:

                  void generateGrid(const vector<ObjectData>& objects) {
    // ... generate flat grid points
    
    for (const auto& obj : objects) {
        double r_s = 2.0 * G * mass / (c * c);
        double dx = worldX - objPos.x;
        double dz = worldZ - objPos.z;
        double dist = sqrt(dx * dx + dz * dz);
        
        if (dist > r_s) {
            double deltaY = 2.0 * sqrt(r_s * (dist - r_s));
            y += static_cast<float>(deltaY) - 3e10f;
        }
    }
}
                

The grid deformation formula approximates the famous "rubber sheet" analogy commonly used to explain general relativity. By watching how the grid bends and warps around the black hole, students can develop intuitive understanding of how mass curves spacetime.

Evaluating Physics Accuracy

What Works Well

After thoroughly analyzing the implementation, several aspects demonstrate impressive physical accuracy. The light bending follows true geodesic paths rather than simplified Newtonian trajectories, which proves crucial for realistic results. The event horizon acts as an absolute boundary: no information or light escapes once crossing this threshold. Conservation of energy and angular momentum remains mathematically rigorous throughout the simulation.

Furthermore, the Schwarzschild radius calculation is implemented correctly, and the resulting gravitational lensing effects match expectations from general relativity theory. These strengths make the simulation valuable both as an educational tool and as a foundation for more advanced work.

Acknowledging the Limitations

However, I should also note some necessary simplifications that limit the simulation's complete realism. Most significantly, the implementation only handles non-rotating black holes. In reality, virtually all astrophysical black holes rotate (described by Kerr geometry), which creates additional frame-dragging effects that would alter the light paths.

Additionally, the simulation assumes vacuum spacetime, meaning it doesn't account for matter density or electromagnetic fields that would affect light propagation in real astrophysical environments. The simulation also uses classical ray tracing without quantum effects like Hawking radiation, and treats the accretion disk as a simple geometric object rather than a complex plasma system with magnetic fields and relativistic motion.

Nevertheless, these limitations don't detract from the educational value. They represent reasonable simplifications that make core physics concepts accessible while keeping the implementation manageable.

Performance Deep Dive

Understanding the Computational Scale

The performance characteristics reveal staggering numbers. The simulation's performance scales as O(W × H × N) where W×H represents screen resolution and N indicates the number of integration steps. For real-time performance at reasonable quality:

• 800×600 resolution requires 480,000 individual rays

• Each ray needs approximately 20,000 RK4 integration steps

• Each RK4 step involves 4 function evaluations

• Total: roughly 38 billion function evaluations per frame

This scale explains why GPU acceleration isn't just helpful. Tt becomes absolutely essential for real-time performance.

Optimization Opportunities

Several areas exist where performance could be improved further. Adaptive step sizing could use smaller steps near the black hole where curvature is extreme, but larger steps in relatively flat regions far away. Early ray termination could stop integration when rays clearly escape to infinity, and hierarchical integration could use different quality levels for different screen regions based on importance.

These optimizations would prove valuable for scaling to higher resolutions or adding more complex physics, though they would add implementation complexity.

Educational Impact and Applications

Teaching Relativity Through Visualization

As someone interested in physics education, I find this implementation serves as an outstanding teaching tool for several reasons. First, it provides visual intuition. People can actually see how spacetime curvature affects light paths rather than just working with abstract equations. The interactive exploration capabilities allow real-time investigation of different perspectives and scenarios.

Most importantly, the code accessibility means that the progression from 2D to 3D makes these complex concepts approachable for those with programming backgrounds. They can modify parameters, experiment with different scenarios, and see immediate visual feedback.

Broader Research Applications

While necessarily simplified for real-time performance, the techniques Kavan demonstrates could extend to more sophisticated applications. Researchers working on gravitational wave visualization could build on this foundation. Scientists comparing theoretical predictions with Event Horizon Telescope observations could use similar ray-tracing approaches. The techniques also have applications in relativistic astrophysics education and computer graphics research in physically-based rendering.

Technical Implementation Excellence

Camera System Design

I was particularly impressed by the orbital camera implementation, which provides intuitive navigation while maintaining physical realism:

                  vec3 position() const {
    float clampedElevation = glm::clamp(elevation, 0.01f, float(M_PI) - 0.01f);
    return vec3(
        radius * sin(clampedElevation) * cos(azimuth),
        radius * cos(clampedElevation),
        radius * sin(clampedElevation) * sin(azimuth)
    );
}
                

The camera always orbits the black hole center, preventing users from accidentally positioning the viewpoint inside the event horizon where the physics would break down. This kind of thoughtful constraint prevents confusing or unphysical scenarios.

Correct Ray Generation

The 3D ray generation properly handles perspective projection and camera orientation through careful coordinate transformations. In the compute shader, rays are generated with:

                  // In compute shader
vec3 rayDir = normalize(
    camera.right * (screenCoord.x * camera.tanHalfFov * camera.aspect) +
    camera.up * (screenCoord.y * camera.tanHalfFov) +
    camera.forward
);
                

This ensures rays originate from the correct camera position with proper field-of-view distortion, maintaining the geometric relationships necessary for realistic rendering.

Final Thoughts

Kavan's black hole simulation represents a remarkable achievement in combining theoretical physics with practical computer graphics programming. The implementation successfully demonstrates how Einstein's field equations can be translated into real-time visual experiences, making abstract concepts tangible through interactive simulation.

The code progresses logically from fundamental concepts to advanced optimizations, serving as both an educational resource and a solid foundation for more sophisticated relativistic simulations. While the physics contains necessary simplifications, the core principles remain sound and the visual results prove genuinely striking.

The project showcases the remarkable power of modern GPU computing in making previously intractable physics simulations accessible to broader audiences. As computational resources continue to advance, implementations like this will become increasingly important tools for science education and public outreach.

For developers interested in physics simulation or educators seeking interactive tools for teaching relativity, this codebase provides an excellent starting point. It demonstrates that even complex astrophysical phenomena can be made accessible through thoughtful programming and clear presentation.

The work serves as a reminder that some of the most rewarding programming projects exist at the intersection of multiple disciplines - in this case, theoretical physics, numerical methods, computer graphics, and user interface design. Kavan has created something that's simultaneously educational, technically impressive, and genuinely engaging to explore.

Full credit goes to Kavan (YouTube: kavan, GitHub: kavan010) for creating this impressive simulation and making the source code available for educational use. His original video and codebase provide the foundation for this entire analysis.

GitHub Repository: kavan010's Black Hole Simulation

Original Video: Simulating Black Holes in C++

Read more at: Lazzerex’s Blog