#
# 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.
#
r"""
This scenario mirrors ``examples/mujoco/scenarioStochasticDrag.py`` but uses
the ``spacecraft.Spacecraft`` dynamics model (``spacecraft.h`` based) together
with ``MeanRevertingNoiseStateEffector``.
The atmospheric density used by drag is:
.. math::
\rho_{\text{stoch}} = \rho_{\text{exp}}(1 + \delta_\rho)
where :math:`\delta_\rho` is an Ornstein-Uhlenbeck process implemented as a
state in the state effector.
Alternatively, running with ``useIgbm=True`` uses :ref:`igbmNoiseStateEffector`
instead: the factor :math:`(1 + \delta_\rho)` then follows an inhomogeneous
geometric Brownian motion (multiplicative noise). Its exact process is positive
and mean-reverting, so it is a natural model for a positive density factor;
:ref:`dragDynamicEffector` clamps the corrected density to be non-negative to
guard against the rare non-positive value an explicit integrator can produce.
Both models are configured with the same stationary standard deviation and time
constant, so the two runs are statistically comparable.
Illustration of Simulation Results
----------------------------------
The following images illustrate a possible simulation result.
The orbit is plotted in the orbital plane:
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_orbit.svg
:align: center
The altitude as a function of time is plotted.
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_altitude.svg
:align: center
The atmospheric density as a function of altitude is plotted in lin-log space.
This shows two lines: the deterministic, exponential density (should appear
linear); and the stochastic density.
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_density.svg
:align: center
The atmospheric density correction, which should have a standard deviation
of 0.15.
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_densityDiff.svg
:align: center
The magnitude of drag force over time is plotted in lin-log space.
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_drag.svg
:align: center
Illustration of Simulation Results with IGBM
--------------------------------------------
The following images illustrate a possible simulation result with ``useIgbm=True``
(figure names gain an ``igbm`` tag). The overall trajectory is statistically similar to
the OU case since both processes share the same stationary standard deviation and time
constant; the qualitative difference is in the density factor, which is a positive
mean-reverting multiplicative factor under IGBM.
The atmospheric density as a function of altitude:
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_igbm_density.svg
:align: center
The atmospheric density correction :math:`X - 1`, which should have a standard deviation
of 0.15 about a positive mean-reverting factor :math:`X`:
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_igbm_densityDiff.svg
:align: center
The magnitude of drag force over time:
.. image:: /_images/Scenarios/scenarioStochasticDragSpacecraft_igbm_drag.svg
:align: center
"""
import os
import numpy as np
import matplotlib.pyplot as plt
from typing import Optional
from Basilisk.simulation import spacecraft
from Basilisk.simulation import svIntegrators
from Basilisk.simulation import exponentialAtmosphere
from Basilisk.simulation import dragDynamicEffector
from Basilisk.simulation import meanRevertingNoiseStateEffector
from Basilisk.simulation import igbmNoiseStateEffector
from Basilisk.simulation import zeroWindModel
from Basilisk.utilities import SimulationBaseClass
from Basilisk.utilities import macros
from Basilisk.utilities import orbitalMotion
from Basilisk.utilities import simIncludeGravBody
from Basilisk.utilities import simSetPlanetEnvironment
fileName = os.path.basename(os.path.splitext(__file__)[0])
[docs]
def run(showPlots: bool = False, rngSeed: Optional[int] = None, useWind: bool = False,
useIgbm: bool = False):
"""
Run the spacecraft-based stochastic drag scenario.
Args:
showPlots: If True, display figures.
rngSeed: Optional stochastic integrator seed for reproducibility.
useWind (bool): If True, link a ``ZeroWindModel`` via SPICE so drag is computed against
the atmosphere-relative velocity. If False (default), the inertial spacecraft
velocity is used directly.
useIgbm (bool): If True, model the density correction factor with an inhomogeneous
geometric Brownian motion (multiplicative noise, strictly positive factor)
instead of the additive Ornstein-Uhlenbeck correction. Defaults to False.
Returns:
Dict of matplotlib figure handles.
"""
initialAlt = 250 # [km]
planet = simIncludeGravBody.BODY_DATA["earth"]
# Match the same initial orbit used by scenarioStochasticDrag.py
oe = orbitalMotion.ClassicElements()
oe.a = planet.radEquator + initialAlt * 1000.0 # [m]
oe.e = 0.0
oe.i = 33.3 * macros.D2R
oe.Omega = 48.2 * macros.D2R
oe.omega = 347.8 * macros.D2R
oe.f = 85.3 * macros.D2R
rN, vN = orbitalMotion.elem2rv(planet.mu, oe)
oe = orbitalMotion.rv2elem(planet.mu, rN, vN)
orbitPeriod = 2.0 * np.pi / np.sqrt(planet.mu / oe.a**3)
dt = 10.0 # [s]
tf = 7.1 * orbitPeriod # [s]
simTaskName = "simTask"
simProcessName = "simProcess"
scSim = SimulationBaseClass.SimBaseClass()
dynProcess = scSim.CreateNewProcess(simProcessName)
dynProcess.addTask(scSim.CreateNewTask(simTaskName, macros.sec2nano(dt)))
# Spacecraft setup (point-mass equivalent to the MuJoCo cannonball body)
scObject = spacecraft.Spacecraft()
scObject.ModelTag = "bsk-Sat"
scObject.hub.mHub = 1.0 # [kg]
scObject.hub.IHubPntBc_B = np.identity(3) # [kg*m^2]
scObject.hub.r_CN_NInit = rN
scObject.hub.v_CN_NInit = vN
scSim.AddModelToTask(simTaskName, scObject)
# Stochastic integrator for states with diffusion. This scenario looks at a single
# sample trajectory, so a STRONG integrator is used (strong order = pathwise
# accuracy; weak order only controls statistics across many samples). SOSRA is the
# recommended strong method for the OU correction (additive noise); the IGBM factor
# has multiplicative noise, which SOSRA does not support, so it uses SRIW1 (strong
# order 1.5 for diagonal/scalar noise).
if useIgbm:
integrator = svIntegrators.svStochasticIntegratorSRIW1(scObject)
else:
integrator = svIntegrators.svStochasticIntegratorSOSRA(scObject)
if rngSeed is not None:
integrator.setRNGSeed(rngSeed)
scObject.setIntegrator(integrator)
# Gravity model (point-mass Earth)
gravFactory = simIncludeGravBody.gravBodyFactory()
earth = gravFactory.createEarth()
earth.isCentralBody = True
gravFactory.addBodiesTo(scObject)
# Deterministic atmosphere
atmo = exponentialAtmosphere.ExponentialAtmosphere()
atmo.ModelTag = "ExpAtmo"
simSetPlanetEnvironment.exponentialAtmosphere(atmo, "earth")
atmo.addSpacecraftToModel(scObject.scStateOutMsg)
scSim.AddModelToTask(simTaskName, atmo)
# Stochastic atmospheric density correction state effector
if useIgbm:
# Multiplicative IGBM: the factor (1 + delta) reverts to 1 and stays positive
stochasticAtmo = igbmNoiseStateEffector.IgbmNoiseStateEffector()
stochasticAtmo.setMean(1.0)
else:
# Additive Ornstein-Uhlenbeck correction centered at zero
stochasticAtmo = meanRevertingNoiseStateEffector.MeanRevertingNoiseStateEffector()
# Both models are configured in stationary form with the same statistics
stochasticAtmo.ModelTag = "StochasticExpAtmo"
stochasticAtmo.setStationaryStd(0.15)
stochasticAtmo.setTimeConstant(1.8 * 60.0) # [s]
scObject.addStateEffector(stochasticAtmo)
# Cannonball drag model
drag = dragDynamicEffector.DragDynamicEffector()
drag.ModelTag = "DragEff"
drag.coreParams.dragCoeff = 2.2
drag.coreParams.projectedArea = 10.0 # [m^2]
drag.atmoDensInMsg.subscribeTo(atmo.envOutMsgs[0])
drag.densityCorrectionStateName = stochasticAtmo.getStateName()
scObject.addDynamicEffector(drag)
scSim.AddModelToTask(simTaskName, drag)
if useWind:
spiceObject = gravFactory.createSpiceInterface(
time="2020 MAY 21 18:28:03 (UTC)",
)
spiceObject.zeroBase = "Earth"
scSim.AddModelToTask(simTaskName, spiceObject, -1)
windModel = zeroWindModel.ZeroWindModel()
windModel.ModelTag = "ZeroWind"
windModel.planetPosInMsg.subscribeTo(spiceObject.planetStateOutMsgs[0])
windModel.addSpacecraftToModel(scObject.scStateOutMsg)
scSim.AddModelToTask(simTaskName, windModel)
drag.windVelInMsg.subscribeTo(windModel.envOutMsgs[0])
# Recorders
stateRecorder = scObject.scStateOutMsg.recorder()
deterministicDensityRecorder = atmo.envOutMsgs[0].recorder()
dragRecorder = drag.logger("forceExternal_B")
scSim.AddModelToTask(simTaskName, stateRecorder)
scSim.AddModelToTask(simTaskName, deterministicDensityRecorder)
scSim.AddModelToTask(simTaskName, dragRecorder)
scSim.InitializeSimulation()
scSim.ConfigureStopTime(macros.sec2nano(tf))
scSim.ExecuteSimulation()
figures = plotOrbits(
timeAxis=stateRecorder.times(),
posData=stateRecorder.r_BN_N,
velData=stateRecorder.v_BN_N,
dragForce=dragRecorder.forceExternal_B,
deterministicDenseData=deterministicDensityRecorder.neutralDensity,
oe=oe,
mu=planet.mu,
planetRadius=planet.radEquator,
dragCoeff=drag.coreParams.dragCoeff,
dragArea=drag.coreParams.projectedArea,
figureTag="igbm" if useIgbm else "",
)
if showPlots:
plt.show()
return figures
[docs]
def plotOrbits(timeAxis, posData, velData, dragForce, deterministicDenseData, oe, mu, planetRadius, dragCoeff,
dragArea, figureTag=""):
"""Plot orbit, altitude, atmosphere, correction, and drag.
``figureTag`` is an optional tag inserted into figure names (e.g. "igbm"), so the
variants get their own documentation images.
"""
figureList = {}
baseName = fileName + (f"_{figureTag}" if figureTag else "")
timeHours = timeAxis * macros.NANO2HOUR
figureList[baseName + "_orbit"], ax = plt.subplots()
ax.axis("equal")
ax.add_artist(plt.Circle((0, 0), planetRadius / 1000.0, color="#008800"))
rData = []
fData = []
for idx in range(len(posData)):
oeData = orbitalMotion.rv2elem(mu, posData[idx], velData[idx])
rData.append(oeData.rmag)
fData.append(oeData.f + oeData.omega - oe.omega)
ax.plot(
np.array(rData) * np.cos(fData) / 1000.0,
np.array(rData) * np.sin(fData) / 1000.0,
color="#aa0000",
linewidth=1.0
)
ax.set_xlabel("$i_e$ Cord. [km]")
ax.set_ylabel("$i_p$ Cord. [km]")
figureList[baseName + "_altitude"], ax = plt.subplots()
ax.ticklabel_format(useOffset=False, style="plain")
alt = (np.array(rData) - planetRadius) / 1000.0
ax.plot(timeHours, alt)
ax.set_xlabel("$t$ [h]")
ax.set_ylabel("Altitude [km]")
velNorm = np.linalg.norm(velData, axis=1)
dragNorm = np.linalg.norm(dragForce, axis=1)
denseData = 2.0 * dragNorm / (dragCoeff * dragArea * np.maximum(velNorm**2, 1e-30))
figureList[baseName + "_density"], ax = plt.subplots()
ax.semilogy(alt[1:], denseData[1:], label="Stochastic")
ax.semilogy(alt[1:], deterministicDenseData[1:], label="Exponential")
ax.legend(loc="upper right")
ax.set_xlabel("Altitude [km]")
ax.set_ylabel("$\\rho$ [kg/m$^3$]")
figureList[baseName + "_densityDiff"], ax = plt.subplots()
ax.plot(timeHours[1:], (denseData / deterministicDenseData)[1:] - 1.0)
ax.set_xlabel("Time [hr]")
ax.set_ylabel(r"$(\rho_{stoch} / \rho_{exp}) - 1$ [-]")
figureList[baseName + "_drag"], ax = plt.subplots()
ax.semilogy(timeHours[1:], np.linalg.norm(dragForce, 2, 1)[1:])
ax.set_xlabel("$t$ [hr]")
ax.set_ylabel("$|F_{drag}|$ [N]")
return figureList
if __name__ == "__main__":
run(True, useWind = False)