1. History

Predicting the position of the planets in the sky was already performed in ancient times. Careful observations and geometrical calculations produced a model of the motion of the Solar System known as the Ptolemaic system, which was based on an Earth-centered system. The parameters of this theory were improved during the Middle Ages by Indian and Islamic astronomers.

The work of Tycho Brahe, Kepler, and Isaac Newton in early modern Europe laid a foundation for a modern heliocentric system. Future planetary positions continued to be predicted by extrapolating past observed positions as late as the 1740 tables of Jacques Cassini.

The problem is that, for example, the Earth is not only gravitationally attracted by the Sun, which would result in a stable and easily predicted elliptical orbit, but also in varying degrees by the Moon, the other planets and any other object in the solar system. These forces cause perturbations to the orbit, which change over time and which cannot be exactly calculated. They can be approximated, but to do that in some manageable way requires advanced mathematics or very powerful computers. It is customary to develop them into periodic series which are a function of time, e.g. (a+bt+ct²+...)×cos(p+qt+rt²+...) and so forth one for each planetary interaction. The factor a in the preceding formula is the main amplitude, the factor q the main angular velocity, which is directly related to a harmonic of the driving force, that is a planetary position. For example: q= 3×(length of Mars) + 2×(length of Jupiter). (The term 'length' in this context refers to the ecliptic longitude, that is the angle over which the planet has progressed in its orbit in unit time, so q is an angle over time too. The time needed for the length to increase over 360° is equal to the revolution period.)

It was Joseph Louis Lagrange in 1781, who carried out the first serious calculations, approximating the solution using a linearization method. Others followed, but it was not until 1897 that George William Hill expanded on the theories by taking second order terms into account. Third order terms had to wait until the 1970s when computers became available and the vast amounts of calculations to be performed in developing a theory finally became manageable.

2. Variations Séculaires Des Orbites Planétaires

2.1. VSOP82

Pierre Bretagnon completed a first phase of this work by 1982 and the results of it are known as VSOP82. But because of the long period variations, his results are expected not to last more than a million years (and much less, maybe 1000 years only on very high accuracy).

A major problem in any theory is that the amplitudes of the perturbations are a function of the masses of the planets (and other factors, but the masses are the bottlenecks). These masses can be determined by observing the periods of the moons of each planet or by observing the gravitational deflection of spacecraft passing near a planet. More observations produce greater accuracy. Short period perturbations (less than a few years) can be quite easily and accurately determined. But long period perturbations (periods of many years up to centuries) are much more difficult, because the timespan over which accurate measurements exist is not long enough, which may make them almost indistinguishable from constant terms. Yet it is these terms which are the most important influence over the millennia.

Notorious examples are the great Venus term and the Jupiter–Saturn great inequality. Looking up the revolution periods of these planets, one may notice that 8×(period of Earth) is almost equal to 13×(period of Venus) and 5×(period of Jupiter) is about 2×(period of Saturn).

A practical problem with the VSOP82 was that since it provided long series only for the orbital elements of the planets, it was not easy to figure out where to truncate the series if full accuracy was not needed. This problem was fixed in VSOP87, which provides series for the positions as well as for the orbital elements of the planets.

2.2. VSOP87

In VSOP87 especially these long period terms were addressed, resulting in much higher accuracy, although the calculation method itself remained similar. VSOP87 guarantees for Mercury, Venus, Earth-Moon barycenter and Mars a precision of 1" for 4000 years before and after the 2000 epoch. The same precision is ensured for Jupiter and Saturn over 2000 years and for Uranus and Neptune over 6000 years before and after J2000.^[1] This, together with its free availability has resulted in VSOP87 being widely used for planetary calculations; for example, it is used in Celestia and Orbiter.

Another major improvement is the use of rectangular coordinates in addition to the elliptical. In traditional perturbation theory it is customary to write the base orbits for the planets down with the following six orbital elements (gravity yields second order differential equations which result in two integration constants, and there is one such equation for each direction in three-dimensional space):

• a semi-major axis
• e eccentricity
• i inclination
• Ω longitude of the ascending node
• ω argument of perihelion (or longitude of perihelion ϖ = ω + Ω)
• T time of perihelion passage (or mean anomaly M)

Without perturbations these elements would be constant and are therefore ideal to base the theories on. With perturbations they slowly change, and one takes as many perturbations in the calculations as possible or desirable. The results are the orbital element at a specific time, which can be used to compute the position in either rectangular coordinates (X,Y,Z) or spherical coordinates: longitude, latitude and heliocentric distance. These heliocentric coordinates can then fairly easily be changed to other viewpoints, e.g. geocentric coordinates. For coordinate transformations, rectangular coordinates (X,Y,Z) are often easier to use: translations (e.g. heliocentric to geocentric coordinates) are performed through vector addition, and rotations (e.g. ecliptic to equatorial coordinates) through matrix multiplication.

VSOP87 comes in six tables:

• VSOP87 Heliocentric ecliptic orbital elements for the equinox J2000.0; the 6 orbital elements, ideal to get an idea of how the orbits are changing over time
• VSOP87A Heliocentric ecliptic rectangular coordinates for the equinox J2000.0; the most useful when converting to geocentric positions and later plot the position on a star chart
• VSOP87B Heliocentric ecliptic spherical coordinates for the equinox J2000.0
• VSOP87C Heliocentric ecliptic rectangular coordinates for the equinox of the day; the most useful when converting to geocentric positions and later compute e.g. rise/set/culmination times, or the altitude and azimuth relative to your local horizon
• VSOP87D Heliocentric ecliptic spherical coordinates for the equinox of the day
• VSOP87E Barycentric ecliptic rectangular coordinates for the equinox J2000.0, relative to the barycentre of the solar system.

The VSOP87 tables are publicly available and can be retrieved from VizieR.^[2]

2.3. VSOP2000

VSOP2000 has an accuracy that is a factor of 10-100 better than its predecessors. The uncertainty for Mercury, Venus and the Earth is reported to be around 0.1 mas (milliarcsecond) for the interval 1900–2000, and that for the other planets a few milliarcseconds.^[3] The publication of and the data for VSOP2000 are publicly available.^[4]

2.4. VSOP2002

Bretagnon's last work was on the implementation of relativistic effects, which was supposed to improve the accuracy with another factor of 10. This version was never finished, and still had weaknesses for Uranus and Neptune.^[5]

2.5. VSOP2010

The VSOP2010 files contain the series of the elliptic elements for the 8 planets Mercury, Venus, Earth-Moon barycenter, Mars, Jupiter, Saturn, Uranus, Neptune and for the dwarf planet Pluto. The VSOP2010 solution is fitted to the DE405 numerical integration over the time interval +1890...+2000.^[6] The numerical precision is 10 times better than VSOP82. Over a greater interval −4000...+8000 a comparison with an internal numerical indicates that the VSOP2010 solutions are about 5 times better than VSOP2000 for the telluric planets and 10 to 50 times better for the outer planets.^[7]

2.6. VSOP2013

The VSOP2013 files contain the series of the elliptic elements for the 8 planets Mercury, Venus, Earth-Moon barycenter, Mars, Jupiter, Saturn, Uranus, and Neptune and for the dwarf planet Pluto of the solution VSOP2013. The planetary solution VSOP2013 is fitted to the numerical integration INPOP10a built at IMCCE, Paris Observatory over the time interval +1890...+2000.^[8]

The precision is of a few 0.1″ for the telluric planets (1.6″ for Mars) over the time interval −4000...+8000.^[9]

2.7. Theory of the Outer Planets

This is an analytical solution for the (spherial and rectangular) positions (rather than orbital elements) of the four planets Jupiter, Saturn, Uranus, and Neptune and the dwarf planet Pluto.

TOP2010

This solution is fitted to the Ephemeris DE405 over the time interval +1890...+2000. The reference system in the solution TOP2010 is defined by the dynamical equinox and ecliptic J2000.0.^[10]

TOP2013

This solution is fitted to the numerical integration INPOP10a built at IMCCE (Paris Observatory) over the time interval +1890...+2000. The reference system in the solution TOP2013 is defined by the dynamical equinox and ecliptic of J2000.0.^[11]

The TOP2013 solution is the best for the motion over the time interval −4000...+8000. Its precision is of a few 0.1″ for the four planets, i.e. a gain of a factor between 1.5 and 15, depending on the planet, compared to VSOP2013. The precision of the theory of Pluto remains valid up to the time span from 0 to +4000.^[12]