Interfacing to SPICE and JPL DE ephs

A pykep.planet is able to represent, with a unique interface, objects moving through space. and we here showcase the case of asteroids/comets and satellites whose motion is not integrated numerically, but rather fitted to observations, simulations or historical data and encapsulated into NAIF SPICE kernels.

Under the hoods, the SPICE C code made available by NAIF, is used and interfaced to the pykep relevant classes.

import pykep as pk
import numpy as np

Let us compute the position of the barycenter of the Jupiter system with respect to the Solar system barycenter. We do this using the JPL DE440 ephemerides which, released in 2022, are considered to be the most accurate ones.

First we must download the binary .bsp file that contains the corresponding data. In the case of the DE440 ephemerides pykep actually includes such a file in its distribution (for convenience). We may thus, for this case, skip this step and get the file path directly:

de440s_kernel = pk.udpla.de440s.kernel_file()
print(de440s_kernel)

/home/runner/local/lib/python3.13/site-packages/pykep/data/de440s.bsp

If we have no idea on what objects are contained in this kernel, we may get their NAIF IDS using one of the utilities in pykep.

naifids = pk.utils.inspect_spice_kernel(de440s_kernel)

and convert into more readible names as follows:

names = [pk.naifid2name(item) for item in naifids]
print(names)

['MERCURY BARYCENTER', 'VENUS BARYCENTER', 'EARTH BARYCENTER', 'MARS BARYCENTER', 'JUPITER BARYCENTER', 'SATURN BARYCENTER', 'URANUS BARYCENTER', 'NEPTUNE BARYCENTER', 'PLUTO BARYCENTER', 'SUN', 'MERCURY', 'VENUS', 'MOON', 'EARTH']

Nice!, inspecting this list we realize that most solar system bodies are contained in this kernel (their barycenters) as well as some of the non barycentric positions. Now, lets proceed in our task to compute the Jupyter barycenter position. The first thing to do (following the common usage of SPICE) is to pre-load in memory the kernel file:

pk.utils.load_spice_kernels(de440s_kernel)

This is done only once and we can now forget about it (unless memory is an issue and we want to unload the data, in which case we provide a corresponding unload_spice_kernels() function).

We now instantiate a planet form the UDPLA spice.

udpla = pk.udpla.spice("JUPITER BARYCENTER", "ECLIPJ2000", "SSB")
jupiter = pk.planet(udpla)
when = pk.epoch("2023-01-02")
r, v = jupiter.eph(when)
print(f"Position (m): {r}\nVelocity (m/s): {v}")

Position (m): [722180808588.1804, 157535374702.5074, -16810696007.16372]
Velocity (m/s): [-2933.2858571285688, 13378.581606366935, 10.115066760074676]

And inspect it:

print(jupiter)

JUPITER BARYCENTER - SPICE

Note how many physical parameters of the body are not defined, as for this particular UDPLA (the pykep.udpla.spice), the ephemerides computations are interpolated from tables and hence do not require any physical parameters. The user can always define them if needed.

The DE 440 JPL Ephemerides

What we have presented above is valid in general and works with any SPICE kernel, be it a comet, a rover a spacecraft. Clearly the backdraw is that we have to first download and position the correct kernel file and then load it into memory for it to work.

Since most of times, or anyway the most common usage pattern, is to query the position of the solar system planets, in pykep we offer a dedicated udpla for that, one that avoid loading kernels in memory. Lets start unloading the current kernel as to not create confusion.

pk.utils.unload_spice_kernels(de440s_kernel)

The UDPLA pykep.udpla.de440s automatically loads the correct kernel shipped with pykep module.

udpla = pk.udpla.de440s("JUPITER BARYCENTER", "ECLIPJ2000", "SSB")
jupiter = pk.planet(udpla)
when = pk.epoch("2025-01-02")
r, v = jupiter.eph(when)
print(f"Position (m): {r}\nVelocity (m/s): {v}")

Position (m): [156005590351.0843, 743270596831.1477, -6573233296.777874]
Velocity (m/s): [-12935.993235030832, 3306.5234815642566, 275.73217606979927]

What if we wanted a different reference frame? Well then we may either instantiate a different pykep.udpla.de440s with a different frame, and/or observer, or perform the transformation ourselves. If we want to perform the transformation ourselves, we need a rotation matrix, and we would then get it as follow:

R = pk.utils.rotation_matrix("ECLIPJ2000", "J2000")

Note that we have not specified an epoch, since these are both inertial frames and their relative orientation does not depend on time. We then perform the transformation.

r_j2000 = np.dot(R,r)
print(f"Position (m): {r_j2000}")

Position (m): [1.56005590e+11 6.84552122e+11 2.89625240e+11]

The same result would be obtained instantiating directly the pykep.planet in the final reference frame.

udpla = pk.udpla.de440s("JUPITER BARYCENTER", "J2000", "SSB")
jupiter_j2000 = pk.planet(udpla)
r, v = jupiter_j2000.eph(when)
print(f"Position (m): {r}")

Position (m): [156005590351.0843, 684552121902.1022, 289625240455.7204]