Open Events in the Horizon: A Black Hole Visualization Project

Ever since the first actual image of a black hole was released by the Event Horizon Telescope team in 2019, I’ve been fascinated by these cosmic objects. As a programmer with a curiosity for astrophysics, I challenged myself to create a physically accurate visualization of a black hole that would run in real-time.

This project, which I’ve named “Open Events in the Horizon” (OEH), combines Python, CUDA, and OpenGL to create an interactive black hole simulation that anyone can run on their computer.

The Science Behind the Visualization

While I’m not an astrophysicist myself, I based my simulation on peer-reviewed research. Specifically, I implemented the magnetically dominated accretion disk model from a 2003 paper by Pariev, Blackman & Boldyrev:

“Extending the Shakura-Sunyaev approach to a strongly magnetized accretion disc model” (Astronomy & Astrophysics, 407, 403-421)

This model accounts for how strong magnetic fields affect the structure and appearance of the disk of matter spiraling into the black hole. The result is a more realistic temperature profile and visual appearance than simpler models.

Technical Components

Creating this simulation involved several technical challenges:

1. Ray Tracing in Curved Spacetime

To visualize a black hole correctly, I had to trace the paths of light rays as they bend around the intense gravitational field. I implemented a modified ray tracer in CUDA that:

Approximates light paths in Schwarzschild geometry
Computes gravitational lensing effects
Calculates relativistic effects like Doppler shifting and redshift


<span>@cuda.jit</span><span>(</span><span>device</span><span>=</span><span>True</span><span>)</span>
<span>def</span> <span>trace_ray_equatorial</span><span>(</span>
    <span>cam_x</span><span>:</span> <span>float</span><span>,</span> <span>cam_y</span><span>:</span> <span>float</span><span>,</span>
    <span>dir_x</span><span>:</span> <span>float</span><span>,</span> <span>dir_y</span><span>:</span> <span>float</span><span>,</span>
    <span>mass_bh_cgs</span><span>:</span> <span>float</span><span>,</span>
    <span>max_steps</span><span>:</span> <span>int</span><span>,</span>
    <span>dt</span><span>:</span> <span>float</span><span>,</span>
    <span>horizon_radius</span><span>:</span> <span>float</span><span>,</span>
    <span>integrator_choice</span><span>:</span> <span>int</span>
<span>)</span> <span>-></span> <span>float</span><span>:</span>
    <span>"""</span><span> Trace a photon in the equatorial plane around a BH with improved physics. </span><span>"""</span>
    <span># Position </span>    <span>x</span><span>,</span> <span>y</span> <span>=</span> <span>cam_x</span><span>,</span> <span>cam_y</span>
    <span># Velocity (normalized to c) </span>    <span>vx</span><span>,</span> <span>dir_x</span> <span>*</span> <span>C</span>
    <span>vy</span> <span>=</span> <span>dir_y</span> <span>*</span> <span>C</span>
    <span># Physics implementation... </span>
<span>@cuda.jit</span><span>(</span><span>device</span><span>=</span><span>True</span><span>)</span>
<span>def</span> <span>trace_ray_equatorial</span><span>(</span>
    <span>cam_x</span><span>:</span> <span>float</span><span>,</span> <span>cam_y</span><span>:</span> <span>float</span><span>,</span>
    <span>dir_x</span><span>:</span> <span>float</span><span>,</span> <span>dir_y</span><span>:</span> <span>float</span><span>,</span>
    <span>mass_bh_cgs</span><span>:</span> <span>float</span><span>,</span>
    <span>max_steps</span><span>:</span> <span>int</span><span>,</span>
    <span>dt</span><span>:</span> <span>float</span><span>,</span>
    <span>horizon_radius</span><span>:</span> <span>float</span><span>,</span>
    <span>integrator_choice</span><span>:</span> <span>int</span>
<span>)</span> <span>-></span> <span>float</span><span>:</span>
    <span>"""</span><span> Trace a photon in the equatorial plane around a BH with improved physics. </span><span>"""</span>
    <span># Position </span>    <span>x</span><span>,</span> <span>y</span> <span>=</span> <span>cam_x</span><span>,</span> <span>cam_y</span>
    <span># Velocity (normalized to c) </span>    <span>vx</span><span>,</span> <span>dir_x</span> <span>*</span> <span>C</span>
    <span>vy</span> <span>=</span> <span>dir_y</span> <span>*</span> <span>C</span>

    <span># Physics implementation... </span>
@cuda.jit(device=True)
def trace_ray_equatorial(
    cam_x: float, cam_y: float,
    dir_x: float, dir_y: float,
    mass_bh_cgs: float,
    max_steps: int,
    dt: float,
    horizon_radius: float,
    integrator_choice: int
) -> float:
    """ Trace a photon in the equatorial plane around a BH with improved physics. """
    # Position     x, y = cam_x, cam_y
    # Velocity (normalized to c)     vx, dir_x * C
    vy = dir_y * C

    # Physics implementation...

Enter fullscreen mode Exit fullscreen mode

2. Physically Accurate Accretion Disk

I implemented the magnetically dominated disk model, which calculates:

Proper temperature profile based on radius and magnetic field strength
Realistic emission spectrum using blackbody radiation
Electron scattering effects that modify the spectrum
Turbulence and plasma physics at different disk radii


<span>@cuda.jit</span><span>(</span><span>device</span><span>=</span><span>True</span><span>)</span>
<span>def</span> <span>disk_radiance</span><span>(</span><span>r_cm</span><span>:</span> <span>float</span><span>,</span> <span>m_bh</span><span>:</span> <span>float</span><span>,</span> <span>nu</span><span>:</span> <span>float</span><span>,</span> <span>b_field_exp</span><span>:</span> <span>float</span><span>)</span> <span>-></span> <span>float</span><span>:</span>
    <span>"""</span><span> Returns total specific intensity I_nu from the disk at radius r_cm for frequency nu. Includes both blackbody and modified blackbody effects depending on the optical depth regime, following the Pariev+2003 model. </span><span>"""</span>
    <span># Get surface temperature at this radius </span>    <span>T</span> <span>=</span> <span>disk_surface_temperature</span><span>(</span><span>r_cm</span><span>,</span> <span>m_bh</span><span>,</span> <span>b_field_exp</span><span>)</span>
    <span># Calculate effective optical depth to determine emission regime </span>    <span>r_g</span> <span>=</span> <span>G</span> <span>*</span> <span>m_bh</span> <span>/</span> <span>(</span><span>C</span><span>*</span><span>C</span><span>)</span>
    <span>ratio</span> <span>=</span> <span>r_cm</span> <span>/</span> <span>(</span><span>10.0</span> <span>*</span> <span>r_g</span><span>)</span>
    <span># Implementation of physics model... </span>
<span>@cuda.jit</span><span>(</span><span>device</span><span>=</span><span>True</span><span>)</span>
<span>def</span> <span>disk_radiance</span><span>(</span><span>r_cm</span><span>:</span> <span>float</span><span>,</span> <span>m_bh</span><span>:</span> <span>float</span><span>,</span> <span>nu</span><span>:</span> <span>float</span><span>,</span> <span>b_field_exp</span><span>:</span> <span>float</span><span>)</span> <span>-></span> <span>float</span><span>:</span>
    <span>"""</span><span> Returns total specific intensity I_nu from the disk at radius r_cm for frequency nu. Includes both blackbody and modified blackbody effects depending on the optical depth regime, following the Pariev+2003 model. </span><span>"""</span>
    <span># Get surface temperature at this radius </span>    <span>T</span> <span>=</span> <span>disk_surface_temperature</span><span>(</span><span>r_cm</span><span>,</span> <span>m_bh</span><span>,</span> <span>b_field_exp</span><span>)</span>

    <span># Calculate effective optical depth to determine emission regime </span>    <span>r_g</span> <span>=</span> <span>G</span> <span>*</span> <span>m_bh</span> <span>/</span> <span>(</span><span>C</span><span>*</span><span>C</span><span>)</span>
    <span>ratio</span> <span>=</span> <span>r_cm</span> <span>/</span> <span>(</span><span>10.0</span> <span>*</span> <span>r_g</span><span>)</span>

    <span># Implementation of physics model... </span>
@cuda.jit(device=True)
def disk_radiance(r_cm: float, m_bh: float, nu: float, b_field_exp: float) -> float:
    """ Returns total specific intensity I_nu from the disk at radius r_cm for frequency nu. Includes both blackbody and modified blackbody effects depending on the optical depth regime, following the Pariev+2003 model. """
    # Get surface temperature at this radius     T = disk_surface_temperature(r_cm, m_bh, b_field_exp)

    # Calculate effective optical depth to determine emission regime     r_g = G * m_bh / (C*C)
    ratio = r_cm / (10.0 * r_g)

    # Implementation of physics model...

Enter fullscreen mode Exit fullscreen mode

3. Real-time Rendering Pipeline

The OpenGL rendering pipeline includes:

GPU-accelerated simulation using CUDA and Numba
Real-time post-processing effects
Interactive camera controls
Shader-based effects for bloom, exposure, contrast, etc.


<span>def</span> <span>render_frame</span><span>(</span><span>self</span><span>):</span>
    <span>"""</span><span>Renders a single frame by running the raytracer and displaying the result.</span><span>"""</span>
    <span># Apply auto-rotation if enabled </span>    <span>self</span><span>.</span><span>_update_camera_for_rotation</span><span>()</span>
    <span># Run the simulation only if not paused </span>    <span>if</span> <span>not</span> <span>self</span><span>.</span><span>paused</span> <span>or</span> <span>self</span><span>.</span><span>last_image</span> <span>is</span> <span>None</span><span>:</span>
        <span>t_start</span> <span>=</span> <span>time</span><span>.</span><span>time</span><span>()</span>
        <span># Run the simulation with the current parameters </span>        <span>image</span> <span>=</span> <span>run_simulation</span><span>(</span>
            <span>custom_camera_position</span><span>=</span><span>self</span><span>.</span><span>camera_position</span><span>,</span>
            <span>custom_black_hole_mass</span><span>=</span><span>self</span><span>.</span><span>black_hole_mass</span><span>,</span>
            <span>custom_fov</span><span>=</span><span>self</span><span>.</span><span>fov</span><span>,</span>
            <span>b_field_exponent</span><span>=</span><span>self</span><span>.</span><span>b_field_exponent</span><span>,</span>
            <span>integrator_choice</span><span>=</span><span>self</span><span>.</span><span>integrator_choice</span>
        <span>)</span>
        <span># Update image cache and timing </span>        <span>self</span><span>.</span><span>last_image</span> <span>=</span> <span>image</span>
        <span>self</span><span>.</span><span>render_time_ms</span> <span>=</span> <span>(</span><span>time</span><span>.</span><span>time</span><span>()</span> <span>-</span> <span>t_start</span><span>)</span> <span>*</span> <span>1000</span>
<span>def</span> <span>render_frame</span><span>(</span><span>self</span><span>):</span>
    <span>"""</span><span>Renders a single frame by running the raytracer and displaying the result.</span><span>"""</span>
    <span># Apply auto-rotation if enabled </span>    <span>self</span><span>.</span><span>_update_camera_for_rotation</span><span>()</span>

    <span># Run the simulation only if not paused </span>    <span>if</span> <span>not</span> <span>self</span><span>.</span><span>paused</span> <span>or</span> <span>self</span><span>.</span><span>last_image</span> <span>is</span> <span>None</span><span>:</span>
        <span>t_start</span> <span>=</span> <span>time</span><span>.</span><span>time</span><span>()</span>

        <span># Run the simulation with the current parameters </span>        <span>image</span> <span>=</span> <span>run_simulation</span><span>(</span>
            <span>custom_camera_position</span><span>=</span><span>self</span><span>.</span><span>camera_position</span><span>,</span>
            <span>custom_black_hole_mass</span><span>=</span><span>self</span><span>.</span><span>black_hole_mass</span><span>,</span>
            <span>custom_fov</span><span>=</span><span>self</span><span>.</span><span>fov</span><span>,</span>
            <span>b_field_exponent</span><span>=</span><span>self</span><span>.</span><span>b_field_exponent</span><span>,</span>
            <span>integrator_choice</span><span>=</span><span>self</span><span>.</span><span>integrator_choice</span>
        <span>)</span>

        <span># Update image cache and timing </span>        <span>self</span><span>.</span><span>last_image</span> <span>=</span> <span>image</span>
        <span>self</span><span>.</span><span>render_time_ms</span> <span>=</span> <span>(</span><span>time</span><span>.</span><span>time</span><span>()</span> <span>-</span> <span>t_start</span><span>)</span> <span>*</span> <span>1000</span>
def render_frame(self):
    """Renders a single frame by running the raytracer and displaying the result."""
    # Apply auto-rotation if enabled     self._update_camera_for_rotation()

    # Run the simulation only if not paused     if not self.paused or self.last_image is None:
        t_start = time.time()

        # Run the simulation with the current parameters         image = run_simulation(
            custom_camera_position=self.camera_position,
            custom_black_hole_mass=self.black_hole_mass,
            custom_fov=self.fov,
            b_field_exponent=self.b_field_exponent,
            integrator_choice=self.integrator_choice
        )

        # Update image cache and timing         self.last_image = image
        self.render_time_ms = (time.time() - t_start) * 1000

Enter fullscreen mode Exit fullscreen mode

Integration Challenges

One of the biggest challenges was integrating the physics, computation, and visualization components. I ran into issues with:

Performance optimization: Finding the balance between physical accuracy and real-time rendering
Memory management: Efficiently transferring data between the CPU, CUDA cores, and OpenGL
Numerical stability: Ensuring the integration methods worked properly in extreme conditions near the event horizon

The Development Process

This project stretched my programming skills and taught me a lot about GPU programming, physics simulation, and scientific visualization. Here’s how I approached it:

Research phase: Understanding the physics papers and translating equations into code
Core simulation: Implementing the basic ray tracing and physics calculations
GPU acceleration: Porting the code to CUDA for massive parallelization
Visualization: Creating the OpenGL-based visualization pipeline
Refinement: Continuous improvement of visuals and physical accuracy

I also leveraged modern AI tools (Claude, Gemini Pro, and GPT) to help with debugging and optimization. These tools were particularly helpful in identifying numerical issues and suggesting code optimizations.

Try It Yourself

The code is open source under the GNU license. You can try it yourself:


<span># Clone the repository</span>
git clone https://github.com/EricsonWillians/OEH.git
<span>cd </span>OEH
<span># Install with UV (preferred package manager)</span>
uv pip <span>install</span> <span>-e</span> <span>.</span>
<span># Run the simulation</span>
python <span>-m</span> oeh.main
<span># Clone the repository</span>
git clone https://github.com/EricsonWillians/OEH.git
<span>cd </span>OEH

<span># Install with UV (preferred package manager)</span>
uv pip <span>install</span> <span>-e</span> <span>.</span>

<span># Run the simulation</span>
python <span>-m</span> oeh.main
# Clone the repository
git clone https://github.com/EricsonWillians/OEH.git
cd OEH

# Install with UV (preferred package manager)
uv pip install -e .

# Run the simulation
python -m oeh.main

Enter fullscreen mode Exit fullscreen mode

You’ll need a CUDA-capable GPU to run the simulation at full speed.

Conclusion

Creating this black hole visualization was both a technical challenge and a fascinating learning experience. It demonstrates how consumer-grade hardware can now simulate complex physics that previously required supercomputers.

Even if you’re not an astrophysicist or specialized graphics programmer, modern tools and libraries make it possible to create impressive scientific visualizations with relatively accessible technology.

If you try out the project or have suggestions for improvement, I’d love to hear from you!

Resources

原文链接：Building a Physically Accurate Black Hole Visualization with Python, CUDA, and OpenGL

展开阅读全文

文章版权声明 1、本网站名称：拾光赋
2、本站永久网址：https://www.blogs.ink
3、本网站的文章部分内容可能来源于网络，仅供大家学习与参考，如有侵权，请联系站长QQ：805375623进行删除处理。
4、本站一切资源不代表本站立场，并不代表本站赞同其观点和对其真实性负责。
5、本站一律禁止以任何方式发布或转载任何违法的相关信息，访客发现请向站长举报
6、本站资源大多存储在云盘，如发现链接失效，请联系我们我们会第一时间更新。

THE END