Cassini Position Example

Below is an interactive example that uses spiceypy to plot the position of the Cassini spacecraft relative to the barycenter of Saturn.

from pyscript import display
import numpy as np
import matplotlib
matplotlib.use("AGG")
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import spiceypy as spice

First import spiceypy and test it out.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import spiceypy as spice
# Print out the toolkit version
print(f'SpiceyPy for {spice.tkvrsn("TOOLKIT")} ready!')

'CSPICE_N0067'

We will need to load some kernels. If you are running this example outside of this browser page, you will need to download the following kernels from the NAIF servers via the links provided. After the kernels have been downloaded, you can define a meta kernel file or simply provide the paths directly to the spiceypy.furnsh method. For more on defining meta kernels in spice, please consult the Kernel Required Reading.

# the kernels with "_s" postfixes are subset from the NAIF hosted versions to make this page load faster
spice.furnsh([
    'kernels/lsk/naif0008.tls',
    'kernels/sclk/cas00084.tsc',
    'kernels/pck/cpck05Mar2004.tpc',
    'kernels/fk/cas_v37.tf',
    'kernels/ik/cas_iss_v09.ti',
    'kernels/spk/020514_SE_SAT105_s.bsp',
    'kernels/spk/981005_PLTEPH-DE405S_s.bsp',
    'kernels/spk/030201AP_SK_SM546_T45_s.bsp'
])
print(f'Loaded {spice.ktotal("ALL")} kernels')

# we are going to get positions between these two dates
etOne = spice.str2et("Jun 20, 2004")
etTwo = spice.str2et("Dec 1, 2005")
print(f"ET One: {etOne}, ET Two: {etTwo}")

ET One: 140961664.18440723, ET Two: 186667264.18308285

# get times for 4000 steps
steps = 4000
times = [x * (etTwo - etOne) / step + etOne for x in range(steps)]

# check first few times:
print(times[0:3])

[140961664.18440723, 140973090.5844069, 140984516.98440656]

# check the documentation on spkpos before continuing
help(spice.spkpos)

the above command will return something like the following:

Help on function spkpos in module spiceypy.spiceypy:

spkpos(
    targ: str,
    et: Union[float, numpy.ndarray],
    ref: str,
    abcorr: str,
    obs: str
) -> Union[Tuple[numpy.ndarray, float], Tuple[numpy.ndarray, numpy.ndarray]]
    Return the position of a target body relative to an observing
    body, optionally corrected for light time (planetary aberration)
    and stellar aberration.

    https://naif.jpl.nasa.gov/pub/naif/misc/toolkit_docs_N0067/C/cspice/spkpos_c.html

    :param targ: Target body name.
    :param et: Observer epoch in seconds past J2000 TDB.
    :param ref: Reference frame of output position vector.
    :param abcorr: Aberration correction flag.
    :param obs: Observing body name.
    :return:
            Position of target in km,
            One way light time between observer and target in seconds.

# Run spkpos as a vectorized function
positions, lightTimes = spice.spkpos(
    "Cassini", times, "J2000", "NONE", "SATURN BARYCENTER"
)

# Positions is a 3xN vector of XYZ positions
print("Positions: ")
print(positions[0])

# Light times is a N vector of time
print("Light Times: ")
print(lightTimes[0])

Positions:
[-5461446.61080924 -4434793.40785864 -1200385.93315424]
Light Times:
23.8062238783

# Clean up the kernel pool, not important for this example, but good practice
spice.kclear()

We will use matplotlib's 3D plotting to visualize Cassini's coordinates. We first convert the positions list to a 2D numpy array for easier indexing in the plot.

positions = (
    positions.T
)  # positions is shaped (4000, 3), let's transpose to (3, 4000) for easier indexing
fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(111, projection="3d")
ax.plot(positions[0], positions[1], positions[2])
plt.title("SpiceyPy Cassini Position Example from Jun 20, 2004 to Dec 1, 2005")
display(fig, target="mpl", append=False)