All the orbit determinations and the observation predictions used in this paper have been performed with the OrbFit software package. OrbFit has been developed by a consortium formed in 1996 by the groups led by myself, by M. Carpino (Obs. Milan/Brera), K. Muinonen (Univ. Helsinki) and Z. Knezevic (Obs. Belgrade). In the development of the portion of OrbFit under my responsibility, I have been supported by my students and former students L. Cattaneo, S. Baccili and A. La Spina, and by M.E. Sansaturio (Univ. Valladolid). This software package is distributed as free software (under a GNU type public licence), and can be obtained either by anonymous ftp (at the address given above) or by contacting any one of the authors. The capabilities of the OrbFit system, whenever they implement innovative algorithms, will be documented in publications by some of the members of the consortium, such as this one. In this appendix I will shortly describe some of the algorithms which have been implemented in OrbFit under my responsibility, but of course the software system would not work without other parts for which I can only thank the other members of the consortium.
To compute accurate orbit computations for natural solar system bodies, two main choices have to be done. First, we need to select a dynamical model, e.g. which bodies have to be included with their gravitational attraction, which masses etc. have to be used. Second, a numerical integration algorithm has to be selected, with all the relevant parameters, such as discretisation stepsize, convergence controls, order; for non self-starting schemes, a starter algorithm has also to be selected. Both sets of choices need to be done in a consistent way, namely there is no point in having a numerical accuracy much better than the model accuracy, and vice versa.
To satisfy the requirement of a dynamical model not only accurate but consistent, we use the well known JPL ephemerides as a source for the positions of all the planets (and the Moon), with all the constants (e.g. masses) extracted from the same source. For main belt asteroids, all the planets Mercury to Neptune are used as perturbing masses; Pluto is added for transneptunian orbits, and the Moon (as a separate mass from the Earth-Moon barycentre) can be optionally added (but it is relevant only for close approaches to the Earth). General relativistic perturbations are handled with the Einstein-Infield-Hoffman equation for the gravitational field of the Sun. The motion of the observer is also fully accounted, taking into account the motion of the Earth around the Earth-Moon barycentre, precession, nutation, and Earth rotation; aberration is also computed.
Perturbations from a few asteroids for which the mass is known can be optionally added, but in this case the source of the perturbing bodies positions is an ephemerides file generated by an auxiliary program.
The software package OrbFit incorporates three numerical integration schemes: (1) an arbitrary order Gauss-Jackson-Cowell multistep, (2) an arbitrary (even) order implicit Runge-Kutta-Gauss, and (3) the 15th order implicit Runge-Kutta-Radau. (1) is used whenever possible, that is for moderate eccentricity, non planet-crossing orbits, because it is by far the fastest (one function evaluation per step); the stepsize is chosen automatically to achieve a prescribed truncation error, according to the theory developed in [Milani and Nobili 1988]. (2) is an arbitrary order symplectic method, which is very accurate but slow, and is used essentially only as starter for (1). The Everhart method implemented in (3) is a well known public domain software [Everhart 1985], and is useful to handle the orbits where variable stepsize is mandatory. A special algorithm has been developed for close approaches, based upon [Milani et al. 1990] and [Baccili 1995]: the integration is slowed down to accurately compute the orbit and to record the closest approach.
The integration method and the dynamical model can be selected automatically, taking into account the different conditions of main belt, Mars crossing, Earth crossing, Jupiter crossing, Trojans, Edgeworth-Kuiper Objects, and Centaurs. The most delicate adjustment has been to select an order, and therefore a stepsize, in such a way that a main belt asteroid orbit can be computed with the multistep (1) even if Mercury is introducing a comparatively high frequency perturbation. A general relativistic correction is automatically introduced for Near Earth orbits only. Manual selection is possible for all the options, but it is very seldom necessary.
The package contains also differential correction algorithms which have been recycled from the satellite geodesy software previously developed by our group [Milani et al. 1995].
To give an idea of the performance, on a computer with a Pentium Pro 200 MHz processor, the determination of the multiple solutions used for Figure 15, including 624 iteration of single arc differential corrections, and the accurate computation of 128 orbits for 31 years, has required 143 seconds of CPU time. This was a main belt case; the CPU times can be significantly longer for an Earth-crossing orbit, requiring a shorter and variable stepsize; however, even the (719) Albert computations used for the Figures 6-7 have required only about one hour.
The package can also generate graphic output, by using the public domain software GNUPLOT, available on essentially all computers; all the Figures of this paper have been generated with OrbFit, though for some of them the output files have been graphically post-processed with Matlab.
Acknowledgements This research has been conducted as part of the overall project ITANET, funded by C.N.R. with the contributions 96.00302.CT02, 97.00026.CT02 and 97.05100.CT12. The data for Figure 17 have been supplied by T. Bowell (Lowell Observatory). The astrometric observations of the asteroids used in the examples of this paper have been obtained from the Minor Planet Center, by using the Extended Computer Service. During the preparation, and revision, of the manuscript I have been greatly helped and encouraged by discussions (either direct or remote, by E-mail) on the most important issues raised by this work with A. Harris, P. Farinella, C. Chapman, D. Morrison and with my coworkers of OrbFit. A. Boattini and M. Cavagna have suggested practical ways to implement recovery strategies. F. Giannessi has helped in the definition of the constrained optimisation algorithm of Section 5.1. The three referees (K. Muinonen, B. Marsden and D. Tholen) have contributed with thoughtful reports; Z. Knezevic and M.E. Sansaturio have also suggested corrections to the manuscript.