Linux DevCenter    
 Published on Linux DevCenter (
 See this if you're having trouble printing code examples

Building a Black Hole With OpenGL

by Chris Halsall

This article and accompanying program are dedicated to Stephen Hawking, one of the greatest explorers of our time. His book A Brief History of Time is a must read for anyone interested in the universe in which we find ourselves.

Related Articles:

Creating Real Time 3D Graphics With OpenGL

OpenGL Rendering and Drawing

Preparing to Build an OpenGL Application

Previous Features

More from the Linux DevCenter

Black holes might be considered the ultimate heroes of gravity. They have incredible density and generate such powerful gravitational fields that not even light can escape their proximity. Black holes were first suggested as early as the late 1700s and were made famous with Einstein's law of general relativity.

A black hole is so compressed that it folds space and time around itself. To the outside observer, it no longer has size and becomes a zero-dimensional singularity. If any physical size is at all meaningful when speaking about a black hole, it would be its event horizon -- the sphere around the singularity from which nothing will ever been seen again.

Invite 1000 particles to play around a black hole, and you just know things are going to get crazy.

Invite 1000 particles to play around a black hole, and you just know things are going to get crazy.

Only as recently as the mid-1970s were the effects best explained by black holes first observed, thus providing direct physical evidence of their existence. But they have been generally accepted by both the public at large and by scientific communities for decades, thereby becoming a favorite resource for both sci-fi story teller and theoretical physicist. Black holes are just so extreme -- literal points of no return.

Because of the incredible gravitational field a black hole generates, the space around it tends to be a very dynamic place. Any particle with mass that enters the orbit of the black hole will quickly accelerate as its orbit decays and the particle moves closer to the black hole itself. The gravity shear, or gradient, is also unique in that the force of gravity exists at different strengths at only slightly different distances from the black hole. Objects are either torn apart immediately, or end up spinning at high speeds and are ripped apart because of centripetal forces.

Even light particles (photons) that pass near the black hole, but outside its event horizon, will still be affected. Light bends around the singularity as a function of the inverse of its distance -- or -- the closer it is, the more it bends. Black holes can be thought of as gigantic lenses. There have been recent astronomical observations of distant starlight bending around what is assumed to be a black hole that happens to be in between ourselves and the more distant star.

You can build a black hole with OpenGL

It has been estimated that a hydrogen bomb made with all the water on the planet would likely cause the matter in the center of the bomb to be compressed into a black hole. While we wait for the accounting department to reject our request, we might as well look for another way to have some fun with black holes that might cause less of a stir.

Drawing the axis lines.

Drawing the axis lines.

Program and Makefile

• bhole.c - The OpenGL Black Hole simulator.

• Makefile - Simple Makefile.

What we can do fairly inexpensively is simulate the physics of a very strong gravity field acting upon a number of particles. By simulating a thousand particles or so, a very pleasing visual display can be created. Using a simple iterative solution, the physics can be calculated quickly, while still resulting in true "Kepler" orbits, within limits.

We'll use OpenGL as the 3D API in order to leverage any 3D hardware that might be available on the machine. Some basic understanding of OpenGL is assumed in this article. For an OpenGL refresher, please see these articles. The accompanying program is similarly structured as well.

First we need a universe

Before we can build our black hole, we need a universe within which to place it. To make the math easier when calculating the particle orbits, we're going to make our universe center on the black hole at (0,0,0), which will also be the only source of gravitational pull.

In the code, main() first calls the usual OpenGL GLUT initialization functions, registers several of our functions for event and drawing callbacks, and then calls ourInit(). This function's job is to allocate and configure everything else needed by our program before passing control off to OpenGL's GLUT for event processing.

ourInit() starts by calling ourBuildTextures(), which builds a single bit-plane texture known as an alpha-channel map, containing a simple dot. The texture is used to control how transparent our particles will be, and is progressively more opaque towards the center. There is a 6% randomness to the alpha channel values, which hides any obvious transitions. Remove the + ourRand(15), and you'll be able to see obvious concentric rings when a particle is nearby.

Of course, 3D motion is required.

Of course, 3D motion is required.

ourBuildStarfield() is called next, with the number of stars desired passed as a parameter. This function builds a display list, which is a way of having OpenGL remember a sequence of calls for use later. The glNewList() function takes an integer "name" as its first parameter; we're calling with STAR_FIELD, which has already been defined as "1." The second parameter can either be GL_COMPILE or GL_COMPILE_AND_EXECUTE. Because we're pre-assembling this list, we'll use the first option.

ourBuildStarfield() takes advantage of the fact that OpenGL display lists record all the parameters made to the OpenGL calls. Thus, we don't have to worry about remembering where our stars are, because OpenGL will do it for us. If we didn't use a display list, we would have had to remember the star locations (and sizes), or use a repeatable pseudo-random number generator.

Another trick the function does is use OpenGL to spin the universe around a random amount while always placing the next star in the same location. This is instead of calculating the star's location ourselves, and demonstrates how rotating the MODELVIEW frame of reference can be used to quickly produce interesting results.

It should be noted that this results in slower run-time rendering speed, however, as OpenGL is having to do a matrix multiply (updating 16 values) for each rotation, which is not always implemented in hardware. A star field with pre-calculated star positions could easily be written, and would be recommended for applications where frame rate is critical.

Returning to ourInit(), calls are made to enable blending, but disabling alpha testing. We also ask for flat shading, since we'll see no benefit from smooth shading in this application, and then we ask for our blend function to be (GL_SRC_ALPHA,GL_ONE).

This is a special blending mode which says that whatever is in the destination buffer should be added to by the contents of the source buffer, according to the source buffer's alpha level. The destination is what has already been drawn, with the source being whatever it is that's currently being drawn.

This blending mode must be used with care in most situations, as it's a rather unnatural form of matter. Objects rendered in this way are completely transparent (as in, they don't block any light from behind), but they also add to the scene any light reflecting off or radiating from the object. No real world objects are truly like this, although fire, laser beams in dusty rooms, and luminescent particles come close.

Since we're modeling the latter, this mode makes sense. Note that overlapping objects quickly add to the RGB values of a pixel, so individual particles should be fairly transparent. A pixel's color values will continue to climb as additional objects are drawn on top of it, with clipping occurring at white. This can be, but is not always, a desirable effect.

ourInit() next calls ourCalcObs() which is used to figure out the X, Y, and Z coordinates based on the observer's angle around, distance from, and height above the origin. Then the cbResizeScene callback function is used in order to set the correct screen perspective matrix. gluNewQuadric() is next called in order to allocate a GLU object, which we'll use for the sphere of our black hole's event horizon.

Lastly, ourAllocParticles() is called to allocate an array of particles of the initial size desired. This function is handy because it can be called with the number of desired particles, and the array will be reallocated as needed. To finish off, the first particle in the array (our "root particle") is introduced to the system with a call to ourFireParticleGun().

At this point, our universe is ready to roll. ourInit() returns, main() prints a small help message to STDOUT, and control is passed to GLUT by calling glutMainLoop().

Hey! We're in orbit!

As cbRenderScene() was registered as both the Idle and Display callbacks, that is where we next receive control back. Here we draw the next frame of our display, and if motion is enabled, request our "physics engine" to update the position of each particle.

Heads Up status display turned on.

Heads Up status display turned on.

First cbRenderScene() clears the screen (including the depth buffer), ensures we're working with the MODELVIEW, and loads the "identity" matrix. Note that the DEPTH_TEST and Depth Mask are both enabled. Next, gluLookAt() is called to look at (0,0,0) from either the location of the observer, or from the root particle, depending on the global variable On_Orbit. This, and many other variables, are controlled from the keyboard, and are detailed below.

Texturing is disabled, and then a truly black sphere is drawn by way of a call to gluSphere(). The stars, particle vectors, and axis are drawn next, as appropriate, and then the depth mask is set to false, which disables writing any new depth information to the depth buffer. This is in preparation for the drawing of the particles themselves, which shouldn't change the depth buffer because they're completely transparent. Note that depth test is still enabled so particles behind the black hole will be obscured.

For the particles themselves, texturing is re-enabled if appropriate, and then the array of particles is iterated through. For each particle that isn't yet "running," there's a 0.03% chance it will be introduced into the system with a call to ourFireParticleGun(). This is a way of slowly introducing particles to the system at the start.

For particles that are already running, we simply set the color as defined by the particle (which is yellow for the root particle), and then draw two intersecting quads centered on the particles position (a "poor man's" billboard). If motion is enabled, we call ourMoveParticle(), the entry point for our physics engine, which will calculate the next position for the particle.

After the particle loop, cbRenderScene() calls ourRenderHeadsUp() to draw the heads-up display, if desired. At this point all drawing is complete, so it swaps the OpenGL display buffers so our just-drawn scene is on-screen. Single-stepping is handled next, then frame-per-second statistics are collected, and finally control is returned to GLU.

"Scotty, I need impulse power or we're all dead!"

Two of the functions registered as callbacks are cbKeyPressed() and cbSpecialKeyPressed(). These are used to adjust various variables for the simulation, and include "impulse engines" to adjust the velocity (and thus the orbit) of the root particle.

Some place you really don't want to go today.

Some place you really don't want to go today.

The controls for the root particle are 'O' for toggling between the observer's position and being "on-orbit" with the root particle. The 'x,' 'y' (or 'c'), and 'z' keys provide the root particle with a small impulse in the positive direction of the respective axis; the shifted capital keys provide a negative impulse. The velocity impulses are too small to directly escape the black hole's gravity field, but over several orbits, escape is possible with carefully planned periods of impulse.

It's also possible to move the root particle into a closer orbit. In fact, playing with different impulse directions at different points in the root particle's orbit can be an entertaining way of learning about orbital mechanics. The short form is -- everything's backwards. You slow down to go closer (and thus faster), and any impulses done on one side of the orbit have an effect on the opposite side.

While changing the root particle's orbit, it's often handy to have the axis lines displayed. This is done with the 'A' key. The 'S' key toggles the stars display, and the 'H' key toggles the "Heads Up" display, which gives several current values including our particle gun's location and injection velocity and the root particle's position and velocity. Such information can be useful when changing the orbit. If you get into trouble, the 'R' key will re-fire the root particle out of the particle gun.

When not on-orbit, the observation point is controlled with the arrow keys to change the angle around the Y axis and the height (above/below the XZ plane). These let you spin around the black hole, getting different views of the orbits. The Page Up and Page Down keys control how close the observation point is from the Y axis. The function ourCalcObs() converts these values into XYZ position coordinates for our drawing routine.

Other controls include the 'B' key for particle brightness, and the '+' and '-' keys to introduce or remove a hundred new particles to the system. Using these keys, you can wipe out all the particles in orbit and then restart, without having to stop the program. Lastly, the numeric keypad (when enabled) controls the position and injection velocities for the particle gun. Changing these values can result in very different and interesting orbits.

The particle gun itself is simply a conceptual point in space from which particles are injected into the system to begin their orbits around the black hole. There is some randomness given to the particles fired from the gun, to avoid "follow-the-leader" type effects. Occasionally the randomness is set to be extreme or to be backwards. This is done to add to the complexity of the scene, and it's always entertaining when a particle on an unusual orbit happens to come whipping by our point of observation.

The gun is fired whenever a particle is eaten by the black hole or escapes, occasionally randomly, and also whenever a particle is not yet running and the 'i' key is pressed. The 'I' key injects all non-running particles at one time, and can be quite amusing to watch. One last key of some importance: 'm' will start single-stepping through the simulation, with 'M' re-starting normal motion.

The physics engine

The term "physics engine" is a black-box wrapper phrase for the function, or set of functions, that implement the physics for a system. In our case, we are going to be calculating the gravitational effects of a massive object on a number of particles; no interactions between particles will be considered.

Making sure the algorithim exhibits Kepler's Laws.

Making sure the algorithim exhibits Kepler's Laws.

The math needed to calculate accurately the motion of objects within a gravity field (in one, two or even three dimensions) is actually quite simple. We'll be using an iterative method to solve what's called the "two body" problem. We simplify the problem even further to assume the smaller particle has no effect on the black hole.

For each moment in time, we update each particle's position based on its current velocity vector. Next, we calculate the square of the distance between our particle and the black hole. We use this as the (possibly scaled) denominator of a gravity constant to find the force of gravity being felt by the particle, thus implementing the inverse-square-law of gravity (and most other forces). The particle's current velocity vector is then updated to include the force of gravity (the vectors are added together). Lastly, we trim a small amount of velocity away from each particle, in order to simulate a slow decay of the orbits.

This simple algorithm produces true orbital dynamics which subscribe to Kepler and Newton's laws, providing the particles don't approach the singularity too closely. As seen in the image above, the orbits are ellipses, as first realized by Johannes Kepler in the early 1600s after studying the orbits of the planets. Newton later came up with his law of universal gravitation, and showed that orbits can also be parabolic or hyperbolic. Both of the latter orbit forms are "open" (the particle escapes), and can appear in our simulation.

If a particle comes too close to the black hole, however, the simulation must iterate at a finer temporal (time) resolution or it will no longer accurately simulate the physics. For our simulation, we simply kill any particles that get too close, assuming they were eaten by the black hole ("thinking it were a carrot"). The reason for the error is worth working through though, in order to understand why the algorithm fails to conserve energy.

If you examine the image below, you'll see a magnified view of a particle that is so close to the black hole that the force of gravity (drawn in red, skewering the black hole) is close to the same magnitude as the particle's velocity (drawn in green). Note that the velocity vector has already been adjusted by the pull of gravity, so at the next step in time the particle will move to the tip of the green line, the gravity vector will then be recalculated from that point, and the velocity vector will be updated for the next step.

Viewing the particle vectors -- one particle is REALLY close.

Viewing the particle vectors -- one particle is REALLY close.

The problem should be fairly obvious from the image: The direction of the vector representing the pull of gravity should be varying by approximately 45 degrees or more during this step's time period. As a result, as the particle passes along the green velocity vector, the vector itself should be bending towards the black hole, bringing the particle closer to its capture.

Instead, because our temporal resolution is too low, the particle will actually get a velocity boost. This can be amusingly demonstrated by removing the test condition to eliminate particles within the EVENT_HORIZON, which results in particles being fired off to infinity at high speed. Use the single-step mode (the 'm' key) with the vectors display (the 'V' key) to watch the vectors of close passing particles. Rotate around the black hole to get a good perspective at each step.

The solution to this problem is fairly simple -- you must break the temporal steps down into a finer resolution, with the force of gravity being divided evenly across each step. This can either be done for all the particles, or only for those that are too close. Keep in mind also that the error is a function of the inverse square of the distance, so the number of steps needed to maintain accuracy will vary considerably for each particle.

Because we're only playing, and are more interested in high frame rates than in perfect simulation, this doesn't matter too much. But for truly accurate predictions, the temporal scale must be quite small for strong gravity fields.

Don't want a 3D point source? Eliminate terms!

This simulation accurately handles frictionless motion within the 3D space around a point-source gravity field. By eliminating terms (by setting constants for some of the equations such that they come out to be 0 or 1, as appropriate, to be able to avoid some calculations), one-dimensional accelerated frames of reference can be simulated that are appropriate for on-planet environments.

Traditionally, such simulations assume the force of gravity will be the same for all objects, and so the calculation of object to gravity source distance can also be eliminated. Thus, the change in velocity (i.e. the acceleration) becomes a constant. Here on earth, the acceleration experienced is 9.8 meters per second per second (m/s2). Note that the mass of the objects doesn't enter the equations -- the mass doesn't have any effect except when dealing with inter-object gravitational influences (or when introducing "drag," simulating atmospheric friction).

Eliminate terms to change the physics for terrestrial environments.

Eliminate terms to change the physics for terrestrial environments.

Shown in the image above is the result of our black hole simulator with the ourMoveParticle() function modified such that gravity is only felt along the Y axis, and it is a constant rather than being a function of altitude. Also, a test is made to see if the particle has moved beneath the "ground" (Y < 0), and if so, the Y component of the velocity vector is reversed. Ta da! The world's largest trampoline! Don't feel so nice? Make the particles crash and take damage instead.

Closing notes

If we were interested, the singularity in this simulation could be eliminated, and a mass value could be added for each particle. By calculating the distance between each particle and the force of gravity each particle's mass exerted upon all the others, accurate n-body predictions could be made. In this way, our simulation would grow from an "m" times 2-body environment to a truly chaotic n-body system.

It is in some ways reassuring to know that the simple algorithm, "foreach (Particle) { P += V; V += GravSum(ParticleList); }", turns out to be the only way to calculate n-body trajectories when "n" grows to be more than (currently) about six; the more direct equations quickly become too complex for we puny humans, or our computers, to solve. In fact, specialized hardware has been built that does nothing but run essentially the n-body version of this algorithm to simulate our solar system at very fine temporal resolution for highly accurate predictions.

In the less abstract, limited term simulations can be used for various effects such as simulated jet exhaust smoke trails, artillery trajectories, or fireworks. By replacing the ourMoveParticle() function, radically different types of simulations can be built such as non-deterministic bird or bee behavioral motion.

Program and Makefile

• bhole.c - The OpenGL black hole simulator.

• Makefile - Simple makefile.

Particle systems can be a wonderful way of adding complexity to an environment with a minimum of programmer effort, providing the hardware can support enough particles to make the aggregate effect convincing. The results can often be mesmerizing, as I found during the construction of this program.

And, of course, no software is finished, only abandoned; there are lots of things that could be added to this program to make it even more interesting. Suggestions include optical warping effects around the black hole of those particles behind it, rendering each particle using a single quad as a true "billboard" facing the observer, and doing some kind of visual effect when a particle is captured. Having a mode where the particle's color is a function of its position or one of its vectors would likely be educational too.

It's your universe; imagine, then create!

Chris Halsall

Chris Halsall is the Managing Director of Ideas 4 Lease (Barbados). Chris is a specialist... at automating information gathering and presentation systems.

Related Articles:

Creating Real Time 3D Graphics With OpenGL

OpenGL Rendering and Drawing

Preparing to Build an OpenGL Application

Discuss this article in the O'Reilly Network Linux Forum.

Return to the Linux DevCenter.


Copyright © 2009 O'Reilly Media, Inc.