Jorge Martinez

Aerospace Engineer and Senior Software Developer


Learning SPICE at ESAC

Jorge Martínez Garrido

April 29, 2023

software astrodynamics


The last SPICE training session was held at ESAC, in Madrid, and I had the pleasure to attend it.

Attendes photo

What is SPICE?

SPICE (Spacecraft Planet Instrument C-matrix Events) is a software toolkit developed by NASA’s Navigation and Ancillary Information Facility (NAIF) for space mission design, navigation, and science data analysis. It provides a set of software libraries and tools for manipulating and visualizing space mission geometry, spacecraft instrument pointing, and other related information.

SPICE consists of a set of subroutines and data files that provide information on the geometric and time-varying relationships between various objects in space, such as planets, satellites, asteroids, and comets. This information includes the positions and orientations of these objects, as well as the pointing directions and field-of-view geometries of spacecraft instruments.

SPICE is used extensively in NASA’s planetary exploration missions, including missions to Mars, Jupiter, Saturn, and other celestial bodies. It is also used in other space-related applications, such as satellite mission design, orbit determination, and space debris tracking.

SPICE logo

Why learning SPICE?

Learning SPICE can be beneficial if you are interested in the field of space exploration, spacecraft design, navigation, and science data analysis. Here are a few reasons why learning SPICE can be useful:

Overall, learning SPICE can enhance your understanding of space exploration, spacecraft design, navigation, and science data analysis, and it can open up exciting career opportunities in the space industry.

A quick example: plotting the evolution of the barycenter of the Solar System

As a personal challenge, I proposed myself to generate a figure showing the evolution of the barycenter of the Solar System (SSB). The reason is that this figure usually appears in literature to demonstrate that the Sun is not the center of our system. However, this figures usually lack of critial information such us the frame of reference and the ephemerides used to create the figure.

Collecting the required kernels

Let us start by downloading all the required kernels for this lesson. Required kernels include the leap seconds kernel, the de430s ephemerides model, and a planetary constants kernel. For the purpose of this tutorial, all kernels are saved in a folder named kernels/.

! wget -O kernels/latest_leapseconds.tls https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/latest_leapseconds.tls
! wget -O kernels/de432s.bsp https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de432s.bsp
! wget -O kernels/pck00011.tpc https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00011.tpc
--2024-03-31 19:54:58--  https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/latest_leapseconds.tls
Resolving naif.jpl.nasa.gov (naif.jpl.nasa.gov)... 137.78.232.95
Connecting to naif.jpl.nasa.gov (naif.jpl.nasa.gov)|137.78.232.95|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 5257 (5,1K) [text/plain]
Saving to: ‘kernels/latest_leapseconds.tls’

          kernels/l   0%[                    ]       0  --.-KB/s               kernels/latest_leap 100%[===================>]   5,13K  --.-KB/s    in 0s      

2024-03-31 19:54:59 (23,7 MB/s) - ‘kernels/latest_leapseconds.tls’ saved [5257/5257]

--2024-03-31 19:54:59--  https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de432s.bsp
Resolving naif.jpl.nasa.gov (naif.jpl.nasa.gov)... 137.78.232.95
Connecting to naif.jpl.nasa.gov (naif.jpl.nasa.gov)|137.78.232.95|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 10895360 (10M)
Saving to: ‘kernels/de432s.bsp’

kernels/de432s.bsp    0%[                    ]       0  --.-KB/s               kernels/de432s.bsp    0%[                    ]  40,00K   145KB/s               kernels/de432s.bsp    2%[                    ] 216,00K   396KB/s               kernels/de432s.bsp    8%[>                   ] 888,00K  1,07MB/s               kernels/de432s.bsp   26%[====>               ]   2,80M  2,77MB/s               kernels/de432s.bsp   45%[========>           ]   4,73M  3,89MB/s               kernels/de432s.bsp   65%[============>       ]   6,82M  4,81MB/s               kernels/de432s.bsp   88%[================>   ]   9,16M  5,64MB/s               kernels/de432s.bsp  100%[===================>]  10,39M  5,99MB/s    in 1,7s    

2024-03-31 19:55:01 (5,99 MB/s) - ‘kernels/de432s.bsp’ saved [10895360/10895360]

--2024-03-31 19:55:02--  https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00011.tpc
Resolving naif.jpl.nasa.gov (naif.jpl.nasa.gov)... 137.78.232.95
Connecting to naif.jpl.nasa.gov (naif.jpl.nasa.gov)|137.78.232.95|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 131226 (128K) [text/plain]
Saving to: ‘kernels/pck00011.tpc’

kernels/pck00011.tp   0%[                    ]       0  --.-KB/s               kernels/pck00011.tp  31%[=====>              ]  40,00K   151KB/s               kernels/pck00011.tp 100%[===================>] 128,15K   321KB/s    in 0,4s    

2024-03-31 19:55:03 (321 KB/s) - ‘kernels/pck00011.tpc’ saved [131226/131226]

Finding the position of the SSB in the ecliptic frame (ECLIPJ2000) as seen from the Sun

The position of a target with respect to an obeserver can be found by using the spkpos function. The code below generates a list of epochs ranging from the desired starting and ending dates. Finally, the position of the SSB is computed in the Ecliptic J2000 frame as seen from the Sun without considering any light abberation correction.

import numpy as np
import spiceypy as spice

# Load all the desired kernels
for file in ["latest_leapseconds.tls", "de432s.bsp", "pck00011.tpc"]:
    spice.furnsh(f"kernels/{file}")
    
# Declare the desired year epochs
start_year, end_year, delta_year, n_points = 1990, 2025, 5, 1_000
year_epochs = [spice.str2et(f"{year}-01-01") for year in range(start_year, end_year+1, delta_year)]
epochs = np.linspace(spice.str2et(f"{start_year}-01-01"), spice.str2et(f"{end_year}-01-01"), n_points)

# Compute the position of the SSB at desired epoch in the Ecliptic J2000 w.r.t. the Sun
frame = "ECLIPJ2000"
sun_pos_at_epochs = [spice.spkpos("SSB", epoch, frame, "NONE", "SUN")[0] for epoch in epochs]
sun_pos_at_years = [spice.spkpos("SSB", epoch, frame, "NONE", "SUN")[0] for epoch in year_epochs]

Visualizing the evolution of the SSB

The following code generates a figure of the SSB as seen from the Sun in the Ecliptic J2000 frame.

import matplotlib
import matplotlib.pyplot as plt

def draw_circle(ax, position_xy, radius, color="white", border_color="black"):
    circle = matplotlib.patches.Circle(
        position_xy,
        radius,
        facecolor=color,
        edgecolor=border_color,
        alpha=0.7
    )
    ax.add_patch(circle)

def draw_sun(ax, position_xy):
    _, (SUN_RADIUS, *_) = spice.bodvrd("SUN", "RADII", 3)
    draw_circle(ax, (0, 0), SUN_RADIUS, color="orange")
    ax.text(-0.12 * SUN_RADIUS, -0.02 * SUN_RADIUS, "Nucleus")
    draw_circle(ax, (0, 0), 0.25 * SUN_RADIUS, color="yellow")
    

# Generate the plot
fig, ax = plt.subplots(figsize=(10,10))
draw_sun(ax, (0, 0))

# Draw the position evolution of the solar system barycenter
cmap = matplotlib.colormaps["gist_heat_r"]
sun_pos_x_at_epochs = [pos[0] for pos in sun_pos_at_epochs]
sun_pos_y_at_epochs = [pos[1] for pos in sun_pos_at_epochs]
for i in range(n_points - 1):
    ax.plot(sun_pos_x_at_epochs[i:i+2], sun_pos_y_at_epochs[i:i+2], color=cmap(i/n_points), linestyle="-")

# Draw the position at each year
sun_pos_x_at_years = [pos[0] for pos in sun_pos_at_years]
sun_pos_y_at_years = [pos[1] for pos in sun_pos_at_years]
for ith_epoch, epoch in enumerate(year_epochs):
    ax.scatter(sun_pos_x_at_years[ith_epoch], sun_pos_y_at_years[ith_epoch], c="k", marker="o")
    ax.text(
        sun_pos_x_at_years[ith_epoch] + 2e4, 
        sun_pos_y_at_years[ith_epoch], 
        spice.et2datetime(epoch).year,
    )

# Add titles and labels
ax.set_title("Evolution of the SSB in the ECLIPJ2000 frame as seen from the SUN")
ax.set_xlabel("Distance [km]")
ax.set_ylabel("Distance [km]")

# Mark the years
ax.set_aspect("equal")