from planetarypy.spice.spicer import SpicerSurface Illumination with Spicer
The Spicer class provides surface illumination calculations for solar system bodies using SPICE kernels. Given a body, location, and time, it computes solar incidence angles, surface flux, solar longitude (season), and local solar time.
No per-body subclasses needed — Spicer("MARS"), Spicer("MOON"), Spicer("ENCELADUS") all work.
Requirements: The NAIF generic kernels (LSK for time conversion, PCK for body constants) are auto-loaded by planetarypy. The set of supported bodies is defined by the PCK kernel (pck00011.tpc) — this includes the Sun, all 8 planets, Pluto, 70+ moons, and a few asteroids (Ceres, Vesta, Lutetia). For emission/phase angles with a specific observer (e.g. MRO), mission-specific kernels are needed.
Body properties
A Spicer object wraps a NAIF body with its radii, reference frame, and SPICE ID.
mars = Spicer("MARS")
print(mars)
print(f"Radii: {mars.radii} km")
print(f"Reference frame: {mars.ref_frame}")
print(f"NAIF ID: {mars.target_id}")# What bodies can Spicer work with? Depends on the loaded PCK kernel.bodies = Spicer.supported_bodies()print(f"{len(bodies)} bodies available (from generic PCK kernel):\n")for code, name, radius in bodies[:5]: print(f" {code:>8d} {name:<20s} R={radius:.0f} km")print(" ...")Solar longitude and season
Solar longitude (L_s) tells you the season. On Mars: 0°=northern spring equinox, 90°=northern summer solstice, 180°=northern autumn equinox, 270°=northern winter solstice (southern summer).
print(f"Mars L_s right now: {mars.solar_longitude():.1f}°")
# L_s at specific times
for date in ["2024-01-01", "2024-04-01", "2024-07-01", "2024-10-01"]:
ls = mars.solar_longitude(date)
print(f" {date}: L_s = {ls:.1f}°")# Sub-solar point: where is the sun directly overhead?
ss_lon, ss_lat = mars.subsolar_point("2024-06-15T12:00:00")
print(f"Sub-solar point: lon={ss_lon:.1f}°, lat={ss_lat:.1f}°")
print(f"Solar constant at Mars distance: {mars.solar_constant('2024-06-15'):.0f} W/m²")Illumination at a surface point
The main function: given lon/lat and time, compute everything about the illumination at that point.
# Gale Crater at local noon
illum = mars.illumination(lon=137.4, lat=-4.6, time="2024-06-15T00:00:00")
print(illum)import matplotlib.pyplot as plt
import numpy as np
# Diurnal flux curve at Gale Crater
hours = np.arange(0, 24, 0.5)
fluxes = []
local_times = []
for h in hours:
hh = int(h)
mm = int((h % 1) * 60)
illum = mars.illumination(lon=137.4, lat=-4.6, time=f"2024-06-15T{hh:02d}:{mm:02d}:00")
fluxes.append(illum.solar_flux)
local_times.append(illum.local_solar_time)
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(hours, fluxes, "b-", lw=2)
ax.set_xlabel("UTC hour")
ax.set_ylabel("Solar flux [W/m²]")
ax.set_title(f"Gale Crater diurnal flux — L_s={mars.solar_longitude('2024-06-15'):.0f}°")
ax.set_xlim(0, 24)
ax.grid(True, alpha=0.3)
# Mark local noon
noon_idx = np.argmax(fluxes)
ax.annotate(f"local {local_times[noon_idx]}",
xy=(hours[noon_idx], fluxes[noon_idx]),
xytext=(hours[noon_idx]+2, fluxes[noon_idx]-50),
arrowprops=dict(arrowstyle="->"), fontsize=10)
plt.show()Slope and aspect effects
Real terrain isn’t flat. The slope and aspect parameters tilt the surface normal to simulate sloped terrain:
- slope: tilt angle in degrees (0 = flat, 90 = vertical)
- aspect: direction the slope faces, in degrees clockwise from north (0 = north-facing, 180 = south-facing)
The flux is computed using the angle between the sun vector and the tilted surface normal.
time = "2024-06-15T00:00:00" # Gale Crater near local noon
flat = mars.illumination(lon=137.4, lat=-4.6, time=time)
print(f"L_s = {flat.l_s:.0f}° (southern summer)")
print(f"Solar incidence: {flat.solar_incidence:.1f}°")
print()
for aspect, label in [(0, "north"), (90, "east"), (180, "south"), (270, "west")]:
illum = mars.illumination(lon=137.4, lat=-4.6, time=time, slope=30, aspect=aspect)
print(f"30° {label}-facing: {illum.solar_flux:.0f} W/m²")
print(f"\nFlat surface: {flat.solar_flux:.0f} W/m²")# How does flux vary with slope angle for different aspects?
slopes = np.arange(0, 75, 5)
fig, ax = plt.subplots(figsize=(10, 5))
for aspect, label, color in [(0, "north-facing", "blue"),
(90, "east-facing", "green"),
(180, "south-facing", "red"),
(270, "west-facing", "orange")]:
flux_vals = [mars.illumination(lon=137.4, lat=-4.6, time=time,
slope=s, aspect=aspect).solar_flux
for s in slopes]
ax.plot(slopes, flux_vals, label=label, color=color, lw=2)
ax.axhline(flat.solar_flux, color="gray", ls="--", label="flat")
ax.set_xlabel("Slope [deg]")
ax.set_ylabel("Solar flux [W/m²]")
ax.set_title(f"Gale Crater flux vs slope — L_s={flat.l_s:.0f}° (sun from south)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()Astropy units support
Set units=True to get astropy Quantities instead of plain floats. Useful for calculations where unit tracking prevents errors.
mars_u = Spicer("MARS", units=True)
illum = mars_u.illumination(lon=137.4, lat=-4.6, time="2024-06-15T00:00:00")
print(f"Solar incidence: {illum.solar_incidence}")
print(f"Solar flux: {illum.solar_flux}")
print(f"L_s: {illum.l_s}")
print(f"\nType: {type(illum.solar_flux)}")Integration with Point class
Spicer works directly with planetarypy.geo.Point objects:
from planetarypy.geo import Point
gale = Point(lon=137.4, lat=-4.6, crs="IAU_2015:49900")
illum = mars.illumination_at(gale, time="2024-06-15T00:00:00")
print(f"Gale Crater: {illum.solar_flux:.0f} W/m², local time {illum.local_solar_time}")
# Compare with Jezero Crater
jezero = Point(lon=77.45, lat=18.44, crs="IAU_2015:49900")
illum_j = mars.illumination_at(jezero, time="2024-06-15T00:00:00")
print(f"Jezero: {illum_j.solar_flux:.0f} W/m², local time {illum_j.local_solar_time}")Sun direction for azimuth calculation
To compute the solar azimuth at a point (e.g. for shadow direction), Spicer can return a nearby surface point in the direction of the sun. Combined with Point.azimuth_to(), this gives the solar azimuth without relying on index metadata.
time = "2024-06-15T00:00:00"sun_lon, sun_lat = mars.sun_direction_at(137.4, -4.6, time=time)print(f"Point toward sun: lon={sun_lon:.4f}°, lat={sun_lat:.4f}°")print("Origin point: lon=137.4000°, lat=-4.6000°")print("\nThe sun direction point is within a pixel of the origin,")print("shifted toward the subsolar point.")Validation: SPICE vs HiRISE indexWe can validate our SPICE-based solar azimuth against the SUB_SOLAR_AZIMUTH stored in the HiRISE RDR index — they should agree closely.
from planetarypy.pds import get_index# Look up ESP_013807_2035_RED in the HiRISE RDR indexrdr = get_index("mro.hirise.rdr", allow_refresh=False)obs = rdr[rdr["PRODUCT_ID"] == "ESP_013807_2035_RED"].squeeze()# Actual image center from corner coordinatescenter_lat = np.mean([obs[f"CORNER{i}_LATITUDE"] for i in range(1, 5)])center_lon = np.mean([obs[f"CORNER{i}_LONGITUDE"] for i in range(1, 5)])obs_time = str(obs["START_TIME"])# SPICE: geographic solar azimuth (CW from north)spice_az = mars.solar_azimuth_at(center_lon, center_lat, time=obs_time)# HiRISE index: CW from 3 o'clock → convert to CW from north# For equirectangular RDR with north=top: geographic = (hirise + 90) % 360hirise_az = obs["SUB_SOLAR_AZIMUTH"]hirise_geographic = (hirise_az + 90) % 360print(f"ESP_013807_2035 at ({center_lon:.1f}°, {center_lat:.1f}°)")print(f" SPICE solar azimuth: {spice_az:.1f}° (CW from north)")print(f" HiRISE index: {hirise_az:.1f}° (CW from right) → {hirise_geographic:.1f}° (CW from north)")print(f" Difference: {abs(spice_az - hirise_geographic):.1f}°")print("\n Agreement within 1 degree — SPICE and index are consistent.")Multi-body comparison
Since Spicer works for any body, we can compare illumination across the solar system:
time = "2024-06-15T12:00:00"
for body, lon, lat, label in [
("MARS", 0, 0, "Mars equator"),
("MOON", 0, 0, "Moon equator"),
("MERCURY", 0, 0, "Mercury equator"),
]:
s = Spicer(body)
sc = s.solar_constant(time)
illum = s.illumination(lon=lon, lat=lat, time=time)
print(f"{label:20s} solar constant={sc:.0f} W/m² flux={illum.solar_flux:.0f} W/m² incidence={illum.solar_incidence:.1f}°")Outer solar system: on-demand ephemerisThe generic SPICE kernels only include ephemeris for the inner solar system (Mercury through Mars + moons). For Jupiter, Saturn, Neptune, and Pluto systems, Spicer automatically downloads the satellite ephemeris SPK on first use. These are large files (100 MB to 1+ GB) but only need to be downloaded once.
# Triton — automatically downloads Neptune satellite SPK (~100 MB) on first calltriton = Spicer("TRITON")illum = triton.illumination(lon=0, lat=0, time="2024-06-15T12:00:00")print(f"Triton: solar constant = {triton.solar_constant('2024-06-15'):.1f} W/m²")print(f" (That far from the Sun, not much light!)")# Also try: Spicer("ENCELADUS"), Spicer("EUROPA"), Spicer("CHARON")Command-Line InterfaceAll of this is also available from the terminal via plp spicer:bash# Current Mars status$ plp spicer Mars MARS ==== Radii: 3396.2 x 3396.2 x 3376.2 km Reference frame: IAU_MARS Solar longitude (L_s): 262.7 Sub-solar point: lon=230.02, lat=-24.98 Solar constant: 711 W/m²# With surface point$ plp spicer Mars --lon 137.4 --lat -4.6 Surface at (137.4, -4.6): Solar incidence: 90.4° Solar flux: 0 W/m² Local solar time: 05:49:40 Solar azimuth: 115.1° (CW from N)# Any body works$ plp spicer Moon$ plp spicer Enceladus$ plp spicer Europa --time 2024-06-15T12:00:00And the quick property mars.Ls gives you the current solar longitude without any arguments:
# Quick check: what's the current Mars season?print(f"Mars L_s right now: {mars.Ls:.1f}°")