Consume wind msgs in all drag effectors

docs/source/Support/bskReleaseNotesSnippets/1325-atmo_rel_vel_all_drag.rst

@@ -0,0 +1,2 @@
+- ``dragDynamicEffector``, ``facetDragDynamicEffector``, and ``cannonballDrag`` now accept an optional ``windVelInMsg`` (:ref:`WindMsgPayload`) to compute atmosphere-relative drag velocity; supersedes the former ``useAtmosphereRelativeVelocity``/``planetOmega_N`` flags.
+- Updated :ref:`scenarioDragDeorbit`, :ref:`scenarioDragRendezvous`, :ref:`scenarioAerocapture`, and :ref:`scenarioStochasticDragSpacecraft` to illustrate optional use of wind velocity input.

examples/scenarioAerocapture.py

@@ -86,6 +86,7 @@
 # Used to get the location of supporting data.
 from Basilisk import __path__
 from Basilisk.simulation import dragDynamicEffector
+from Basilisk.simulation import zeroWindModel
 
 # import simulation related support
 from Basilisk.simulation import spacecraft
@@ -155,13 +156,16 @@ def sph2rv(xxsph):
     return rvec_N, uvec_N
 
 
-def run(show_plots, planetCase):
+def run(show_plots, planetCase, useWind=False):
     """
     The scenarios can be run with the followings setups parameters:
 
     Args:
         show_plots (bool): Determines if the script should display plots
-        planetCase (string): Specify if a `Mars` or `Earth` arrival is simulated
+        planetCase (string): Specify if a ``Mars`` or ``Earth`` arrival is simulated
+        useWind (bool): If True and ``planetCase == 'Earth'``, link a ``ZeroWindModel`` via SPICE
+            so drag is computed against the atmosphere-relative velocity.  If False (default),
+            or for Mars, the inertial spacecraft velocity is used directly.
 
     """
 
@@ -255,6 +259,20 @@ def run(show_plots, planetCase):
     # attach gravity model to spacecraft
     gravFactory.addBodiesTo(scObject)
 
+    if useWind and planetCase == "Earth":
+        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)
+        dragEffector.windVelInMsg.subscribeTo(windModel.envOutMsgs[0])
+
     if planetCase == "Earth":
         r = 6503 * 1000
         u = 11.2 * 1000
@@ -398,4 +416,4 @@ def run(show_plots, planetCase):
 
 
 if __name__ == "__main__":
-    run(True, "Mars")  # planet arrival case, can be Earth or Mars
+    run(True, "Mars", False)  # planet arrival case, can be Earth or Mars

examples/scenarioDragDeorbit.py

@@ -26,6 +26,12 @@
 :ref:`dragDynamicEffector` dynamics module. The simulation is executed until the altitude falls below some threshold,
 using a terminal event handler.
 
+An optional wind model (:ref:`zeroWindModel`) can be enabled via the ``useWind`` argument (default: ``False``).
+When active, SPICE is loaded with the default ``IAU_EARTH`` frame and drag is computed against the
+atmosphere-relative velocity rather than the inertial spacecraft velocity.  This is a more accurate treatment
+for scenarios where Earth's rotation is relevant, but adds a SPICE dependency.  When ``useWind=False``, the
+inertial spacecraft velocity is used directly, which is the conventional approximation for simple deorbit studies.
+
 The script is found in the folder ``basilisk/examples`` and executed by using::
 
       python3 scenarioDragDeorbit.py
@@ -72,6 +78,24 @@
 
     dragEffector.atmoDensInMsg.subscribeTo(atmo.envOutMsgs[0])
 
+Optionally (``useWind=True``), a :ref:`zeroWindModel` can be linked so drag is computed against the
+atmosphere-relative velocity.  SPICE is loaded with the default ``IAU_EARTH`` frame so that
+``J20002Pfix_dot`` is populated; :ref:`windBase` reads those derivatives to derive the Earth rotation
+rate, giving a physically consistent co-rotation::
+
+    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)
+    dragEffector.windVelInMsg.subscribeTo(windModel.envOutMsgs[0])
+
 Illustration of Simulation Results
 ----------------------------------
 
@@ -142,6 +166,7 @@
     exponentialAtmosphere,
     msisAtmosphere,
     spacecraft,
+    zeroWindModel,
 )
 from Basilisk.utilities import (
     SimulationBaseClass,
@@ -157,7 +182,7 @@
 fileName = os.path.basename(os.path.splitext(__file__)[0])
 
 
-def run(show_plots, initialAlt=250, deorbitAlt=100, model="exponential"):
+def run(show_plots, initialAlt=250, deorbitAlt=100, model="exponential", useWind=False):
     """
     Initialize a satellite with drag and propagate until it falls below a deorbit altitude. Note that an excessively
     low deorbit_alt can lead to intersection with the Earth prior to deorbit being detected, causing some terms to blow
@@ -168,6 +193,8 @@ def run(show_plots, initialAlt=250, deorbitAlt=100, model="exponential"):
         initialAlt (float): Starting altitude in km
         deorbitAlt (float): Terminal altitude in km
         model (str): ["exponential", "msis"]
+        useWind (bool): If True, link a :ref:`zeroWindModel` so drag is computed against the atmosphere-relative
+            velocity.  If False (default), the inertial spacecraft velocity is used directly.
 
     Returns:
         Dictionary of figure handles
@@ -238,6 +265,25 @@ def run(show_plots, initialAlt=250, deorbitAlt=100, model="exponential"):
     mu = planet.mu
     gravFactory.addBodiesTo(scObject)
 
+    # Optionally link a zero-wind model so drag is computed against the atmosphere-relative velocity.
+    # SPICE is loaded with the default IAU_EARTH frame so that J20002Pfix_dot is populated and
+    # WindBase can derive the true Earth rotation rate from it, keeping the co-rotation consistent
+    # with the actual planet orientation.  Omitting this (useWind=False) falls back to using the
+    # inertial spacecraft velocity directly.
+    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)
+        dragEffector.windVelInMsg.subscribeTo(windModel.envOutMsgs[0])
+
     # Set up a circular orbit using classical orbit elements
     oe = orbitalMotion.ClassicElements()
     oe.a = planet.radEquator + initialAlt * 1000  # meters
@@ -383,5 +429,6 @@ def register_fig(i):
         show_plots=True,
         initialAlt=250,
         deorbitAlt=100,
-        model="msis"   # "msis", "exponential"
+        model="msis",   # "msis", "exponential"
+        useWind=False
     )

examples/scenarioDragRendezvous.py

@@ -85,6 +85,7 @@
     facetDragDynamicEffector,
     simpleNav,
     exponentialAtmosphere,
+    zeroWindModel,
 )
 from Basilisk.utilities import RigidBodyKinematics as rbk
 from Basilisk.utilities import SimulationBaseClass
@@ -255,7 +256,7 @@ def setup_spacecraft_plant(rN, vN, modelName):
 
 
 def drag_simulator(
-    altOffset, trueAnomOffset, densMultiplier, ctrlType="lqr", useJ2=False
+    altOffset, trueAnomOffset, densMultiplier, ctrlType="lqr", useJ2=False, useWind=False
 ):
     """
     Basilisk simulation of a two-spacecraft rendezvous using relative-attitude driven differential drag. Includes
@@ -265,6 +266,11 @@ def drag_simulator(
     Args:
         altOffset - double - deputy altitude offset from the chief ('x' hill direction), meters
         trueAnomOffset - double - deputy true anomaly difference from the chief ('y' direction), degrees
+        densMultiplier - double - multiplicative scale factor applied to the exponential atmosphere base density
+        ctrlType (str): Control law type; ``'lqr'`` for static LQR or ``'dtv'`` for desensitized time-varying gain
+        useJ2 (bool): If True, include J2 spherical harmonic gravity perturbation
+        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.
     """
 
     startTime = time.time()
@@ -341,6 +347,22 @@ def drag_simulator(
     atmosphere.addSpacecraftToModel(chiefSc.scStateOutMsg)
     chiefDrag.atmoDensInMsg.subscribeTo(atmosphere.envOutMsgs[-1])
 
+    if useWind:
+        spiceObject = gravFactory.createSpiceInterface(
+            time="2020 MAY 21 18:28:03 (UTC)",
+        )
+        spiceObject.zeroBase = "Earth"
+        scSim.AddModelToTask(dynTaskName, spiceObject, 999)
+
+        windModel = zeroWindModel.ZeroWindModel()
+        windModel.ModelTag = "ZeroWind"
+        windModel.planetPosInMsg.subscribeTo(spiceObject.planetStateOutMsgs[0])
+        windModel.addSpacecraftToModel(depSc.scStateOutMsg)
+        windModel.addSpacecraftToModel(chiefSc.scStateOutMsg)
+        scSim.AddModelToTask(dynTaskName, windModel, 919)
+        depDrag.windVelInMsg.subscribeTo(windModel.envOutMsgs[0])
+        chiefDrag.windVelInMsg.subscribeTo(windModel.envOutMsgs[1])
+
     #   Add all dynamics stuff to dynamics task
     scSim.AddModelToTask(dynTaskName, atmosphere, 920)
     # scSim.AddModelToTask(dynTaskName, ephemConverter, 921)
@@ -470,10 +492,10 @@ def drag_simulator(
 
 
 def run(
-    show_plots, altOffset, trueAnomOffset, densMultiplier, ctrlType="lqr", useJ2=False
+    show_plots, altOffset, trueAnomOffset, densMultiplier, ctrlType="lqr", useJ2=False, useWind=False
 ):
     results = drag_simulator(
-        altOffset, trueAnomOffset, densMultiplier, ctrlType=ctrlType, useJ2=useJ2
+        altOffset, trueAnomOffset, densMultiplier, ctrlType=ctrlType, useJ2=useJ2, useWind=useWind
     )
 
     timeData = results["dynTimeData"]
@@ -598,4 +620,5 @@ def run(
         1,  #    Density multiplier (nondimensional)
         ctrlType="lqr",
         useJ2=False,
+        useWind=False,
     )

examples/scenarioStochasticDragSpacecraft.py

@@ -72,6 +72,7 @@
 from Basilisk.simulation import exponentialAtmosphere
 from Basilisk.simulation import dragDynamicEffector
 from Basilisk.simulation import meanRevertingNoiseStateEffector
+from Basilisk.simulation import zeroWindModel
 from Basilisk.utilities import SimulationBaseClass
 from Basilisk.utilities import macros
 from Basilisk.utilities import orbitalMotion
@@ -80,13 +81,16 @@
 
 fileName = os.path.basename(os.path.splitext(__file__)[0])
 
-def run(showPlots: bool = False, rngSeed: Optional[int] = None):
+def run(showPlots: bool = False, rngSeed: Optional[int] = None, useWind: 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.
     Returns:
         Dict of matplotlib figure handles.
     """
@@ -159,6 +163,20 @@ def run(showPlots: bool = False, rngSeed: Optional[int] = None):
     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()
@@ -247,4 +265,4 @@ def plotOrbits(timeAxis, posData, velData, dragForce, deterministicDenseData, oe
 
 
 if __name__ == "__main__":
-    run(True)
+    run(True, useWind = False)