Planetary Ephemerides¶

This tutorial demonstrates satkit's planetary ephemeris capabilities, covering both high-precision JPL ephemerides and faster low-precision analytical models.

Overview¶

Satkit provides two approaches for computing solar system body positions:

• sk.jplephem — High-precision positions from NASA/JPL Development Ephemerides (DE440/DE441). These are interpolated from pre-computed Chebyshev polynomial coefficients and deliver sub-kilometer accuracy across centuries.
• sk.sun / sk.moon — Low-precision analytical models from Vallado, useful when speed matters more than accuracy. The Sun model is accurate to ~0.01 degrees (1950--2050), and the Moon model to ~0.3 degrees in ecliptic longitude and ~1275 km in range.

All positions are returned in meters in the GCRF (Geocentric Celestial Reference Frame) unless otherwise noted.

In [1]:

import satkit as sk
import numpy as np
import matplotlib.pyplot as plt
import scienceplots  # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']

import satkit as sk import numpy as np import matplotlib.pyplot as plt import scienceplots # noqa: F401 plt.style.use(["science", "no-latex", "../satkit.mplstyle"]) %config InlineBackend.figure_formats = ['svg']

Getting Planet Positions¶

The sk.jplephem module provides geocentric and barycentric positions for all major solar system bodies. Bodies are specified using the sk.solarsystem enum.

Geocentric Positions¶

sk.jplephem.geocentric_pos returns the position of a body relative to the Earth center, in the GCRF frame.

In [2]:

t = sk.time(2026, 3, 20)

# Sun position relative to Earth center
sun_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t)
print(f"Sun distance from Earth: {np.linalg.norm(sun_pos)/1e9:.3f} million km")
print(f"Sun position (GCRF): {sun_pos / 1e3} km")

# Moon position relative to Earth center
moon_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t)
print(f"\nMoon distance from Earth: {np.linalg.norm(moon_pos)/1e3:.0f} km")
print(f"Moon position (GCRF): {moon_pos / 1e3} km")

t = sk.time(2026, 3, 20) # Sun position relative to Earth center sun_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t) print(f"Sun distance from Earth: {np.linalg.norm(sun_pos)/1e9:.3f} million km") print(f"Sun position (GCRF): {sun_pos / 1e3} km") # Moon position relative to Earth center moon_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t) print(f"\nMoon distance from Earth: {np.linalg.norm(moon_pos)/1e3:.0f} km") print(f"Moon position (GCRF): {moon_pos / 1e3} km")

Sun distance from Earth: 148.962 million km
Sun position (GCRF): [ 1.48940181e+08 -2.32276697e+06 -1.00753548e+06] km

Moon distance from Earth: 370245 km
Moon position (GCRF): [362541.91468222  59535.30489146  45829.30373204] km

In [3]:

# Positions of other planets
planets = [
    ("Mercury", sk.solarsystem.Mercury),
    ("Venus",   sk.solarsystem.Venus),
    ("Mars",    sk.solarsystem.Mars),
    ("Jupiter", sk.solarsystem.Jupiter),
    ("Saturn",  sk.solarsystem.Saturn),
]

for name, body in planets:
    pos = sk.jplephem.geocentric_pos(body, t)
    dist_au = np.linalg.norm(pos) / 1.496e11  # convert meters to AU
    print(f"{name:>8s}: {dist_au:.3f} AU from Earth")

# Positions of other planets planets = [ ("Mercury", sk.solarsystem.Mercury), ("Venus", sk.solarsystem.Venus), ("Mars", sk.solarsystem.Mars), ("Jupiter", sk.solarsystem.Jupiter), ("Saturn", sk.solarsystem.Saturn), ] for name, body in planets: pos = sk.jplephem.geocentric_pos(body, t) dist_au = np.linalg.norm(pos) / 1.496e11 # convert meters to AU print(f"{name:>8s}: {dist_au:.3f} AU from Earth")

 Mercury: 0.680 AU from Earth
   Venus: 1.609 AU from Earth
    Mars: 2.313 AU from Earth
 Jupiter: 4.880 AU from Earth
  Saturn: 10.486 AU from Earth

Barycentric Positions¶

sk.jplephem.barycentric_pos returns positions relative to the solar system barycenter. This is useful for interplanetary trajectory work or for understanding the geometry of the solar system.

Note that the Sun's barycentric position is close to the origin but not exactly zero, since the barycenter shifts due to the gravitational influence of the giant planets.

In [4]:

# Barycentric positions
sun_bary = sk.jplephem.barycentric_pos(sk.solarsystem.Sun, t)
earth_bary = sk.jplephem.barycentric_pos(sk.solarsystem.EMB, t)
jupiter_bary = sk.jplephem.barycentric_pos(sk.solarsystem.Jupiter, t)

print(f"Sun offset from barycenter: {np.linalg.norm(sun_bary)/1e3:.0f} km")
print(f"Earth-Moon barycenter:      {np.linalg.norm(earth_bary)/1e11:.4f} x 10^8 km")
print(f"Jupiter:                    {np.linalg.norm(jupiter_bary)/1e11:.4f} x 10^8 km")

# Barycentric positions sun_bary = sk.jplephem.barycentric_pos(sk.solarsystem.Sun, t) earth_bary = sk.jplephem.barycentric_pos(sk.solarsystem.EMB, t) jupiter_bary = sk.jplephem.barycentric_pos(sk.solarsystem.Jupiter, t) print(f"Sun offset from barycenter: {np.linalg.norm(sun_bary)/1e3:.0f} km") print(f"Earth-Moon barycenter: {np.linalg.norm(earth_bary)/1e11:.4f} x 10^8 km") print(f"Jupiter: {np.linalg.norm(jupiter_bary)/1e11:.4f} x 10^8 km")

Sun offset from barycenter: 903585 km
Earth-Moon barycenter:      1.4932 x 10^8 km
Jupiter:                    7.8332 x 10^8 km

Vectorized Queries¶

All sk.jplephem functions accept arrays of times, returning an Nx3 array of positions. This is much faster than calling the function in a Python loop.

In [5]:

# Create an array of times spanning one year
t0 = sk.time(2026, 1, 1)
times = np.array([t0 + sk.duration(days=d) for d in range(365)])

# Get Mars geocentric position over the year (returns 365x3 array)
mars_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Mars, times)
mars_dist = np.linalg.norm(mars_pos, axis=1) / 1.496e11  # AU

print(f"Mars position array shape: {mars_pos.shape}")
print(f"Mars distance range: {mars_dist.min():.2f} to {mars_dist.max():.2f} AU")

# Create an array of times spanning one year t0 = sk.time(2026, 1, 1) times = np.array([t0 + sk.duration(days=d) for d in range(365)]) # Get Mars geocentric position over the year (returns 365x3 array) mars_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Mars, times) mars_dist = np.linalg.norm(mars_pos, axis=1) / 1.496e11 # AU print(f"Mars position array shape: {mars_pos.shape}") print(f"Mars distance range: {mars_dist.min():.2f} to {mars_dist.max():.2f} AU")

Mars position array shape: (365, 3)
Mars distance range: 0.92 to 2.41 AU

State Vectors¶

sk.jplephem.geocentric_state returns both position and velocity, which is essential for orbit determination and trajectory planning.

In [6]:

t = sk.time(2026, 3, 20)

# Moon state vector
pos, vel = sk.jplephem.geocentric_state(sk.solarsystem.Moon, t)
print("Moon state vector (GCRF):")
print(f"  Position: [{pos[0]/1e3:.1f}, {pos[1]/1e3:.1f}, {pos[2]/1e3:.1f}] km")
print(f"  Velocity: [{vel[0]:.3f}, {vel[1]:.3f}, {vel[2]:.3f}] m/s")
print(f"  Speed:    {np.linalg.norm(vel):.1f} m/s")

# Sun state vector
pos, vel = sk.jplephem.geocentric_state(sk.solarsystem.Sun, t)
print("\nSun state vector (GCRF):")
print(f"  Position: [{pos[0]/1e9:.4f}, {pos[1]/1e9:.4f}, {pos[2]/1e9:.4f}] x 10^6 km")
print(f"  Velocity: [{vel[0]/1e3:.3f}, {vel[1]/1e3:.3f}, {vel[2]/1e3:.3f}] km/s")
print(f"  Speed:    {np.linalg.norm(vel)/1e3:.2f} km/s")

t = sk.time(2026, 3, 20) # Moon state vector pos, vel = sk.jplephem.geocentric_state(sk.solarsystem.Moon, t) print("Moon state vector (GCRF):") print(f" Position: [{pos[0]/1e3:.1f}, {pos[1]/1e3:.1f}, {pos[2]/1e3:.1f}] km") print(f" Velocity: [{vel[0]:.3f}, {vel[1]:.3f}, {vel[2]:.3f}] m/s") print(f" Speed: {np.linalg.norm(vel):.1f} m/s") # Sun state vector pos, vel = sk.jplephem.geocentric_state(sk.solarsystem.Sun, t) print("\nSun state vector (GCRF):") print(f" Position: [{pos[0]/1e9:.4f}, {pos[1]/1e9:.4f}, {pos[2]/1e9:.4f}] x 10^6 km") print(f" Velocity: [{vel[0]/1e3:.3f}, {vel[1]/1e3:.3f}, {vel[2]/1e3:.3f}] km/s") print(f" Speed: {np.linalg.norm(vel)/1e3:.2f} km/s")

Moon state vector (GCRF):
  Position: [362541.9, 59535.3, 45829.3] km
  Velocity: [-244.528, 916.786, 484.697] m/s
  Speed:    1065.5 m/s

Sun state vector (GCRF):
  Position: [148.9402, -2.3228, -1.0075] x 10^6 km
  Velocity: [0.987, 27.441, 11.896] km/s
  Speed:    29.92 km/s

Low-Precision Models¶

For applications where speed matters more than accuracy (e.g., quick visibility checks, rough force models), satkit provides analytical low-precision models:

• sk.sun.pos_gcrf — Sun position accurate to ~0.01 degrees (Vallado Algorithm 29)
• sk.moon.pos_gcrf — Moon position accurate to ~0.3 degrees in ecliptic longitude and ~1275 km in range (Vallado Algorithm 31)

These do not require JPL ephemeris data files and are faster to compute.

In [7]:

t = sk.time(2026, 3, 20)

# Low-precision Sun position
sun_lp = sk.sun.pos_gcrf(t)
print(f"Sun distance (low-precision): {np.linalg.norm(sun_lp)/1e9:.3f} million km")

# Low-precision Moon position
moon_lp = sk.moon.pos_gcrf(t)
print(f"Moon distance (low-precision): {np.linalg.norm(moon_lp)/1e3:.0f} km")

# Moon phase information
illum = sk.moon.illumination(t)
phase = sk.moon.phase_name(t)
print(f"\nMoon illumination: {illum*100:.1f}%")
print(f"Moon phase: {phase}")

t = sk.time(2026, 3, 20) # Low-precision Sun position sun_lp = sk.sun.pos_gcrf(t) print(f"Sun distance (low-precision): {np.linalg.norm(sun_lp)/1e9:.3f} million km") # Low-precision Moon position moon_lp = sk.moon.pos_gcrf(t) print(f"Moon distance (low-precision): {np.linalg.norm(moon_lp)/1e3:.0f} km") # Moon phase information illum = sk.moon.illumination(t) phase = sk.moon.phase_name(t) print(f"\nMoon illumination: {illum*100:.1f}%") print(f"Moon phase: {phase}")

Sun distance (low-precision): 148.954 million km
Moon distance (low-precision): 370539 km

Moon illumination: 1.2%
Moon phase: moonphase.NewMoon

Comparing JPL vs Low-Precision Models¶

How much accuracy do you sacrifice by using the analytical models? Let's compare Sun and Moon positions from both methods over an entire year and plot the position error.

In [8]:

# Generate times over one year at daily intervals
t0 = sk.time(2026, 1, 1)
times = np.array([t0 + sk.duration(days=d) for d in range(365)])
days = np.arange(365)

# JPL high-precision positions
sun_jpl = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, times)
moon_jpl = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times)

# Low-precision positions
sun_lp = sk.sun.pos_gcrf(times)
moon_lp = sk.moon.pos_gcrf(times)

# Compute position errors
sun_err = np.linalg.norm(sun_jpl - sun_lp, axis=1) / 1e3   # km
moon_err = np.linalg.norm(moon_jpl - moon_lp, axis=1) / 1e3  # km

print(f"Sun position error  — mean: {sun_err.mean():.0f} km, max: {sun_err.max():.0f} km")
print(f"Moon position error — mean: {moon_err.mean():.0f} km, max: {moon_err.max():.0f} km")

# Generate times over one year at daily intervals t0 = sk.time(2026, 1, 1) times = np.array([t0 + sk.duration(days=d) for d in range(365)]) days = np.arange(365) # JPL high-precision positions sun_jpl = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, times) moon_jpl = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times) # Low-precision positions sun_lp = sk.sun.pos_gcrf(times) moon_lp = sk.moon.pos_gcrf(times) # Compute position errors sun_err = np.linalg.norm(sun_jpl - sun_lp, axis=1) / 1e3 # km moon_err = np.linalg.norm(moon_jpl - moon_lp, axis=1) / 1e3 # km print(f"Sun position error — mean: {sun_err.mean():.0f} km, max: {sun_err.max():.0f} km") print(f"Moon position error — mean: {moon_err.mean():.0f} km, max: {moon_err.max():.0f} km")

Sun position error  — mean: 13114 km, max: 28262 km
Moon position error — mean: 2518 km, max: 3859 km

In [9]:

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)

ax1.plot(days, sun_err)
ax1.set_ylabel("Position Error (km)")
ax1.set_title("Sun: Low-Precision vs JPL DE440")
ax1.grid(True, alpha=0.3)

ax2.plot(days, moon_err)
ax2.set_xlabel("Day of Year 2026")
ax2.set_ylabel("Position Error (km)")
ax2.set_title("Moon: Low-Precision vs JPL DE440")
ax2.grid(True, alpha=0.3)

fig.tight_layout()
plt.show()

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True) ax1.plot(days, sun_err) ax1.set_ylabel("Position Error (km)") ax1.set_title("Sun: Low-Precision vs JPL DE440") ax1.grid(True, alpha=0.3) ax2.plot(days, moon_err) ax2.set_xlabel("Day of Year 2026") ax2.set_ylabel("Position Error (km)") ax2.set_title("Moon: Low-Precision vs JPL DE440") ax2.grid(True, alpha=0.3) fig.tight_layout() plt.show()

Practical Example: Moon Distance Over a Month¶

The Moon's orbit is noticeably elliptical, with perigee (closest approach) around 363,000 km and apogee (farthest point) around 405,000 km. Let's compute and visualize the Moon's distance over a month to identify perigee and apogee passages.

In [10]:

# Compute Moon distance at hourly intervals over one month
t0 = sk.time(2026, 3, 1)
hours = np.arange(0, 30 * 24)  # 30 days in hours
times = np.array([t0 + sk.duration(hours=float(h)) for h in hours])

moon_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times)
moon_dist = np.linalg.norm(moon_pos, axis=1) / 1e3  # km

# Find perigee and apogee
perigee_idx = np.argmin(moon_dist)
apogee_idx = np.argmax(moon_dist)

print(f"Perigee: {moon_dist[perigee_idx]:.0f} km on {times[perigee_idx].datetime()}")
print(f"Apogee:  {moon_dist[apogee_idx]:.0f} km on {times[apogee_idx].datetime()}")

# Compute Moon distance at hourly intervals over one month t0 = sk.time(2026, 3, 1) hours = np.arange(0, 30 * 24) # 30 days in hours times = np.array([t0 + sk.duration(hours=float(h)) for h in hours]) moon_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times) moon_dist = np.linalg.norm(moon_pos, axis=1) / 1e3 # km # Find perigee and apogee perigee_idx = np.argmin(moon_dist) apogee_idx = np.argmax(moon_dist) print(f"Perigee: {moon_dist[perigee_idx]:.0f} km on {times[perigee_idx].datetime()}") print(f"Apogee: {moon_dist[apogee_idx]:.0f} km on {times[apogee_idx].datetime()}")

Perigee: 366857 km on 2026-03-22 12:00:00+00:00
Apogee:  404384 km on 2026-03-10 14:00:00+00:00

In [11]:

days_elapsed = hours / 24.0

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(days_elapsed, moon_dist / 1e3, linewidth=1.5)
ax.plot(days_elapsed[perigee_idx], moon_dist[perigee_idx] / 1e3, 'v',
        markersize=10, color='green', label=f'Perigee: {moon_dist[perigee_idx]:.0f} km')
ax.plot(days_elapsed[apogee_idx], moon_dist[apogee_idx] / 1e3, '^',
        markersize=10, color='red', label=f'Apogee: {moon_dist[apogee_idx]:.0f} km')

ax.set_xlabel("Days since March 1, 2026")
ax.set_ylabel("Distance (x 1000 km)")
ax.set_title("Earth-Moon Distance — March 2026")
ax.legend()
ax.grid(True, alpha=0.3)

fig.tight_layout()
plt.show()

days_elapsed = hours / 24.0 fig, ax = plt.subplots(figsize=(10, 5)) ax.plot(days_elapsed, moon_dist / 1e3, linewidth=1.5) ax.plot(days_elapsed[perigee_idx], moon_dist[perigee_idx] / 1e3, 'v', markersize=10, color='green', label=f'Perigee: {moon_dist[perigee_idx]:.0f} km') ax.plot(days_elapsed[apogee_idx], moon_dist[apogee_idx] / 1e3, '^', markersize=10, color='red', label=f'Apogee: {moon_dist[apogee_idx]:.0f} km') ax.set_xlabel("Days since March 1, 2026") ax.set_ylabel("Distance (x 1000 km)") ax.set_title("Earth-Moon Distance — March 2026") ax.legend() ax.grid(True, alpha=0.3) fig.tight_layout() plt.show()