Accuracy test of the Orbiter integration engine

[Keithth G]

OK, a quick summary and update:

So far, I've gone through an exercise of building an off-line numerical integrator with which to test the accuracy of Orbiter's integration engine. To make a like-for-like comparison, I have:

a. adjusted Orbiter's configuration files so that the Solar System has only two gravitating bodies (the Sun and Mars);

b. built interpolation functions based on Orbiter output so that my integrator has, to within a few metres, the Sun and Mars following the same trajectories as used by Orbiter;

c. checked to make sure that I was using the same masses (and other parameters) as Orbiter.

In principle, then, both Orbiter's integration engine and mine should integrate a trajectory in much the same way. To test this, I worked out the initial conditions needed to send a ship on a near parabolic orbit around the Sun on a 111 day and 4 AU journey from a point on one side of the Solar System near Earth's orbit to the other side and a (high-speed) rendezvous with the centre of Phobos. I worked out that my firing solution should be accurate to within a metre or so, and yet, when the same initial conditions were run through Orbiter's integration engine, the ship arrived at Phobos around 1.7 seconds too late (and, consequently, about 8.1 km 'off target'). This discrepancy, though small, had me stumped. Basically, there didn't seem to be any control parameters left to explain the difference.

[Just to put this in perspective, and to highlight how pedantic this is, one can liken the exercise to lobbing a stone across a garden and fretting anxiously that the stone hit the ground 'off-target' by the width of a few atoms. Not the sort of stuff likely to cause empires to crumble.]

Bu, anyway, a number of helpful comments were forthcoming but one in particular from 'martins' (a.k.a., presumably, Martin Schweiger, author of Orbiter) relating to the fact that Orbiter was a bit picky about which contributions to the gravitational fields applied to trajectory determination at any point in time. As it turned out, I quickly worked out that the contribution to the gravitational field from Mars 'switched' on at some point between MJD 51949.7 and MJD 51949.8 - roughly 32.5 days before arrival at Phobos - i.e. when one is still a very, very long way from Mars.

Now, this 'switching on' of the gravitation field isn't, strictly speaking, the way that Newtonian gravity works. Normally one would sum the components of the gravitational field from all of the gravitating bodies (in this case, the Sun and Mars) no matter how far away one is. Notionally, then, we are all affected by the slight gravitational tug of Jupiter and Saturn, say. So, why throw away the contribution from Mars if one is a far away from it? Well, basically for the same reason that we don't include the gravitation contributions from Alpha Centauri in our calculations: the magnitude of those contributions is generally deemed to be too small to matter. But calculating those contributions is computationally costly so that to ensure that Orbiter doesn't spend an inordinate amount of time calculating gravitational contributions of little consequence, the very long tails of the gravitational potentials is simply cut off. This simplifies things enormously and allows Orbiter to devote more resources to displaying breathtaking imagery without much by way of a significant deviation from real physics.

And generally this is all very reasonable - except, that is, when one is attempting to test the accuracy of Orbiter's integration engine using precision integrators. In this case, and rather an exceptional case at that, Orbiter's treatment of long-range gravitational forces do matter and one needs to take account of the gravitational cut-off if one is to explain the 1.7 second delay.

So, then, what is the explanation for the delay. Well, when calculating the initial conditions for my trajectory using my integrator, I assumed that the tail of Mars' gravitation field was there. To counter this additional tug from Mars, the initial conditions included a very, very small amount of additional 'outward' velocity. However, when run through Orbiter, this pull from Mars wasn't there so that the additional 'outward' velocity sent the ship on a slightly longer (and slower) trajectory to Mars - hence the 1.7 second delay in the Phobos 'rendezvous'.

An interesting theory, but how do things work out in practice. Well, in my integrator I can add rule that says that if the date is before MJD 51949.75 the trajectory should be calculated as if Mars was not there; and after that date, the trajectory should include Martian gravitational field. This should bring my model of the gravitational field encountered by the ship (a.k.a. 'test particle') into much closer alignment with that used by Orbiter.

And now, we should be comparing apples with apples.

Using my integrator, I calculate a revised set of initial conditions based upon this modified gravitational potential and then run these initial conditions through Orbiter and, sure enough, the ship arrives at the centre of Phobos to within about 0.01 seconds of the scheduled arrival date; and now just 300m from the notional centre of Phobos. In one fell swoop, all but 4% of the anomaly has been explained. Mystery solved.

What of the residual 300m. Let's put this in context. We are now talking about an anomaly that is approximately 5 x 10^-10 of the total arc length of the trajectory. It isn't much - and probably too small now to worry much about (even for me). But I can hazard a guess as to where the bulk of the residual difference lies: i am using a very slightly different interpolated trajectory for Mars. The difference is just a few metres here and there - but over time, these differences add up. There may also be some issues relating to my rather cavalier treatment of fixed time-steps, but I would put my money on the interpolation issue being the principal cause.

Anyway, having lobbed my stone across the garden and having had it safely land on target (to within the width of an atomic nuclei at any road), I'm happy.

 

[Jumper_001]

I really enjoyed this thread. One interest I have had lately is if I learned to use GMAT well enough, could I use it as a navigational tool for planning missions in Orbiter, and get the data "back into" one of the simple Delta Velocity programs. I was trying to fly BrianJ's excellent LADEE add-on "by the numbers" but planning the 3 TCMs for the phasing loops is difficult with either IMFD or TransX. So this comparison of Orbiter's integration to something else is interesting as I was thinking of doing something similar with GMAT.

And I agree with the poster above that I have had years of enjoyment out of this "applied math" project. It was a great inspiration to me while I was studying for my degree. I am lucky that some hardware that I have worked on since has made it's way to Mars for real - but I made the journey dozens of times in Orbiter before I made it that far in my career. And I vaguely recall that flying to the outer planets when the time acceleration only went to 10,000x (was that Orbiter 2003 or 2005?) took several days of leaving Orbiter running 24/7!

 

[Keithth G]

Hi, Jumper_001

A few comments about GMAT.

In many ways, I think, GMAT does many of the same things that Orbiter does, but GMAT talks High German, whereas Orbiter talks colloquial French. There is a translation issue that needs to be addressed but the basic idea of doing 'offline' mission planning in GMAT and then executing that plan in Orbiter is, in principle, plausible - albeit one that is likely to entail a few complications.

Here are some of the issues that you will need to address:


Time measurement: Because Orbiter uses the popular VSOP87 ephemeris, its measure of time is a thing known as Barycentric Dynamical Time (TDB). GMAT, being a Nasa sponsored widget, uses a JPL ephemeris. I'm not quite sure which ephemeris that is, but it is probably the successor to the DE405/6 ephemeris. This ephemeris' basic measure of time is a relativistic coordinate time, known as Ephemeris Time (ET) which differs from the commonly used Terrestrial Time (TT) which differs from TDB. Sounds complicated? It can be. Getting your measure of time right can be an ordeal in itself.

Reference Frame: The VSOP87 ephemeris uses a coordinate system that isn't quite inertial. If I remember correctly, it rotates slowly with respect to FK5 and ICRF inertial reference frames. One or other of the latter inertial reference frames is probably used by GMAT - I don't know which - so you also need to make sure that you rotate coordinates from FK5/IRCF to the dynamical reference frame used by VSOP87. Another little nuisance with which to complicate life.

Ephemeris: In Orbiter, as in GMAT, the centre of masses of each of the planets and their satellites run on rails. These rails are determined by an ephemeris engine. For Orbiter, the ephemeris of choice is VSOP87. This is a semi-analytical theory developed by the Bureau de Longitudes. It is now a little outdated but remains popular due, primarily, to its compactness. More modern ephemeris are those produced by JPL - such as DE405/6. These are based on direct numerical co-integration of all of the planets (given a set of initial conditions). Although close, the two ephemeris have the planets running on slightly different tracks. And, as seen, small differences in planetary positions can add up over time.

Gravitational model: as 'martins' has advised earlier in this thread, Orbiter chops the tails of the gravitational fields of some of the planets - depending upon the location of the spacecraft in the gravitational potential. I'm pretty sure that GMAT does not. In addition, GMAT may take into account a number of other gravitational perturbations that it would be difficult to replicate in Orbiter.


So, it can be done; it will entail complications; and it is likely that some differences will remain that may be significant depending on the level of accuracy that you require from your mission planning.


The other way forward that achieve much the same thing is to reverse-engineer Orbiter's integration engine and use that knowledge to build an offline mission planning tool that very closely replicates the behaviour of Orbiter's own integration engine. ("It's physics, Jim, but not quite as we know it.") On the whole, I suspect that this is likely to lead to more satisfactory results - assuming, that is, that one wishes to execute mission plans exclusively within in Orbiter rather than in the 'real world'.

 

[martins]

Excellent news about finding the cause for the bulk of the discrepancy. I'll add an option in the next version to allow users to set the threshold for dropping gravitational sources from the update calculation (or disabling that cutoff altogether). For a not too densely populated solar system the performance penalty may not be very high, but in cases with both a large number of gravity sources M and dynamically updated vessels N, the computational cost will be O(M*N), whereas at the moment it's probably more like O(sqrt(M)*N).

 

[Keithth G]

I think adding "an option in the next version to allow users to set the threshold for dropping gravitational sources from the update calculation (or disabling that cutoff altogether)" is a good idea. This should keep the Puritans (of which I may be one) happy.

Having said that, I can fully appreciate the performance reasons for setting a threshold in the first place: computing power is not infinite and there is a non-negligible cost for computing all of the Newtonian forces.

Aside from issues of 'purity', from a mission planning perspective the issue is not so much that all of the forces from all of the gravitation bodies is calculated all of the time. Rather, the issue is that the 'rules of the game' are understood. Mission planning is really an exercise in emulation - i.e., "if I execute a series of manoeuvres in my planning tool, will Orbiter (or Nature) execute those manoeuvres in the same way." So long as one understands the principles used by Orbiter to drop gravitational sources from its force calculations, it is then quite straightforward to emulate that behaviour outside of Orbiter. This is basically what I did - albeit in 'hot-wired' form but to good effect - in my simple test. Another way forward, then, and one that may be equally useful, is to make those rules more transparent.

On another issue: in an earlier comment relating to making public integrators you mentioned that "this would come in handy in various places, e.g. for mission planning tools, or for creating our own power series ephemeris solutions for celestial bodies (I'm still looking for a solution for Hyperion)." What by way of ephemeris solutions are you looking for?

 

[martins]

On another issue: in an earlier comment relating to making public integrators you mentioned that "this would come in handy in various places, e.g. for mission planning tools, or for creating our own power series ephemeris solutions for celestial bodies (I'm still looking for a solution for Hyperion)." What by way of ephemeris solutions are you looking for?

Essentially I was thinking about a machinery that can generate a series expansion of a perturbed orbit from a numerical solution over a reasonably long period centered at the present - essentially what VSOP does for the major planets. In the case of Hyperion, this would be referenced to the true position of Saturn, or barycenter position of Saturn and its moons (and ideally in the same frame of reference as VSOP, to avoid the frame transformations you mentioned earlier).

As far as I know, Hyperion's may be a particularly difficult orbit to model, because of resonances with Titan and high eccentricity, so a series expansion diverges either very quickly or requires a large number of terms. But since you have already created a precise numerical gravity simulator, it might be interesting to use it to see just how well it can be approximated with a series solution.

 

[Keithth G]

I will have a think about that.

 

[boogabooga]

IMHO, an error of ~2 seconds in 111 days and 8km in 4AU is already a success. But I get that it was an exercise in mathematical pureism.

From a practical point that I have made before, Orbiter is not necessarily for professional planning, and small performance trade-offs seem quite logical. The types of errors you are talking about are outside of our abilities to actually execute at those tolerances using existing MFDs.

Even a real mission to Phobos is not so precise at first and would include several MCCs.

Orbiter planning in GMAT is a compelling idea. However, I don't think that the complications that Keithth G lists will be so severe for missions in the near future/past; again, errors on the order of seconds per month should not be so critical in Orbiter.

The propagators in GMAT are user defined from a buffet of options. You can select as many gravitational influences as you wish from the planets, the Moon, and the Sun. However, even in official tutorials they don't turn on everything all at one. They usually suggest a zonal approach with for example Earth, Sun, Luna (Moon) while under Earth SOI, then a deep-space propagator with inner planets and perhaps Jupiter, etc.

GMAT is quite scalable in complexity. For example, you could approximate maneuvers as impulsive, if you wish.

As a practical matter, I would worry more about GMAT's different coordinate systems. I've found the default GMAT "VNB" system to be quite similar to the IMFD "P30 LVLH" mode (albeit with the orbit normal directions reversed). I believe that this is also similar to TransX's coordinate system. However, AFAIK GMAT doesn't replicate IMFD's "Vel. Frame" mode.

 

[Keithth G]

Essentially I was thinking about a machinery that can generate a series expansion of a perturbed orbit from a numerical solution over a reasonably long period centered at the present - essentially what VSOP does for the major planets. In the case of Hyperion, this would be referenced to the true position of Saturn, or barycenter position of Saturn and its moons (and ideally in the same frame of reference as VSOP, to avoid the frame transformations you mentioned earlier).

As far as I know, Hyperion's may be a particularly difficult orbit to model, because of resonances with Titan and high eccentricity, so a series expansion diverges either very quickly or requires a large number of terms. But since you have already created a precise numerical gravity simulator, it might be interesting to use it to see just how well it can be approximated with a series solution.

The machinery to do this exists. I will do a 'proof of concept' exercise to see how well it works out in practice. It may take a while before I can report back.

 

[meson800]

This might be better in a new thread, but it seems related:

How does one go about building the numerical integrator? I understand how to numerically integrate using Simpson's rule and friends, so I imagine it would go something like the following:

1. Calculate total acceleration due to gravity using gravity equations and locations of ship+planets
2. Use Simpson's rule to calculate dV for a certain dT
3. Using that dV, calculate a dX(and dY and dZ) for that dT
4. Do it again!

Looking at the Runge-Kutta methods, it looks an awful lot like Simpson's rule applied to differential equations. Does the higher-order Runge-Kutta simply integrate using a cubic/quintic/etc function analogue instead of the quadratic Simpson's uses?

 

[Keithth G]

In the broad, you've got it right.

Range-Kutta methods are a general class of methods for integrating systems of ordinary differential equations. There may be others on this forum that know more about this class of integrator, but the all of the coefficients and weights are carefully chosen to ensure that error terms up to a particular order (e.g., 4, 6 or 8) are exactly zero. There isn't, I suspect, an assumption that between the start and end of each time-step that the object of study follows a path of a certain polynomial order - although it may be possible to interpret a Runge-Kutta integration scheme in that way.

However, for many physical systems, Runge-Kutta methods tend to have the disturbing property that the total energy of the system isn't bound. For example, if you take the equations of, say, a simple pendulum and then integrate them using straightforward Runge-Kutta techniques, the amplitude of the swing of the pendulum will either increase without limit, or tend to zero. In the absence of friction of other dissipative forces, this isn't what is supposed to happen and so physicists tend to frown a little on this kind of mathematical behaviour.

This is where symplectic integrators come in. Physical systems generally have a property of being symplectic. No need to worry too much about what it is except to note that it exists and physical systems have it. If one chooses the weights of the Runge-Kutta so as to kill all terms up to a certain order AND preserve the symplectic character of the physical system THEN the energy of the system tends to be bound. For our simple pendulum this means that once set in motion with a particular amplitude then, within some tolerance, the amplitude remains the same pretty much indefinitely thereafter. In other words, they tend to be much more stable and more representative of the true behaviour of the physical system being modelled. Physicists tend to like this.

But then there other kinds of integrators that are also symplectic integrators - there are, for example, variational integrators (based on Gauss-Lobatto quadrature which do model the path of the 'object' as following a polynomial path of some specified order). These are also symplectic. And then there are integrators based on splitting of the Hamiltonian (the function used to set up the dynamics of the physical system). Again, these are usually constructed so as to be symplectic.

In short-term, a symplectic integrator is generally deemed a 'good thing'. Being a Runge-Kutta integrator is fine, but for integrators probably not as important as it being symplectic.

Rather than go into too much detail right now, what I think I will do is put together a short note on constructing some symplectic integrators and post it as a separate thread.

 

[samplereturn]

2015 PDC and Orbiter accuracy

Thank you for the interesting and timely thread about the precision of ORBITER's integration engine. Keithth G's post #6 answered my question of whether ORBITER's engine uses Sun-Centered or Solar System Barycenter elements/coordinates. I tried a much more humble experiment of loading asteroid data from JPL's HORIZONS database into a scenario and was incredibly impressed with how closely ORBITER'S results agreed.

Does the community know about the IAA Planetary Defense Conference 2015 that met just last week, or the "2015 PDC Hypothetical Asteroid Impact Scenario" discussed there?

Go To --> http://neo.jpl.nasa.gov/pdc15

In addition to the interesting online app to play with (toward bottom of that page), orbital data for the fictitious asteroid "2015 PDC" is available in the HORIZONS database, just like any other real moon, comet, asteroid or listed spacecraft. I've loaded asteroid elements from this before and have flown to them, but this exercise offered a test: loaded into ORBITER, would 2015 PDC's encounter with Earth accurately match the scenario described at the conference ?

To make a long story short, I requested one set of osculating elements for 10 Sept 2021 (almost a year before impact predicted on 3 Sept 2022), used Sun Centered coordinates (assuming ORBITER would make the correction to barycenter internally), made sure that the MJD of ORBITER corresponded to the Coordinate Time output of HORIZONS, modeled the asteroid as a massive vessel so it would be perturbed by all the major planets out to Saturn, included the light-pressure parameter and ran the scenario mostly at 100,000 time compression.

Success! Earth encounter was within the conference's (whole-asteroid, Day 4) predicted impact footprint in the South China Sea; in my case: north of the western end of Palawan Island, west of Mindoro Island, of the Philippines. My 2015 PDC asteroid entered the atmosphere (200 km alt, DNP = 0.03 Pa) on 3 Sept 2022 at 03:52:09 UT. Conference impact prediction (fragment, Day 5): about 03:50 UTC.

I Love ORBITER! Thank you, Martin! Thank you, everybody!