Initial Thoughts
Speed
The output of the motion integration system needs to be fast enough for realtime physics calculation, on limited hardware (phones/tablets). It also needs to be fast enough for realtime orbit prediction, with some form of larger/adaptive timestep.Stable Orbits
Some bodies require very stable orbits, and only have forces applied on them by the body they are orbiting. For example, planets, space stations, moons, etc which could all be simplified to two body gravity problems to provide a form of artificial stability in the simulation, so long as the integration method is accurate enough to reproduce the underlying stability of the orbit.Unstable/Variable Orbits
Some bodies can have many different forces applied to them, for example ships, asteroids, debris.Conservation of Energy
If an orbit trajectory calculated in a two body problem is unstable, this might be slightly offputting at worst, and could even be put down instability of the underlying system (nearby moons, or the sun, changing the orbit), but if there is little conservation of energy, the body in orbit will soon crash into the surface.Summary
The perfect motion integration algorithm would be one which is very fast, very accurate in both the two body problem, and the n body problem, and finally fairly simple and easy to extend (with atmospheric, propulsion, collisions, explosions and other interesting effects, not just gravity).Testing Results
With that in mind, I set out to find out what method (or combination of methods) I should use. I did this by writing a python program as a demonstration and benchmark for the different popular motion integration algorithms. For the initial parameters/situation of this benchmark, I used planet earth, with the ISS orbiting on a slightly unusual and catastrophic elliptical orbit, like something out of the movie "Gravity". Credit goes to matplotlib for helping me plot these pretty graphs.For all but the Kepler example, the orbital trajectories have been calculated and plotted over a period of 100000 seconds, with the same timestep of 100 seconds.
Euler
Verlet
Time Corrected Verlet
Runge Kutta 4 (RK4)
Kepler's Formula
Performance Benchmarks
| Method Name | Results |
|---|---|
| Euler | 450884 function calls in 0.366 seconds |
| Verlet | 454897 function calls in 0.387 seconds |
| Time Corrected Verlet | 1835832 function calls in 1.638 seconds |
| Runge Kutta 4 | 1799879 function calls in 1.548 seconds |
| Kepler's Formula | 106710 function calls in 0.1112 seconds |
Discussion
Speed
As expected, both the Euler and Verlet numeric integration methods are much faster than the Runge Kutta 4 method. Kepler's formula requires an expensive precalculation step for each change in the orbit parameters. It is the fastest here because the benchmark is only for plotting out the orbit given the current orbit parameters. It is around 14 times faster than RK4 the only numeric method with comparable accuracy.
Accuracy and Stability
For the two body problem, there is little doubt that Kepler's formula will be the most accurate, because it solves the problem algebraically, and there is no propagation of error between orbits.
The attempt at time corrected Verlet integration method which I used, is obviously the worst, and with the fact that this one takes the longest, I see little reason to use it, although there is a good chance that I've made a mistake in its implementation.
RK4 as shown here is very stable, and quite accurate, but is 4 times slower than Euler and Verlet. Available literature suggests RK4 is still more efficient speed/accuracy (otherwise why would anybody use it, rather than just decreasing the timestep on either Euler or Verlet).
Conservation of Energy
They all seem to do a reasonable job at this I think (apart from time corrected Verlet), reasonable enough for the purposes of the game. Unfortunately, while this works out well for gameplay, the projection of orbits seems heavily affected by the accuracy. Stable orbits look a lot neater, and are a lot easier to comprehend. If unstable (or inaccurately calculated) orbits are to be commonplace, then a method for perhaps colouring the projected path by time, or allowing the user to select how far ahead they want to predict the orbit will be necessary.
Conclusion
For static scenery objects such as planets or moons, which can be treated simple two body problems, the obvious choice is to use Kepler's Formula. This requires an expensive precalculation step, but seeing as the orbital parameters are unlikely to change during gameplay, this shouldn't be a problem.
Static scenery objects are expected to be on extremely stable orbits, which makes Kepler another good match. There could be many static scenery objects, yet another reason to use Kepler's formula to quickly plot and calculate the motion of these objects.
For dynamic objects such as ships, debris, asteroids, etc, it is have a much harder problem. Not only is this potentially an N body problem, or at least a two body gravity, with extra applied forces (from thrusters, collisions, friction etc), but the orbit predictions are likely to be unstable, no matter which numeric method I choose to calculate them. This means that unless I use something like patched conics (which I'd like to avoid if possible), it's impossible to use kepler's formula to predict the orbit of this body when it is in freefall (only gravity forces being applied).
Calculating orbit trajectories using larger than realtime timesteps is required not only for showing the player where their ship, or other bodies are headed, but also as a form of precalculation for bodies in freefall, whether that the body be away from the player, or the player is using warp drive on their ship.
Large timesteps can be used to plot a course over a large period of time, but some form of interpolation will be required between the calculated points for a continous position of the object over time. Linear interpolation might be good enough for positioning the object on the screen, but it probably won't be good enough for the initial velocity required if calculation is resumed (object enters player's influence, or player exits time warp)
One possible solution is to use some form of polynomial/bezier interpolation between the points and apply it to velocity as well. Another solution, which I think might be better, is that when a body exits interpolated motion, it recalculates it's current position and velocity starting from the previous plotted point, but using a much smaller timestep.
So based on these results I've decided I'd like to use Kepler's Formular for static scenery objects, and some form of Numeric Integration, whether that be RK4, Euler or Verlet remains to be seen, but if performance isn't an issue, then preferably RK4.
In order to deal with the numeric calculations efficiently, each body will have a required timestep accuracy, which is based upon the player's position, whether the player is using warp drive, and various other influences. This timestep is used to precalculate an interpolation buffer, which is reset if the required timestep shrinks, as described in the previous paragraph.
Each set of dynamic objects which can interact gravitationally with each other should be self contained, and it should be possible to put each set into a seperate thread. Most dynamic objects don't interact with each other, and so it should fairly simple to make this suitable for parallel processing.
Most modern phones and tablets these days are quad core. A single thread/Core will be dedicated will be dedicated to input, bullet physics, graphics, and game logic. This leaves three other threads to be split over motion integration, and other system simulations.
Making the motion integration algorithm play ball with a general purpose physics engine such as Bullet will be the subject of future posts, but I'll record here a few thoughts about what I've been planning to try and do. The main thread will feed values out of the interpolation buffer into bullet, calculating the required linear velocity to reach the next interpolation point at the correct time. If a collision occurs, interpolation buffer is reset, and the resulting change in velocity is fed back to the motion integrator to be used in the initial parameters for buffer recalculation.
Scheduling will be tricky. When the main thread goes through and looks for values from the interpolation buffer, if any buffers are empty, thread/s are spawned, and they are immediately calculated to a reasonable distance (a couple of timesteps perhaps), pausing all the continuous threads operating on the buffer. This is because if the zero value of objects in the buffer are required but not available, this will block the main thread, and thus takes priority over all other calculations.
Anyway, that's enough for now. Next post I'm hoping to have a simple 2d touch example on Android using Verlet or RK4 in realtime, ready to show off.
