Looking for advice on solar system integrator

[Deadstar]

I've been building a solar system integrator using C++, OpenGL for graphics, Qt for UI and I would like some info about the choice of integrator. The project is at a pretty advanced stage now and I'm currently adding all moons since graphically it's doing well, I can move about the whole system freely and I've got several integrators available to use mid-simulation (as well as the option to add in bodies during the simulation and control over time steps etc). The purpose of this project is still sort of undecided. I want it to run in real time, have decent graphics, have the ability to move throughout the solar system and be able to add bodies when you like. I'm not trying to get a perfect simulator but would like it be as accurate as possible while maintaining the above conditions (obviously violated by the add body feature but besides that). All of this (except the graphics aren't great yet) is currently done but as I want to ask some questions before I advance things any further.

I am currently using an Adams-Bashforth stage 5 integrator although I can also use stage 3 and 4 Adams (easily extendable to higher orders), Runge-Kutta, Euler integrators and have a decent framework set up so that I can add more pretty easily. This is fine and the results are really good but adding moons makes things slightly more awkward in that step sizes have to decrease a bit for the orbits to be anywhere near accurate for a decent length of time. I use data from the NASA HORIZONS site btw.

So, question time.

I was wondering if there are any integrators (surely must be) that are built for this sort of problem and use keplers laws as a basis and adds perturbations? It would have to deal with the effects of extra planets that I randomly throw in so a pure keplerian integrator wouldn't do it as it wouldn't take into account any other body, just the one it orbits.The way I currently set the integration up is to calculate gravitational effects using the centre of mass of planet-moon system and only use the effect from individual moons when an object is within the planets hill sphere. Does this seem like a decent idea and will it be accurate enough for short term (<1000 year) simulations? I figure I currently have scaled units such that 1 unit of time is a day, 1 unit of distance is 1AU and 1 unit of mass is the mass of the earth. Does this seem reasonable or will it become problematic with precision when including smaller moons (since also radius of moons are scaled by 1AU, they can get pretty small). Would perhaps using mercury mass and semi major as base units help this or not? Quite difficult to find a balance between scale when dealing with 100km moons and billion kilometer semi major axis.Again, I don't exactly have clear defined goals yet. It's sort of an exercise in integration, programming, graphics and that sort of thing. I do have a certain idea I want to implement but that would be after everything else is done... I suppose I would be working towards something like Universe Sandbox (although I can't get it to run on my laptop, I've seen images of it and looks very interesting). I think that uses Euler integration though with tiny time steps and as I'm very much interested in numerical integration I'd prefer to use this as a chance to implement several different ones. I've also got a symplectic Runge Kutta written up somewhere but I imagine it would drop the frame rate too much if used...Thanks in advance for any info or advice.

 

[gneill]

You might want to have a look here for some ideas about integration methods. The site author has done some nifty simulations, particularly for Klemperer Rosettes.

 

[D H]

I was wondering if there are any integrators (surely must be) that are built for this sort of problem and use keplers laws as a basis and adds perturbations? It would have to deal with the effects of extra planets that I randomly throw in so a pure keplerian integrator wouldn't do it as it wouldn't take into account any other body, just the one it orbits.

Yes, there are. This is essentially what NASA used back in the 1960s to calculate precision orbits of the Apollo missions (Encke-Nystrom). This approach used a baseline elliptical trajectory propagated via Kepler's laws and a deviation from this baseline with perturbations caused by other bodies. Occasionally the baseline had to be re-anchored when the difference between calculated and baseline grew too large.

Another approach is to use Lagrange's planetary equations (variants: Gauss' planetary equations, Hill's planetary equations, ...). These give the partial derivatives of some set of orbital elements with respect to disturbing forces. Which set of orbital elements? That's why there are some many different sets of planetary equations. Another factor that comes into play is whether the forces are conservative (gravitation) or non-conservative (drag).

Yet another approach that falls into this camp are those that use the KL (Karhunen-Loève) transform. You may have heard the phrase "Euler angles are evil." So are orbital elements. In fact, the longitude of ascending node, the inclination, and the argument of periapsis are a ZXZ Euler sequence. One way around using Euler angles to describe rotations is to use quaternions. Can quaternions be used as an alternative to orbital elements? The answer is yes. The KL transform uses a pair of quaternions to get around those Euler angle type problems that are associated with orbital elements.

One problem with all of these techniques is that the states being integrated are not cartesian position and velocity. For an N-body problem, this means there is a lot of back and forth between the state as modeled by the integrator and the cartesian positions needed to compute planet-to-planet interactions. Another problem is that it is mean anomaly (or something very similar) rather than time that is the independent variable. This again creates problems with the N-body problem.

That said, there is a lot of renewed interest in these techniques as of late. They are potentially very useful for solving the orbital debris problem. Here, body-to-body gravitational interactions are essentially non-existent. The problem separates nicely into a bunch of independent orbit propagators of objects that follow perturbed ellipses. Not so nicely, the amount of debris is growing, rapidly. Compounding this is the desire to track ever smaller chunks of debris. The old techniques can't keep up, even with the help of Moore's law.

The way I currently set the integration up is to calculate gravitational effects using the centre of mass of planet-moon system and only use the effect from individual moons when an object is within the planets hill sphere. Does this seem like a decent idea and will it be accurate enough for short term (<1000 year) simulations?

Consider the Earth's Moon. It is essential to model the perturbations by the Sun on the Moon's orbit if you want to have even short term accuracy for its orbit about the Earth. Some of the solar perturbations were known to the ancients (obviously not the cause, but the effect certainly was). If you want to have long term accuracy (100s of years) you also need to model perturbations by Jupiter and Venus.

You may also need to model moon-to-moon interactions. For example, Io, Europa, and Ganymede are locked in a 1:2:4 orbital resonance. Accounting for these inter-moon interactions is essential for long term accuracy. (This pertains to the perturbations by the large moons. The perturbations by smaller moons don't matter.)

Bottom line: A Hill sphere approach won't cut it if you want long term accuracy.

I currently have scaled units such that 1 unit of time is a day, 1 unit of distance is 1AU and 1 unit of mass is the mass of the earth. Does this seem reasonable or will it become problematic with precision when including smaller moons (since also radius of moons are scaled by 1AU, they can get pretty small).

A bigger concern is your timestep. You need to have r/(vΔt) and v/(aΔt) be the "right size". Too small and numerical errors will kill your accuracy. Too large and deficiencies inherent in the integration technique will kill your accuracy. What that "right size" is depends on your integration technique. You can use adaptive step size techniques to overcome this to some extent, but now you have a new problem.

That new problem is stiffness. This is a stiff problem. Adaptive step size techniques tend to have a hard time with stiff problems. In fact, this stiffness is a bit pervasive. You need multiple integration rates to get Mercury and Neptune right. Jupiter is a fairly big perturber of Mercury's orbit, so once again, stiffness rears it's ugly head. It rears its ugly head yet again with the moons of the gas giants.
One last thing: You mentioned you are using NASA's Horizons to obtain initial conditions. Just something to beware of: The velocities most likely will kill your long term accuracy. The underlying data are sets of polynomial coefficients (Chebyshev polynomials) that yield position as a function of time. The coefficients are tuned to yield accurate positions and to make position piecewise continuous across set boundaries. While velocity is an easy calculation (it's just the time derivative of the polynomial that yields position), the coefficients are not tuned to yield good velocities.

 

[Deadstar]

Hey D_H, late reply here. Thanks for your post.

Should probably clear up one or two things.

With regards to time step and units used. I suppose I was meaning in terms of finite computer precision and whether certain values 'work better' than others at not accumulating errors.

With regards to perturbations. I think I explained it wrong. Using any moon of Jupiter as an example... I would include perturbations from all other moons of Jupiter, the sun, and the centre of mass of all individual planet-moon systems. If I was say, flying a satellite through my solar system sim, I would calculate effects from all planet-moon centre of masses except if the satellite was inside the hill sphere of the planet and then I would use individual moons effects instead of the centre of mass approach. Sort of a barnes-hut tree approach using the planets hill sphere as a region.

I'm needing to properly learn about quaternions so I'll maybe use this as a chance.