Planetary Positions with VSOP87

Written by User 1

Last Updated: August 25, 2019, 03:05 pm (UTC)
Originally created on July 25, 2012

VSOP87 provides a method for computing the positions of the 8 planets (and the Sun) efficiently and accurately without the major headaches that astronomers from past centuries had to deal with.

Ever since ancient times, humans have observed and tracked the movements of the 5 brightest planets -- Mercury, Venus, Mars, Jupiter and Saturn -- in addition to the Sun and Moon. Today, we know of two additional planets -- Uranus and Neptune -- in addition to a wealth of other mobile celestial bodies including asteroids and comets. As technology improved, so did our ability to predict the movements of these objects. Today, using advanced algorithms like VSOP87, we can predict the positions of the planets to within a ten-thousandth of a degree thousands of years into the future.

The VSOP87 theory and solutions, by Pierre Bretagnon and Gerard Francou of Bureau des Longitudes, are likely the most used and most accurate algorithms available today for determining the positions of the planets without using interpolation. It consists of a large number of periodic terms that are then added up together in a special way to produce the 3-dimensional heliocentric coordinates of any planet at any moment in time for thousands of years into the future and the past. These coordinates can then be converted through a few transformations into geocentric coordinates which can be used to show their position as seen from Earth. As the name might suggest, VSOP87 was developed in 1987 -- a long time ago by technological standards and is based off of the DE200 ephemerides which have since been obsoleted. Nonetheless, for all practical purposes, the error originating in VSOP87 is negligible and unnoticeable by all but the most powerful telescopes.

Despite this apparent usefulness, there is little documentation or examples available online besides a few user implementations and the provided Fortran code. Looking at these, it is hardly obvious what they are doing unless you are either familiar with VSOP87 already or have significant experience with Fortran. For novice programmers, the apparent complexity of the code may be enough for frighten them away for good. That's why for the past several decades, the go-to resource for learning to use VSOP87 has been Jean Meeus' renowned book, Astronomical Algorithms. Even well into the age of the internet, it effectively remains the only resource available for learning to use VSOP87. Despite this, the algorithm is a relatively simple one -- one that does not require an entire book to understand if it is all that is needed. This page will attempt to clarify any confusing aspects of the algorithm so that anyone with a computer can quickly and easily implement it.

Getting Started

The VSOP87 data files are available from the VSOP87 FTP site and are organized by solution and series. There are several VSOP87 solutions designed for various purposes. They are:

  • VSOP87 - orbital elements, of J2000.0

  • VSOP87A - rectangular, of J2000.0

  • VSOP87B - spherical, of J2000.0

  • VSOP87C - rectangular, equinox of date

  • VSOP87D - spherical, equinox of date

  • VSOP87E - rectangular, barycentric, of J2000.0

All of the solutions are centered on the Sun except VSOP87E which is centered on the solar system barycenter.

For determining the positions of the planets in the sky, it is most convenient to use VSOP87C which produces rectangular coordinates that are easy to transform (Meeus uses a truncated version of VSOP87D in his book). As a result, this page will concentrate on using this particular solution though it will briefly touch on applying the same concepts to the other solutions.

To get started, download the necessary data files from the site -- it is necessary to have a copy of the data for every planet whose position is to be computed for the selected solution. In addition, if calculating the apparent position of the planet from Earth is wanted, then the data for Earth is also needed.

Understanding the Data

The data for each planet is divided into three sections labeled by three different letters: X, Y and Z for VSOP87C, and L, B and r for VSOP87D. Each of these sections are computed individually to obtain the coordinate represented by that section. For VSOP87C, the X, Y and Z form the heliocentric rectangular (Cartesian) coordinates of each planet, while for VSOP87D, the L (ecliptic longitude), B (ecliptic latitude), r (distance) form the heliocentric spherical coordinates of each planet. In the data files, each of these are labelled as "Variable 1" (X or L), "Variable 2" (Y or B) and "Variable 3" (Z or r).

Within each section, the data is further subdivided into individual series of terms generally referred to as X0, X1, X2,... (or X0, X1, X2,...) though are labelled in the data as *T**0, *T**1, *T**2,... Each of these series consists of a list of sets of 3 coefficients (the 3 rightmost numbers in the data) generally referred to as A, B and C. Each set of coefficients is used to calculate one term of the series using a particular formula. The terms of each series are then summed together. These series sums are then put in another formula to obtain the actual value of the coordinate. This is then repeated for each coordinate of the planet to obtain the heliocentric position of the data.

Mathematically, this is just a very big set of parametric equations that are functions of time although it usually is not written as such due readability concerns stemming from its tremendous size.

Using the Data: VSOP87C

This section will focus on VSOP87C for Venus and its X-coordinate. The formulas presented here are equally applicable for any coordinate of every other solution of VSOP87, except just with different numbers.

In the data file (VSOP87C.ven), the first set of 3 coefficients the X0 series are 0.72268045621, 3.17614669179 and 10213.52936369450. To find the term this set represents, use the following formula:

Term = A * cos(B + C * T)

Here, A, B and C are the coefficients above and T is the number of Julian millennia since J2000.0. It can be found with the following formula:

T = (JDE - 2451545) / 365250

JDE is the Julian Ephemeris Date (the Julian Date corresponding to Terrestrial Time) of the moment for which the coordinate is to be computed for.

To determine the first term for December 21, 2012 at 0:00 TT (JDE 2456282.5), we find:

T = (2456282.5 - 2451545) / 365250 = 0.012970568104
Term = 0.72268045621 * cos(3.17614669179 + 10213.52936369450 * 0.012970568104)
= -0.611163163207

Note that cos() should be set to accept radians for these terms to work.

Repeat this calculation using the same value of T with the remaining coefficients of the series to obtain the other terms of the series X0. When complete, take the sum of the terms and save this value as X0. Doing so, we find that X0 = -0.60495199432. Repeat this calculation for the remaining series: X1, X2, X3 and X4.

The results (for December 21, 2100 at 0:00 TT, as above) should be:

X0 = -0.604951994318
X1 = -0.000475966265
X2 = +0.000208557920
X3 = +0.000000953670
X4 = -0.000000348009

These values can then be used to determine the actual value of X with the following formula:

X = X0 + X1 * T + X2 * T2 + X3 * T3 + X4 * T4

This formula explains the reasoning behind the labels for each of these series (*T**0, *T**1, *T**2,...) as they indicate what each series must be multiplied by to obtain the actual value for each series. In other words, to obtain the term that includes X4 (as calculated above), labeled as *T**4, use X4 * T**4. The formula above is effectively the sum of each of these terms. (Note that "term" here refers to the terms in the formula above, not the terms that form each series as was used previously)

Plugging in the values of T and X0 through X4 into the formula, we obtain X = -0.604958132783. Units for VSOP87C is the astronomical unit (AU) so if the position of Venus is plotted in a Cartesian space centered on the Sun, its X coordinate would be at -0.604958132783 AU.

To find the full set of coordinates, repeat this with the Y and Z series using the same value of T obtained above.

Notes on Using Spherical Coordinates

The procedure illustrated above is identical when using the VSOP87A solution which also returns rectangular coordinates. However, with VSOP87B, VSOP87D and VSOP87E which use spherical coordinates, there are a few considerations to take note of that arise from the use of the different coordinate system.

The first important note is that result of L (ecliptic longitude) and B (ecliptic latitude) obtained through the method outlined above will output values in radians. In astrodynamics, units of radians are rarely (if ever) used, so the values should first be converted to degrees by multiplying them by 180 / Pi.

In addition, the values should be also modulated to within the proper range before being used -- longitude is from 0 to 360 degrees, latitude from -90 to 90 degrees. Otherwise, the algorithm may output values on the order of several millions of degrees which is indecipherable to humans and may cause issues when given to computers.

Beyond these added considerations, the procedures for obtaining L, B and r are identical to the one used to obtain the X coordinate for Venus above.

Finding Geocentric Coordinates

The results of VSOP87C above are heliocentric coordinates, centered on the Sun. Since none of us actually live on the Sun, it might be more useful to obtain geocentric coordinates, centered on the Earth. Using rectangular coordinates, a simple transformation is all that is needed to convert.

X = Xobj - XEarth
Y = Yobj - YEarth
Z = Zobj - ZEarth

Determine the heliocentric rectangular coordinates of the planet to locate (Xobj, Yobj and Zobj) as well as those of Earth (XEarth, YEarth and ZEarth) and subtract. This will now set the position of every planet in a space where the origin is the position of Earth. To verify this, set Xobj, Yobj and Zobj to XEarth, YEarth and ZEarth and, as expected, the new position of Earth is (0, 0, 0), the origin.

To determine the geocentric position of the Sun, set Xobj, Yobj and Zobj all to zero (since it was at the origin in the heliocentric system). The formula for the Sun can then be simplified to:

X = -XEarth
Y = -YEarth
Z = -ZEarth

To transform this into spherical Ecliptic coordinates, use:

tan(Lon) = Y / X
tan(Lat) = Z / sqrt(X2 + Y2)

Be sure to obtain these to the correct quadrants either manually or by using atan2(). These can then be quickly converted to equatorial coordinates using the typical ecliptic-to-equatorial formulas.

View other pages under Computations

Copyrighted © 2007-2022, The Caglow Project.
Material is available under AL.