Orbit determination from observations

Last updated 2014 August 8

Italian version of this page

French version of this page

Find_Orb users group

There is an e-mail group for Find_Orb users at http://groups.yahoo.com/group/find_orb/. Announcements of new versions will be made there; questions about the program or about orbit determination in general, or requests for new features, are welcome.

What is Find_Orb?

Find_Orb can take a set of observations of an asteroid or comet, given in the MPC (Minor Planet Center) format or the NEODyS or AstDyS formats, and find the corresponding orbit. The MPC 80-column format is a "standard" used by most astrometry software. You can click here for details of the MPC 80-column format.

Find_Orb can determine orbits of artificial Earth satellites and for satellites of other planets. It exists as a Windows GUI program, and as C/C++ source code that can be compiled to make Windows or Linux or Mac console applications (one interactive, one non-interactive). In what follows, I'll just refer to "Find_Orb", without making a distinction between the various flavors.

Click here to download the Windows version (about 660 KBytes).

Click here to download the C/C++ source code (about 240 KBytes).

Find_Orb is a user-friendly program that handles the initial determination of the orbit using the methods of Gauss, Herget, and Väisälä. Given some more observations, it can find a "best fit" orbit using the method of least squares . The physical model includes most significant effects on comets and asteroids.

Why does Find_Orb exist?

There are several situations wherein one might want to use Find_Orb:

• If you're tracking a near-earth object (NEO), you may not be able to wait for a revised orbit from the MPC. With Find_Orb, you can take any available astrometry, then generate an orbit and/or ephemeris. At least in some cases, if you had to wait for MPC to update the orbit, you could lose the object. (This is increasingly rare; MPC has done a lot of work in this area.)

• You can use Find_Orb to test the quality of astrometric observations. If you feed your astrometry through Find_Orb, along with some comparison data from others, and your positions line up nicely with theirs, resulting in an orbit with residuals of around an arcsecond, you can be reasonably confident about the quality of your astrometry. So astrometrists in general may have an interest in the software.

• Inquiring Minds are sometimes curious about the entire "black art" of turning astrometric data into orbital elements. This was my main motivation for writing Find_Orb. (And it's turned out to be useful for instructional purposes; I occasionally hear from people using Find_Orb as part of their coursework.) You can find out what effect Jupiter has on the orbit of an object of interest, or how certain observations worked to change the knowledge of the orbit.

• It can be interesting, just as an educational exercise, to take your own astrometry and generate an orbit from it. If you've set up a telescope, gathered astrometry for an object, and then find its orbit, you have essentially done all of the steps involved in asteroid discovery.

• Occasionally, I hear from people doing more esoteric things with Find_Orb, such as attempting to link the orbit of a comet with possible archived data from the past. Some people have used it to do impact probabilities for objects such as the Shoemaker-Levy 9 fragments that hit Jupiter in 1994, or the object 2007 WD5 that briefly seemed likely to hit Mars. Most notably, my announcement that 2008 TC3 was going to hit the earth was based on feeding the data through Find_Orb, and some people were able to use it to determine the impact point and time. (Since 2008, the code has been changed to make it quite easy to compute short-term impacts of this sort.)

Licensing

Find_Orb is distributed under the GPL (GNU General Public License), version 2, which basically means you can do what you want with it (and with the source code), and can use/distribute it without having to pay for it. Use in or with commercial software is forbidden without written permission from Project Pluto:

Project Pluto
168 Ridge Road
Bowdoinham ME 04008
(UNITED STATES)
tel (207) 666 5750
plutó@projectplutô.cóm
http://www.projectpluto.com

If you use it for something interesting, though, please let me know!

As is required by the GPL, source code is provided.

Other orbit determination software of interest

Perhaps the most widely-used publicly available software comparable to Find_Orb, allowing you to include perturbing objects and to do both initial and full least square orbit fitting, is the OrbFit software. It's also free software (with restrictions similar to those for Find_Orb). It was put together by about half a dozen professional astronomers, and can run in Windows or Unix. FORTRAN source code is provided. (You can also click here for information about the C/C++ code for Find_Orb.)

OrbFit also includes code to do some quite advanced tasks, such as figuring out proper orbital elements, impact probabilities, use radar observations, compute error values for almost all parameters, etc. In most ways, it's much more sophisticated than Find_Orb. (Which sometimes translates into "not as easy to use".)

Pasquale Tricarico has informed me of the development of the ORSA (Orbit Reconstruction, Simulation and Analysis) project. This will provide a core library and a set of tools for orbital determination tasks.

Jim Baer's CODES software (http://home.earthlink.net/~jimbaer1/), "performs comet/asteroid initial and least-squares orbit determination, collision analysis, ephemeris calculation, and object identification. CODES utilizes integrated n-body mechanics that account for cometary thrusting (if applicable) and gravitational perturbations (including relativistic effects) from the Sun, nine planets, Earth's Moon, and up to 300 asteroids." It also does some things in the area of computing covariance matrices and sigmas that are not (yet, anyway) implemented in Find_Orb. Since it's written in Java, it's somewhat more cross-platform than Find_Orb.

Aldo Vitagliano's SOLEX package includes an orbit-determination program, FindOrb (same as my software but without the underscore). SOLEX also has some impressive capabilities (not available anywhere else, as far as I know) for integrating the entire solar system forward/backward in time. So it can compute where planets were tens of thousands of years ago or from now. (Most software, my own included, is accurate over only a few thousand years.)

John Rogers' CAA (Computer-Aided Astrometry) software has an orbit-determination part, including the ability to determine ephemeris uncertainty via the Monte Carlo method. Unfortunately, it appears to be no longer available.

If you know of other orbit determination software not listed here, please e-mail me (after removing the accent from the 'ó' in the e-mail address!)

Where does one get observations to feed to Find_Orb?

Ideally, one should get a CCD camera and telescope, gather images and positions for some asteroids, and report those positions to the Minor Planet Center. As a side benefit, one can use those positions in Find_Orb.

Originally, I tested Find_Orb using observations reported in Minor Planet Electronic Circulars (MPECs). The observations are reported in MPC's 80-column format, of course, so you can feed the MPECs directly into Find_Orb, without even bothering to eliminate the other text associated with an MPEC; Find_Orb will recognize which lines contain observations and which contain other data.

Also, one can get observations for most asteroids from NEODyS and AstDyS, in their .rwo format. This format contains all the data that can be found in the standard MPC format, except (regrettably) the five-character reference code. It also includes uncertainties for both positions and magnitudes, "bias" data allowing for a rough correction of systematic biases in star catalogs, and it specifies whether the observation was included in NEODyS or AstDyS' orbit solution.

You can use the MPC's MPCOBS service to get observations for any given object. (This used to be a subscription-only service, but now, anyone can use it.) Or you can click here for access to MPC's databases of astrometric data. Note that the files with data for all numbered or unnumbered minor planets are enormous, and Find_Orb can't actually load all of them; you need to edit the files and extract observations of interest.

Data for artificial satellites is tougher to come by, but the MPC's DASO (Distant Artificial Satellite Observation) Circulars provide some data for these objects. (Note that a link to a file containing all astrometry from all DASOs is provided at the bottom of the page.)

A file containing "observations" for two totally imaginary asteroids, a (real) satellite of Jupiter, and simulated observations of the Mir space station is supplied with Find_Orb. The asteroid cases are the most practical; they ensure that you will have two "known" cases to work with before trying your own observations.

Starting Find_Orb

For Windows, you need only download Find_Orb from this page, put it in a folder of your choosing (usually a new one), and unZIP it. (Linux and OS/X users will have to download the source code and compile the program.) Windows users can then add a desktop shortcut to find_o32.exe; you'll see a rather ugly icon showing a yellow sun with an object orbiting it.

When you start it up, you'll get a big dialog box. There is an "Open..." button in the upper left corner; click on it and select the example file from the Open File dialog box.

Right below that "Open..." button is a small list box that will show all the objects Find_Orb found in that file. In the case of this file, it will find observations for two imaginary objects, 1996 XX1 and 1997 ZZ99, the Mir space station, and Sinope (Jupiter IX, a faint outer satellite of Jupiter).

Running Find_Orb with 1996 XX1

We'll try the easy one first. Double-click on 1996 XX1, and Find_Orb will immediately find the following pretty good orbit:

Orbital elements:
1996 XX1
   Perihelion 1998 Mar 8.969361 TT = 23:15:52 (JD 2450881.469361)
Epoch 1997 Oct 23.0 TT = JDT 2450744.5   Earth MOID: 0.0615   Ma: 0.0049
M 324.90486              (2000.0)            P               Q
n   0.25622617     Peri.   72.47300     -0.50277705     -0.86441415
a   2.45501256     Node    47.71103      0.79213685     -0.46159042
e   0.5741561      Incl.    0.14296      0.34602664     -0.19930493
P   3.85           H   15.8     G   0.15   q 1.04545211  Q 3.86457302
From 13 observations 1997 Oct. 12-22;   RMS error 0.485 arcseconds

In this instance (and in most instances), Find_Orb can get something close to the right orbit without you having to adjust it. The residual errors are in the proper ballpark (under an arcsecond or so.) The orbit has a large enough MOID with the Earth so that it's not a danger to us, though Martians might worry about the Mars MOID of .0049 AU. Again, you'll usually get an answer that indicates a sensible, non-threatening object.

In this case, you are basically ready to go; you could compute an ephemeris or feed the orbital elements to another program without further work. This won't be the case for shorter arcs, because the orbit won't be very well determined. It also won't be the case for longer arcs or objects coming near a planet (usually the Earth), where perturbations are a factor. It also wouldn't be the case if there are some erroneous observations that ought to be excluded from the orbit fit. An example of the perturbation problem follows.

Running Find_Orb with example asteroid 1997 ZZ99

The second "example" asteroid, 1997 ZZ99, is mostly similar, but it does throw some problems into the system.

Double-click on 1997 ZZ99, and right away, you'll see the following unusual idea of an "orbit".

Orbital elements:
1997 ZZ99
   Perigee 1997 Apr 22.935609 TT = 22:27:16 (JD 2450561.435609)
Epoch 1997 Apr 23.0 TT = JDT 2450561.5                  Find_Orb
q 14941.6690km           (2000.0)            P               Q
H   22.9  G 0.15   Peri.  356.23652      0.62657228     -0.76960718
                   Node   312.51135     -0.71062719     -0.49939125
e   5.8196325      Incl.    9.59986     -0.32002525     -0.39788586
From 35 observations 1997 Apr. 21-22;   RMS error 1.128 arcseconds

In this case, Find_Orb has automatically latched on to the correct orbit, and has detected that the earth and moon are significant perturbing forces. The check-boxes for both objects in the "Perturbers" section are turned on. Note that Find_Orb isn't always bright enough (yet) to detect this sort of situation.

Also, you'll notice that the perigee distance is comfortably greater than the 6378 kilometer radius of the earth. So we won't get hit this time.

Find_Orb has detected that at the epoch, 1997 April 23, the object is so close to the earth that it makes more sense to talk about its orbit relative to us, rather than its orbit relative to the sun. That is to say, the object is within our "sphere of influence". If you'd rather see the heliocentric orbit, you can use the Settings... Heliocentric Orbits Only check-box. Or, you can set the epoch to be a few days ahead of the encounter on April 23, or a few days after it. (The elements will be quite different, since the asteroid is obviously seriously perturbed by the encounter.)

Somewhat oddly, this object actually comes within the very small sphere of influence of the moon, and if you set the epoch to April 22.5 and do a Full Step, Find_Orb will show you a selenocentric orbit. (1997 ZZ99 is a carefully contrived example. In reality, the only cases for which I've seen selenocentric orbits have been artificial satellites such as LCROSS and STEREO-A and STEREO-B, which were specifically targeted to impact or pass the moon.)

Running Find_Orb with S/2000 S 11, a satellite of Saturn

At one time, I was at somewhat of a loss for astrometry to provide for irregular satellites of gas giants. Fortunately, this is no longer the case. The ideas discussed below will be quite general, but the example used will be the satellite S/2000 S 11, the discovery of which was announced on 19 December 2000.

Click here to download the astrometry for this object. This file is the Minor Planet Electronic Circular (MPEC) giving the original discovery astrometry. Load it into Find_Orb, and you'll get a Centaur orbit with not-very-good residuals by default. To find a Saturnicentric orbit, Find_Orb will need a little coaching.

Find_Orb is bright enough to realize that its default orbit puts the object near to Saturn, and therefore sets the "Saturn" checkbox in the Perturbers group.

The data arc spans 9 November 2000 to 17 December 2000. Consulting with an ephemeris program (or with the JPL Horizons system) tells you that, on those dates, Saturn was about 8.148 and 8.257 AU from the Earth. Since a natural satellite basically is just dragged around by a planet, you can be confident that the distance to the satellite was about (8.148 + delta) and (8.257 + delta) AU on those dates (assuming only a fraction of an orbit was completed).

If you look at the orbits of known irregular satellites of Saturn, you'll see that they have semimajor axes in the ballpark of a=.1 AU. I usually use initial guesses, therefore, of delta = -.1 AU (a guess assuming the object is on the near side of the planet) and delta = .1 AU (assuming the far side of the planet). Going with the first guess, set R1 = 8.048 and R2 = 8.157, click on "Herget Step" a few times, then on "Full Step" a few times. Rather quickly, the program will settle down on an orbit resembling the following:

Orbital elements:
S/2000 S 11
   Perisaturn 2001 Jul 6.675722 TT
Epoch 2000 Dec  2.0 TT = JDT 2451880.5
M 272.38363              (2000.0)            P               Q
n   0.40436636     Peri.    5.58450     -0.36000669     -0.07200746
a   0.1193108      Node   110.56628      0.81289929     -0.51345927
e   0.7685389      Incl.   83.45621      0.45780992      0.85508743
P   2.44           H   11.0           G   0.15      q 0.0276158
From 15 observations 2000 Nov. 9-Dec. 17;   RMS error 0.552 arcseconds

Well, this has reasonably low residuals, but the inclination and eccentricity look a little bizarre. In any case, let's now try the alternative assumption of delta = +.1. Set R1 = 8.248, R2 = 8.357, repeat the "Herget step" and "Full Steps", and you get:

Orbital elements:
S/2000 S 11
   Perisaturn 2002 Apr 1.933825 TT
Epoch 2001 Apr  1.0 TT = JDT 2452000.5
M 211.67402              (2000.0)            P               Q
n   0.40533552     Peri.   73.30298     -0.83559992      0.05564627
a   0.1191205      Node   107.06023     -0.17734258     -0.96891430
e   0.3806275      Incl.   34.86662      0.51992536     -0.24105718
P   2.43           H   10.9           G   0.15      q 0.0737800
From 15 observations 2000 Nov. 9-Dec. 17;   RMS error 0.474 arcseconds

If you look at the file of observations you downloaded, you'll see that this is essentially the solution the MPC published.

In most cases, this approach will get you an orbit. In a few cases, you need to play around with delta a bit. In a few cases, especially if the arc is very short, all you get is diverging solutions with absurd eccentricities and residuals; I've got some special, pretty weird code to handle such situations.

This sort of "double solution" is a common pattern. After you get a few more observations, you can usually discard one solution because its RMS errors start to grow. I (and everybody else who looked at this object when it was first found) went with the second solution, not just because of its lower errors, but because the first orbit was so unlike the orbits of the other satellites. And indeed, the object was followed up and got the permanent designation Saturn XXVI, and the name Albiorix; and using most of the data, one gets this orbit:

Orbital elements:
Saturn XXVI = Albiorix = S/2000 S 11
   Perisaturn 2001 Dec 27.043157 TT =  1:02:08 (JD 2452270.543157)
Epoch 2001 Apr  1.0 TT = JDT 2452000.5                  Find_Orb
M 235.3085               (2000.0)            P               Q
n   0.4617462      Peri.   60.1568      -0.8204340      -0.0720332
a   0.1092016      Node   110.2299       0.0006200      -0.9921431
e   0.472264       Incl.   37.1909       0.5717407      -0.1022900
P   2.13/779.54d   H   11.1     G   0.15   q 0.0576295  Q 0.1607736
From 114 observations 2000 Nov. 9-2009 Feb. 20;   RMS error 0.552 arcseconds

Running Find_Orb with (simulated) astrometry for the Mir space station

It is possible, though not necessarily very useful, to run Find_Orb to determine the orbit of an artificial Earth satellite. I didn't have any actual test data to use, unfortunately, so I generated an example. I did this by starting up Guide, setting my lat/lon to that of an MPC station, and generating an ephemeris for a pass of the Mir space station. If you run Find_Orb, open the EXAMPLE file, and double-click on "Mir", you'll load up the simulated observations in question.

Again, a decent guess for R1 and R2 matters here. The default guess of 1 AU is obviously nonsense for an earth-orbiting object. Reset these to .00001 AU. This corresponds to about 1500 km, which is a good "low-earth orbit" distance.

Check the "Earth" box and run a few Herget steps. You'll quickly see convergence to these values:

Orbital elements:
Mir
   Perigee 1998 Jun 20.810218 TT
Epoch 1998 Jun 20.8 TT = JDT 2450985.3
M 171.29447              (2000.0)            P               Q
n5660.54255256     Peri.  132.79498     -0.81371242     -0.44835486
a    6720.118 km   Node    28.13785     -0.15470880     -0.44641400
e   0.0057829      Incl.   51.66791      0.56030106     -0.77439813
P  91.58m           H   24.2           G   0.15      q 6681.256 km
From 10 observations 1998 Jun. 20-20;   RMS error 158.892 arcseconds

The RMS errors are admittedly a bit shocking. The reason for their size will be explained further on. You'll notice that the period of 91 minutes is decent for a low-earth object. The perigee distance still keeps Mir from plowing into the Earth (remembering that the radius of the Earth at the equator is 6378 km).

Before you click on the "Full Step" button, I must suggest that you take an unusual step. By default, when you full-stepped, Find_Orb would shift the epoch back to Jun 16.0. The problem is that attempting to extrapolate this short arc back by about five days (about 700 times its own length) would be unstable; you need an epoch that is closer to the observations. So reset the epoch to 2450985.34. This is right in the middle of the actual observation set, and will be quite stable.

Having done that, a few Full Steps get you this orbit.

Orbital elements:
Mir
   Perigee 1998 Jun 20.810114 TT
Epoch 1998 Jun 20.8 TT = JDT 2450985.3
M 168.93682              (2000.0)            P               Q
n5652.74897459     Peri.  132.46347     -0.81106791     -0.45313414
a    6726.294 km   Node    28.13381     -0.15214408     -0.44720528
e   0.0048460      Incl.   51.67437      0.56481946     -0.77115296
P  91.71m           H   24.2           G   0.15      q 6693.698 km
From 10 observations 1998 Jun. 20-20;   RMS error 39.292 arcseconds

For a longer arc, or a higher satellite, the Moon might come into play. But if you check that box and do more Full Steps, you won't see much of a difference.

The reason for those huge residuals is a simple one. The MPC format allows for a maximum precision in time of .0864 seconds = 10-6 days. In that seemingly brief period of time, Mir can move about 700 meters. When it's 1000 km away (and for some of these observations, it was much nearer), 700 meters = 140 arcseconds. Until Find_Orb can read other (non-MPC) formats, this will be a problem.

Using Find_Orb to predict impact locations

There are several situations wherein one might want to know the location on the earth, or another planet or moon, where an asteroid or space junk is going to impact. Find_Orb has been used to compute the impact location of 2008 TC3, which impacted in northern Sudan on 2010 October 7; LCROSS, which impacted on the Moon on 2009 October; Genesis, which impacted in Utah in 2004; several of the Shoemaker-Levy 9 objects, which impacted Jupiter in 1994; and to produce impact predictions for 2007 WD5, which briefly had a small chance of hitting Mars. It has also been used to find the impact point of the Japanese Haybusa probe, but unfortunately, not all the astrometry for that case is publicly available yet.

The example case we will use will be 2008 TC3. First, go to MPCOBS and request data for 2008 TC3. (You can alternatively get the astrometry here, from NEODyS.)

You can then save that data to a file, and run Find_Orb, and click on Open, and load the data. I usually use MPCOBS or NEODyS to display the data; then I copy everything to the clipboard, run Find_Orb, and click with the right mouse button on the orbital elements area; then I select "Open from Clipboard". Either way, Find_Orb will pause a bit before loading up the data and producing these elements (note that the 'heliocentric elements only' check-box in Settings should be un-checked!) :

2008 TC3
   Perigee 2008 Oct 7.118796 TT =  2:51:03 (JD 2454746.618796)
Epoch 2008 Oct  7.0 TT = JDT 2454746.5                  Find_Orb
q  5849.7208km      (2000.0)            P               Q
H   30.6  G 0.15   Peri.  148.3640      -0.4988773       0.8437959
                   Node   330.1500       0.6888919       0.2475819
e   1.591833       Incl.   23.4175       0.5258794       0.4761423
From 859 observations 2008 Oct. 6-7 (18.8 hr);   RMS error 1.748 arcseconds
IMPACT at   7 Oct 2008  2:45:54.16 lat +20.60295 lon E33.11633

It may appear as if Find_Orb has just done everything for you, as if by magic: here we have the exact place and time of impact! This is slightly misleading. If you click on "Worst Obs", you'll see that there are some strong outliers in the data (mostly cases where the clock wasn't set correctly; you'll see that these have dT values of several seconds, indicating clock errors at that level, but the "cross" residuals are small.) So I'd recommend going into Settings, and setting a "Max Filtered Residual" value of two arcseconds, and clicking OK.

If you do that, and then click on "Filter Observations" several times, you'll eventually get a "No changes made!" message. The elements will then be:

 2008 TC3
   Perigee 2008 Oct 7.118799 TT =  2:51:04 (JD 2454746.618799)
Epoch 2008 Oct  7.0 TT = JDT 2454746.5                  Find_Orb
q  5849.8084km      (2000.0)            P               Q
H   30.6  G 0.15   Peri.  148.3645      -0.4988776       0.8437960
                   Node   330.1495       0.6888962       0.2475869
e   1.591828       Incl.   23.4170       0.5258736       0.4761395
From 718 observations 2008 Oct. 6-7 (18.8 hr);   RMS error 0.872 arcseconds
IMPACT at   7 Oct 2008  2:45:54.50 lat +20.60191 lon E33.11760 

But we still aren't quite done yet. First off, I'd turn on all of the perturbers. (As it happens, it's not particularly relevant in this case. The perturbations between the time 2008 TC3 was found and the time it impacted scarcely mattered. But it's good practice.)

Also, the epoch of the elements is at midnight on 7 October. We need to reset the epoch to be somewhat closer to that of impact. Change the 'epoch' field to read 2008 Oct 7.12, click on "Full Step" again, and you'll get

2008 TC3
   Perigee 2008 Oct 7.118801 TT =  2:51:04 (JD 2454746.618801)
Epoch 2008 Oct  7.12 TT = JDT 2454746.62                Find_Orb
q  5848.4013km      (2000.0)            P               Q
H   30.6  G 0.15   Peri.  148.4778      -0.4994245       0.8433497
                   Node   330.0722       0.6889788       0.2478319
e   1.592770       Incl.   23.4243       0.5252459       0.4768024
From 718 observations 2008 Oct. 6-7 (18.8 hr);   RMS error 0.872 arcseconds
IMPACT at   7 Oct 2008  2:45:54.43 lat +20.59288 lon E33.12104 

That leaves us with a couple of fine points to consider. First, we really ought to account for the fact that a small object such as this will hit the atmosphere and explode a few dozen kilometers above the earth's surface. To account for this, before doing the above steps, you should edit the file 'environ.dat' and add a line such as

COLLISION_ALTITUDE=50000 

to specify the altitude, in meters, for which a collision lat/lon should be computed.

At some point, I'd like to have the program include the effects of atmospheric drag. But at present, drag isn't included in the physical model.

The other fine point to consider is that of uncertainty. If you hit the "Monte Carlo" or "Statistical Ranging" buttons, Find_Orb will start to compute alternative virtual asteroids that match the input data. The resulting "virtual impact" locations can then be shown in Guide. (Which I realize not all Find_Orb users have. But Guide has some very good geographic display capabilities that I didn't want to attempt to replicate in Find_Orb. One possible workaround may be to have Find_Orb produce .kmz files of the sort that could be displayed in Google Earth.)

To make such a chart of impact locations, let the "Monte Carlo" or "Statistical Ranging" process run long enough to generate a few dozen or more virtual asteroids. Their orbital elements and impact locations will be stored in the file virtual.txt. Put this file in the same folder as Guide. Then, download this file:

http://www.projectpluto.com/impact.tdf

to your Guide folder, and run Guide. In Guide, hit the ':' button; this will cause Guide to show a solar eclipse path on the earth. But also, as you zoom in, Guide will put symbols on the map of the Earth showing where the virtual impactors landed.

Such charts can be very illustrative of how the uncertainty area shrinks with additional data. For example, when I first announced that 2008 TC3 was a likely impactor...


http://tech.groups.yahoo.com/group/mpml/message/21070

...we had only the initial data from (G96) and (854), and about 10% of the virtual asteroids didn't even hit the earth. Here's what the impact area looked like with that very short arc. (Note that at the time of this impact, Find_Orb wasn't nearly so well equipped to handle the situation. About all I knew was that the nominal perigee, and that of about 90% of the virtual asteroids, was less than the earth's radius; the code to compute impact lat/lons and to create these charts was added in response to 2008 TC3.)

Fortunately, Cristovao Jacques and Gordon Garradd ignored my statement that the object probably couldn't be recovered. With their data from (E12) and (D90), the uncertainty area was narrowed down to a strip in northern Sudan. With still more data, the uncertainty area shrank to an area about one kilometer long.

Selecting perturbers

Usually, Find_Orb does a good job of deciding which perturbers are actually important for a particular case. As you may have noticed in the above examples, it could recognize that the Earth and Moon are important perturbers for objects that come near the earth, and that Saturn should be used for an object coming near to that planet.

However, its selections aren't always perfect. It errs on the side of leaving out perturbers, because including them can make the program run somewhat slowly. (You'll probably not notice this unless you've an object with a longer orbital arc -- usually some years -- or a slower machine than is usual these days.) You can, therefore, select whichever perturbers you'd like using check-boxes in the main Find_Orb dialog. The next time you do a 'full step' or Herget step, or compute an ephemeris, those perturbers will be used.

There is also a "Toggle Perturbers On" button. When clicked, this causes Find_Orb to use all major perturbers, currently defined to be Mercury through Pluto, with the Earth and Moon handled separately. Asteroids are not toggled on with this button, since they almost never matter very much and slow the program down (admittedly, not by a lot in most cases). Click here for more information on asteroid perturbers.

I expect that people will have different ideas as to what 'All Perturbers' should mean. One can edit the file environ.dat in a text editor, and add a line such as

DEFAULT_PERTURBERS=7007fe

to cause 'All Perturbers' to include Pluto and asteroids. Alternatives are:

DEFAULT_PERTURBERS=7fe (all perturbers except asteroids)
DEFAULT_PERTURBERS=3fe (Mercury-Pluto, with the Earth and Moon treated as a single object at their barycenter)
DEFAULT_PERTURBERS=1fe (same as above, but with Pluto excluded)
DEFAULT_PERTURBERS=7005fe (all objects except Pluto)

Or, if you dislike what 'All Perturbers' provides, you can just toggle objects on manually.

Be advised that I've yet to see a case where Pluto mattered. You have to search a bit to find cases where asteroids matter. Usually, it's because the object passed close by to a large asteroid, at a slow enough speed that the smaller object was perturbed. When this happens, however, it can be used to determine the mass of the asteroid causing those perturbations.

Physical model:

Find_Orb can include the perturbations due to the planets and Earth's moon, those of Jupiter's four largest satellites, several of Saturn's satellites, and those of up to 300 perturbing asteroids. By default, no perturbations are considered (a two-body model) unless you have checked some of the various boxes for perturbers. However, if the object comes within a planet's Hill sphere (the volume where the planet's perturbations are particularly significant), perturbations will be automatically turned on as long as the object remains close to the planet.

Satellites of Jupiter and Saturn are automatically included as the object gets close to those planets. Otherwise, their masses are "thrown into" the planet.

Also, as objects get close to the Earth, Mars, or gas giants, the effects of oblateness are included. Specifically, zone terms J2, J3, and J4 are used. It probably wouldn't be particularly difficult to add in such terms for the Moon, and should be feasible to include non-zonal terms (these would be important for certain lower-orbiting artificial satellites of a sort not yet encountered by Find_Orb.)

In the Settings dialog, there are options to have the effects of solar radiation pressure, or those of comet non-gravitational effects.

The effects of general relativity are always included; see a discussion of GR in orbit determination here.

At present, atmospheric drag is not handled (though I'd like to add it; it would improve handling of meteor observations). Also, in a better world, the masses of inner planets would be "thrown into" the sun unless their check-boxes were explicitly set. (The definition of which planets are considered to be "inner" would depend on the type of object; for a TNO, for example, Jupiter and Saturn are "inner objects".)

The "Auto-Solve" feature

Ideally, it would be helpful if you could load up a set of observations, sit back, and watch while Find_Orb solved for the orbit. This is not quite possible yet, but the "Auto-Solve" button can now do this in most cases.

The idea is that, after loading observations with "Open", you click "Auto-Solve". There is generally a bit of flickering as the automatic solver tries assorted solutions and converges (we hope) on the correct orbit. Finally, the hourglass turns back to a cursor and you're done.

This almost always works for "normal" objects. Near-Earth objects sometimes give it trouble, especially if it turns out that the object is headed almost straight toward or straight away from us. Satellites, both natural and artificial, usually don't work. (I'm working on both limitations. One thing to try is to shorten the initial arc, or shut off the current arc and select some other part of the full arc of observations to start with.) It's always worth giving this function a try, though; it has sometimes surprised me, finding orbits with initial guesses so stupid that I never expected them to lead to solutions.

The "Simplex" method

The main Find_Orb dialog has a "Simplex" button. This refers to the "downhill simplex" (DS) method of function minimization. Much of what I know about DS comes from this section of Numerical Recipes in C; I'd refer you to that for mathematical details, but a capsule explanation follows.

Click on this button, and Find_Orb will use the downhill simplex method in an attempt to find the "best" orbit fitting the current observations. Unlike other methods, this takes into account the actual nature of asteroid orbits; ordinary main-belt orbits will be chosen in preference to exotic near-light-speed ones. You may find that you have to click this a couple of times before it converges somewhere near that point. You can click here for more information about the downhill simplex method. (Note that at least thus far, I've not added the "superplex" method to Windows Find_Orb. It is available in the console version, but doesn't seem to offer any tremendous advantage over the standard simplex method.)

Finding Väisälä orbits

When you have a very short arc of observations (say, a newly-found object that you've seen over the last week or so), the orbit is usually so poorly defined that neither you nor Find_Orb have much of an idea as to how far away the object really is. In such cases, if you use the "Herget Step" or "Full Step" buttons, the orbit may diverge and become totally ridiculous.

Still, you really need some sort of orbit so you can make predictions as to where the object might be in a few nights. The most common way to do this is with a Väisälä orbit, where you assume that the object is near perihelion, and make a guess as to the perihelion distance. In Find_Orb, this is done by entering your guess in the R1 field and clicking "Väisälä". R2 is ignored.

The "classical" sort of Väisälä orbit uses only two observations (usually the first and last) and ignores all others. The method implemented in Guide does take advantage of more than two observations, using the others to minimize errors. You don't really need to know how or why this is done, but if you're curious, you can click here for the mathematical details.

Finding orbits with the method of Gauss

The method of Gauss is one of the most widely used initial orbit determination methods. It was most famously used when the first asteroid, (1) Ceres, was discovered shortly before conjunction with the sun; Gauss was able to determine a good enough orbit to acquire the object after it came out from behind the sun.

The method is described a bit in Methods of Orbit Determination for the Micro Computer and in Fundamentals of Celestial Mechanics. The method of Gauss has been refined and adjusted by swarms of mathematicians, and therefore comes in endless flavors; I chose that described in the former book.

The method of Gauss works with three observations. Click on the "Gauss" button in Find_Orb, and the first and last unexcluded observations in the arc, along with whatever unexcluded observation comes closest to the midpoint of these two, will be used. The resulting orbit will fit those three observations with near-zero residuals, and the other observations with passably small residuals.

Be warned: this method was ideally suited to the pre-computer era. It can get you a good preliminary orbit with relatively little effort. But it does have some quirks, some inherent in the method itself, some inherent in how I implemented that method:

• It's best used on an arc of weeks to a few months. Use a longer or shorter arc, and you may not end up with a reasonable orbit. (Though sometimes, you will get a perfectly good orbit with longer or shorter arcs... and it may fail with intermediate arcs, especially if the object's motion is along an essentially straight line over that time.)

• There can be zero to three solutions returned by this method. If no valid solution was found, you'll get a message to that effect. Otherwise, clicking on "Gauss" again will cycle through the valid solutions. In certain cases, the solution will only be nominally "valid", and you'll see wacky residuals. Usually, that means you're looking at the wrong solution, and clicking on "Gauss" again will get you the right one.

In some cases, you'll get two or three solutions with decent residuals. You may be able to toss out one or two because they correspond to weird orbits beyond Alpha Centauri or objects three kilometers away from Earth. Or, all the solutions may look reasonable, in which case all you can do is wait for more data to find out which solution is "right".

• As currently implemented, this method is intended solely for objects in heliocentric orbit. I am a little tempted to make suitable adjustments to allow for planetocentric orbits, thereby creating an additional tool for handling artificial and natural satellites. But I've not done that yet, and doubt that it would lead to better orbit determination (Find_Orb generally does quite well on such objects already.)

The Settings dialog

Click on the "Settings..." button, and you are provided with control over most of the inner workings of Find_Orb:

Constraints: Allows you to specify that you want, say, a circular or parabolic orbit, or one with a=3, etc. This is quite a powerful (and therefore involved) feature, and therefore has its own section.

Reference: By default, the "reference" shown in the orbital elements is to Find_Orb, but you can change it to give your name. This text will appear among the orbital elements on-screen, and in any "pseudo-MPEC"s you generate.

Monte Carlo noise: One can set the amount of "noise" used in the Monte Carlo and statistical ranging functions; the default amount is .5 arcseconds.

Element precision: By default, this is set to 5, meaning that the angular elements are shown to five digits; other quantities have their precisions set accordingly. You can reset it to any number of digits from 1 to 10.

Note that it can be tempting to say that, if the uncertainties shown indicate that the elements are only good to (say) 0.01 degree, you should set the precision to two or three digits. Oddly, this is really not a good idea. The problem is that if someone attempts to re-create an ephemeris with truncated elements, they'll get serious errors when compared to the original observations. So keeping at least four digits is advisable.

"Use sigma file" check-box: By default, observations are all given the same uncertainties (sigmas) in position, magnitude, and timing. If they are in .rwo format, then the sigmas given in that file will be used. Turn on this check-box, and sigmas from #sigma tags or from sigma.txt will be used.

Alternative orbit format: Normally, you would leave this box checked. Find_Orb will then show MOID information for most planets, and uncertainties in orbital elements if a "full step" is performed. However, this format does not exactly match the eight-line format used by MPC. Un-check this box if you need MPC format (for example, when supplying data to a planetarium program.)

Heliocentric-only: By default, Find_Orb will show heliocentric or planet-centric orbits based on the object with dominant influence at the epoch. If an object is within a planet's "sphere of influence" (the point where planet-centered osculating orbital elements are a better approximation to the object's motion than heliocentric elements), Find_Orb will usually switch from heliocentric to planet-centered elements. You will usually be perfectly happy with this, but if you really need heliocentric elements for some purpose, you can click on "Heliocentric Only" and elements will then be shown in that manner.

• Precise residual check-box: Checking this will cause residuals to be show with an extra digit. This option has long been available in console Find_Orb; I've found it handy mostly in ensuring that Find_Orb was replicating results from NEODyS and AstDyS, right down to the milliarcsecond. (In the normal course of things, the only observations I can think of that might be good to such a level would be HST and Hipparcos ones.)

Comet magnitudes section: The choice here is between "total" or "nuclear" comet magnitudes. (The choice is irrelevant for anything except comets.) This controls the type of magnitude shown in the orbital elements and ephemerides. In each case, only total or nuclear magnitudes from the observations will be considered in computing the absolute magnitude M, and the orbital element display will show M(T) or M(N). (These parameters are also known as M1 and M2. But I can never remember which is which, and assume others may have the same problem.)

There are multiple problems associated with this change. If you've got data with some nuclear magnitudes and no total magnitudes, and tell Guide to use total magnitudes, the orbital elements and ephemerides will show no magnitudes at all. The MPC's optical astrometry format lets you specify a comet magnitude as total or nuclear, but doesn't let you specify a band for it. Observers differ wildly in the validity of their magnitude data, and that's especially true for comet magnitudes; ideally, we'd have some uncertainties (sigmas) for the photometry. (Find_Orb allows one to set magnitude uncertainties, but it's not of much practical use, since one rarely knows what those uncertainties are!)

Physical model: The "standard" model is the simple gravitational one (and the only one most people will ever use): sun and any selected perturbers, general relativity, oblateness of planets. In this model, a "full step" leads to Find_Orb fitting six parameters: the position and velocity of the object at the time of epoch. Select "Include SRP", and the effects of solar radiation pressure will be considered, leading to a seventh parameter being fitted. This is frequently useful when dealing with astrometry for artificial satellites when the arc gets long enough. As the arc grows, you first notice trends in the residuals; when these grow unacceptably, you'll often find that if you turn on the use of SRP and do a "full step", the residuals will become reasonable, and Find_Orb will provide a decent estimate of the ratio of the object's area and mass. (This has also been used for a very few small asteroids that lingered in Earth's neighborhood long enough for SRP to be noticeable.)

Comets occasionally show more complex non-gravitational effects, due to vaporization as they approach the sun. There is a standard comet model that approximates this behavior, and encapsulates it in parameters A1 and A2. Select "Comet Non-Grav", and the non-gravitational parameters A1 and A2 will be fitted. This only makes sense for comets, and only if the arc is pretty long.

Note that A1 and A2 are in units of AU/day^2. MPC uses km/day^2. I may revise things so that you can select which system is used. I've stuck with AU/day^2, because those units are used by Patrick Rocher at the IMCCE, allowing me to verify that Find_Orb is doing all this correctly.

Filtering section: The "filtering" section determines what happens when you click on "Filter Obs". The default is to filter as before, in the manner probably used by almost all orbit computers: observations greater than a "Max residual" limit are excluded, and those within that limit are included. The limit is specified in sigmas. By default, Find_Orb assumes identical uncertainties for all observations, though this can be adjusted by setting uncertainties for observations in a variety of ways.

Alternatively, one can assign a "blunder probability". This new way of filtering observations is described here. The basic idea is that observations far from prediction are included to the degree we think they are not blunders. (By default, I'm assuming that 2% of all observations are blunders. If we knew, a priori, which observations were in that 2%, we'd just eliminate them... but we don't.) An observation close to prediction will be included at nearly full weight; one ten sigmas off prediction will have almost zero effect on the solution. An observation two or three sigmas will be included, but to a lesser degree than would otherwise be the case.

Over-observing parameters: The "over-observing parameters" section addresses situations where one observer gets so much data that her results dominate the orbit solution. For example: on occasion, you'll find that an object was observed three times by one observer, four by another, three by yet another... and then there's somebody who made 57 measurements on one night. In the normal course of things, the guy who made 57 measurements will dominate the solution, even if he made 57 bad observations. (The usual cause of over-observing is data gathered for a light curve, but there can be other reasons.)

I've tried to come up with a good solution to this issue, and have implemented one based on ideas kindly suggested independently by Marco Micheli and Alan Harris, currently called the "t-and-k scheme". The idea is that if, over a t-day time period, the number of observations n_obs made is greater than a constant k, then the observations should be reweighted downward by k/n_obs. Thus, the weight given to the entire arc will be as if only k observations had been made. This effectively caps how much influence one observer can have over a given time span t.

In the Settings dialog, one can set t ("time range", defaulting to one day) and k ("observation ceiling", defaulting to five observations.) With these parameters, if Find_Orb notices that a given observation is one of twenty made from the same observatory during the same day, the observation will be given half the weight it would otherwise receive.

Thus, the total effect of the 57 observations would be as if five had been made. All 57 will still be used; this helps to remove random noise that would have been apparent had we, say, just taken five observations and tossed the rest. But that observatory will no longer dominate the orbit.

Set the "time range" to be zero, and it becomes impossible to over-observe an object, and this scheme is effectively disabled.

Getting constrained orbits

Usually, clicking on "full step" results in a general solution in which all six orbital parameters are adjusted. But under some circumstances, it can be helpful to impose a "constraint" on the orbit. The most common examples are specifying that e=1 (a parabolic orbit) or that the object is at perihelion at the time of observation (a Väisälä orbit.)

It's possible to use Find_Orb to specify much more general constraints. For example, let's say you have loaded data for an object into Find_Orb and have found the value a=4.5 for the semimajor axis. If you're familiar with minor planet groups, you'll know that such a value would be possible, but very unlikely. It's more likely that the object is a Hilda (the next group in, with a=4.2 on its outer edges) or a Jupiter Trojan (the next group out, with a=5.05 on its inner edge.)

You can test the first possibility by clicking on the Settings button, and entering "a=4.2" into the Constraints box. Click OK, and then click on Full Step. Find_Orb will look for an orbit with that semimajor axis. You may have to hit Full Step several times before it converges on an a=4.2 solution, and when it does, you may find that the RMS error grows to an unacceptably high value. If so, the possibility that the object is a Hilda will have to be rejected.

Similarly, you can set "a=5.05" in that edit box and see if that converges, to confirm or dismiss the possibility of a Jupiter Trojan orbit.

You can also set constraints on the eccentricity, perihelion distance, orbital period, and inclination, and even combine constraints, such as:

e=.3        (constrain eccentricity to be .3)
e=0         (force a circular orbit)
q=1         (force an orbit whose perihelion brushes the earth's orbit)
i=15        (force an inclination of 15 degrees)
P=248       (force orbital period of 248 years,  that of a Plutino)
a=5.1,e=0   (combined forcing of a 5.1 AU semimajor axis and circular orbit)
e=1,i=144   (constraint for a Kreutz sungrazer type of orbit)

Usually, when you have a very short arc of observations (or a very distant object), you'll find that forcing odd constraints does not necessarily raise the RMS error very much. This simply means that the range of possible values for each parameter is still very large. As more astrometry is gathered, the range drops sharply, meaning that the parameters are better-determined than was once the case. (Try to force an orbit too far from its valid range, and Find_Orb will often fail to converge to a solution.)

It's also possible to constrain the magnitude slope parameters. For comets, this defaults to K=10; for asteroids, to G=.15. However, you can set constraints such as K=18 or G=.12. This isn't necessarily a good idea (it would be better if Find_Orb had the ability to determine K or G from the observed photometry), but one can do it if desired.

The "filter obs" button

Click on the "Filter Obs" button, and Find_Orb will perform observational filtering using the scheme selected in the filtering section of the Settings dialog. If this filtering doesn't result in any changes, it will report that fact. Otherwise, it will do a "full step". Applying this button a few times will usually result in convergence on a good solution, and automatically rejects outliers.

Determining uncertainties with the Monte Carlo and statistical ranging methods

Often, Find_Orb breaks a cardinal tenet of physics: it gives all sorts of numbers without giving you any idea as to how accurate they are. You can give it a few observations made over a one-week period in 1973 and generate ephemerides for the present, and it will happily crunch out positions to .1-arcsecond precision, totally ignoring the fact that said positions are probably nowhere near the right part of the sky. (You can usually say that a short arc will lead to huge errors when extrapolated in this way, and that a long arc will have high accuracy, but it's surprisingly hard to be sure unless you actually compute the errors.)

If you're an observer, you have a real need to know just how exact the ephemeris is. If the positions are good to within a few arcseconds, you may decide that there is no point in re-observing the object; the current knowledge of its orbit is good enough. If the uncertainty is a few arcminutes, you know that taking an image or two will nail down the orbit nicely. And if the uncertainty is huge, you'd at least like to know what part of the sky should be searched.

Find_Orb can now generate data to be used in plots such as this. Here, you can see the uncertainty region (which has elongated into a line and developed a bit of a banana-like curve), showing it as it was at the time of the first recovery observation... and that observation (yellow symbol) is embedded in the lower part of that banana.

There are several ways of computing ephemeris uncertainties. Jim Baer's CODES software and the JPL Horizons system both use a covariance matrix method; this is very fast and works beautifully for small/medium uncertainties. (It's the method Find_Orb uses when sigmas are shown on-screen for orbital elements.) Another common method is the "observational Monte Carlo" one, used by John Rogers' CCD Astrometry and by the OrbFit software (OrbFit also supports several other methods of computing uncertainties). Steve Chesley described six different ways of going at the uncertainty problem in an MPML post.

Usually, the method to use in Find_Orb is the "Monte Carlo" approach. In this scheme, you take the initial observations and add some Gaussian "noise" to them, and solve for the orbit that fits these slightly-mangled points best. This creates a "virtual asteroid", representing a possible orbit that fits your actual observations passably well. Then you repeat.

Eventually, you have a swarm of "virtual asteroids" that, it's hoped, combine to illustrate the region in which the real asteroid might actually be. That region usually starts out as an ellipse, and when plotted on the sky, the virtual asteroids look like a squashed globular cluster.

As time passes by, the "virtual asteroids" spread out in one direction (usually along the line of variation), and the uncertainty region resembles a long, thin line. The line will develop some kinks and curves over time, especially if part of the line passes close to a planet. Sometimes, virtual asteroids will cluster at one or more points on the line, indicating areas where the object is most likely to be found.

To use this method in Find_Orb, proceed as follows. Solve for the orbit as you normally would. Set the "epoch" to be passably near the time span over which you want to know the positional uncertainty, and click on "Monte Carlo".

Find_Orb will pause while it computes the first virtual asteroid. You'll see the element and residuals display flash with a revised orbit, one with somewhat higher residuals (usually) than you originally had. Then it will flash again and again, as new virtual asteroids are computed. The "Monte Carlo" button will change to count the number of virtual asteroids created thus far. You can click on that button to stop computing virtual asteroids, and click on it again to resume computing them. (You can also switch to other tasks while virtual asteroids are computed in the background.)

As the virtual asteroids are computed, their orbital elements are added to the file mpcorb.dat. This file has the format of the MPC's MPCORB database, which can be read by most desktop planetarium programs.

One key point here: you may well want to set the amount of Gaussian noise added in Monte Carlo computations. If you click on "Settings...", you'll see an option to set the "Monte Carlo noise" level; this defaults to .5 arcseconds. Exactly what this value really ought to be is a matter of some debate, but it ought to reflect the uncertainties in the observations. That is, if you figure they're good to one arcsecond, a noise level of 1 makes sense. I've used the method with older photographic-era measurements good to about four arcseconds, and used a noise level of four, and have used a smaller noise level with more precise astrometry.

Use your planetarium software to display the mpcorb.dat generated by Find_Orb, and you will see a slew of numbered asteroids illustrating the current uncertainty area of the "real" asteroid. (I chose this approach because it meant people could view the uncertainty area in the program of their choice; it also makes it easy to see how that area changes over time.)

If your arc is quite short, the Monte Carlo orbits will usually "blow up"; you'll get orbits light-years away that make no physical sense. If that happens, the thing to do is to stop the MC process and click instead on "Statistical Ranging".

SR is, from the standpoint of the user, not much different from Monte Carlo. The advantage of SR is that you can use it even if you have only two observations, and/or a really short arc. So why isn't it the preferred method? Simple: if the arc is much longer (long enough so that "full steps" converge), SR becomes extremely slow. So only use it if MC fails.

SR works by taking a guess at the distance and radial velocity of the first observation. (It has some logic in it that enables it to set some limits on both distance and radial velocity, so it won't waste time considering hyperbolic orbits or orbits that are light-years away.) It then finds the orbit that would bring an object to the last observation. If the residuals are low enough that the orbit can be considered "reasonable", then the orbit is written out to mpcorb.dat and virtual.txt (just as would happen in Monte Carlo). Hence my comment that, from the viewpoint of the user, SR doesn't "look" all that different from MC. They're just different ways of finding alternative orbits that fit whatever data you have.

To be added: when you generate an ephemeris, all the variants should be integrated so that scatter-plots can be shown of "where the asteroid might be at time X", as well as just showing that "the uncertainty at time X is Y arcseconds at position angle Z". Also, uncertainties can often be computed just using the covariance matrix; and uncertainties in the orbital elements from any of these methods really should be shown to the user.

Saving orbital elements and residuals to a file

At any point, if you want to save the results, click on "Save Elements" and you'll be able to create an ASCII file containing the text in the "orbital elements" box. You can also click on "Save Residuals" and create a text file containing all the data in the observations/residuals list box.

Also, whatever elements are shown in the dialog are also saved in the file elements.txt (in the "standard" eight-line MPC format, same format as shown in the dialog), and in the single-line MPC format in the file mpc_fmt.txt. This latter format is the one used in the MPCORB.DAT and similar files on the MPC's ftp site and the mirrors of that ftp site.

Be advised that elements.txt, and any similar file you make using "Save Elements", contains data beyond the eight-line MPC formatted data. The file will also contain a state vector, MOIDs for eight planets, Tisserand criteria with respect to the Earth, Jupiter, and Neptune, and some "housekeeping" data such as when the elements were written, which perturbers were included, and the version date for Find_Orb. You'll probably just ignore most of this data, but it can come in handy every now and then.

By default, the format of the "Save Residuals" file is essentially that of the list box. However, if you specify a .res extension, then Find_Orb will realize that you want the MPC format for residuals. In this abbreviated format, only the date, observatory code, and residuals are stored for each observation, packed so that there are three observations/line, and a brief blurb is given for each observatory. Thus, something like

 1  3 24.30998   704 12 10 22.03  -08 50 15.1  -0.06  -0.04   0.562  1.556
 1  3 24.33638   704 12 10 16.94  -08 49 50.4  -0.25   0.20   0.562  1.556
 1  3 24.34986   704 12 10 14.38  -08 49 38.8   0.19  -0.72   0.562  1.556
 1  3 25.19163   859 12 07 38.08  -08 36 23.3   0.30   1.86   0.565  1.560
 1  3 25.19471   859 12 07 37.48  -08 36 22.1   0.07   0.18   0.566  1.560
 1  3 25.40278   474 12 07  0.19  -08 33  3.6   0.05  -0.32   0.566  1.561
 1  3 25.40405   474 12 06 59.96  -08 33  2.3   0.16  -0.19   0.566  1.561
 1  3 25.40925   474 12 06 58.97  -08 32 57.9  -0.10  -0.58   0.566  1.561

is transformed to

010324 704  .06-  .04-    010325 859  .30+  1.9+    010325 474  .16+  .19-
010324 704  .25-  .20+    010325 859  .07+  .18+    010325 474  .10-  .58-
010324 704  .19+  .72-    010325 474  .05+  .32-

Station data:
(474) Mount John Observatory, Lake Tekapo  (S43.9876 E170.4650)
(649) Powell Observatory, Louisburg   (N38.6464 W94.6997)
(704) Lincoln Laboratory ETS, New Mexico  (N33.8185 W106.6591)
(859) Wykrota Observatory-CEAMIG  (S19.8240 W43.6903)
(918) Badlands Observatory, Quinn   (N43.9908 W102.1306)

Creating and saving an ephemeris

Click on the "Make Ephemeris" option, and Find_Orb will bring up a small dialog box with a starting time, a number of steps, and a proposed step size in days, and a lat/lon. Edit these and click on "Go", and it will generate an ephemeris in the list box. Click on "Save", and you can write the results to a file.

The starting time for the ephemerides can be input in a variety of formats; click here for details.

In computing the ephemeris, Find_Orb will include the effects of any perturbers you selected in the "main" dialog.

You can enter a three-character MPC observatory code in the latitude box. Do this, and Find_Orb will ignore the longitude box and just set your position to that of the MPC code. Use an MPC code of 'Sun' for heliocentric ephemerides, 'Mer' for Mercuricentric ephemerides, 'Ven' for Venus-centered (Cytherian? Venerian?), 'Mar' for Mars, and so on, up to 'Plu' for Pluto and 'Lun' for the Moon. 'Ast1' will generate Ceres-centric ephemerides, and 'Ast15' will generate (15) Eunomia centered ephemerides. Note that only the 300 asteroids used in Find_Orb are available.

You can also add your own 'MPC code' for topocentric sites on the Earth, Moon, Sun, or planets; see rovers.txt for details. (The most common use of this has been to support observations of artificial satellites from people who don't have MPC codes. But it can also be used for generating ephemerides from the surface of another planet. I've also used it, for example, to specify the exact locations of all telescopes on Mauna Kea, instead of leaving them lumped together under the (568) code.)

There are six radio buttons to choose from six different types of ephemerides. The default is "Observables", which provides data such as RA/dec, alt/az, and similar quantities. If you choose to make an "Observables" ephemeris, you can use check-boxes for items such as motion details (how fast the object is moving, in arcminutes per hour, and the position angle at which it is moving), alt/az, radial velocity, and so forth. (Of course, the alt/az checkbox doesn't accomplish anything for geocentric ephemerides.)

Also available is a "State Vectors" radio button. Click on this, and each ephemeris line will give a Julian Day followed by the position and velocity in Cartesian coordinates, in AU and AU/day, relative to the observer. (Thus, one can set the "observer" to be station code 500, center of the earth; or Sun, center of the sun; or Lun, center of the moon, and so on, entering the station code in the "Latitude" box.)

The "Cartesian coords" results in the same sort of output, except that only the positions are given (i.e., the velocities are omitted).

Occasionally, it may be useful to have a table of the orbital elements, to show how they change over time. The "MPCORB elements" and "8-line elements" options will accomplish this, resulting in the orbital elements at each ephememis step being output.

Finally, the "Close Approaches" radio button will cause only the dates/times of closest approach over the ephemeris span to be displayed. For this, it's best to select geocentric (or heliocentric or planetocentric) coordinates; a short step size would be needed for topocentric coordinates, to account for the relatively high-frequency rotation of the earth.

Observable ephemerides use the UTC (Coordinated Universal Time) scale. Others use the TT (Terrestrial Time) timescale.

Making a "pseudo-MPEC"

Whenever you generate an ephemeris, Find_Orb creates a "pseudo-MPEC" with the filename mpec.htm. Clicking on the "Pseudo-MPEC" button will bring it up in your browser. This looks a bit like the standard Minor Planet Center MPECs, listing the observations, stations supplying the observations, orbital elements, and ephemerides (with the time span determined by that given in the "Make Ephemeris" dialog.) It differs from the MPC version in that it is heavily linked; you can click on an observation to see its residuals, or on a residual to see its observation, or on a three-character observatory station code to find out which observatory that is, or on the observatory's lat/lon to get a map of that observatory, or (usually) on an observatory name to get that observatory's home page. Click here for some examples.

A few hints to customize/improve the "pseudo-MPECs" you generate:

Getting data about the observations

Click on any observation in the list box at the bottom of the dialog box, and you'll get a few lines of text just below the list box giving data about that observation. Here's an example:

Elong 179.3    Phase   0.7    RA vel -4.28'/hr   decvel  9.15'/hr   dT=-1.49 sec
ang vel 10.10'/hr at PA 334.9   radial vel -1.150 km/s  cross -0.22  9.6 days ag
Delta= .00807  r= 0.9918  mag=18.40  mag (computed)=18.34   2009 Jan 16  6:38:01
MPEC ????-B14  Obj alt 73.5 az 127.9  Sun alt -73.2 az 310.0
(G96) Mt. Lemmon Survey  (N32.44283 W110.78872)
     K09B00D* C2009 01 16.27641 07 53 19.03 +21 35 31.9          18.4 V EB014G96

The first line gives the elongation and phase at the time of the observation, followed by the motion of the object in RA and declination. (The motions are in units of arcminutes per hours, which are the same as arcseconds per minute.) Next, the "residual in time", dT, is given. This describes the difference between the computed time of the observation and that given in the observation data.

The "residual in time" is especially helpful if you think the observer may have had a clock error. You will sometimes see observations of very fast moving objects (VFMOs) with large, but consistent, residuals. When that happens, you should check the dT values and the "cross" value (given just below the dT value). The "cross" value indicates how far off the observation was from the predicted track.

If you see a series of observations wherein the "cross" values are all normal residual values (say, under an arcsecond), but the residuals in RA and dec are abnormally large, you can suspect a timing error, and dT will tell you how large that timing error is.

The second line tells you the rate of apparent motion of the object, and the position angle at which it is moving. This is essentially the same data as "RA vel" and "dec vel", expressed in a different form. It is followed by a radial velocity, in km/second, and the cross-track residual. If the observation was reasonably recent, you'll also get an indication of how many days (or hours or minutes) ago it was made.

The third line gives the distance to the target, delta; the distance between the target and the Sun, r; the observed magnitude (this may be blank); the computed magnitude (this may also be blank, if no magnitude data was reported for the target); and the time of observation, expressed in "conventional" calendar/hh:mm:ss notation.

The fourth line shows the altitude/azimuth of the object and of the sun at the time of the observation. That can be a useful sanity check at times; if the object is below the horizon, or the sun above it, you should probably disregard the observation. In some cases, there will be a reference preceding that data, to an MPEC or MPC IAUC or other publication.

The fifth line gives the observatory name, three-character code, latitude, and longitude.

The sixth and final line shows the observation data in the original MPC 80-byte form.

Setting coordinate epochs other than J2000

By default, the input from an MPC-formatted or .rwo- formatted file is assumed to be in J2000. However, you can insert #coord epoch lines to specify a given epoch. For example,


#coord epoch 1905.0
     J96X01X  C1997 10 18.97256 00 30 26.10 +03 07 12.3                      691
     J96X01X  C1997 10 19.65954 00 29 27.92 +03 01 07.2                      691
#coord epoch 0
     J96X01X  C1997 10 21.03352 00 27 34.62 +02 49 01.9                      709
     J96X01X  C1997 10 21.72050 00 26 37.40 +02 43 00.9          17.3        709
#coord epoch 2000
     J96X01X  C1997 10 22.06400 00 26 09.76 +02 40 01.9                      709
     J96X01X  C1997 10 22.75098 00 25 13.19 +02 34 04.0                      709

The first two observations would be assumed to be 1905.0 coordinates. The next two, with "epoch 0", would be assumed to be mean coordinates of date. (Properly speaking, there should be a provision for apparent coordinates of date, which see common use in older publications.) And for the last two, the default epoch would return to J2000.

Note that the coordinates will be read in, then precessed immediately. So the coordinates shown in Guide will be in J2000.

Epoch of elements box

Find_Orb will initially set a "reasonable" epoch for the orbital elements of the object, rounded to an integer Julian Day (noon UT of a particular day). Normally, you will not bother to re-set this. However, you can reset the epoch to any desired value, by entering it in year/month/day form (such as "2005 4 13" for 2005 April 13); or in Julian Day form (such as "2453400.5" for JD 2453400.5 = 2005 Jan 30.) In truth, you can input the epoch in a variety of formats; click here for details.

Some may want all orbits to be on a fixed epoch. (MPC orbits, for example, are usually at 200-day epoch intervals.) One can edit the file environ.dat and add a line such as

EPOCH=2455600.5 

to have all orbits start out with epoch JD 2455600.5 = 2011 Feb 8. (One can then override that fixed epoch with the epoch edit box, of course.) One can specify the epoch in environ.dat with any of the time formats used in the epoch box.

Asteroid perturbers

If you have occasion to turn on the "asteroids" perturber box, you may wonder exactly which asteroids are included. Commonly, this would be (1) Ceres, (2) Pallas, and (4) Vesta (the "CPV" model), and indeed, Find_Orb used to include only these.

However, Find_Orb can now include perturbations from the 300 most massive asteroids. Click here for details. (As you'll see, it's now required that you download a file of orbital elements in order to get this capability.)

When Find_Orb doesn't find an orbit

In the above examples, Find_Orb converged (eventually) to a solution. This is not always the case, which is why the program can't (yet) just read the MPC observations and spit out an orbit, unattended. Sometimes it will latch on to a solution involving a Warp 10 orbit around Alpha Centauri; sometimes the Herget method collapses toward r1=r2=0 (that is, the object is floating just in front of the telescope). Such "problem cases" have gotten rarer over the years, but they still happen sometimes.

One way to avoid this is to give Find_Orb a better handle on the actual distance to the object. You can do that by setting R1 and R2 to a better guesstimate, and then clicking on "Herget" a few times. As you have seen, this was absolutely essential for the S/2000 S 11 and Mir cases.

Also, in the case of some comets, I've found it helpful to restrict the output to a cometary orbit. If you've run through a 'full solution' and the eccentricity is pretty close to 1, try clicking on the Settings button and entering "e=1" in the Constraints box. Then do a few more iterations; you'll probably see that the residuals don't change very much. (This is especially important when figuring orbits for SOHO and STEREO sun-grazing comets. These can give the software a serious headache! For these, I usually set r1=r2=1, then use a constraint of "e=1,i=144", corresponding to a Kreutz sun-grazer. A few Herget steps usually get a passable answer, with the full steps polishing it off. If the arc is long enough, I can then switch the constraint to just "e=1" and do full steps successfully.)

For non-comets, use of other orbital constraints can be helpful. You might try setting the eccentricity and/or the perihelion distance to a likely-seeming value, for example; if it works, you can then experiment with other nearby values to see if the residuals drop.

Sometimes, if only a few observations are available, you'll get wild orbits that still have small residuals. This just means that the orbit was very poorly determined. The resulting orbit will generally still give decent apparent positions for a few days, even if it does look rather odd. After those few days, you can then observe the target again and get more positions, which may be enough to get a "realistic" orbit.

I've yet to see the program fail completely, though sometimes a little playing with parameters was required before it behaved itself properly and didn't generate "RMS error 345184.9 arcseconds"-like nonsense.

Excluding observations

Quite often, you'll find that, while most of the observations "fit" the orbit well, with residuals of under an arcsecond, there are a few observations with suspiciously high residuals, that make you wonder if perhaps that observation was of poor quality. If so, double-click on the observation in question; Find_Orb will re-fit the orbit, with that observation excluded. Double-click on the observation again, and you can toggle it back into consideration.

You can also select multiple observations, then click on the "Toggle Obs" button. All the observations will then be turned on or off.

You can also add lines to the input file of astrometry to indicate that particular observations are shut off. For example, if you loaded the following data:

     J97Z99Z  C1997 04 21.12432 13 56 34.48 -09 11 06.7         x13.5 V      413
     J97Z99Z  C1997 04 21.16887 13 56 49.70 -09 11 08.1          13.5 V      413
     J97Z99Z  C1997 04 21.21342 13 57 03.97 -09 11 16.6                      413
#suppress_obs
     J97Z99Z  C1997 04 21.25797 13 57 16.44 -09 11 31.6                      413
     J97Z99Z  C1997 04 21.30252 13 57 26.35 -09 11 51.4          13.4 V      413
     J97Z99Z  C1997 04 21.34707 13 57 33.13 -09 12 14.3                      413
#include_obs
     J97Z99Z  C1997 04 21.39162 13 57 36.42 -09 12 37.9                      413
     J97Z99Z  C1997 04 21.43617 13 57 36.13 -09 12 59.7                      413
     J97Z99Z  C1997 04 21.48072 13 57 32.47 -09 13 17.1                      413

...the middle three observations would be turned off when the data were first loaded. Also, the first observation would be toggled off by the 'x' in column 65 (just before the magnitude value, not otherwise used in an MPC report).

You should consider the possibility that the observation really is good. Find_Orb follows the classical method of finding an orbit that reduces the "overall" error, as expressed by the sum of the squares of the residuals. In some unusual circumstances, this can cause a perfectly decent observation to have artificially high residuals. For example, if you have a "cluster" of observations and add in one observation from the distant past or future, that distant observation can have anomalously large residuals.

Under no circumstances should high residuals keep you from reporting the astrometric observations in question to the MPC! The people at MPC will be quite capable of deciding whether to ignore a given observation, and (in part because they will have more observations) will be far better equipped to judge which observations are really problem cases.

Setting the uncertainties for observations

Not all observations are created equal. Back in the era of photographic astrometry, results good to a few arcseconds were perfectly okay. In the beginning era of CCD astrometry, results were often good to an arcsecond (still all-too-often the case, I'm afraid... which is a pity; getting much better astrometry isn't necessarily difficult.) People using more modern star catalogues, such as the UCAC, Tycho-2, and CMC-14, are routinely getting results good to a tenth of an arcsecond or so.

In theory, every observation should come with an uncertainty (sigma) that gives you an idea as to how much it should be trusted. We would then know that the data that is good to 0.01 arcsecond should have a much greater impact on the orbit solution than data that's only good to ten arcseconds. We should also know the uncertainty in the magnitude and in the observation time.

That's the theory. Unfortunately, we rarely know what the sigmas actually are. The 80-column MPC format allows one to report the time of observation, RA/dec, and magnitude, but not the uncertainties in those quantities. So we have to resort to the following crude methods.

At the very lowest level, the precision of the input sets a lower limit on the uncertainties. If the declination is given to 0.1 arcseconds (the usual value), the sigma on the position will be at least 0.1 arcseconds. If the magnitude is given to an integer magnitude, then the sigma on the magnitude will be at least one magnitude. Usually, the uncertainty will actually be worse than this (and will be determined using one of the following methods), but looking at the number of digits in the input data at least lets us avoid the error of assuming that "2014 January 28.14" represents time to the millisecond level.

If your data is in the .rwo format used by AstDyS and NEODyS, the sigmas on magnitude and position are given along with the data. (AstDyS and NEODyS use a complicated model due to Chesley, Baer, and Monet to estimate astrometric and photometric uncertainties.) .rwo data does not (at least right now) include timing uncertainties (it has an "accuracy" determined from the number of places given in the input data), so those are determined as they would be for MPC's 80-column astrometric data format, as described below.

The 80-column format doesn't allow you to specify any sigmas whatsoever for optical data. So for data in that format, Find_Orb makes some guesses based on the history of some stations as to the uncertainties in time, date, and magnitude. These guesses (essentially rules such as "(J95) really gets good timing, down to about 0.1 seconds" and "astrometry before 1993 is usually photographic and not good to better than three arcseconds") are given in the file sigma.txt, which replaces the now-obsolete weights.txt. You can edit sigma.txt with any text editor or word processor. (Click here to see the default version. It's quite heavily commented and should be easy to follow.)

This is a really barbaric way of assigning uncertainties, but is much better than assuming all observations to be equal.

If you wish, you can force all observations to be treated equally, by un-checking the "Use Weights" check-box in the Settings dialog. If this is turned off, all astrometry is assumed good to 0.5"; all magnitudes to 0.5 magnitudes; and all times to one second (except radar observations; those are always assumed to be "perfectly" timed, which is pretty close to reality.)

Note that you won't normally want to do this. About the only situation where you would use uniform uncertainties for all observations would be to allow you to compare Find_Orb's results to those from a different program.

If you want to change these 0.5"/0.5 magnitude/one-second default uncertainties, edit the file sigma.txt and look for the "catch-all" line near the bottom, just above the paragraph about "how to assign sigmas".

You can set sigmas within the input file of astrometry, by adding lines such as

(One or more records 80-column astrometric data; group A)
#Sigma 0.3
#Mag sigma 0.01
(One or more records 80-column astrometric data; group B)
#Time sigma 15
#Sigma 0
(One or more records 80-column astrometric data; group C)
#Mag sigma 0
(One or more records 80-column astrometric data; group D)

Here, uncertainties for the observations in group A are assigned "normally", using the sigma.txt rules. Observations in group B will all be given positional uncertainties of 0.3 arcseconds and photometric uncertainties of 0.01 magnitude.

In group C, we set a timing uncertainty of 15 seconds. #Sigma 0 tells Find_Orb to revert to the default sigma.txt rules for astrometric uncertainty; the photometric uncertainty is left at 0.01 magnitude. In group D, we go back to sigma.txt for photometric sigmas, too.

Finally, you can specify uncertainties by selecting one or more observations, then clicking on the "Sigmas" button in the Find_Orb dialog. You'll be prompted for the positional uncertainty for the observation(s), in arcseconds. This makes it look as if only the uncertainty in position can be reset in this manner. In truth, there's a little hidden trickery here: you can enter (for example) m.3 to reset the magnitude uncertainty for the observation(s) to 0.3 mags, or t.01 to set the uncertainty in timing to 0.01 seconds.

Weaknesses in the error model: Find_Orb has, I think, the most advanced modelling of errors of any orbit determination program. In addition to allowing one to set uncertainties in position, time, and magnitude, it has blunder management to mitigate the effects of dubious or completely wrong data, and "over-observing" parameters to ensure that an observatory that makes huge numbers of measurements can't dominate the result too greatly. However, aside from the fact that determining uncertainties is, in itself, an uncertain business, it still has the following weaknesses.

Correlations between observations: Find_Orb mostly assumes that errors between successive observations are uncorrelated. In reality, it's common for people to get a few observations on a given night where the timing is consistently wrong (i.e., if you start out with your timing wrong by ten seconds, it'll probably still be wrong by ten seconds at the end of the night.) If your photometry setup isn't too good and you're measuring magnitudes that are .3 mags too faint, they'll probably all be .3 mags too faint. And so on. In this regard, Find_Orb's simple model in which errors are random and (mostly) Gaussian falls flat. (The deviation from Gaussian is mostly handled by the use of blunder management.)

I've only vague ideas as to how this might be fixed. The currently implemented "over-observing" strategy does help get around the most noticeable aspect of this problem: namely, the fact that if one observatory gets a zillion observations on one night, and they're all consistently off, it can overwhelm the (correct) observations made elsewhere and dominate the solution.

Positional uncertainties are "circular": Find_Orb assumes that the uncertainties in RA and dec are equal. In reality, the uncertainty in position is often at least somewhat elliptical. If a formal uncertainty is determined by the astrometric fit, it can be very elliptical in some cases. The most likely real-world case of which I can think involves sungrazer astrometry. The images from SOHO and STEREO sometimes give you a few reference stars, with the comet not especially near those reference stars. In that situation, where you're extrapolating from the area of "known stars" to a relatively distant point, the uncertainty area may get distorted into an ellipse.

Another situation would be astrometry done in the micrometer era. This was done in a variety of ways, but many of them resulted in the right ascension being measured to greater or lesser accuracy than the declination.

The .rwo format used by AstDyS and NEODyS does allow the RA and dec to have different sigmas, which would handle the micrometer era observations pretty well. It wouldn't help very much with the uncertainty ellipses from SOHO or STEREO or similar cases, where the uncertainty ellipses aren't aligned with RA and dec.

I've stuck with the assumption that the RA and dec uncertainties are identical for several reasons. First, they usually are close enough to equal anyway; most of the time, people get enough stars around all sides of the target to allow the uncertainty to approach being a circle. Second, it's rare to get an uncertainty ellipse (or, equivalently, the correlation between RA/dec or a covariance matrix) with an observation; I've never seen it happen. In reality, getting any sort of uncertainty data at all is rare; that's why I've gone through such extreme exercises to estimate uncertainties from other sources.

Why timing uncertainty matters: Traditionally, orbit solutions are based on the assumption that the observer measured the time exactly. However, not everybody does get timing right. With a fast-moving object, you can often see than one or more observers failed to set the timing correctly, and their positions fall a few seconds "ahead" or "behind" everybody else's. In Find_Orb, you can see this by looking at the information shown when you click on such an observation: you'll see that the "cross-track" error is reasonable, even though the overall residuals may be quite high. In short, the object really did pass by that RA/dec; it just didn't do so at the time specified in the observation.

Most objects are slow enough so that timing errors are not too troublesome. The problems arise when the object is moving quickly. For example, suppose your astrometry is good to 0.7 arcseconds, and you observe an object moving 0.4 arcseconds per second, and your timing error is six seconds. That is to say, your timing error could add an error of 2.4 arcseconds along the object's line of motion. (The perpendicular component -- the "cross-residual" -- would be unaffected by the timing error, and still good to 0.7 arcseconds.)

The total error along the direction of motion would be sqrt( .72+2.42)=2.5 arcseconds. Essentially, the object's uncertainty ellipse is no longer a circle of radius 0.7 arcsec; it has been stretched out to have a major axis of 2.5 arcsec along the direction of motion, but still has a 0.7 arcsec minor axis perpendicular to that.

I think people haven't worried about this because it didn't matter so much for a long time. Few fast-movers were observed, and the positional uncertainty wasn't as low as can be achieved with current CCDs. But as astrometry has gotten better and faster objects are observed, the effects of poor timing have become more noticeable.

In some cases, things might be a lot worse. Somebody may have gotten an image of a meteor, for example, from which we can say: "The object can be seen at such-and-such RA/dec in the image, but we have only a vague idea as to when the image was taken." The Chelyabinsk meteor, caught on many dashboard cameras, springs to mind as an example.

This is still useful data. It still constrains the orbit (though not as well as a properly-timed observation). However, to make use of it, we have to include the uncertainty in the timing. I know of no other orbit software that includes timing uncertainty, though I very much doubt that doing so is original with me.

Ideally, we'd "fix" the problem by getting observers to time their observations correctly, but that's not as easy as it might seem. I think a lot of observers just make sure their computer's clock is set to within a few milliseconds, but remain unaware of timing issues caused by the fact that the "image time" is for the beginning or end of the exposure instead of its midpoint; or because there is some latency in the shutter; or maybe the "image time" is the time all the data finished coming in from the camera. For some professional setups, things just get worse: the shutter is big enough and takes long enough to open/close that the time is a function of where you are in the image. About the only way around this is to observe fast-moving objects with "known" ephemerides, and to check that your positions match. There are a few artificial satellites that fill the bill; I'll add comments about this later.

Note that even if we did get everybody to fix this problem, we would still have the issues of archived observations that were not properly timed, plus problems arising from meteor camera observations.

As mentioned above, one very serious weakness in how timing uncertainty is handled remains, and that is that we've assumed random timing errors. In reality, the error is usually strongly correlated over a night's observations: for example, the computer's clock might have been three seconds fast all night. (Or the software might store the time of the beginning or end of the exposure, or when the data was received, as described above.)

It's been suggested by several people that the "clock error" of a particular site could be included as an extra parameter to be solved for, much as one can solve for (say) asteroid masses or non-gravitational parameters as extra parameters. This is tempting, though I see possible problems with trying to solve for too many parameters with too little data.

Setting precise and non-standard-format times and RA/decs for observations

The current MPC format for optical astrometry allows one to specify the time to a precision of .000001 = 1e-6 day, or 86.4 milliseconds. Usually, that's sufficient. (In fact, most observations are given to only five decimal places, or .864 second precision). But certain NEOs have been moving so fast that it would be best to specify the time more precisely, especially for video-based observations.

Similarly, positions are usually given to .01 second in RA and 0.1 arcsecond in declination, with the option of using another digit if you're suitably confident in your data. On rare occasion, though, you may want more (or less) precision than this.

Find_Orb provides a couple of ways to get around these problems. For timing, I'd suggest using any of several alternative time formats. Essentially, instead of using the MPC's YYYY MM DD.ddddd format, you can specify the time in the same columns using different formats, thereby allowing you to give more digits and/or specify the time in HH:MM:SS.ssss form, or in JD or MJD form.

Alternatively, you can put a line of the following sort in front of the observation in question:

#time 2009 Mar 3.1415926535897
#time 2009 3 3 14:15:09.265358
#time 2454894.31415926 

Basically, one can reset the time to arbitrary precision, and with a somewhat arbitrary format (any of the formats used when specifying epochs or ephemeris start dates can be used; click here for details.) The subsequent observation will be assumed to have been made at the time specified, and the time stored in the MPC observation line will be ignored.

Yet another alternative is use of the .rwo format, which allows time to be specified to 10-10 day, or 8.64 nanoseconds.

A warning, though: the double-precision values (64-bit floats) used in Find_Orb have 53 bits of precision. Find_Orb handles times in the Julian Day system, so its precision for current-day events is about 25 microseconds. Anything beyond that just won't get handled correctly. And as you approach that level of accuracy, roundoff error will become a big issue. I don't really know at what point one would actually have to worry about all of this, but it might be well before approaching the 25 microsecond level. (There are a couple of ways of getting around this, if it ever becomes an issue. Find_Orb could switch to "long doubles" with 64 bits of precision, giving us a 2048-fold improvement. A precision of 12 nanoseconds should be enough for anybody. Times could be stored in days since 1 January 2000, instead of in the Julian Day system of days since noon on 1 January -4712; this would give us some added precision for dates/times in the current epoch.)

RA/decs can be stored in a variety of ways, to allow for input in decimal hours or decimal degrees, or for extra digits. Columns 33 to 56 of 80-column MPC-formatted data can contain any of the following. (The first two versions are the "usual" MPC format, and the "more precise" MPC format. The others are not official MPC formats, though the low-precision flavors have shown up for some ancient comet observations and for SOHO/STEREO data, which tend to be of coarse quality.) Here, 'z' = 'h'ours or 'd'egrees.

hh mm ss.sss            (MPC uses this for 'precise' RAs)
zz mm ss.ss             (MPC uses for most RAs & 'precise' decs)
zz mm ss.s              (MPC uses for most decs & low-precision RA)
zz mm ss                (Used by MPC, only rarely)
zz mm                   (Used _very_ rarely by MPC)
zz mm.m                 (Maybe used by MPC once or twice)
zz mm.mm
zz mm.mmm
zz mm.mmmm
zz mm.mmmmm
zz mm.mmmmmm
zz.                     (degrees or hours)
zz.z
zz.zz
zz.zzz
zz.zzzz              ... can go up to nine places in RA or
                        eigt places for declination
ddd.                   (used for RA only)
ddd.d
ddd.dd               ... can go up to eight places (again, RA only)
HHMMSSs               (RA to .1 second)
ddmmSSs               (dec to .1 arcsecond)
HHMMSSss              (RA to .01 second)
ddmmSSss              (dec to .01 arcsecond)
        ... and so forth until...
ddmmSSsssss           (dec to 10 microarcseconds)
HHMMSSssssss          (RA to one microsecond)

Find_Orb supports all these formats because I've had input data from sources such as Mars Express (astrometry for Deimos and Phobos) and VLBA observations of Saturn. That last is extremely precise, requiring the last two formats shown above (click here to see the result for the Saturn VLBA data). This precision will also be needed, I hope, for data from the Gaia mission. Gaia is supposed to reach 20 micro-arcsecond accuracy at magnitude 15.

Extending the time range for observations and ephemerides

By default, if you attempt to use Find_Orb to generate ephemerides for dates before the year -1000 or after the year +3000, you'll get an error message, and observations are assumed to be valid only for the years +1000 to 2020. Usually, dates outside those ranges indicate some sort of a blunder, and the user should be warned about them. But if you're actually doing something that requires larger time ranges, you should edit the text file environ.dat and look for these lines:

TIME_RANGE=-1000,3000 TIME_RANGE_OBS=1000,2020

Change those numbers, and you can revise the accepted time ranges. Note that Find_Orb is not really set up for long-range integrations; if you want to see what happens to an object over hundreds of thousands of years, you're pretty much out of luck. However, if you use the JPL DE-431 ephemeris (see the next section on how to do that), you can get good results over the entire range of that ephemeris, which is roughly -13000 to +17000.

Speeding Find_Orb up by using JPL DE ephemerides

By default, Find_Orb computes planetary and lunar positions using the PS-1996 and ELP series. These don't require much file space, but they can be quite slow. You'll never notice with most of these examples. But in some cases (very near-earth objects and such), you'll see that Find_Orb takes a long time to do its work. Most of that time is spent in computing planetary positions.

To speed things up, you can switch to use of the JPL DE ephemerides. These are the fundamental ephemerides from which series approximations such as PS-1996, ELP, and VSOP were derived. Use of DE is only marginally more precise (not many people will notice the difference), but frequently a lot faster, especially with longer orbital arcs that require more planetary positions to be computed.

You can learn a bit about JPL DE ephemerides here, and find out how to purchase or download (free) DE data here. (You can also find a DE file on the second Guide 8.0 CD-ROM; the Guide 9.0 DVD contains the file /jpl_eph/pl_eph.422 covering years -3000 to +3000). You will then have a file with an extension such as '.200', '.405', or '.4xx' (corresponding to the DE-200, DE-405, or DE-4xx ephemeris). In most cases, you can just put that file in your Find_Orb folder, and the program will automatically detect it and make use of it. The only change you are apt to notice is that the program will run faster; and when you start the program, there will be a line at the bottom of the dialog saying "Using DE-xxx", and the years covered by the ephemeris file.

The automatic detection relies on the file name being one of those given in jpl_eph.txt, which will usually be the case for any JPL ephemeris file you've downloaded or copied from the Willmann-Bell CD or from a Guide 8.0 or 9.0 disk. However, suppose you've downloaded the source files for a DE ephemeris and built it yourself, giving it a new name. Or suppose you've put the JPL ephemeris file in some other folder, perhaps because some other program is also using it. In such cases, you would edit environ.dat and add a line such as

JPL_FILENAME=d:\unix.405

Note: Linux users would instead, or also, add something like:

LINUX_JPL_FILENAME=/mnt/cdrom/unix.405

Incidentally, if you wind up solving orbits that run outside the span of the DE file you've provided, Find_Orb will fall back on the PS-1996 and ELP solutions.

C/C++ source code for Find_Orb

You can click here to download the source code for a console-mode version of Find_Orb (about 241 KBytes.) This has been compiled and run successfully on Windows (with the Visual C/C++, OpenWatcom, and MinGW (gcc for Windows) compilers), and on Linux, and Mac OS/X using the GNU gcc compiler. You'll also need to get the source code for basic astronomical computations, also on this site; and you'll have to have downloaded the Windows version at the top of this page, since the DOS and Linux versions use many of the same files for tasks such as computing planetary positions, figuring out which observatory is which, and so on.

If you're compiling for DOS/Windows, the OpenWatcom version is simplest: build the basic astronomical functions library with wmake -f watmake and then use wmake -f watfind.mak to build find_orb.exe. Similarly, with the gcc compiler in Linux or Mac, one can use make -f linlunar.mak and make -f linmake. With the MinGW compiler, one would run make -f lunar.mng, then make -f makefile.mng.

On all platforms, two programs are compiled: find_orb, an interactive program which works much as the Windows software does, allowing you to modify the orbit and create ephemerides as you wish; and fo, a non-interactive program that reads an input file of astrometry, then computes orbital elements for each object without human intervention.

find_orb uses the curses library, which is either already included with most Linux distributions or easily installed. A simplified (but sufficient) version of the library for the OpenWatcom compiler is provided with find_sou.zip. For Microsoft C/C++ or MinGW, you'll need to get the PDCurses (Public Domain Curses) library, and build it. (Which is easy to do.)

For the Microsoft compiler, build the library of basic functions with nmake lunar.mak. Then run nmake dos_find.mak. You'll then be able to run, say, find_orb example or fo example to process the objects in the 'example' file.

Mac version: Eric Christensen did some work to get find_orb to compile on a Mac. The process is essentially that described for Linux. However, if you're running on an older, non-Intel Mac, you must be sure to set up the program to use JPL ephemerides. The default method used for planetary/lunar ephemerides in Find_Orb is byte-order sensitive (at least right now) and won't work on Mac (or other non-Intel-endian machines.) However, it's recommended that you use JPL ephemerides anyway, because it makes the program immensely faster.

Running find_orb: Once built, you can simply run, for example, find_orb example (in DOS) or ./find_orb example in Linux. The objects in the file you've specified will be listed, and you can use the cursor keys to select the object you're interested in. (If there's only one object, it's automatically selected.)

The overall function of the console version is very similar to that of the Windows version; it's just that things are done in text mode instead of graphically. find_orb uses the same underlying functionality as its Windows cousin; only the interface is different. Hit '?' (or any other key not understood by find_orb) and you'll get a list of keyboard commands, straight from the dos_help.txt file.

At some point, I'll probably post the Windows UI code. But posting source code for Charon will be a higher priority. If you'd like to see the Windows source "as is", e-mail me, and I'll send it out in all its current lack of glory (Windows GUI code is never glorious.)

Non-interactive (batch-processing) mode:

The second program compiled from the source code , fo, can read in all observations from a file of input astrometry, puzzle out which objects are present, and spit out orbital elements for them. This somewhat experimental feature works as follows. Suppose you have a file of astrometry, astrometry.txt . You would run the following command on the DOS or Linux or Mac command line:

./fo astrometry.txt 

fo would process each object, telling you its name and summary data, then finding its orbit. A very basic summary of its semimajor axis, eccentricity, and inclination would be shown on the console. The full elements would be written out in the one-line MPCORB format to the file mpc_fmt.txt. They would simultaneously be written out in the eight-line MPC format to the file elements.txt.

That's about all there is to it. This is still somewhat experimental, and there's still work to be done: for example, the code just puzzles out the initial solution, and if the arc is very long, it won't go through the sort of process the "Auto-Solve" button uses to get the full arc and to exclude outliers. Still, it's useful if one just has a slew of short-arc data to examine.

Some more things Find_Orb will eventually do

As is always the case, I have a wish list of things I'd like to see Find_Orb do someday:

Definitions of a few terms:

Cometary orbit. Many comets have orbits with eccentricities that are extremely close to 1, meaning their orbits are highly elongated ellipses (Hale-Bopp, for example) or closely resemble parabolas. Go into the Settings dialog and set the "constraint" of e=1, and Find_Orb will search for a cometary orbit. This can be especially helpful if you have only a few observations to work with; in such cases, telling Find_Orb that the orbit must be parabolic (cometary) makes it much more likely that it will converge to a solution. (This is usually a good idea for a preliminary orbit, where there are few observations. When you find that both a cometary and non-cometary orbit converge, but the non-cometary one has a significantly smaller RMS error, switch to the non-cometary.)

Elongation. The angle Sun-Earth-Target. When close to zero degrees, the object is in conjunction with the sun. When close to 180 degrees, the object is almost exactly opposite the sun. See also phase.

Epoch. In the case of an orbit that includes the effects of perturbing objects, the orbital elements will change slowly with time. Thus, elements eventually become "outdated", and you need to replace them with new elements with an epoch closer to the time of observation. That's why Find_Orb has a small edit box to select the desired epoch, as a Julian day .

Full step. Once you have found a good approximation to an object's orbit using (for example) the Herget method , you can find the "best fit" orbit (the one that best matches all the observations of the object) by taking a few full steps. The "full step" uses a mathematical tool known as the method of least squares.

Herget method. This is one pretty good method for computing an approximate orbit from a set of observations. It works by assuming that two observations (usually the first and last observations) were at known distances, and then working to improve the quality of those distances to decrease the residuals. Eventually, the "steps" in the Herget method lead to a minimum set of residuals, and you have a fair approximation to the object's true orbit. Once that happens, you can switch to taking full steps.

The more mathematically adventurous may be interested in this description of some of the math underlying the method of Herget , with particular emphasis on some of the ways it is both 'simplified' and made more powerful in Find_Orb.

Julian Day. (Not to be confused with the Julian calendar, which is a very different thing!) The Julian Day (JD) system is a means of counting off time as the number of days elapsed since noon, 1 Jan -4713. By simply using the number of days (including fractional days) elapsed from that instant, we avoid the complexities associated with the traditional day/month/year, hour/minute/second system. Also, by putting the starting instant in the distant past, over 6700 years ago, we avoid negative numbers.

Irregular satellite. Natural satellites tend to fall in one of two groups. The first, "regular" sort orbits close to a planet, in the plane of the equator of the planet, in the same direction the planet rotates (a prograde orbit). All known regular satellites have had their rotation tidally locked to their primary, so that (as with the earth's moon) they do not appear to rotate as seen from the primary. The assumption is that they were formed at the same time as the planet.

Irregular satellites, on the other hand, have orbits independent of the equator or rotation of the planet, and can be either prograde or retrograde. They also are quite distant from the primary. All four gas giants have such satellites.

Least squares. The method of least squares is a common mathematical tool, originated by C. F. Gauss sometime around 1800. It starts with the assumption that you have gathered some observations (in our case, of positions of a celestial object) and have a theory, or "model", to describe them (in our case, Newton's laws of motion). The differences between observed and computed values are called residuals; the method of least squares states that the "best" model will find a minimum of the sum of the squares of the residuals. (There's some hedging to account for the fact that some observations will be more precise than others, and an additional assumption that the distribution of errors resembles a normal distribution, also known as a "Gaussian distribution" or, more popularly, a "bell curve".)

Magnitude parameters. The magnitudes of asteroids and comets are predicted using quantities for each object known as "magnitude parameters". These are H and G for asteroids, with G (the "slope parameter") usually left set to .15, and M(T) or M(N) and K for comets. In this situation, K defaults to 10, and one can choose between M(T) (for total magnitudes) or M(N) (for nuclear magnitudes) in the Settings menu. You can also reset the K and G magnitude parameters if you wish, though this is probably not a good idea in general.

MOID. The Minimum Orbit Intersection Distance, or MOID, is the minimum distance between two orbits, neglecting perturbations. It's computed by assuming that two objects (usually the earth and an asteroid or comet) will continue swinging around in their orbits, and looking for the place where the distance between those orbits would be a minimum.

Sometimes, that minimum distance will be so small that it looks as if the two objects will collide, but that's not necessarily the case for two reasons. First, when object A reaches that point, object B will usually be somewhere else. Second, the computation ignores perturbations. So it's a good first estimate as to how close the object can get to you, but it takes a much more detailed analysis to determine if the object is going to hit you. A large MOID usually means the object is no hazard at all; a small one means you should keep an eye on it and find out if it really is a hazard.

Find_Orb shows the MOID not only between the asteroid/comet and the Earth (if it's less than one AU), but also the MOIDs between the object and all planets (if those MOIDs are small enough to be interesting.)

Perturbations. Objects orbiting the Sun would, in theory, trace paths of perfect ellipses, parabolas, or hyperbolas. However, the gravitational forces (a.k.a. perturbations) exerted by all the planets, satellites, and (less significantly) asteroids in the solar system make the actual motion much more complex. Find_Orb has a set of check-boxes that allow you to include (or, by default, not include) the perturbing effects of these objects.

Phase. The angle Sun-Target-Earth. When close to zero degrees, the object is almost fully illuminated by the sun. When close to 180 degrees, the object is about to transit the sun, and is almost totally un-illuminated. See also elongation.

P, Q, R vectors. MPC orbital elements, and the elements shown in Find_Orb, give the P and Q vectors in columns to the right of the 'standard' elements. The vectors are frequently used in computing positions from orbital elements. P is a unit vector giving the direction from the sun's center to the perihelion point of the object's orbit. Q is a unit vector at right angles to P, but in the plane of the object's orbit. R is a unit vector at right angles to both of these (and therefore, at right angles to the plane of the orbit). Thus, Q = R x P, R = P x Q, and P = Q x R, where 'x' is the vector cross-product operation. The benefit of all this is that, if we know the object is at true anomaly v and distance from the primary r, we can say that its position is

r * (P cos v + Q sin v) 

'R', the normal vector to the plane of the orbit, is almost never shown, since it does not appear in the above formula.

Residuals. When you attempt to determine an orbit for an object, you will always find some difference between the "observed" positions (measured from a CCD image), and the "computed" positions (computed using the orbit you have determined). These differences are known variously as "observed minus computed", or "O-C", or as "residual errors", or "residuals". Ideally, they should represent the random errors in observation that inevitably happen in the real world.

Retrograde. As viewed from a point far "above" the ecliptic, almost all planets and satellites orbit in a counterclockwise direction, also known as "prograde". A few satellites of the gas giants, however, go in the opposite direction, also known as "retrograde". All of these objects are quite small and faint.

RMS error. "RMS" = "Root mean square"; it's a modified version of the average of the residuals. (To be mathematically precise, it is the square root of the average value of the square of the residuals. If you took the squares of all the residuals, averaged them, and took the square root of the result, you would have the RMS error.) With CCD-based observations and with an accurate orbit, the RMS error should be about an arcsecond or better. If the errors are much greater than that, you either should continue working on the quality of the orbit, or check to see if there are problems with the observations.

Sphere of Influence. For an object in orbit around a planet, it makes sense to use planet-centric orbital elements. For an object far from any planet, it makes sense to use heliocentric elements. There is an almost spherical region around each planet where planet-centric elements are a better choice: that is, a two-body solution centered on the planet would better represent the object's motion than a two-body solution centered on the sun. That region is the planet's "sphere of influence."

You can override Find_Orb's judgment and insist on only heliocentric orbits; there is a check-box for that.

Danby shows that for a planet of mass Mp at distance R from the sun (of mass Ms), this sphere has radius r where

r = R (Mp/Ms)(2/5)

Following are the radii of the sphere of influence for most of the planets. (For the moon, we have to switch to r=R(Mmoon/Mearth)(2/5), where R is the earth-moon distance. So far, I've yet to see a "real" object within the moon's sphere of influence, and therefore haven't seen Find_Orb generate a selenocentric orbit... though note that one can be generated for the example case 1997 ZZ99.)

Planet    r(AU)      r(km)
Mercury   0.001     112000
Venus     0.004     614000
Earth     0.006     921000
Mars      0.004     575000
Jupiter   0.322   48000000
Saturn    0.365   54300000
Uranus    0.346   51600000
Neptune   0.579   86000000
Pluto     0.022    3300000
Moon      0.00046    69000

Tisserand criterion: The Tisserand criterion is sometimes used to see if an object might have been captured from parabolic orbit, usually by Neptune or Jupiter.

Väisälä orbit. A newly-discovered object may have only two known observations, or may have been observed over so short a period that its orbit is not at all well-defined (that is, you can't tell if it's a nearby object moving slowly or a more-distant object moving rapidly, and any possibility in that range seems equally likely.) In fact, this is usually the case for a newly-found object. When this happens, it's common to determine a Väisälä orbit: one in which it's assumed that the object is at a particular distance from the sun, and that it is moving at right angles to the sun (that is, it's either at perihelion or aphelion). Given that assumption, the orbit is precisely defined, and you can generate an ephemeris that will (usually) be good enough to recover the object a few nights later. Click here for details on how to determine a Väisälä orbit.