[#1490] Add zeroBase origin support to planetEphemeris

schaubh · schaubh · commit ddec43c550db · 2026-07-20T17:48:58.000-06:00
Add an optional zeroBase setting to PlanetEphemeris so generated planet
states can be expressed relative to another body propagated by the same
module instance. Preserve the existing heliocentric behavior when
zeroBase is left empty.

Resolve the selected body during Reset using case-insensitive name
matching and report an initialization error when the requested body is
not configured in the module. Compute all heliocentric states before
subtracting the selected body's position and velocity, ensuring the
translation is independent of body ordering. Leave planet names,
timestamps, attitudes, attitude rates, and orientation flags unchanged.

Expand the PlanetEphemeris unit tests to verify:

- the default empty zeroBase value
- position and velocity translation at multiple simulation times
- zero position and velocity for the selected base body
- preservation of attitude and message metadata
- case-insensitive body-name matching
- rejection of bodies from outside the module instance

Document zeroBase configuration, default behavior, and the requirement
that the reference body be propagated by the same PlanetEphemeris
instance. Clarify that changing the output origin does not introduce
hierarchical or mutually coupled orbital propagation.

Update scenarioCustomGravBody to use Itokawa as the output origin,
placing the planet and spacecraft states in a common asteroid-centered
translational frame.
diff --git a/examples/scenarioCustomGravBody.py b/examples/scenarioCustomGravBody.py
@@ -53,9 +53,12 @@
 be customized.  The gravity body ephemeris states are connected to the :ref:`planetEphemeris` planet
 state output messages.
 
-Finally, the recorded states will all be relative to the inertial origin at the sun.  :ref:`planetEphemeris` does not
-have the ``zeroBase`` capability as :ref:`spiceInterface` has.  This script also records the asteroid
-states so that the plot is done of the spacecraft motion relative to the asteroid.
+The ephemeris outputs remain heliocentric because the Sun gravity body is not
+connected to an ephemeris message.  Its default zero state therefore represents
+the Sun at the heliocentric origin.  Setting ``planetEphemeris.zeroBase`` to
+Itokawa would incorrectly place both the Sun and Itokawa at the origin.  To plot
+the spacecraft motion relative to the asteroid, the script instead records the
+Itokawa state and subtracts it from the spacecraft state during post-processing.
 
 The simulation executes and shows a plot of the spacecraft motion relative to the asteroid.
 
@@ -148,6 +151,9 @@ def run(show_plots):
 
     # specify celestial object orbit
     gravBodyEphem.planetElements = planetEphemeris.classicElementVector([oeAsteroid, oeEarth])
+    # Keep these outputs heliocentric. The unlinked Sun gravity body below uses
+    # its default zero state to represent the Sun at the heliocentric origin.
+    # The asteroid state is subtracted from the spacecraft state when plotting.
     # specify celestial object orientation
     gravBodyEphem.rightAscension = planetEphemeris.DoubleVector([0.0 * macros.D2R, 0.0 * macros.D2R])
     gravBodyEphem.declination = planetEphemeris.DoubleVector([0.0 * macros.D2R, 0.0 * macros.D2R])
diff --git a/src/simulation/environment/planetEphemeris/_UnitTest/test_planetEphemeris.py b/src/simulation/environment/planetEphemeris/_UnitTest/test_planetEphemeris.py
@@ -46,6 +46,189 @@
 from Basilisk.architecture.bskLogging import BasiliskError
 
 
+def _make_zero_base_test_elements():
+    """Create two heliocentric element sets for the zero-base tests.
+
+    Returns:
+        planetEphemeris.classicElementVector: Independent Earth and Moon-like
+        heliocentric element sets.
+    """
+    earthElements = planetEphemeris.ClassicElements()
+    earthElements.a = orbitalMotion.SMA_EARTH * 1000.0  # [m]
+    earthElements.e = 0.01  # [-]
+    earthElements.i = 10.0 * macros.D2R  # [deg] -> [rad]
+    earthElements.Omega = 30.0 * macros.D2R  # [deg] -> [rad]
+    earthElements.omega = 20.0 * macros.D2R  # [deg] -> [rad]
+    earthElements.f = 90.0 * macros.D2R  # [deg] -> [rad]
+
+    moonElements = planetEphemeris.ClassicElements()
+    moonElements.a = 1.01 * orbitalMotion.SMA_EARTH * 1000.0  # [m]
+    moonElements.e = 0.02  # [-]
+    moonElements.i = 5.0 * macros.D2R  # [deg] -> [rad]
+    moonElements.Omega = 110.0 * macros.D2R  # [deg] -> [rad]
+    moonElements.omega = 220.0 * macros.D2R  # [deg] -> [rad]
+    moonElements.f = 180.0 * macros.D2R  # [deg] -> [rad]
+
+    return planetEphemeris.classicElementVector(
+        [earthElements, moonElements]
+    )
+
+
+def _configure_zero_base_test_module(module):
+    """Configure a ``PlanetEphemeris`` module with two oriented bodies.
+
+    Args:
+        module (planetEphemeris.PlanetEphemeris): Module to configure.
+    """
+    module.setPlanetNames(
+        planetEphemeris.StringVector(["earth", "moon"])
+    )
+    module.planetElements = _make_zero_base_test_elements()
+    module.rightAscension = planetEphemeris.DoubleVector(
+        [0.0, 10.0 * macros.D2R]
+    )  # [rad]
+    module.declination = planetEphemeris.DoubleVector(
+        [90.0 * macros.D2R, 80.0 * macros.D2R]
+    )  # [rad]
+    module.lst0 = planetEphemeris.DoubleVector(
+        [0.0, 30.0 * macros.D2R]
+    )  # [rad]
+    module.rotRate = planetEphemeris.DoubleVector(
+        [1.0e-5, 2.0e-5]
+    )  # [rad/s]
+
+
+def test_zero_base_translation():
+    """Verify that ``zeroBase`` translates positions and velocities only.
+
+    Validation Test Description
+    ---------------------------
+    Two identical modules propagate the same oriented bodies.  One retains the
+    default heliocentric origin while the other selects Earth as ``zeroBase``.
+    The relative outputs are compared with differences of the heliocentric
+    outputs at three simulation times.
+
+    Test Parameter Discussion
+    -------------------------
+    Mixed-case ``"EaRtH"`` verifies case-insensitive body-name matching.  The
+    bodies use distinct orbits and attitudes so position, velocity, and
+    orientation behavior are all observable.
+
+    Description of Variables Being Tested
+    ---------------------------------------
+    The test checks ``PositionVector``, ``VelocityVector``, ``J20002Pfix``,
+    ``J20002Pfix_dot``, ``computeOrient``, and ``PlanetName``.
+    """
+    simulation = SimulationBaseClass.SimBaseClass()
+    process = simulation.CreateNewProcess("zeroBaseProcess")
+    taskStep = 0.5  # [s]
+    process.addTask(
+        simulation.CreateNewTask(
+            "zeroBaseTask",
+            macros.sec2nano(taskStep),
+        )
+    )
+
+    heliocentricModule = planetEphemeris.PlanetEphemeris()
+    heliocentricModule.ModelTag = "heliocentricPlanetEphemeris"
+    assert heliocentricModule.zeroBase == ""
+    _configure_zero_base_test_module(heliocentricModule)
+
+    relativeModule = planetEphemeris.PlanetEphemeris()
+    relativeModule.ModelTag = "relativePlanetEphemeris"
+    _configure_zero_base_test_module(relativeModule)
+    relativeModule.zeroBase = "EaRtH"
+
+    simulation.AddModelToTask("zeroBaseTask", heliocentricModule, 10)
+    simulation.AddModelToTask("zeroBaseTask", relativeModule, 10)
+
+    heliocentricRecorders = [
+        message.recorder() for message in heliocentricModule.planetOutMsgs
+    ]
+    relativeRecorders = [
+        message.recorder() for message in relativeModule.planetOutMsgs
+    ]
+    for recorder in heliocentricRecorders + relativeRecorders:
+        simulation.AddModelToTask("zeroBaseTask", recorder, 0)
+
+    simulation.InitializeSimulation()
+    simulationDuration = 1.0  # [s]
+    simulation.ConfigureStopTime(macros.sec2nano(simulationDuration))
+    simulation.ExecuteSimulation()
+
+    basePosition = heliocentricRecorders[0].PositionVector
+    baseVelocity = heliocentricRecorders[0].VelocityVector
+    positionTolerance = 1.0e-6  # [m]
+    velocityTolerance = 1.0e-10  # [m/s]
+    attitudeTolerance = 1.0e-14  # [-]
+    for bodyIndex in range(2):
+        expectedPosition = (
+            heliocentricRecorders[bodyIndex].PositionVector - basePosition
+        )
+        expectedVelocity = (
+            heliocentricRecorders[bodyIndex].VelocityVector - baseVelocity
+        )
+        np.testing.assert_allclose(
+            relativeRecorders[bodyIndex].PositionVector,
+            expectedPosition,
+            rtol=0.0,
+            atol=positionTolerance,
+        )
+        np.testing.assert_allclose(
+            relativeRecorders[bodyIndex].VelocityVector,
+            expectedVelocity,
+            rtol=0.0,
+            atol=velocityTolerance,
+        )
+        np.testing.assert_allclose(
+            relativeRecorders[bodyIndex].J20002Pfix,
+            heliocentricRecorders[bodyIndex].J20002Pfix,
+            rtol=0.0,
+            atol=attitudeTolerance,
+        )
+        np.testing.assert_allclose(
+            relativeRecorders[bodyIndex].J20002Pfix_dot,
+            heliocentricRecorders[bodyIndex].J20002Pfix_dot,
+            rtol=0.0,
+            atol=attitudeTolerance,
+        )
+        np.testing.assert_array_equal(
+            relativeRecorders[bodyIndex].computeOrient,
+            heliocentricRecorders[bodyIndex].computeOrient,
+        )
+
+    assert relativeModule.planetOutMsgs[0].read().PlanetName == "earth"
+    assert relativeModule.planetOutMsgs[1].read().PlanetName == "moon"
+
+
+def test_zero_base_requires_local_body():
+    """Verify that ``zeroBase`` rejects bodies owned by another instance.
+
+    Validation Test Description
+    ---------------------------
+    A module configured with Earth and Moon selects Mars as its translational
+    origin.  Initialization must stop with a ``BasiliskError`` because Mars is
+    not propagated by that module instance.
+    """
+    simulation = SimulationBaseClass.SimBaseClass()
+    process = simulation.CreateNewProcess("invalidZeroBaseProcess")
+    taskStep = 0.5  # [s]
+    process.addTask(
+        simulation.CreateNewTask(
+            "invalidZeroBaseTask",
+            macros.sec2nano(taskStep),
+        )
+    )
+
+    module = planetEphemeris.PlanetEphemeris()
+    _configure_zero_base_test_module(module)
+    module.zeroBase = "mars"
+    simulation.AddModelToTask("invalidZeroBaseTask", module)
+
+    with pytest.raises(BasiliskError, match="not one of the bodies configured"):
+        simulation.InitializeSimulation()
+
+
 # Uncomment this line is this test is to be skipped in the global unit test run, adjust message as needed.
 # @pytest.mark.skipif(conditionstring)
 # Uncomment this line if this test has an expected failure, adjust message as needed.
@@ -274,6 +457,8 @@ def planetEphemerisTest(show_plots, setRAN, setDEC, setLST, setRate):
 # stand-along python script
 #
 if __name__ == "__main__":
+    test_zero_base_translation()
+    test_zero_base_requires_local_body()
     test_module(
                  False,           # show plots flag
                  True,           # setRAN
diff --git a/src/simulation/environment/planetEphemeris/planetEphemeris.cpp b/src/simulation/environment/planetEphemeris/planetEphemeris.cpp
@@ -17,13 +17,27 @@
 
  */
 #include "simulation/environment/planetEphemeris/planetEphemeris.h"
+#include <algorithm>
+#include <cctype>
 #include <cstdio>
 #include <iostream>
 #include "architecture/utilities/astroConstants.h"
 #include "architecture/utilities/linearAlgebra.h"
 #include "architecture/utilities/macroDefinitions.h"
 #include "architecture/utilities/rigidBodyKinematics.h"
 
+namespace {
+/*! Compare two planet names without regard to letter case. */
+bool planetNamesMatch(const std::string& firstName, const std::string& secondName)
+{
+    return firstName.size() == secondName.size()
+           && std::equal(firstName.begin(), firstName.end(), secondName.begin(),
+                         [](unsigned char firstCharacter, unsigned char secondCharacter) {
+                             return std::tolower(firstCharacter) == std::tolower(secondCharacter);
+                         });
+}
+} // namespace
+
 
 /*! This constructor initializes the variables.
  */
@@ -64,14 +78,30 @@ void PlanetEphemeris::setPlanetNames(std::vector<std::string> names)
 
 
 
-void PlanetEphemeris::Reset(uint64_t CurrenSimNanos)
+void PlanetEphemeris::Reset(uint64_t CurrentSimNanos)
 {
-    this->epochTime = CurrenSimNanos*NANO2SEC;
+    this->epochTime = CurrentSimNanos*NANO2SEC;
 
     /*! - do sanity checks that the vector arrays for planet names and ephemeris have the same length */
     if (this->planetElements.size() != this->planetNames.size()) {
-        bskLogger.bskError("Only %lu planet element sets provided, but %lu plane names are present.",
-                  this->planetElements.size(), this->planetNames.size());
+        bskLogger.bskError("Only %zu planet element sets provided, but %zu planet names are present.",
+                           this->planetElements.size(), this->planetNames.size());
+    }
+
+    /*! - Resolve the optional translational origin from the bodies owned by this module. */
+    this->zeroBaseIndex = -1;
+    if (!this->zeroBase.empty()) {
+        for (size_t c = 0; c < this->planetNames.size(); c++) {
+            if (planetNamesMatch(this->zeroBase, this->planetNames.at(c))) {
+                this->zeroBaseIndex = static_cast<int>(c);
+                break;
+            }
+        }
+        if (this->zeroBaseIndex < 0) {
+            bskLogger.bskError(
+                "PlanetEphemeris zeroBase '%s' is not one of the bodies configured in this module instance.",
+                this->zeroBase.c_str());
+        }
     }
 
     /*! - See if planet orientation information is set */
@@ -84,28 +114,32 @@ void PlanetEphemeris::Reset(uint64_t CurrenSimNanos)
     }
 
     if (computeAttitudeFlag) {
-        /*! - check that the right number of planet local sideral time angles are provided */
+        /*! - check that the right number of planet local sidereal time angles are provided */
         if (this->lst0.size() != this->planetNames.size()) {
-            bskLogger.bskError("Only %lu planet initial principal rotation angles provided, but %lu planet names are present.",
-                      this->lst0.size(), this->planetNames.size());
+            bskLogger.bskError(
+                "Only %zu planet initial principal rotation angles provided, but %zu planet names are present.",
+                this->lst0.size(), this->planetNames.size());
         }
 
         /*! - check that the right number of planet polar axis right ascension angles are provided */
         if (this->rightAscension.size() != this->planetNames.size()) {
-            bskLogger.bskError("Only %lu planet right ascension angles provided, but %lu planet names are present.",
-                      this->rightAscension.size(), this->planetNames.size());
+            bskLogger.bskError(
+                "Only %zu planet right ascension angles provided, but %zu planet names are present.",
+                this->rightAscension.size(), this->planetNames.size());
         }
 
         /*! - check that the right number of planet polar axis declination angles are provided */
         if (this->declination.size() != this->planetNames.size()) {
-            bskLogger.bskError("Only %lu planet declination angles provided, but %lu planet names are present.",
-                      this->declination.size(), this->planetNames.size());
+            bskLogger.bskError(
+                "Only %zu planet declination angles provided, but %zu planet names are present.",
+                this->declination.size(), this->planetNames.size());
         }
 
         /*! - check that the right number of planet polar rotation rates are provided */
         if (this->rotRate.size() != this->planetNames.size()) {
-            bskLogger.bskError("Only %lu planet rotation rates provided, but %lu planet names are present.",
-                      this->rotRate.size(), this->planetNames.size());
+            bskLogger.bskError(
+                "Only %zu planet rotation rates provided, but %zu planet names are present.",
+                this->rotRate.size(), this->planetNames.size());
         }
     }
 
@@ -133,6 +167,7 @@ void PlanetEphemeris::UpdateState(uint64_t CurrentSimNanos)
     double omega_NP_P[3];                   // [r/s] angular velocity of inertial frame relative to planet frame in planet frame components
     double tilde[3][3];                     // [] skew-symmetric matrix
     double theta_PN[3];                     // [rad] 3-2-3 planet orientation coordinates
+    std::vector<SpicePlanetStateMsgPayload> planetStates(this->planetNames.size());
 
 
 
@@ -146,7 +181,7 @@ void PlanetEphemeris::UpdateState(uint64_t CurrentSimNanos)
     for(it = this->planetNames.begin(), c=0; it != this->planetNames.end(); it++, c++)
     {
         //! - Create new planet output message copy
-        SpicePlanetStateMsgPayload newPlanet;
+        SpicePlanetStateMsgPayload& newPlanet = planetStates.at(c);
         newPlanet = this->planetOutMsgs.at(c)->zeroMsgPayload;
         //! - specify planet name in output message
         std::snprintf(newPlanet.PlanetName, sizeof(newPlanet.PlanetName), "%s", it->c_str());
@@ -183,9 +218,25 @@ void PlanetEphemeris::UpdateState(uint64_t CurrentSimNanos)
             m33SetIdentity(newPlanet.J20002Pfix); // set default zero orientation
         }
 
-        //! - write output message
-        this->planetOutMsgs.at(c)->write(&newPlanet, this->moduleID, CurrentSimNanos);
+    }
+
+    /*! - Translate every position and velocity to the selected zero-base origin. */
+    if (this->zeroBaseIndex >= 0) {
+        double zeroBasePosition_N[3];         // [m] zero-base heliocentric position
+        double zeroBaseVelocity_N[3];         // [m/s] zero-base heliocentric velocity
+        v3Copy(planetStates.at(this->zeroBaseIndex).PositionVector, zeroBasePosition_N);
+        v3Copy(planetStates.at(this->zeroBaseIndex).VelocityVector, zeroBaseVelocity_N);
+
+        for (auto& planetState : planetStates) {
+            v3Subtract(planetState.PositionVector, zeroBasePosition_N, planetState.PositionVector);
+            v3Subtract(planetState.VelocityVector, zeroBaseVelocity_N, planetState.VelocityVector);
+        }
+    }
 
+    /*! - Write the translated planet states after every absolute state is available. */
+    for (size_t planetIndex = 0; planetIndex < planetStates.size(); planetIndex++) {
+        this->planetOutMsgs.at(planetIndex)->write(
+            &planetStates.at(planetIndex), this->moduleID, CurrentSimNanos);
     }
     return;
 }
diff --git a/src/simulation/environment/planetEphemeris/planetEphemeris.h b/src/simulation/environment/planetEphemeris/planetEphemeris.h
@@ -20,6 +20,7 @@
 #ifndef planetEphemeris_H
 #define planetEphemeris_H
 
+#include <string>
 #include <vector>
 #include "architecture/_GeneralModuleFiles/sys_model.h"
 
@@ -52,12 +53,15 @@ class PlanetEphemeris: public SysModel {
 
     std::vector<double> rotRate;                //!< [r/s] planet rotation rate
 
+    std::string zeroBase{};                     //!< -- Output origin body; empty preserves heliocentric states
+
     BSKLogger bskLogger;                      //!< -- BSK Logging
 
 private:
     std::vector<std::string> planetNames;       //!< -- Vector of planet names
     double epochTime;                           //!< [s] time of provided planet ephemeris epoch
-    int computeAttitudeFlag;                    //!< -- flag indicating if the planet orienation information is provided
+    int computeAttitudeFlag;                    //!< -- flag indicating whether planet orientation is provided
+    int zeroBaseIndex = -1;                     //!< -- index of the body selected as the translational origin
 };
 
 
diff --git a/src/simulation/environment/planetEphemeris/planetEphemeris.rst b/src/simulation/environment/planetEphemeris/planetEphemeris.rst
