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.
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.
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.
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.