The earliest modern scientific model considered only the gravitational attraction between the Sun and each planet, with the resulting orbits being unvarying Keplerian ellipses.
In reality, all the planets exert slight forces on each other, causing slow changes in the shape and orientation of these ellipses.
Increasingly complex analytical models have been made of these deviations, as well as efficient and accurate numerical approximation methods.
VSOP was developed and is maintained (updated with the latest data) by the scientists at the Bureau des Longitudes in Paris.
An updated version, VSOP87, computed the positions of the planets directly at any moment, as well as their orbital elements with improved accuracy.
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.
Third order terms had to wait until the 1970s when computers became available and the vast numbers of calculations to be performed in developing a theory finally became manageable.
In VSOP87 especially these long period terms were addressed, resulting in much higher accuracy, although the calculation method itself remained similar.
[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.
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): Without perturbations these elements would be constant and are therefore ideal to base the theories on.
The planetary solution VSOP2013 is fitted to the numerical integration INPOP10a built at IMCCE, Paris Observatory over the time interval +1890...+2000.
Masses multiplied by the gravitational constant of the Sun, the planets and the five big asteroids are used values from INPOP10a.
[10] This solution is fitted to the numerical integration INPOP10a built at IMCCE (Paris Observatory) over the time interval +1890...+2000.