Notebook

Introduction to the VSOP2013 planetary theory¶

Starting from version 0.15, heyoka.py provides an implementation of the VSOP2013 planetary theory. VSOP2013 consists of a set of (large) analytical formulae which calculate the heliocentric orbital elements of the 8 planets of the Solar System (plus the dwarf planet Pluto) as a function of time. Because the VSOP2013 theory is expressed as a set of analytical, differentiable formulae, it can be implemented in heyoka.py's expression system and it can thus be used in the formulation of differential equations.

API overview¶

heyoka.py provides three functions to generate the analytical formulae for the VSOP2013 theory.

The first one, vsop2013_elliptic(), returns the formulae for the heliocentric elliptic orbital elements for a given planet. The 6 elliptic orbital elements are, in order:

the semi-major axis $a$ (in astronomical units),
the mean longitude $\lambda$ (in radians),
$k=e\cos\varpi$ and $h=e\sin\varpi$, where $e$ is the eccentricity and $\varpi=\Omega+\omega$ is the longitude of the perihelion,
$q=\sin\frac{i}{2}\cos\Omega$ and $p=\sin\frac{i}{2}\sin\Omega$, where $i$ is the inclination and $\Omega$ is the longitude of the ascending node.

This function requires as input the planet index (from 1 to 9) and the orbital element index (from 1 to 6). The orbital elements are referred to the inertial frame defined by the dynamical equinox and ecliptic J2000. Note that the VSOP2013 solution for the Earth (planet index 3) actually refers to the Earth-Moon barycentre.

The second function, vsop2013_cartesian(), returns the heliocentric Cartesian state vector $\left( x, y, z, v_x, v_y, v_z \right)$ for a given planet. The Cartesian state is computed via a direct conversion of the orbital elements returned by vsop2013_elliptic(), and thus it also refers to the inertial frame defined by the dynamical equinox and ecliptic J2000.

The third function, vsop2013_cartesian_icrf(), returns the heliocentric Cartesian state vector for a given planet in the International Celestial Reference Frame (ICRF). The ICRF state vector is computed by applying a rotation to the output of vsop2013_cartesian().

Precision and truncation threshold¶

The VSOP2013 solution is formulated as large Fourier/Poisson series whose numerical coefficients decay exponentially. When using the full VSOP2013 solution, the precision with respect to high-accuracy numerical integrations is roughly sub-meter for the inner planets (except Mars) and kilometers for the outer planets (except Pluto). By dropping the smaller terms from the series, the precision of the theory degrades but computational efficiency greatly increases. See the VSOP2013 readme for more details and estimates of the precision at various truncation levels.

By default, all vsop2013_*() functions filter out from the solution terms whose numerical coefficients have a magnitude smaller than $10^{-9}$. The truncation threshold can be changed via the thresh keyword argument:

In [1]:

import heyoka as hy

for thr in [1e-7, 1e-9, 1e-11]:
    print("Size of Mars' semi-major axis solution @ {}: {}".format(thr, len(hy.vsop2013_elliptic(4, 1, thresh = thr))))

Size of Mars' semi-major axis solution @ 1e-07: 4916
Size of Mars' semi-major axis solution @ 1e-09: 48574
Size of Mars' semi-major axis solution @ 1e-11: 467781

Time coordinate¶

The precision of the VSOP2013 solution is highest betwen $1890\,\mathrm{CE}$ and $2000\,\mathrm{CE}$, and it slowly degrades outside this time range. See the VSOP2013 readme for estimates of the precision degradation up to a few thousand years in the past/future.

VSOP2013 uses barycentric dynamical time (TDB), expressed in thousands of Julian years from J2000 (JD2451545.0). Note than one Julian year = 365250 days.

By default, all vsop2013_*() functions use hy.time to represent the time variable in the VSOP2013 formulae. This means that, by default, when the VSOP2013 solution is used in an ODE system, time is assumed to be measured in thousands of Julian years of TDB and $t = 0$ corresponds to the Julian date $2451545.0$.

It is possible to change the expression used to represent the time variable in the VSOP2013 solution via the time keyword argument. This allows, for instance, to rescale and change the origin of time in the VSOP2013 formulae. Let us see a few examples:

In [2]:

print("Mars' semi-major axis solution, threshold = 6e-5, default time expression:\n{}\n".format(hy.vsop2013_elliptic(4, 1, thresh = 6e-5)))
print("Mars' semi-major axis solution, threshold = 6e-5, time represented by variable 'x':\n{}\n".format(hy.vsop2013_elliptic(4, 1, time=hy.expression("x"), thresh = 6e-5)))
print("Mars' semi-major axis solution, threshold = 6e-5, time rescaled by a factor of 100:\n{}".format(hy.vsop2013_elliptic(4, 1, time=hy.time / 100., thresh = 6e-5)))

Mars' semi-major axis solution, threshold = 6e-5, default time expression:
(((2.1017753733317330e-07 * sin(((2.0000000000000000 * (6.2035000141410004 + (3340.6124341454570 * t))) + (-2.0000000000000000 * (0.59954610703500000 + (529.69096156232501 * t)))))) + (6.6017042784613703e-05 * cos(((2.0000000000000000 * (6.2035000141410004 + (3340.6124341454570 * t))) + (-2.0000000000000000 * (0.59954610703500000 + (529.69096156232501 * t))))))) + 1.5236793402339999)

Mars' semi-major axis solution, threshold = 6e-5, time represented by variable 'x':
(((2.1017753733317330e-07 * sin(((2.0000000000000000 * (6.2035000141410004 + (3340.6124341454570 * x))) + (-2.0000000000000000 * (0.59954610703500000 + (529.69096156232501 * x)))))) + (6.6017042784613703e-05 * cos(((2.0000000000000000 * (6.2035000141410004 + (3340.6124341454570 * x))) + (-2.0000000000000000 * (0.59954610703500000 + (529.69096156232501 * x))))))) + 1.5236793402339999)

Mars' semi-major axis solution, threshold = 6e-5, time rescaled by a factor of 100:
(((2.1017753733317330e-07 * sin(((2.0000000000000000 * (6.2035000141410004 + (3340.6124341454570 * (t / 100.00000000000000)))) + (-2.0000000000000000 * (0.59954610703500000 + (529.69096156232501 * (t / 100.00000000000000))))))) + (6.6017042784613703e-05 * cos(((2.0000000000000000 * (6.2035000141410004 + (3340.6124341454570 * (t / 100.00000000000000)))) + (-2.0000000000000000 * (0.59954610703500000 + (529.69096156232501 * (t / 100.00000000000000)))))))) + 1.5236793402339999)

Checking the implementation¶

In order to check the correctness of heyoka.py's implementation of the VSOP2013 theory, we are going to use the NASA HORIZONS service via the Astroquery package. Our objective is to compare the planetary positions computed by VSOP2013 at different threshold levels with those computed by HORIZONS. On the heyoka.py side, we will be using compiled functions to evaluate the VSOP2013 series.

In this example, we will be using Venus as a target body. Because the ecliptical coordinate system used by VSOP2013 does not match exactly any of the coordinate systems available in HORIZONS, we will be using the vsop2013_cartesian_icrf() function instead, which outputs heliocentric Cartesian coordinates in the standard ICRF (which is available in HORIZONS). Let's take a look at the code and at the results:

In [3]:

%matplotlib inline
from matplotlib.pylab import plt
import numpy as np
from astroquery.jplhorizons import Horizons

# Setup the plot.
fig = plt.figure(figsize=(9, 6))
plt.xscale('log')
plt.yscale('log')

# Create the heyoka variables.
x, y, z, tm = hy.make_vars("x", "y", "z", "tm")

# Threshold levels at which we will be computing
# the VSOP2013 solution.
thr_values = [1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10]

# We will perform 10 evaluations of the
# position of Venus in the 10 years following
# J2000.
dates = np.linspace(2451545.0, 2451545.0 + 365 * 10, 10)

for thr in thr_values:
    # Build the VSOP2013 formulae for the heliocentric Cartesian
    # position of Venus (planet index 2) in the ICRF. Replace the
    # default time variable with the "tm" symbolic variable and set
    # a custom threshold level.
    venus_x, venus_y, venus_z = hy.vsop2013_cartesian_icrf(2, time=tm, thresh=thr)[:3]

    # Compile the function for the evaluation of venus_x/y/z.
    venus_cf = hy.make_cfunc([venus_x, venus_y, venus_z])
    
    # Run the evaluation.
    # NOTE: the default time coordinate
    # for VSOP2013 is thousands of Julian
    # years from J2000.
    vsop_states = venus_cf((dates.reshape((1,-1)) - 2451545.0) / 365250)

    # HORIZONS query:
    # - 299 is the ID of Venus,
    # - 500@10 is the observer's location,
    #   the centre of the Sun.
    q = Horizons('299', location='500@10', epochs=dates)
    # Convert the HORIZONS query to the ICRF
    # and extract the cartesian position vector.
    tab = q.vectors(refplane='earth')
    hor_x = tab['x'].quantity[:].value
    hor_y = tab['y'].quantity[:].value
    hor_z = tab['z'].quantity[:].value
    
    hor_states = np.array([hor_x, hor_y, hor_z])

    # Plot the difference (in km) between the positions computed
    # by VSOP2013 and HORIZONS.
    plt.scatter([thr]*len(dates), np.linalg.norm(vsop_states - hor_states, axis=0) * 1.496e8,
               alpha=.5, label="Threshold = {}".format(thr))

# Finish setting up the plot.
plt.legend()
plt.xlabel("Threshold value")
plt.ylabel("Error wrt HORIZONS (km)")
plt.tight_layout();

We can see how VSOP2013 matches up well with respect to HORIZONS, and how lower truncation levels indeed increase the precision of the VSOP2013 solution (up to an error of hundreds of meters for a threshold of $10^{-10}$).

Limitations and caveats¶

A possible issue in astrodynamical applications is that VSOP2013 does not directly provide formulae for the state vector of the Sun. The barycentric position of the Sun can be estimated via the heliocentric planetary formulae after imposing that the system's barycentre is fixed in the origin. It is not clear at this time if this approach is 100% accurate (as VSOP2013 includes perturbations on the planets by asteroids, whose state vectors are not provided in the solution), but it should be adequate for most practical purposes.

A better alternative is to formulate the dynamical equations of motion (e.g., of a spacecraft or a small body) in a heliocentric reference frame (i.e., in the framework of the so-called (N+1)-body problem).