#
# ISC License
#
# Copyright (c) 2026, Autonomous Vehicle Systems Lab, University of Colorado at Boulder
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
#
#
# Basilisk Scenario Script and Integrated Test
#
# Purpose: Verify orbit-propagation consistency for a spherical-harmonics gravity
# field and demonstrate the effect of the planet rotation on the
# tesseral/sectoral terms (issue #1352).
# Author: robotrocketscience (https://github.com/robotrocketscience)
# Creation Date: Jun. 25, 2026
#
r"""
Overview
--------
This scenario verifies that Basilisk's orbit propagation with a spherical-harmonic
gravity field is consistent with independent propagators, and demonstrates *why*
the planet rotation must be supplied when the gravity field contains
tesseral/sectoral terms (spherical-harmonic order :math:`m \ge 1`).
It propagates a near-circular, sun-synchronous LEO orbit (TerraSAR-X-like) under
the GGM03S Earth field, twice, with identical initial conditions, integrator, and
step size. The only difference is whether a planet-orientation message is
connected:
#. **Earth rotating** -- a SPICE interface is attached via the gravity body
factory's ``createSpiceInterface`` method, so the field is evaluated in the
rotating planet-fixed frame.
#. **Earth non-rotating** -- no orientation message is connected, so the gravity
effector falls back to an identity planet orientation and evaluates the field
in the inertial frame.
Background
----------
The spherical-harmonic field is evaluated in the planet-fixed frame. Zonal terms
(order :math:`m = 0`, e.g. :math:`J_2`) depend only on latitude and are unaffected
by the planet rotation. Tesseral and sectoral terms (order :math:`m \ge 1`, e.g.
:math:`C_{22}`) depend on the body longitude. For a rotating planet these terms
average out over each day and produce only small short-period effects. If the
planet is (incorrectly) held fixed in inertial space, those terms no longer
average out and instead force large spurious secular/long-period drift in the
eccentricity vector and the inclination, while the semi-major axis and the
:math:`J_2`-driven RAAN precession remain essentially unchanged.
This is the inconsistency reported in issue #1352, where the eccentricity vector
and inclination of a Basilisk LEO propagation diverged from Orekit, GMAT, and
SpOCK results (which mutually agree). Those tools all model Earth rotation; the
Basilisk script that produced the discrepancy did not connect a planet-orientation
message. With Earth rotation connected, as in the *rotating* case here, the
eccentricity and inclination remain bounded, recovering consistency with the
external propagators.
When a spherical-harmonic body that has order :math:`\ge 1` terms is used without
an orientation message, :ref:`gravityEffector` now logs a ``BSK_WARNING`` at reset.
Running the script
------------------
The script is found in the folder ``basilisk/examples`` and executed by using::
python3 scenarioOrbitConsistencyVerification.py
Illustration of Simulation Results
----------------------------------
The eccentricity and inclination of the non-rotating case drift away secularly,
whereas the rotating case stays bounded -- the bounded behavior is what matches
the externally validated propagators.
The figures below are generated by the integrated test with a short, low-degree,
coarse-step run (``simDays = 12.0, gravDeg = 8, stepSeconds = 30.0``); a longer
higher-degree run (the script default) makes the drift even more pronounced.
.. image:: /_images/Scenarios/scenarioOrbitConsistencyVerification_ecc.svg
:align: center
.. image:: /_images/Scenarios/scenarioOrbitConsistencyVerification_inc.svg
:align: center
"""
import numpy as np
from Basilisk.architecture import bskLogging
from Basilisk.simulation import spacecraft
from Basilisk.utilities import SimulationBaseClass, macros, orbitalMotion, simIncludeGravBody
from Basilisk.utilities.supportDataTools.dataFetcher import get_path, DataFile
try:
from matplotlib import pyplot as plt
except ImportError:
plt = None
# TerraSAR-X-like initial conditions (J2000/ECI, near-polar sun-synchronous LEO).
R_N = [-138514.048884, -888553.853077, 6818674.278964] # [m]
V_N = [-7306.6601, 2126.208457, 128.787296] # [m/s]
EPOCH = "2020 DEC 04 00:58:00"
[docs]
def propagate(useEarthRotation, simDays, gravDeg, sampleMinutes, stepSeconds=10.0):
"""Propagate the reference orbit and return arrays of (days, a, e, inc[deg]).
Args:
useEarthRotation (bool): connect a SPICE planet-orientation message.
simDays (float): propagation duration in days.
gravDeg (int): spherical-harmonic degree/order of the GGM03S field.
sampleMinutes (float): output sampling period in minutes.
stepSeconds (float): integration time step in seconds.
"""
scSim = SimulationBaseClass.SimBaseClass()
proc = scSim.CreateNewProcess("p")
proc.addTask(scSim.CreateNewTask("t", macros.sec2nano(stepSeconds)))
scObject = spacecraft.Spacecraft()
scObject.ModelTag = "tsx"
scSim.AddModelToTask("t", scObject)
gravFactory = simIncludeGravBody.gravBodyFactory()
earth = gravFactory.createEarth()
earth.isCentralBody = True
ggm03s = get_path(DataFile.LocalGravData.GGM03S)
earth.useSphericalHarmonicsGravityModel(str(ggm03s), gravDeg)
mu = earth.mu
if useEarthRotation:
gravFactory.createSpiceInterface(time=EPOCH)
gravFactory.spiceObject.zeroBase = "Earth"
scSim.AddModelToTask("t", gravFactory.spiceObject, 100)
else:
# The non-rotating case deliberately omits the planet-orientation message
# to demonstrate the resulting drift, so gravityEffector's reset-time
# "missing planet rotation" BSK_WARNING (issue #1352) is expected here.
# Raise this module's log level to BSK_ERROR so that expected warning is
# not printed when running the scenario.
quietLogger = bskLogging.BSKLogger()
quietLogger.setLogLevel(bskLogging.BSK_ERROR)
scObject.gravField.bskLogger = quietLogger
scObject.gravField.gravBodies = spacecraft.GravBodyVector(
list(gravFactory.gravBodies.values())
)
scObject.hub.r_CN_NInit = R_N
scObject.hub.v_CN_NInit = V_N
rec = scObject.scStateOutMsg.recorder(macros.sec2nano(sampleMinutes * 60.0))
scSim.AddModelToTask("t", rec, 1)
scSim.InitializeSimulation()
scSim.ConfigureStopTime(macros.sec2nano(simDays * 86400.0))
scSim.ExecuteSimulation()
times = np.array(rec.times()) * macros.NANO2SEC / 86400.0
R = np.array(rec.r_BN_N)
V = np.array(rec.v_BN_N)
a = np.zeros(len(times))
e = np.zeros(len(times))
inc = np.zeros(len(times))
for k in range(len(times)):
oe = orbitalMotion.rv2elem(mu, R[k], V[k])
a[k] = oe.a
e[k] = oe.e
inc[k] = np.degrees(oe.i)
if useEarthRotation:
gravFactory.unloadSpiceKernels()
return times, a, e, inc
[docs]
def run(show_plots, simDays=30.0, gravDeg=20, stepSeconds=10.0):
"""Propagate the reference orbit with and without Earth rotation and compare.
Returns the ratio of the non-rotating to rotating eccentricity excursion
(which should be large, demonstrating the spurious drift) and the figure list.
Args:
show_plots (bool): show the matplotlib plots.
simDays (float): propagation duration in days.
gravDeg (int): spherical-harmonic degree/order of the GGM03S field.
stepSeconds (float): integration time step in seconds.
"""
sampleMinutes = 60.0 # [min] output sampling period
tR, aR, eR, iR = propagate(True, simDays, gravDeg, sampleMinutes, stepSeconds)
tN, aN, eN, iN = propagate(False, simDays, gravDeg, sampleMinutes, stepSeconds)
# Measure the *secular* drift of the eccentricity, i.e. the net change of its
# daily mean from the first to the last day. Averaging over a day removes the
# short-period (orbital) oscillation, which is present in both cases and would
# otherwise mask the secular trend over short runs. The non-rotating case
# drifts secularly; the rotating (physically correct) case stays bounded.
windowDays = 1.0 # [day] averaging window for the daily mean
def secularDrift(t, series):
first = series[t <= (t[0] + windowDays)].mean()
last = series[t >= (t[-1] - windowDays)].mean()
return abs(last - first)
eDriftRot = secularDrift(tR, eR)
eDriftNonRot = secularDrift(tN, eN)
iDriftRot = secularDrift(tR, iR)
iDriftNonRot = secularDrift(tN, iN)
driftRatio = eDriftNonRot / eDriftRot if eDriftRot > 0 else np.inf
print(f"eccentricity excursion: rotating={eDriftRot:.3e}, non-rotating={eDriftNonRot:.3e}")
print(f"inclination excursion [deg]: rotating={iDriftRot:.4f}, non-rotating={iDriftNonRot:.4f}")
figureList = {}
if plt is not None:
plt.close("all")
plt.figure(1)
plt.plot(tR, eR, color="tab:blue", label="Earth rotating (consistent)")
plt.plot(tN, eN, color="tab:red", linestyle="--", label="Earth non-rotating")
plt.xlabel("time [days]")
plt.ylabel("eccentricity [-]")
plt.legend(loc="upper left")
plt.grid(True)
figureList["scenarioOrbitConsistencyVerification_ecc"] = plt.figure(1)
plt.figure(2)
plt.plot(tR, iR, color="tab:blue", label="Earth rotating (consistent)")
plt.plot(tN, iN, color="tab:red", linestyle="--", label="Earth non-rotating")
plt.xlabel("time [days]")
plt.ylabel("inclination [deg]")
plt.legend(loc="upper left")
plt.grid(True)
figureList["scenarioOrbitConsistencyVerification_inc"] = plt.figure(2)
if show_plots:
plt.show()
plt.close("all")
return driftRatio, figureList
if __name__ == "__main__":
run(True, simDays=90.0, gravDeg=20)