So, I had reached a point with my game prototype where I was interested in getting the environment working. Since I’m working on a game that heavily involves space, figuring out how to represent and simulate celestial entities (from asteroids to galaxies) would seem to be an important matter.

I have to pick a solution from the spectrum of solutions which ranges from the practical to the accurate. Most space games invariably choose the practical solutions, which is, typically, having an unchanging, logical and scaled-down representation of the universe. Scales aren’t represented properly; entities may have orbits, but these are fixed; and, finally, entities have a very logical, typically parent-child, relationship, where one would typically enforce planets to be associated with a star, stars with a galaxy, etc. Such models only represent the universe in its most general terms. I want to be able to show true scales, irregular celestial structures (rogue planets, binary/multiple star systems and so on) and a greater degree of dynamics (changing orbits, colliding entities), but I can accept a loss of precision as long as it looks believable.

Scene Representation

To represent the systems in graphics space, I’m experimenting with using a tree of irregular coordinate spaces, where each child node has to be a subdivision of its parent. Each node is then its own coordinate space which is a scaled up space of its parent. It has a defined boundary for each axis (a box) and can be freely positioned and oriented. Each celestial body would then occupy a coordinate space that is appropriate for it. The camera would also occupy a node in the tree, but at the scale where 1 mm can be accurately represented. Some rudimentary testing I did a while ago indicated that I could manage this requirement using a coordinate space that’s bounded to 10,000 units (equiv. to 10 km, if 1 unit = 1 m). To be on the safe side, I’ve bounded it to 1,000 units. The camera then asks the tree to move it to a node that has a boundary that contains its position and has the required scale. Obviously, there will be missing nodes in-between with our requirements (we have to divide several times if we are to go from a node with galaxy-wide boundary down to one that has boundaries acceptable for a camera). These are created dynamically in a grid-like manner and kept track of using a hash map (indexed by its integer division “position”) in the parent node. When rendering a node, we simply move upwards in the coordinate spaces (from the camera node) to the common node and then down to the celestial entity of interest.

Logical Representation

Each celestial entity is represented by a CelestialEntity object, which contains rendering details (models, etc), physical details and a set of references to child entities. This set is used to generate the local OrbitalPhysics object which determines the physical behavior of the system. A CelestialView provides the scene graph with the minimum amount of data needed to accurately render the universe for a specified camera and settings. Some of these settings relate to scale compression, and I’m experimenting with lerping to logarithmic scale to provide better overview of a solar system.

Simulation

As mentioned before, the OrbitalPhysics object determines how the celestial entities move. Early on, I tried to use PhysX to simulate orbits and had a lot of skepticism about whether it would work well. After scaling down all the values to make them tolerable to PhysX, it turned out this skepticism was warranted. As I recall, orbits decayed very fast and it had nowhere the stability I needed. I’m fairly sure PhysX uses a low-order integrator. So, I looked into how Celestia, which is the first program of the top of my head that performs scale-accurate rendering and animation of the solar system, handles these physics and as I understand it they use analytical solutions limited to our solar system (VSOP87) for determining positions. No help there, then. Then I found a good paper on solving N-Body problems using multi-step methods. I decided to go with the 9th order position-acceleration method which is the one he has used himself for his simulations, and am currently using leapfrog to get the initial steps. This had satisfactory performance.

For orbits, I sort all the orbital periods, evaluate celestial states within that period using appropriate time step size given by a division of the lowest orbital period, then reset to the optimal time step size for the second lowest orbital period and evaluate celestial states from the current time index, and so on. The entities will then have a collection of points in their orbits with equal or increasing distance in time between them compared to the previous, but never more than the optimal time step for each body within a single orbit.

Caveats

The scene representation has some issues with floating point precision which is exacerbated by the fact that matrices can be oriented freely, so there are occurrences where a node might not fit into a child node that’s located at the edge of its parent, but when transforming it into parent space, it fits into the parent’s boundary, so, logically, it must fit in the child and the whole thing barfs. I currently use some hacks (fuzzy boundaries, etc) to get around it, but I still want to figure out a better solution.

I am struggling somewhat with rendering accurate orbits, because I want to use greater time steps for the orbits and in some cases I’m getting wildly changing orbit data from frame to frame. I’m testing them on binary stars moving elliptically, and suspect it’s because of the very tight and steep gradients when they’re closest. The 9th order method handles it fine, but it seems like the leapfrog is giving very inaccurate values. So, when the stars are close and I have to recompute the orbit, it starts getting the first 8 values using leapfrog and the 9th order method doesn’t get accurate enough values to begin with. I might try some more advanced single-step methods instead of leapfrog for initialization, but I’m leaning towards trying to reuse some of the previous states instead which seems somewhat nontrivial.

Another issue is estimating orbital periods to determine optimal time step size for orbit computation. Currently, I’m simply assuming a circular orbit (using current velocity and acceleration),  calculating the period of that and adjusting it based on previous orbit data, but it only works as a very rough estimate.