A program for numerical integration of asteroid and comet orbits
Last updated 10 February 2009
Purpose of this program: The most commonly used datasets of asteroid orbital elements are the Minor Planet Center's MPCORB and the Lowell Observatory ASTORB files. These both provide elements for particular standard epochs, at 200-day intervals. So the epoch can be up to 100 days from the present date.
The MPC provides a companion comet.dat file with a different format, but a similar standard-epoch method.
I am posting an open-source program, integrat, to integrate any of these files to any desired epoch.
Usually, the "standard" epoch isn't a big deal. But it does mean that if you compute an asteroid or comet position for "right now" using the simplest two-body methods, you may be ignoring up to 100 days of perturbations. For main-belt objects, that usually results in errors of an arcsecond or so. For near-earth objects that have just recently passed us by, it may cause a few arcminutes or degrees of error. And if you want to find out where the asteroids were, say, fifty years ago (in hopes of locating them on old plates), you'd be completely out of luck; the current orbital elements would be just about useless.
One solution would be to have your software numerically integrate the orbits on demand. I've considered this for my own software, and will probably do it eventually... but the result will not be without limitations in terms of speed and ease of use; it will still be helpful to be able to integrate a large element set once to a desired epoch, then leave it be.
In my commercial software, I adopted the solution of taking the Lowell dataset and running it through an integration routine, generating elements at 200-day intervals covering a span of a few decades. Ask for asteroid elements at a given date/time, and the program figures out which element set is "closest" to that date/time. It's a nice solution, except that it leads to enormous files (even when stored in binary form). With the current number of asteroids, I'm stuck storing about 200 MBytes of compressed asteroid orbital data on the CD. As the number of asteroids goes up, I'll have to cut down the time span or toss something else off the CDs. And people are stuck with whatever data is available at the time the CDs are made.
Using this software: First, you will have to click here to download the software (about 152 KBytes). Create a directory for it and unZIP it, and you can then run INTEGRAT.EXE from a DOS shell. Do so without command-line arguments, and you'll get the following:
INTEGRAT takes as command-line arguments the name of an input file of the MPCORB.DAT or COMET.DAT type; the name of the output file that is to be created; and the epoch of that output file, as either a JD or in YYYYMMDD form. For example, either: integrat mpcorb.dat 2452600.mpc 2452600.5 integrat mpcorb.dat 2452600.mpc 20021122 would read in the 'mpcorb.dat' file, and create a new file updated to the epoch JD 2452600.5 = 22 Nov 2002. One can also specify the date as "today", or "today+7" (i.e., this day next week) or "today-30" (a month ago).
You can then swap in the new 2452600.mpc file in place of MPCORB.DAT in a variety of planetarium programs (including my own) and have that program run with elements in the revised epoch.
A warning: doing this integration can take a while. You'll get a running estimate of time elapsed/time to completion. If you're just taking the current MPCORB and bringing it up to date (that is, you're integrating over about a hundred days or less), it will take about 15 minutes on a 500-MHz machine. The time required is proportional to the time span, though: if you integrate back a few decades, be prepared for an overnight run.
If you are integrating one of the other mpcorb-like files, such as nea.dat, this problem will mostly vanish, simply because the other files contain only a small fraction of the orbits in mpcorb.dat.
The run time is also proportional to the number of asteroids in the input file. Chop out all but a few dozen asteroids (as happens with nea.dat, for example), and you can get perfectly reasonable run times over long time spans. (In such cases, most of the time is spent doing the initial computation of planetary positions; using JPL ephemerides can speed this up nicely.)
C/C++ source code: integrat.cpp and integrat.mak are included with the above ZIP file. You'll also need the basic astronomical functions library and the JPL ephemeris source code to get the whole thing compiled and running.
Using JPL ephemerides for the planetary positions: By default, INTEGRAT computes planetary/lunar positions using a truncated version of the VSOP87 method. This leads to a very small error, not normally significant. Its greatest drawback is that if you want to integrate orbits a long time back or forward, there is a long pause when the program starts while it computes all the positions required to cover that time span.
For both reasons, it can be helpful to switch from VSOP to the JPL ephemerides. To do this, you need to get JPL DE files; these can be downloaded over the Internet, or copied from the second Guide 8.0 CD-ROM, or purchased on CD-ROM from Willmann-Bell. All of these options are discussed at the above link.
Once you have a DE file, you can instruct INTEGRAT to use it with the '-f' switch. For example:
integrat mpcorb.dat 2452600.mpc 2452600.5 -fd:\unix.405
would tell INTEGRAT to use the DE file d:\unix.405 as a source for ephemerides.
Methods used in this program: In almost all cases, INTEGRAT uses a fixed-step Runge-Kutta-Fehlberg method. Using a fixed step has the advantage of allowing the program to set up a table of planetary positions once, rather than having to consume time in computing positions at every step. Also, the method of Encke is used: instead of integrating the orbit "directly", only the difference between the actual motion and elliptical, two-body motion is computed. This raised the complexity of the program a great deal, but means that a much larger step size can be used.
Using RKF means that it has a rough error estimate at each step. If the error estimate is large (usually meaning the asteroid has come close to a planet), the program can compute a smaller step size and try that instead. (Sometimes, it turns out that even the smaller step would lead to an excessive error, and the program has to try a still smaller step size.) Eventually, it covers the "problem step". With any luck, the next step will not be problematic and the program can return to using its pre-computed table of planetary positions. In any case, the error tolerances and step sizes are set up such that "problem steps" won't happen very often. (When they do happen, they'll be handled accurately... they just won't be as fast as "normal" steps.)
In addition to the above RKF method, it would help if the program had a multistep method (probably Gauss-Jackson). Suppose we use an incarnation that requires data from four previous steps; in such a case, we would use RKF to generate four steps, then be off and running with Gauss-Jackson. From time to time, Gauss-Jackson would lead to excessive error; in such cases, the step would be taken with RKF instead. All of this would let us use Gauss-Jackson most of the time (becoming much faster), but still rest secure in the knowledge that, if the asteroid runs near a planet and an error buildup is detected, we can switch to a method that will retain accuracy.
(10 February 2009) The current version can handle epochs that aren't exactly on 0 UT. This required a bit of fudging, though, because the mpcorb.dat format stores epochs with one-day precision; you can't specify fractions of a day. So for some time, I told people looking for an epoch such as 2008 January 30.46 that they were out of luck.
I have "solved" this as follows. If you ask Integrat to integrate to 2008 January 30.46, it will produce elements that give the exact state vector (position and velocity) for that instant. The epoch will be rounded (to 2008 January 30.46) and the mean anomaly adjusted for .46 days of motion.
Hence, you can load the elements into almost any desktop planetarium software and get exact results for 30.46 January. A truly elegant solution would involve scrapping the mpcorb.dat format in favor of one that allows you to specify epochs exactly, instead of to the nearest day (and one that would handle hyperbolic and parabolic orbits, and not lose precision for nearly-parabolic orbits, and store uncertainty values and a covariance matrix, and some more digits, and planetocentric orbits as well as heliocentric ones).
(10 February 2009) There was a bug when integrating the orbit of Pluto when using JPL ephemerides. With JPL ephemerides, Pluto is included as a perturber, and the program would attempt to include Pluto's perturbations on itself! The program would grind to a halt. This problem is now fixed: if the program sees (134340) Pluto in the input mpcorb.dat file, it temporarily disables Pluto as a perturber.
(10 February 2009) The date/time entry is more flexible now; one can use any of the time formats described here. For example,
integrat mpcorb.dat 30jan.dat 20080130.46 integrat mpcorb.dat 30jan.dat 2008 jan 30.46 integrat mpcorb.dat 30jan.dat 2454495.96 integrat mpcorb.dat 30jan.dat 30/1/2008 11:02:24
would all have the same result. (Previously available time specifiers, such as today or today-3, will still work.)
(3 November 2005) The current version now can integrate comet.dat files, as well as mpcorb.dat files.
(3 May 2004) Recently, there have been problems with the NEA orbits provided in mpcorb format on the MPC and mirror ftp sites. Specifically, neatod.dat, which is supposed to provide NEA orbits with today's epoch, has often been missing some recent objects.
To get around that, people have been downloading nea.dat, with its epoch that can be a hundred days from "today", and integrating it to today's epoch to create their own neatod.dat. For example:
integrat nea.dat neatod.dat 20040503
I've simplified this process by having integrat interpret the solitary command-line argument "today" to mean, "integrate nea.dat to today's epoch, storing the result as neato.dat." That is, you need simply run
to get the desired result. Also, one can use "today" in place of the epoch's Julian Day, or the YYYYMMDD, and can use expressions such as "today-7" (days) or "today-30" to get files with an epoch of a week or month ago.
And finally, integrat now adds a few lines to the header of the output file, mentioning that the file has been integrated from such-and-such starting date to its current epoch, and the version date for integrat. Not a big deal at all, but it does mean that if you've forgotten where a given mpcorb-formatted file came from or what steps have been performed upon it, you can find out by checking its header. (I may be the only person who has a need to do that, though!)
(6 February 2004) Posted an update in reply to an inquiry from Cristovao Jacques. The newly-posted version does a better job of handling asteroid perturbations. Previously, integrat assumed that the first, second, and fourth asteroids in mpcorb.dat were (1) Ceres, (2) Pallas, and (4) Vesta, the perturbers considered when integrating asteroid orbits. If you're using the "usual" mpcorb.dat file, this is true; but the files of NEOs and other objects provided by MPC, as well as the mpcorb.dat-format files produced by Find_Orb's Monte Carlo routine, do not have these three objects. integrat now looks for the actual object names, instead of making wrongheaded conclusions based on where an asteroid might be in the file.
Known limitations: Reiner Stoss pointed out that the program assumes that the input objects all have the same epoch. When I wrote the program, I realized that it would be a lot easier to do if every object was integrated over the same time span. I asked Gareth Williams at the MPC, and was told that every perturbed orbit in mpcorb.dat would have a consistent epoch. So I didn't worry about this limitation.
However, Reiner suggested that it would be nice if the unperturbed objects also got integrated. True, this can be a little dodgy, but ought to at least result in an improvement over leaving them in their original epochs. At some point, I may make this change (especially since I now have some code, written in support of Find_Orb, that makes it easier to do.)
Future/planned improvements: The main drawback of this program is that it's not user-friendly. I'll probably focus on that aspect of it before doing anything else.
Once it's user-friendly, with a nice Windows interface, I may work on adding a multi-step integration method. This will cause no difference in the operation of the program or in its accuracy, but will make it run a lot faster (I think.)