Earth to Jupiter distance with Skyfield

Question

I'm trying to use Skyfield to plot the distance in au from Earth to solar system planets as a function of time. This is super simple and is even given in the front page of the package homepage. However while this works perfectly well for mercury, venus and mars, it doesn't work for other planets. I'm not familiar with JPL ephemeris files but it seems that for example Jupiter has no key entry in the file de421.bsp which would explain the issue.

Here is a minimal example (the one from the homepage):

from skyfield.api import load, now

planets = load('de421.bsp')
earth, planet = planets['earth'], planets['jupiter']

jd = now()
position = earth.at(jd).observe(planet)
ra, dec, distance = position.radec()

print(distance)

The error is the below. Note that if you replace 'jupiter' by 'mars' in the code above, it doesn't crash.

---->  earth, planet = planets['earth'], planets['jupiter']
KeyError: "kernel 'de421.bsp' is missing 'JUPITER' - the targets it supports are:
SOLAR SYSTEM BARYCENTER, MERCURY BARYCENTER, VENUS BARYCENTER, EARTH BARYCENTER, 
MARS BARYCENTER, JUPITER BARYCENTER, SATURN BARYCENTER, URANUS BARYCENTER, 
NEPTUNE BARYCENTER, PLUTO BARYCENTER, SUN, MERCURY, VENUS, MOON, EARTH, MARS"

Am I using the ephemeris file in the wrong way (wrong barycenter ?) or is this just a limitation of the de421.bsp file ? I read the description of the ephemeris file on the Skyfield website (here) but not sure I fully understood it.

Any suggestion of how to perform this simple calculation of Earth-Jupiter distance with Skyfield ?

Thanks !

Eyeglasses

Have you tried planets['jupiter barycenter']? – Daladier

Answer 1

Like the error says, you need to use JUPITER BARYCENTER instead of jupiter.

Long

Thanks for the quick answer, it solved the problem. Just out of curiosity is there a difference between MARS BARYCENTER and MARS in the file ? or is it just an alias ? The distance value returned are identical. – Eyeglasses

@Eyeglasses The actual calculation, tabulation and interpolation of the JPL ephemerides is a little complicated, but after reading some of the links (next comment), I think the earth-moon barycenter (center of mass) is treaded as a "planet" in the solar system calculations, and then the relative motion between the earth and moon is treated as a somewhat separate calculation. For earth-moon they are also broken out separately for us, but for the Jupiter barycenter (the planet plus all it's moons) it's not broken out separately in Skyfield. Jupiter is so big that the difference is very small. – Jimenez

See for example this, this, and this – Jimenez

So the case of Pluto and Charon is particularly interesting - the barycenter is in space between the two bodies. – Jimenez

Answer 2

This is just supplemental if it's helpful - the accepted answer solved the problem.

I wanted to show that the since the positions are in barycentric coordinates, that the ['solar system barycenter'] would remain at the origin. But I was foiled, because it returns a single value of zero, instead of a vector (or None). Anyway

import matplotlib.pyplot as plt
from skyfield.api import load, JulianDate

data = load('de421.bsp')
sun  = data['sun']
bary = data['solar system barycenter']

years = [1975+i for i in range(51)]
sunpos, barypos = [], []

for year in years:
    jd = JulianDate(utc=(year, 1, 1))
    sunpos.append(sun.at(jd).position.km)
    barypos.append(bary.at(jd).position.km)

plt.figure()
x, y, z = zip(*sunpos)
plt.plot(years, x)
plt.plot(years, y)
plt.plot(years, z)
# x, y, z = zip(*barypos)
# plt.plot(years, x, '-k')
# plt.plot(years, y, '-k')
# plt.plot(years, z, '-k')

plt.title('suns motion in barycentric frame')
plt.savefig('bary one')
plt.show()

The bottom two plots (below) show the motion of the earth and moon relative to the earth-moon barycenter, called ['earth barycenter'] in Skyfield:

import matplotlib.pyplot as plt
import numpy as np
from skyfield.api import load, JulianDate

data  = load('de421.bsp')
earth = data['earth']
moon  = data['moon']
bary  = data['earth barycenter']

days = range(0, 366, 5)
earthpos, moonpos, barypos = [], [], []
for day in days:
    jd = JulianDate(utc=(2016, 1, day))  # seems to work
    earthpos.append(earth.at(jd).position.km)
    moonpos.append(moon.at(jd).position.km)
    barypos.append(bary.at(jd).position.km)
ep = np.array(earthpos).T
mp = np.array(moonpos).T
bp = np.array(barypos).T

plt.figure(figsize=[9,9])
plt.subplot(5,1,1)
for thing in ep:
    plt.plot(days, thing)
plt.subplot(5,1,2)
for thing in mp:
    plt.plot(days, thing)
plt.subplot(5,1,3)
for thing in bp:
    plt.plot(days, thing)
plt.subplot(5,1,4)
for thing in (ep-bp):
    plt.plot(days, thing)
plt.subplot(5,1,5)
for thing in (mp-bp):
    plt.plot(days, thing)
plt.savefig('bary two')
plt.show()

Jimenez

Thanks for this detailed insight ! I now understand better the 'barycenter' in the JPL files. – Eyeglasses