From 2338b5e3dcd36675719d3ea2c424cae8acdab11a Mon Sep 17 00:00:00 2001 From: Mason Drumwright Date: Mon, 5 Oct 2026 19:16:07 -0600 Subject: [PATCH 1/9] [1596] Enable LinearSpringMassDamper branching on prescribed motion --- .../linearSpringMassDamper.cpp | 87 +++++++++++++++++-- .../linearSpringMassDamper.h | 7 ++ 2 files changed, 89 insertions(+), 5 deletions(-) diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp index 0d3691b2f1c..e5728190b83 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp @@ -78,16 +78,33 @@ LinearSpringMassDamper::~LinearSpringMassDamper() */ void LinearSpringMassDamper::linkInStates(DynParamManager& states) { - // - Grab access to gravity - this->g_N = states.getPropertyReference(this->propName_vehicleGravity); - // - Grab access to c_B and cPrime_B - this->c_B = states.getPropertyReference(this->propName_centerOfMassSC); - this->cPrime_B = states.getPropertyReference(this->propName_centerOfMassPrimeSC); + // Get access to the hub's states needed for dynamic coupling + this->hubOmega = states.getStateObject("hubOmega"); + + const std::string gravName = this->propName_vehicleGravity.empty() ? "g_N" : this->propName_vehicleGravity; + this->g_N = states.getPropertyReference(gravName); + + if (!this->propName_centerOfMassSC.empty()) { + this->c_B = states.getPropertyReference(this->propName_centerOfMassSC); + } + if (!this->propName_centerOfMassPrimeSC.empty()) { + this->cPrime_B = states.getPropertyReference(this->propName_centerOfMassPrimeSC); + } return; } +void LinearSpringMassDamper::linkInPrescribedMotionProperties(DynParamManager& states) +{ + this->prescribedPositionProperty = states.getPropertyReference(this->propName_prescribedPosition); + this->prescribedVelocityProperty = states.getPropertyReference(this->propName_prescribedVelocity); + this->prescribedAccelerationProperty = states.getPropertyReference(this->propName_prescribedAcceleration); + this->prescribedAttitudeProperty = states.getPropertyReference(this->propName_prescribedAttitude); + this->prescribedAngVelocityProperty = states.getPropertyReference(this->propName_prescribedAngVelocity); + this->prescribedAngAccelerationProperty = states.getPropertyReference(this->propName_prescribedAngAcceleration); +} + /*! This is the method for the spring mass damper particle to register its states: rho and rhoDot * * @param[in,out] states Dynamic parameter manager used to register states or properties. @@ -115,6 +132,8 @@ void LinearSpringMassDamper::registerStates(DynParamManager& states) return; } +// Create method linkInPrescribedMotionProperties + /*! This is the method for the SMD to add its contributions to the mass props and mass prop rates of the vehicle * * @param[in] integTime [s] Current integration time. @@ -281,6 +300,64 @@ void LinearSpringMassDamper::computeDerivatives(double integTime [[maybe_unused] return; } +void LinearSpringMassDamper::addPrescribedMotionCouplingContributions(BackSubMatrices& backSubContr) +{ + // SMD frame is denoted with T and prescribed frame is denoted with P within the scope of this method + + // Access prescribed motion properties + Eigen::Vector3d r_PB_B = (Eigen::Vector3d) * this->prescribedPositionProperty; + Eigen::Vector3d rPrime_PB_B = (Eigen::Vector3d) * this->prescribedVelocityProperty; + Eigen::Vector3d rPrimePrime_PB_B = (Eigen::Vector3d) * this->prescribedAccelerationProperty; + Eigen::MRPd sigma_PB(this->prescribedAttitudeProperty->data()); + Eigen::Vector3d omega_PB_P = (Eigen::Vector3d) * this->prescribedAngVelocityProperty; + Eigen::Vector3d omegaPrime_PB_P = (Eigen::Vector3d) * this->prescribedAngAccelerationProperty; + Eigen::Matrix3d dcm_PB = sigma_PB.toRotationMatrix().transpose(); + + // Collect hub states + Eigen::Vector3d omega_BN_B = this->hubOmega->stateView(); + Eigen::Vector3d omega_BN_P = dcm_PB * omega_BN_B; + + // Prescribed motion coupling contributions + Eigen::Vector3d tHat_P = this->pHat_B; + Eigen::Vector3d r_PB_P = dcm_PB * r_PB_B; + Eigen::Matrix3d rTilde_PB_P = eigenTilde(r_PB_P); + backSubContr.matrixB += -this->massSMD * tHat_P * this->aRho.transpose() * rTilde_PB_P; + + Eigen::Matrix3d omegaTilde_PB_P = eigenTilde(omega_PB_P); + Eigen::Vector3d rPPrime_TP_P = this->rPrime_PcB_B; + Eigen::Matrix3d omegaPrimeTilde_PB_P = eigenTilde(omegaPrime_PB_P); + Eigen::Vector3d r_TP_P = this->r_PcB_B; + Eigen::Vector3d rPrimePrime_PB_P = dcm_PB * rPrimePrime_PB_B; + Eigen::Matrix3d omegaTilde_BN_P = eigenTilde(omega_BN_P); + Eigen::Vector3d rPrime_PB_P = dcm_PB * rPrime_PB_B; + Eigen::Vector3d term1 = 2.0 * omegaTilde_PB_P * rPPrime_TP_P + omegaPrimeTilde_PB_P * r_TP_P + + omegaTilde_PB_P * omegaTilde_PB_P * r_TP_P + rPrimePrime_PB_P; + Eigen::Vector3d term2 = + rPrimePrime_PB_P + 2.0 * omegaTilde_BN_P * rPrime_PB_P + omegaTilde_BN_P * omegaTilde_BN_P * r_PB_P; + Eigen::Vector3d term3 = omegaPrime_PB_P + omegaTilde_BN_P * omega_PB_P; + backSubContr.vecTrans += -this->massSMD * term1 - this->massSMD * this->aRho.transpose() * term2 * tHat_P - + this->massSMD * this->bRho.transpose() * term3 * tHat_P; + + // Prescribed motion rotation coupling contributions + backSubContr.matrixC += this->massSMD * rTilde_PB_P * tHat_P * this->aRho.transpose(); + + Eigen::Vector3d r_TB_P = r_TP_P + r_PB_P; + Eigen::Matrix3d rTilde_TB_P = eigenTilde(r_TB_P); + backSubContr.matrixD += +this->massSMD * rTilde_PB_P * tHat_P * this->bRho.transpose() - + this->massSMD * rTilde_TB_P * tHat_P * this->aRho.transpose() * rTilde_PB_P; + + Eigen::Matrix3d rTilde_TP_P = eigenTilde(r_TP_P); + + Eigen::Vector3d vecRotTerm1 = -this->massSMD * rTilde_TB_P * term1; + Eigen::Vector3d vecRotTerm2 = + -this->massSMD * (omegaTilde_BN_P * rTilde_PB_P - omegaTilde_PB_P * rTilde_TP_P) * rPPrime_TP_P - + this->massSMD * omegaTilde_BN_P * rTilde_TB_P * (omegaTilde_PB_P * r_TP_P + rPrime_PB_P); + Eigen::Vector3d vecRotTerm3 = -this->massSMD * this->cRho * rTilde_PB_P * tHat_P; + Eigen::Vector3d vecRotTerm4 = -this->massSMD * rTilde_TB_P * tHat_P * (this->aRho.transpose() * term2) - + this->massSMD * rTilde_TB_P * tHat_P * (this->bRho.transpose() * term3); + backSubContr.vecRot += vecRotTerm1 + vecRotTerm2 + vecRotTerm3 + vecRotTerm4; +} + /*! This method is for the SMD to add its contributions to energy and momentum * * @param[in] integTime [s] Current integration time. diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h index fb426aa453d..c0cf14cbd66 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h @@ -76,6 +76,8 @@ class LinearSpringMassDamper final void Reset(uint64_t CurrentSimNanos) override; void registerStates(DynParamManager& states) override; //!< Method for SMD to register its states void linkInStates(DynParamManager& states) override; //!< Method for SMD to get access of other states + void linkInPrescribedMotionProperties( + DynParamManager& states) override; //!< Method for SMD to get access to prescribed motion properties void retrieveMassValue(double integTime) override; void calcForceTorqueOnBody(double integTime, Eigen::Vector3d omega_BN_B) override; //!< Force and torque on s/c due to linear spring mass damper @@ -93,6 +95,11 @@ class LinearSpringMassDamper final Eigen::Vector3d rDDot_BN_N, Eigen::Vector3d omegaDot_BN_B, Eigen::MRPd sigma_BN) override; //!< Method for each stateEffector to calculate derivatives + void addPrescribedMotionCouplingContributions(BackSubMatrices& backSubContr) + override; //!< Method for adding coupling contributions for state effector branching on prescribed motion + + // Properties required for prescribed motion branching/attachment + StateData* hubOmega; //!< [rad/s] hub inertial angular velocity vector }; #endif /* LINEAR_SPRING_MASS_DAMPER_H */ From e1b520f4445623e2ab0bd50960a3281b45abb862 Mon Sep 17 00:00:00 2001 From: Mason Drumwright Date: Sun, 8 Nov 2026 11:54:16 -0700 Subject: [PATCH 2/9] [1596] Add LinearSpringMassDamper unit test --- ..._PrescribedMotionStateEffectorBranching.py | 119 ++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py index c132aaa20b2..81e26bc811c 100644 --- a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py +++ b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py @@ -35,6 +35,7 @@ from Basilisk.simulation import prescribedLinearTranslation from Basilisk.simulation import prescribedRotation1DOF from Basilisk.simulation import linearTranslationOneDOFStateEffector +from Basilisk.simulation import linearSpringMassDamper from Basilisk.simulation import spinningBodyOneDOFStateEffector from Basilisk.simulation import spinningBodyTwoDOFStateEffector from Basilisk.architecture import messaging @@ -609,6 +610,123 @@ def test_prescribed_branching_linear_translation_one_dof(show_plots): unit_test_check_verification(orb_ang_mom_n, orb_energy, rot_ang_mom_n) +def test_prescribed_branching_linear_spring_mass_damper(show_plots): + r""" + **Verification Test Description** + + This unit test checks the branching capability of the prescribed motion module by connecting a + :ref:`linearSpringMassDamper` to the prescribed motion state effector. Both the + :ref:`prescribedRotation1DOF` and :ref:`prescribedLinearTranslation` modules are configured so that + the prescribed body simultaneously translates and rotates relative to the spacecraft hub. The particle + is set up to have free (undamped) spring motion along its direction of translation. + + **Description of Variables Being Tested** + + The test checks that the conservation quantities of the spacecraft orbital energy, orbital angular momentum, + and rotational angular momentum remain constant over the duration of the simulation. + """ + + task_name = "unitTask" + process_name = "testProcess" + sim = SimulationBaseClass.SimBaseClass() + process_rate = macros.sec2nano(0.001) + process = sim.CreateNewProcess(process_name) + process.addTask(sim.CreateNewTask(task_name, process_rate)) + + # Create the spacecraft + sc_object = create_spacecraft_hub() + sim.AddModelToTask(task_name, sc_object) + + # Create prescribed motion object + prescribed_platform = create_prescribed_body(r_P0M_M, sigma_P0M) + sim.AddModelToTask(task_name, prescribed_platform) + sc_object.addStateEffector(prescribed_platform) + + # Create rotational motion profiler + one_dof_rotation_profiler = prescribedRotation1DOF.PrescribedRotation1DOF() + one_dof_rotation_profiler.ModelTag = "prescribedRotation1DOF" + one_dof_rotation_profiler.setRotHat_M(prescribed_rot_axis_M) + one_dof_rotation_profiler.setThetaDDotMax(prescribed_ang_accel_max) + one_dof_rotation_profiler.setThetaInit(prescribed_theta_init) + one_dof_rotation_profiler.setCoastOptionBangDuration(1.0) + one_dof_rotation_profiler.setSmoothingDuration(1.0) + sim.AddModelToTask(task_name, one_dof_rotation_profiler) + + prescribed_rotation_msg_data = messaging.HingedRigidBodyMsgPayload() + prescribed_rotation_msg_data.theta = prescribed_theta_ref + prescribed_rotation_msg_data.thetaDot = 0.0 + prescribed_rotation_msg = messaging.HingedRigidBodyMsg().write(prescribed_rotation_msg_data) + one_dof_rotation_profiler.spinningBodyInMsg.subscribeTo(prescribed_rotation_msg) + prescribed_platform.prescribedRotationInMsg.subscribeTo(one_dof_rotation_profiler.prescribedRotationOutMsg) + + # Create translational motion profiler + one_dof_translation_profiler = prescribedLinearTranslation.PrescribedLinearTranslation() + one_dof_translation_profiler.ModelTag = "prescribedLinearTranslation" + one_dof_translation_profiler.setTransHat_M(prescribed_trans_axis_M) + one_dof_translation_profiler.setTransAccelMax(prescribed_trans_accel_max) + one_dof_translation_profiler.setTransPosInit(prescribed_pos_init) + one_dof_translation_profiler.setCoastOptionBangDuration(1.0) + one_dof_translation_profiler.setSmoothingDuration(1.0) + sim.AddModelToTask(task_name, one_dof_translation_profiler) + + prescribed_translation_msg_data = messaging.LinearTranslationRigidBodyMsgPayload() + prescribed_translation_msg_data.rho = prescribed_pos_ref + prescribed_translation_msg_data.rhoDot = 0.0 + prescribed_translation_msg = messaging.LinearTranslationRigidBodyMsg().write(prescribed_translation_msg_data) + one_dof_translation_profiler.linearTranslationRigidBodyInMsg.subscribeTo(prescribed_translation_msg) + prescribed_platform.prescribedTranslationInMsg.subscribeTo(one_dof_translation_profiler.prescribedTranslationOutMsg) + + # Create the spring-mass-damper particle + smd = linearSpringMassDamper.LinearSpringMassDamper() + smd.ModelTag = "springMassDamper" + smd.massInit = 20.0 # [kg] + smd.k = 100.0 # [N/m] + smd.c = 0.0 # [N s/m] undamped + smd.rhoInit = 0.5 # [m] + smd.rhoDotInit = 0.0 # [m/s] + smd.pHat_B = [[1.0], [0.0], [0.0]] # Translation axis (x) + smd.r_PB_B = [[1.0], [0.0], [0.0]] # [m] + sim.AddModelToTask(task_name, smd) + + # Attach the particle to the prescribed motion platform + prescribed_platform.addStateEffector(smd) + + # Add Earth gravity to the simulation + earth_gravity = gravityEffector.GravBodyData() + earth_gravity.planetName = "earth_planet_data" + earth_gravity.mu = 0.3986004415E+15 # meters + earth_gravity.isCentralBody = True + sc_object.gravField.gravBodies = spacecraft.GravBodyVector([earth_gravity]) + + # Set up data logging + conservation_data_log = sc_object.logger(["totOrbAngMomPntN_N", "totRotAngMomPntC_N", "totOrbEnergy", "totRotEnergy"]) + prescribed_rotation_data_log = prescribed_platform.prescribedRotationOutMsg.recorder() + sim.AddModelToTask(task_name, conservation_data_log) + sim.AddModelToTask(task_name, prescribed_rotation_data_log) + + # Run the simulation + sim.InitializeSimulation() + sim.ConfigureStopTime(macros.sec2nano(sim_time)) + sim.ExecuteSimulation() + + # Extract the logged variables + orb_energy = conservation_data_log.totOrbEnergy + orb_ang_mom_n = conservation_data_log.totOrbAngMomPntN_N + rot_ang_mom_n = conservation_data_log.totRotAngMomPntC_N + rot_energy = conservation_data_log.totRotEnergy + timespan = prescribed_rotation_data_log.times() * macros.NANO2SEC # [s] + + # Plot results + plot_conservation(timespan, orb_ang_mom_n, orb_energy, rot_ang_mom_n, rot_energy) + + if show_plots: + plt.show() + plt.close("all") + + # Test check + unit_test_check_verification(orb_ang_mom_n, orb_energy, rot_ang_mom_n) + + def create_spacecraft_hub(): mass_hub = 800 # [kg] length_hub = 1.0 # [m] @@ -788,3 +906,4 @@ def unit_test_check_verification(orb_ang_mom_n, orb_energy, rot_ang_mom_n): test_prescribed_branching_spinning_body_one_dof(True) test_prescribed_branching_spinning_body_two_dof(True) test_prescribed_branching_linear_translation_one_dof(True) + test_prescribed_branching_linear_spring_mass_damper(True) From 9085d8b6e677b35b3aed9e7bd8cee0bcb5a00dee Mon Sep 17 00:00:00 2001 From: Leah Kiner Date: Tue, 27 Oct 2026 12:13:30 -0600 Subject: [PATCH 3/9] [1596] Remove private gravity and initialize pointers to nullptr --- .../linearSpringMassDamper.cpp | 12 ++----- .../linearSpringMassDamper.h | 33 +++++++++---------- 2 files changed, 19 insertions(+), 26 deletions(-) diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp index e5728190b83..bf3936025a9 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp @@ -72,7 +72,7 @@ LinearSpringMassDamper::~LinearSpringMassDamper() return; } -/*! Method for spring mass damper particle to access the states that it needs. It needs gravity and the hub states +/*! Method for spring mass damper particle to access the states that it needs. It needs the hub states * * @param[in] states Dynamic parameter manager containing the required states. */ @@ -82,9 +82,6 @@ void LinearSpringMassDamper::linkInStates(DynParamManager& states) // Get access to the hub's states needed for dynamic coupling this->hubOmega = states.getStateObject("hubOmega"); - const std::string gravName = this->propName_vehicleGravity.empty() ? "g_N" : this->propName_vehicleGravity; - this->g_N = states.getPropertyReference(gravName); - if (!this->propName_centerOfMassSC.empty()) { this->c_B = states.getPropertyReference(this->propName_centerOfMassSC); } @@ -209,7 +206,7 @@ void LinearSpringMassDamper::updateContributions(double integTime [[maybe_unused BackSubMatrices& backSubContr, Eigen::MRPd sigma_BN, Eigen::Vector3d omega_BN_B, - Eigen::Vector3d g_N [[maybe_unused]]) + Eigen::Vector3d g_N) { if (this->massSMD <= 0.0) { this->aRho.setZero(); @@ -233,10 +230,7 @@ void LinearSpringMassDamper::updateContributions(double integTime [[maybe_unused dcm_BN = dcm_NB.transpose(); // - Map gravity to body frame - Eigen::Vector3d gLocal_N; - Eigen::Vector3d g_B; - gLocal_N = *this->g_N; - g_B = dcm_BN * gLocal_N; + Eigen::Vector3d g_B = dcm_BN * g_N; // - Define aRho this->aRho = -this->pHat_B; diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h index c0cf14cbd66..d20f749c475 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h @@ -52,23 +52,22 @@ class LinearSpringMassDamper final BSKLogger bskLogger; //!< BSK Logging private: - void validateConfiguration(); //!< Validate the initial particle mass without changing its current state - double cRho; //!< Term needed for back-sub method - double rho; //!< [m] spring mass damper displacement from equilibrium - double rhoDot; //!< [m/s] time derivative of displacement from equilibrium - double massSMD; //!< [kg] mass of spring mass damper particle - Eigen::Vector3d r_PcB_B; //!< [m] position vector form B to center of mass location of particle - Eigen::Matrix3d rTilde_PcB_B; //!< [m] tilde matrix of r_Pc_B - Eigen::Vector3d rPrime_PcB_B; //!< [m/s] Body time derivative of r_Pc_B - Eigen::Matrix3d rPrimeTilde_PcB_B; //!< [m/s] Tilde matrix of rPrime_PcB_B - Eigen::Vector3d aRho; //!< Term needed for back-sub method - Eigen::Vector3d bRho; //!< Term needed for back-sub method - Eigen::MatrixXd* g_N; //!< [m/s^2] Gravitational acceleration in N frame components - StateData* rhoState; //!< state data for spring mass damper displacement from equilibrium - Eigen::MatrixXd* c_B; //!< [m] Vector from point B to CoM of s/c in B frame components - Eigen::MatrixXd* cPrime_B; //!< [m/s] Body time derivative of vector c_B in B frame components - StateData* rhoDotState; //!< state data for time derivative of rho; - static uint64_t effectorID; //!< [] ID number of this panel + void validateConfiguration(); //!< Validate the initial particle mass without changing its current state + double cRho; //!< Term needed for back-sub method + double rho; //!< [m] spring mass damper displacement from equilibrium + double rhoDot; //!< [m/s] time derivative of displacement from equilibrium + double massSMD; //!< [kg] mass of spring mass damper particle + Eigen::Vector3d r_PcB_B; //!< [m] position vector form B to center of mass location of particle + Eigen::Matrix3d rTilde_PcB_B; //!< [m] tilde matrix of r_Pc_B + Eigen::Vector3d rPrime_PcB_B; //!< [m/s] Body time derivative of r_Pc_B + Eigen::Matrix3d rPrimeTilde_PcB_B; //!< [m/s] Tilde matrix of rPrime_PcB_B + Eigen::Vector3d aRho; //!< Term needed for back-sub method + Eigen::Vector3d bRho; //!< Term needed for back-sub method + StateData* rhoState; //!< state data for spring mass damper displacement from equilibrium + Eigen::MatrixXd* c_B = nullptr; //!< [m] Vector from point B to CoM of s/c in B frame components + Eigen::MatrixXd* cPrime_B = nullptr; //!< [m/s] Body time derivative of vector c_B in B frame components + StateData* rhoDotState; //!< state data for time derivative of rho; + static uint64_t effectorID; //!< [] ID number of this panel public: LinearSpringMassDamper(); //!< Constructor From b46e275307e95b49af86ebee8b2829bfe53592b2 Mon Sep 17 00:00:00 2001 From: Mason Drumwright Date: Sun, 8 Nov 2026 12:18:58 -0700 Subject: [PATCH 4/9] [1596] Add release note snippet --- .../1596-configure-SMD-for-prescribed-motion-branching.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/source/Support/bskReleaseNotesSnippets/1596-configure-SMD-for-prescribed-motion-branching.rst diff --git a/docs/source/Support/bskReleaseNotesSnippets/1596-configure-SMD-for-prescribed-motion-branching.rst b/docs/source/Support/bskReleaseNotesSnippets/1596-configure-SMD-for-prescribed-motion-branching.rst new file mode 100644 index 00000000000..a4b481ea12b --- /dev/null +++ b/docs/source/Support/bskReleaseNotesSnippets/1596-configure-SMD-for-prescribed-motion-branching.rst @@ -0,0 +1 @@ +- Enabled :ref:`linearSpringMassDamper` to be attached to a prescribed motion state effector. From 39005b18887a7d66837492145b36e337c5241f39 Mon Sep 17 00:00:00 2001 From: Hanspeter Schaub Date: Tue, 10 Nov 2026 11:38:20 -0700 Subject: [PATCH 5/9] [#1596] Remove hard-coded hub state dependency from spring branching --- .../linearSpringMassDamper.cpp | 13 ++++++------- .../LinearSpringMassDamper/linearSpringMassDamper.h | 4 +--- .../test_PrescribedMotionStateEffectorBranching.py | 10 ++++++++-- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp index bf3936025a9..5b8649b1f6b 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp @@ -72,16 +72,13 @@ LinearSpringMassDamper::~LinearSpringMassDamper() return; } -/*! Method for spring mass damper particle to access the states that it needs. It needs the hub states +/*! @brief Link the spacecraft center-of-mass properties when attached directly to the hub. * * @param[in] states Dynamic parameter manager containing the required states. */ void LinearSpringMassDamper::linkInStates(DynParamManager& states) { - // Get access to the hub's states needed for dynamic coupling - this->hubOmega = states.getStateObject("hubOmega"); - if (!this->propName_centerOfMassSC.empty()) { this->c_B = states.getPropertyReference(this->propName_centerOfMassSC); } @@ -208,6 +205,9 @@ void LinearSpringMassDamper::updateContributions(double integTime [[maybe_unused Eigen::Vector3d omega_BN_B, Eigen::Vector3d g_N) { + // The callback supplies the attached parent's inertial rate, in that parent's frame. + this->omega_PN_P = omega_BN_B; + if (this->massSMD <= 0.0) { this->aRho.setZero(); this->bRho.setZero(); @@ -307,9 +307,8 @@ void LinearSpringMassDamper::addPrescribedMotionCouplingContributions(BackSubMat Eigen::Vector3d omegaPrime_PB_P = (Eigen::Vector3d) * this->prescribedAngAccelerationProperty; Eigen::Matrix3d dcm_PB = sigma_PB.toRotationMatrix().transpose(); - // Collect hub states - Eigen::Vector3d omega_BN_B = this->hubOmega->stateView(); - Eigen::Vector3d omega_BN_P = dcm_PB * omega_BN_B; + // Recover the hub rate from the parent rate supplied to updateContributions(). + Eigen::Vector3d omega_BN_P = this->omega_PN_P - omega_PB_P; // Prescribed motion coupling contributions Eigen::Vector3d tHat_P = this->pHat_B; diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h index d20f749c475..c0951f247b7 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h @@ -63,6 +63,7 @@ class LinearSpringMassDamper final Eigen::Matrix3d rPrimeTilde_PcB_B; //!< [m/s] Tilde matrix of rPrime_PcB_B Eigen::Vector3d aRho; //!< Term needed for back-sub method Eigen::Vector3d bRho; //!< Term needed for back-sub method + Eigen::Vector3d omega_PN_P{}; //!< [rad/s] Parent inertial angular velocity in parent components StateData* rhoState; //!< state data for spring mass damper displacement from equilibrium Eigen::MatrixXd* c_B = nullptr; //!< [m] Vector from point B to CoM of s/c in B frame components Eigen::MatrixXd* cPrime_B = nullptr; //!< [m/s] Body time derivative of vector c_B in B frame components @@ -96,9 +97,6 @@ class LinearSpringMassDamper final Eigen::MRPd sigma_BN) override; //!< Method for each stateEffector to calculate derivatives void addPrescribedMotionCouplingContributions(BackSubMatrices& backSubContr) override; //!< Method for adding coupling contributions for state effector branching on prescribed motion - - // Properties required for prescribed motion branching/attachment - StateData* hubOmega; //!< [rad/s] hub inertial angular velocity vector }; #endif /* LINEAR_SPRING_MASS_DAMPER_H */ diff --git a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py index 81e26bc811c..7e76676c708 100644 --- a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py +++ b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py @@ -610,7 +610,8 @@ def test_prescribed_branching_linear_translation_one_dof(show_plots): unit_test_check_verification(orb_ang_mom_n, orb_energy, rot_ang_mom_n) -def test_prescribed_branching_linear_spring_mass_damper(show_plots): +@pytest.mark.parametrize("hub_omega_name", ["hubOmega", "customHubOmega"]) +def test_prescribed_branching_linear_spring_mass_damper(show_plots, hub_omega_name): r""" **Verification Test Description** @@ -620,6 +621,10 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots): the prescribed body simultaneously translates and rotates relative to the spacecraft hub. The particle is set up to have free (undamped) spring motion along its direction of translation. + **Test Parameter Discussion** + + ``hub_omega_name`` exercises both the default and a custom hub angular-velocity state name. + **Description of Variables Being Tested** The test checks that the conservation quantities of the spacecraft orbital energy, orbital angular momentum, @@ -635,6 +640,7 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots): # Create the spacecraft sc_object = create_spacecraft_hub() + sc_object.hub.nameOfHubOmega = hub_omega_name sim.AddModelToTask(task_name, sc_object) # Create prescribed motion object @@ -906,4 +912,4 @@ def unit_test_check_verification(orb_ang_mom_n, orb_energy, rot_ang_mom_n): test_prescribed_branching_spinning_body_one_dof(True) test_prescribed_branching_spinning_body_two_dof(True) test_prescribed_branching_linear_translation_one_dof(True) - test_prescribed_branching_linear_spring_mass_damper(True) + test_prescribed_branching_linear_spring_mass_damper(True, "customHubOmega") From 1ddd819e7fe65dabb674a7171686db5b7419b561 Mon Sep 17 00:00:00 2001 From: Hanspeter Schaub Date: Tue, 10 Nov 2026 11:41:54 -0700 Subject: [PATCH 6/9] [#1596] Verify spring branching against independent dynamics references --- .../test_linearSpringMassDamperPrescribed.py | 263 ++++++++++++++++++ ..._PrescribedMotionStateEffectorBranching.py | 46 ++- 2 files changed, 297 insertions(+), 12 deletions(-) create mode 100644 src/simulation/dynamics/LinearSpringMassDamper/_UnitTest/test_linearSpringMassDamperPrescribed.py diff --git a/src/simulation/dynamics/LinearSpringMassDamper/_UnitTest/test_linearSpringMassDamperPrescribed.py b/src/simulation/dynamics/LinearSpringMassDamper/_UnitTest/test_linearSpringMassDamperPrescribed.py new file mode 100644 index 00000000000..216a21fd540 --- /dev/null +++ b/src/simulation/dynamics/LinearSpringMassDamper/_UnitTest/test_linearSpringMassDamperPrescribed.py @@ -0,0 +1,263 @@ +# 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. + +import numpy as np +import pytest + +from Basilisk.architecture import messaging +from Basilisk.simulation import linearSpringMassDamper, prescribedMotionStateEffector, spacecraft +from Basilisk.utilities import SimulationBaseClass, macros, pythonVariableLogger +from Basilisk.utilities import RigidBodyKinematics as rbk + + +def create_simulation(): + """Create a simulation with a 1 ms dynamics task.""" + sim = SimulationBaseClass.SimBaseClass() + process = sim.CreateNewProcess("process") + process.addTask(sim.CreateNewTask("task", macros.sec2nano(0.001))) # [s] + return sim + + +@pytest.mark.parametrize("damping", [0.0, 4.0]) # [N s/m] +def test_prescribed_spring_analytical_motion(damping): + r"""Compare the particle and hub motion with the two-mass oscillator solution. + + **Validation Test Description** + + A stationary prescribed body is centered on the hub. The spring passes through their + common center of mass, so the hub and prescribed body recoil without rotating. The + independent oscillator solution uses the reduced mass of the particle and its parent + assembly. A nonzero mount attitude checks the prescribed-frame interpretation of the axis. + + **Test Parameter Discussion** + + ``damping`` selects an undamped and an underdamped oscillator. There is no gravity or + external force. Both default and renamed hub states are covered by the integration test + in :ref:`test_PrescribedMotionStateEffectorBranching`. + + **Description of Variables Being Tested** + + Compare every recorded displacement and velocity sample to the analytical solution, and + check the hub recoil from conservation of the system center of mass. + """ + sim = create_simulation() + hub = spacecraft.Spacecraft() + hub.hub.mHub = 800.0 # [kg] + hub.hub.IHubPntBc_B = np.diag([150.0, 180.0, 210.0]).tolist() # [kg m^2] + parent = prescribedMotionStateEffector.PrescribedMotionStateEffector() + parent.setMass(100.0) # [kg] + parent.setIPntPc_P(np.diag([20.0, 25.0, 30.0]).tolist()) # [kg m^2] + sigma_MB = np.array([0.1, -0.2, 0.15]) # [-] MRPs + parent.setSigma_MB(sigma_MB.tolist()) + particle = linearSpringMassDamper.LinearSpringMassDamper() + particle.massInit = 20.0 # [kg] + particle.k = 100.0 # [N/m] + particle.c = damping + particle.rhoInit = 0.5 # [m] + particle.rhoDotInit = 0.15 # [m/s] + axis_P = np.array([0.0, 1.0, 0.0]) # [-] + particle.pHat_B = axis_P.tolist() + parent.addStateEffector(particle) + hub.addStateEffector(parent) + sim.AddModelToTask("task", hub) + sim.AddModelToTask("task", parent) + particle_log = pythonVariableLogger.PythonVariableLogger({ + "rho": lambda _: rho_state.getState(), + "rhoDot": lambda _: rho_dot_state.getState(), + }) + hub_log = hub.scStateOutMsg.recorder() + sim.AddModelToTask("task", particle_log) + sim.AddModelToTask("task", hub_log) + sim.InitializeSimulation() + rho_state = hub.dynManager.getStateObject(particle.nameOfRhoState) + rho_dot_state = hub.dynManager.getStateObject(particle.nameOfRhoDotState) + assert rho_state is not None and rho_dot_state is not None + duration = 4.0 # [s] + sim.ConfigureStopTime(macros.sec2nano(duration)) + sim.ExecuteSimulation() + + parent_mass = hub.hub.mHub + parent.getMass() + reduced_mass = 1 / (1 / particle.massInit + 1 / parent_mass) + decay_rate = damping / (2 * reduced_mass) + damped_frequency = np.sqrt(particle.k / reduced_mass - decay_rate**2) + times = particle_log.times() * macros.NANO2SEC + cosine = np.cos(damped_frequency * times) + sine = np.sin(damped_frequency * times) + sine_amplitude = (particle.rhoDotInit + decay_rate * particle.rhoInit) / damped_frequency + expected_rho = np.exp(-decay_rate * times) * (particle.rhoInit * cosine + sine_amplitude * sine) + expected_rho_dot = np.exp(-decay_rate * times) * ( + (-decay_rate * particle.rhoInit + damped_frequency * sine_amplitude) * cosine + - (decay_rate * sine_amplitude + damped_frequency * particle.rhoInit) * sine + ) + position_tolerance = 1e-10 # [m] + velocity_tolerance = 1e-10 # [m/s] + np.testing.assert_allclose(particle_log.rho, expected_rho, rtol=0, atol=position_tolerance) + np.testing.assert_allclose(particle_log.rhoDot, expected_rho_dot, rtol=0, atol=velocity_tolerance) + axis_N = rbk.MRP2C(sigma_MB).T @ axis_P + expected_recoil = -particle.massInit / (parent_mass + particle.massInit) * np.outer( + expected_rho - particle.rhoInit, axis_N + ) + np.testing.assert_allclose( + hub_log.r_BN_N - hub_log.r_BN_N[0], expected_recoil, rtol=0, atol=position_tolerance + ) + + +@pytest.mark.parametrize("seed", [1596, 1634, 2026]) +def test_prescribed_spring_newton_euler_balance(seed): + r"""Check the coupled accelerations against independent Newton/Euler equations. + + **Validation Test Description** + + Use non-collinear offsets, nonidentity mount and parent attitudes, simultaneous prescribed + translation and rotation, and nonzero spring velocity and damping. Reconstruct the inertial + accelerations from rigid-body kinematics and the integrated state derivatives. In the absence + of external forces, total force and torque must vanish, and the particle acceleration along + the spring must satisfy its scalar force equation. + + **Test Parameter Discussion** + + ``seed`` chooses three reproducible arbitrary geometries and prescribed kinematic states. + The hub angular-velocity state is renamed to exercise callback-based rate access. + + **Description of Variables Being Tested** + + Check hub translational/angular acceleration and particle acceleration through total force, + total torque, and axial spring-force residuals, without using the backsubstitution formulas. + """ + rng = np.random.default_rng(seed) + sim = create_simulation() + hub = spacecraft.Spacecraft() + hub.hub.nameOfHubOmega = "customHubOmega" + hub.hub.mHub = 800.0 # [kg] + inertia_hub_B = np.diag([150.0, 180.0, 210.0]) # [kg m^2] + hub.hub.IHubPntBc_B = inertia_hub_B.tolist() + omega_BN_B = rng.normal(size=3) * 0.2 # [rad/s] + hub.hub.omega_BN_BInit = omega_BN_B.tolist() + parent = prescribedMotionStateEffector.PrescribedMotionStateEffector() + parent.setMass(100.0) # [kg] + inertia_parent_P = np.diag([20.0, 25.0, 30.0]) # [kg m^2] + parent.setIPntPc_P(inertia_parent_P.tolist()) + r_PM_M = rng.normal(size=3) # [m] + r_MB_B = rng.normal(size=3) # [m] + r_prime_PM_M = rng.normal(size=3) * 0.2 # [m/s] + r_prime_prime_PM_M = rng.normal(size=3) * 0.3 # [m/s^2] + omega_PM_P = rng.normal(size=3) * 0.3 # [rad/s] + omega_prime_PM_P = rng.normal(size=3) * 0.2 # [rad/s^2] + sigma_PM = rng.normal(size=3) * 0.1 # [-] MRPs + sigma_MB = rng.normal(size=3) * 0.1 # [-] MRPs + r_PcP_P = rng.normal(size=3) * 0.3 # [m] + parent.setR_MB_B(r_MB_B.tolist()) + parent.setSigma_MB(sigma_MB.tolist()) + parent.setSigma_PM(sigma_PM.tolist()) + parent.setR_PcP_P(r_PcP_P.tolist()) + translation = messaging.PrescribedTranslationMsgPayload() + translation.r_PM_M = r_PM_M.tolist() + translation.rPrime_PM_M = r_prime_PM_M.tolist() + translation.rPrimePrime_PM_M = r_prime_prime_PM_M.tolist() + translation_msg = messaging.PrescribedTranslationMsg().write(translation) + parent.prescribedTranslationInMsg.subscribeTo(translation_msg) + rotation = messaging.PrescribedRotationMsgPayload() + rotation.sigma_PM = sigma_PM.tolist() + rotation.omega_PM_P = omega_PM_P.tolist() + rotation.omegaPrime_PM_P = omega_prime_PM_P.tolist() + rotation_msg = messaging.PrescribedRotationMsg().write(rotation) + parent.prescribedRotationInMsg.subscribeTo(rotation_msg) + particle = linearSpringMassDamper.LinearSpringMassDamper() + particle.massInit = 20.0 # [kg] + particle.k = 100.0 # [N/m] + particle.c = 2.0 # [N s/m] + axis_P = rng.normal(size=3) + axis_P /= np.linalg.norm(axis_P) + equilibrium_P = rng.normal(size=3) * 0.5 # [m] + particle.pHat_B = axis_P.tolist() + particle.r_PB_B = equilibrium_P.tolist() + particle.rhoInit = 0.3 # [m] + particle.rhoDotInit = 0.1 # [m/s] + parent.addStateEffector(particle) + hub.addStateEffector(parent) + sim.AddModelToTask("task", hub) + sim.AddModelToTask("task", parent) + sim.InitializeSimulation() + parent.UpdateState(0) + hub.equationsOfMotion(0.0, 0.001) # [s] Initial time and integration step + + def state_derivative(name): + """Return one registered state derivative as a flat array.""" + return np.array(hub.dynManager.getStateObject(name).getStateDeriv()).reshape(-1) + + # The initial hub attitude is identity, so inertial and hub components agree. + acceleration_BN_B = state_derivative(hub.hub.nameOfHubVelocity) + alpha_BN_B = state_derivative(hub.hub.nameOfHubOmega) + rho_ddot = state_derivative(particle.nameOfRhoDotState)[0] + dcm_BM = rbk.MRP2C(sigma_MB).T + dcm_BP = dcm_BM @ rbk.MRP2C(sigma_PM).T + omega_PB_B = dcm_BP @ omega_PM_P + omega_prime_PB_B = dcm_BP @ omega_prime_PM_P + r_PB_B = r_MB_B + dcm_BM @ r_PM_M + velocity_PB_B = dcm_BM @ r_prime_PM_M + acceleration_PB_B = dcm_BM @ r_prime_prime_PM_M + + def point_kinematics(position_P, velocity_P, acceleration_P): + """Return hub-relative position and inertial acceleration, both in hub components.""" + arm_B = dcm_BP @ position_P + arm_velocity_B = dcm_BP @ velocity_P + position_B = r_PB_B + arm_B + velocity_B = velocity_PB_B + np.cross(omega_PB_B, arm_B) + arm_velocity_B + relative_acceleration_B = ( + acceleration_PB_B + np.cross(omega_prime_PB_B, arm_B) + + np.cross(omega_PB_B, np.cross(omega_PB_B, arm_B)) + + 2 * np.cross(omega_PB_B, arm_velocity_B) + dcm_BP @ acceleration_P + ) + inertial_acceleration_B = ( + acceleration_BN_B + np.cross(alpha_BN_B, position_B) + + np.cross(omega_BN_B, np.cross(omega_BN_B, position_B)) + + 2 * np.cross(omega_BN_B, velocity_B) + relative_acceleration_B + ) + return position_B, inertial_acceleration_B + + r_parent_B, acceleration_parent_B = point_kinematics(r_PcP_P, np.zeros(3), np.zeros(3)) + r_particle_B, acceleration_particle_B = point_kinematics( + equilibrium_P + particle.rhoInit * axis_P, particle.rhoDotInit * axis_P, rho_ddot * axis_P + ) + force_residual = ( + hub.hub.mHub * acceleration_BN_B + parent.getMass() * acceleration_parent_B + + particle.massInit * acceleration_particle_B + ) + inertia_parent_B = dcm_BP @ inertia_parent_P @ dcm_BP.T + omega_PN_B = omega_BN_B + omega_PB_B + alpha_PN_B = alpha_BN_B + omega_prime_PB_B + np.cross(omega_BN_B, omega_PB_B) + torque_residual = ( + inertia_hub_B @ alpha_BN_B + np.cross(omega_BN_B, inertia_hub_B @ omega_BN_B) + + inertia_parent_B @ alpha_PN_B + np.cross(omega_PN_B, inertia_parent_B @ omega_PN_B) + + np.cross(r_parent_B, parent.getMass() * acceleration_parent_B) + + np.cross(r_particle_B, particle.massInit * acceleration_particle_B) + ) + spring_residual = ( + particle.massInit * np.dot(dcm_BP @ axis_P, acceleration_particle_B) + + particle.k * particle.rhoInit + particle.c * particle.rhoDotInit + ) + force_tolerance = 1e-10 # [N] + torque_tolerance = 1e-10 # [N m] + np.testing.assert_allclose(force_residual, np.zeros(3), rtol=0, atol=force_tolerance) + np.testing.assert_allclose(torque_residual, np.zeros(3), rtol=0, atol=torque_tolerance) + np.testing.assert_allclose(spring_residual, 0, rtol=0, atol=force_tolerance) + + +if __name__ == "__main__": + for damping in (0.0, 4.0): # [N s/m] + test_prescribed_spring_analytical_motion(damping) + for seed in (1596, 1634, 2026): + test_prescribed_spring_newton_euler_balance(seed) diff --git a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py index 7e76676c708..bc498d1083d 100644 --- a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py +++ b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py @@ -615,26 +615,38 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots, hub_omega_na r""" **Verification Test Description** - This unit test checks the branching capability of the prescribed motion module by connecting a - :ref:`linearSpringMassDamper` to the prescribed motion state effector. Both the - :ref:`prescribedRotation1DOF` and :ref:`prescribedLinearTranslation` modules are configured so that - the prescribed body simultaneously translates and rotates relative to the spacecraft hub. The particle - is set up to have free (undamped) spring motion along its direction of translation. + Connect a :ref:`linearSpringMassDamper` to a prescribed body undergoing simultaneous + translation and rotation. The spring axis is transverse to the rotation axis, its equilibrium + point is offset from that axis, and the mount attitude is nonzero, exercising rotational + coupling. Require both particle states to be registered and check conservation at every sample. **Test Parameter Discussion** ``hub_omega_name`` exercises both the default and a custom hub angular-velocity state name. + Run at 1 ms and 0.5 ms. Sampled profiler commands introduce first-order momentum errors at + the task boundaries: bound the relative rotational angular-momentum error by ``1e-4`` and + require it to decrease by at least 40 percent when halving the step. Orbital energy and + angular momentum must remain within ``1e-8`` throughout both runs. **Description of Variables Being Tested** - The test checks that the conservation quantities of the spacecraft orbital energy, orbital angular momentum, - and rotational angular momentum remain constant over the duration of the simulation. + Check orbital energy, orbital angular momentum, and rotational angular momentum. Rotational + energy need not be constant because the prescribed actuators do work. Independent spring + displacement/velocity and Newton/Euler checks are in + :ref:`test_linearSpringMassDamperPrescribed`. """ + coarse_step = 0.001 # [s] + coarse_error = run_prescribed_spring_branch(show_plots, hub_omega_name, coarse_step) + fine_error = run_prescribed_spring_branch(False, hub_omega_name, coarse_step / 2) + assert fine_error < 0.6 * coarse_error + 1e-10 + +def run_prescribed_spring_branch(show_plots, hub_omega_name, time_step): + """Run the profiler-driven branching case and return its maximum relative momentum error.""" task_name = "unitTask" process_name = "testProcess" sim = SimulationBaseClass.SimBaseClass() - process_rate = macros.sec2nano(0.001) + process_rate = macros.sec2nano(time_step) process = sim.CreateNewProcess(process_name) process.addTask(sim.CreateNewTask(task_name, process_rate)) @@ -645,6 +657,7 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots, hub_omega_na # Create prescribed motion object prescribed_platform = create_prescribed_body(r_P0M_M, sigma_P0M) + prescribed_platform.setSigma_MB([0.1, -0.2, 0.15]) # [-] MRPs sim.AddModelToTask(task_name, prescribed_platform) sc_object.addStateEffector(prescribed_platform) @@ -690,8 +703,8 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots, hub_omega_na smd.c = 0.0 # [N s/m] undamped smd.rhoInit = 0.5 # [m] smd.rhoDotInit = 0.0 # [m/s] - smd.pHat_B = [[1.0], [0.0], [0.0]] # Translation axis (x) - smd.r_PB_B = [[1.0], [0.0], [0.0]] # [m] + smd.pHat_B = [[0.0], [1.0], [0.0]] # [-] Axis transverse to prescribed rotation + smd.r_PB_B = [[0.4], [-0.2], [0.6]] # [m] sim.AddModelToTask(task_name, smd) # Attach the particle to the prescribed motion platform @@ -712,6 +725,8 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots, hub_omega_na # Run the simulation sim.InitializeSimulation() + assert sc_object.dynManager.getStateObject(smd.nameOfRhoState) is not None + assert sc_object.dynManager.getStateObject(smd.nameOfRhoDotState) is not None sim.ConfigureStopTime(macros.sec2nano(sim_time)) sim.ExecuteSimulation() @@ -729,8 +744,15 @@ def test_prescribed_branching_linear_spring_mass_damper(show_plots, hub_omega_na plt.show() plt.close("all") - # Test check - unit_test_check_verification(orb_ang_mom_n, orb_energy, rot_ang_mom_n) + # Check conservation at every sample. + max_errors = [] + for history in (orb_ang_mom_n, orb_energy, rot_ang_mom_n): + relative_error = (history - history[0]) / np.linalg.norm(history[0]) + max_errors.append(np.max(np.abs(relative_error))) + assert np.all(np.isfinite(max_errors)) + np.testing.assert_array_less(max_errors[:2], 1e-8) + assert max_errors[2] < 1e-4 + return max_errors[2] def create_spacecraft_hub(): From 9b0800bfb706cd78f83520b530f600b41be9ad48 Mon Sep 17 00:00:00 2001 From: Hanspeter Schaub Date: Tue, 10 Nov 2026 11:44:12 -0700 Subject: [PATCH 7/9] [#1596] Document prescribed spring attachments and validation --- .../linearSpringMassDamper.rst | 88 ++++++++++++++++++- .../prescribedMotionStateEffector.rst | 29 ++++-- 2 files changed, 107 insertions(+), 10 deletions(-) diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.rst b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.rst index d09eff3911e..478b3feb59d 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.rst +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.rst @@ -2,7 +2,8 @@ Executive Summary ----------------- -This class in an instantiation of the state effector class and implements an effector representing a oscillating particle +This state effector models a particle constrained to move along a linear spring and damper. +It can attach directly to a spacecraft hub or to a :ref:`prescribedMotionStateEffector`. The module :download:`PDF Description ` @@ -10,6 +11,17 @@ contains further information on this module's function, how to run it, as well as testing. +Module Assumptions and Limitations +---------------------------------- +The particle has one translational degree of freedom along the unit vector ``pHat_B``. +``rho`` measures displacement from the spring equilibrium point, and the spring and damping +forces are ``-k * rho`` and ``-c * rhoDot`` along that axis. + +Prescribed-motion attachments support constant particle mass. Nested variable-mass effectors +are not supported: a particle whose mass is depleted through :ref:`fuelTank` must be attached +directly to :ref:`spacecraft`. + + Message Connection Descriptions ------------------------------- This state effector does not have any input or output messages. @@ -26,3 +38,77 @@ task. Invalid input raises ``BasiliskError``. integrated displacement, velocity, or mass. In particular, it does not restore depleted mass from ``massInit``. Initial values are applied during state registration. See :ref:`effectorInitialization`. + + +User Guide +---------- +For direct hub attachment, express ``pHat_B`` and ``r_PB_B`` in hub-frame components and call +``scObject.addStateEffector(particle)``. The PDF description above covers this configuration. + +.. _linearSpringMassDamperPrescribed: + +Attachment to a Prescribed Body +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The public variable names retain their direct-attachment spelling. For a prescribed-motion +attachment, their reference frame changes to the prescribed body: + +- ``pHat_B`` is the spring's unit direction expressed in prescribed-body components. +- ``r_PB_B`` is the vector from the prescribed body's origin to the particle's equilibrium + point, expressed in prescribed-body components. +- ``rhoInit`` and ``rhoDotInit`` are the initial displacement and displacement rate relative + to that equilibrium point along ``pHat_B``. + +Given an existing spacecraft ``scObject`` and simulation task ``unitTaskName``, the following +example creates a stationary prescribed parent and attaches an oscillating particle to it. +Add each effector to one parent only, before ``InitializeSimulation()``. + +.. code-block:: python + + from Basilisk.simulation import linearSpringMassDamper, prescribedMotionStateEffector + + prescribed_body = prescribedMotionStateEffector.PrescribedMotionStateEffector() + prescribed_body.ModelTag = "springPlatform" + prescribed_body.setMass(100.0) # [kg] + prescribed_body.setIPntPc_P([[20.0, 0.0, 0.0], + [0.0, 25.0, 0.0], + [0.0, 0.0, 30.0]]) # [kg m^2] + prescribed_body.setR_MB_B([0.5, 0.0, 0.0]) # [m] Mount offset in hub components + prescribed_body.setSigma_MB([0.1, -0.2, 0.15]) # [-] Mount attitude as MRPs + + particle = linearSpringMassDamper.LinearSpringMassDamper() + particle.ModelTag = "springParticle" + particle.massInit = 20.0 # [kg] + particle.k = 100.0 # [N/m] + particle.c = 2.0 # [N s/m] + particle.rhoInit = 0.5 # [m] + particle.rhoDotInit = 0.0 # [m/s] + particle.pHat_B = [0.0, 1.0, 0.0] # [-] Direction in prescribed-body components + particle.r_PB_B = [0.4, -0.2, 0.6] # [m] Equilibrium offset from prescribed origin + + prescribed_body.addStateEffector(particle) + scObject.addStateEffector(prescribed_body) + unitTestSim.AddModelToTask(unitTaskName, prescribed_body) + unitTestSim.AddModelToTask(unitTaskName, particle) + +To move the parent, connect kinematic profiler messages to its ``prescribedTranslationInMsg`` +and/or ``prescribedRotationInMsg`` as described in :ref:`prescribedMotionStateEffector`. +The particle's spring and damper motion remains dynamically coupled to the spacecraft while +the prescribed parent translates and rotates. The particle itself has no message ports. + +After initialization, access the particle states through ``scObject.dynManager`` using +``particle.nameOfRhoState`` and ``particle.nameOfRhoDotState``. See +:ref:`test_linearSpringMassDamperPrescribed` for a displacement and velocity logging example. + + +Module Testing +-------------- +:ref:`test_linearSpringMassDamperPrescribed` compares undamped and damped displacement, +velocity, and hub recoil against the analytical two-mass oscillator solution with a fixed +prescribed parent. It also checks independent Newton/Euler force and torque balances for +three oblique configurations with simultaneous prescribed translation and rotation. + +The spring case in :ref:`test_PrescribedMotionStateEffectorBranching` exercises the kinematic +profiler modules, default and renamed hub states, and non-collinear attachment geometry. +It checks orbital energy and orbital and rotational angular momentum over the full simulation history and +checks convergence of the sampled-profiler momentum error when the timestep is halved. +Rotational energy is not required to be constant while the prescribed actuators do work. diff --git a/src/simulation/dynamics/prescribedMotion/prescribedMotionStateEffector.rst b/src/simulation/dynamics/prescribedMotion/prescribedMotionStateEffector.rst index bda176799bb..2c280896839 100644 --- a/src/simulation/dynamics/prescribedMotion/prescribedMotionStateEffector.rst +++ b/src/simulation/dynamics/prescribedMotion/prescribedMotionStateEffector.rst @@ -87,9 +87,11 @@ capability is provided in the following conference paper L. Kiner, and H. Schaub, `“Backsubstitution Method For Prescribed Motion Actuators With Attached Dynamic Sub-Components” `_, AAS Astrodynamics Specialist Conference, Boston, Massachusetts, August 10–14, 2025 -Currently, this branching capability has only been configured for the spinningBodyOneDOF, spinningBodyTwoDOF, and -linearTranslationOneDOF state effectors. See the User Guide section for how to connect these effectors to the -prescribed motion effector. The scenario scripts :ref:`scenarioPrescribedMotionWithRotationBranching` +This branching capability supports :ref:`spinningBodyOneDOFStateEffector`, +:ref:`spinningBodyTwoDOFStateEffector`, :ref:`linearTranslationOneDOFStateEffector`, and +:ref:`linearSpringMassDamper`. See the User Guide section for attachment instructions and +:ref:`linearSpringMassDamperPrescribed` for a spring-particle example and its frame conventions. +The scenario scripts :ref:`scenarioPrescribedMotionWithRotationBranching` and :ref:`scenarioPrescribedMotionWithTranslationBranching` are more complex examples demonstrating how to connect multiple state effectors to the prescribed motion effector. Nested variable-mass state effectors are not supported and must instead be attached @@ -112,12 +114,21 @@ are conserved for the duration of the simulation to verify the module dynamics. Branching Test ^^^^^^^^^^^^^^ -The branching test script for this module has three separate test functions. All tests configure an identical spacecraft -hub an actuating prescribed motion component. The first test connects a free :ref:`spinningBodyOneDOFStateEffector` -to the prescribed component. The second test connects a free :ref:`spinningBodyTwoDOFStateEffector` to the -prescribed component. The third test connects a free :ref:`linearTranslationOneDOFStateEffector` to the prescribed -component. All tests check that the quantities of spacecraft orbital energy, orbital angular momentum, rotational -angular momentum are conserved for the duration of the simulation. +The branching test script has four test functions, one for each supported attached effector. +The first three attach a free :ref:`spinningBodyOneDOFStateEffector`, +:ref:`spinningBodyTwoDOFStateEffector`, or :ref:`linearTranslationOneDOFStateEffector` to an +actuating prescribed component and compare initial and final orbital energy, orbital angular +momentum, and rotational angular momentum. + +The fourth attaches a :ref:`linearSpringMassDamper` with an offset equilibrium point, +a nonzero mount attitude, and a spring axis transverse to the prescribed rotation axis. +It requires the particle states to be registered, tests default and renamed hub angular-velocity +states, and checks conservation throughout the simulation. Relative orbital energy and angular +momentum errors must remain below ``1e-8``. The relative rotational angular-momentum error must +remain below ``1e-4`` and decrease by at least 40 percent when the timestep is halved from +1 ms to 0.5 ms, accounting for sampled profiler commands. Rotational energy may change because +the prescribed actuators do work. Independent analytical particle-motion and Newton/Euler +balance tests are provided in :ref:`test_linearSpringMassDamperPrescribed`. User Guide From 79e3973b57112d5bfc766be9b4a3890768035e83 Mon Sep 17 00:00:00 2001 From: Hanspeter Schaub Date: Tue, 10 Nov 2026 11:45:30 -0700 Subject: [PATCH 8/9] [#1596] Document prescribed coupling callback contracts --- .../linearSpringMassDamper.cpp | 2 -- .../linearSpringMassDamper.h | 22 +++++++++++++++---- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp index 5b8649b1f6b..909a254a1b7 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.cpp @@ -126,8 +126,6 @@ void LinearSpringMassDamper::registerStates(DynParamManager& states) return; } -// Create method linkInPrescribedMotionProperties - /*! This is the method for the SMD to add its contributions to the mass props and mass prop rates of the vehicle * * @param[in] integTime [s] Current integration time. diff --git a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h index c0951f247b7..1db54c72f0d 100644 --- a/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h +++ b/src/simulation/dynamics/LinearSpringMassDamper/linearSpringMassDamper.h @@ -76,8 +76,14 @@ class LinearSpringMassDamper final void Reset(uint64_t CurrentSimNanos) override; void registerStates(DynParamManager& states) override; //!< Method for SMD to register its states void linkInStates(DynParamManager& states) override; //!< Method for SMD to get access of other states - void linkInPrescribedMotionProperties( - DynParamManager& states) override; //!< Method for SMD to get access to prescribed motion properties + /*! @brief Link the prescribed parent's motion relative to the spacecraft hub. + * @param[in] states Dynamics manager owning the prescribed position, velocity, acceleration, + * attitude, angular velocity, and angular acceleration properties. + * @note The prescribed parent assigns the property names during attachment and registers the + * properties before calling this method. The retained references provide the current integration-stage + * values. This callback is only used for prescribed-motion attachments. + */ + void linkInPrescribedMotionProperties(DynParamManager& states) override; void retrieveMassValue(double integTime) override; void calcForceTorqueOnBody(double integTime, Eigen::Vector3d omega_BN_B) override; //!< Force and torque on s/c due to linear spring mass damper @@ -95,8 +101,16 @@ class LinearSpringMassDamper final Eigen::Vector3d rDDot_BN_N, Eigen::Vector3d omegaDot_BN_B, Eigen::MRPd sigma_BN) override; //!< Method for each stateEffector to calculate derivatives - void addPrescribedMotionCouplingContributions(BackSubMatrices& backSubContr) - override; //!< Method for adding coupling contributions for state effector branching on prescribed motion + /*! @brief Add the coupling terms for a particle attached to a prescribed body. + * @param[in,out] backSubContr Particle backsubstitution contributions from updateContributions(), + * expressed in prescribed-parent frame P components. Adds the terms required to couple the particle + * to the spacecraft hub at point B. + * @note Call linkInPrescribedMotionProperties() during initialization, then updateEffectorMassProps() + * and updateContributions() for the current integration stage before this callback. The latter supplies + * the prescribed parent's inertial angular velocity in P components. The prescribed parent transforms + * the resulting contributions into hub-frame components after this callback returns. + */ + void addPrescribedMotionCouplingContributions(BackSubMatrices& backSubContr) override; }; #endif /* LINEAR_SPRING_MASS_DAMPER_H */ From 20565b1380c056362fe132ad00460b531e00f2e4 Mon Sep 17 00:00:00 2001 From: Hanspeter Schaub Date: Tue, 10 Nov 2026 11:46:07 -0700 Subject: [PATCH 9/9] [#1596] Correct units in prescribed spring branching test --- ...est_PrescribedMotionStateEffectorBranching.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py index bc498d1083d..896818aa976 100644 --- a/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py +++ b/src/simulation/dynamics/prescribedMotion/_UnitTest/test_PrescribedMotionStateEffectorBranching.py @@ -646,7 +646,7 @@ def run_prescribed_spring_branch(show_plots, hub_omega_name, time_step): task_name = "unitTask" process_name = "testProcess" sim = SimulationBaseClass.SimBaseClass() - process_rate = macros.sec2nano(time_step) + process_rate = macros.sec2nano(time_step) # [ns] process = sim.CreateNewProcess(process_name) process.addTask(sim.CreateNewTask(task_name, process_rate)) @@ -667,13 +667,13 @@ def run_prescribed_spring_branch(show_plots, hub_omega_name, time_step): one_dof_rotation_profiler.setRotHat_M(prescribed_rot_axis_M) one_dof_rotation_profiler.setThetaDDotMax(prescribed_ang_accel_max) one_dof_rotation_profiler.setThetaInit(prescribed_theta_init) - one_dof_rotation_profiler.setCoastOptionBangDuration(1.0) - one_dof_rotation_profiler.setSmoothingDuration(1.0) + one_dof_rotation_profiler.setCoastOptionBangDuration(1.0) # [s] + one_dof_rotation_profiler.setSmoothingDuration(1.0) # [s] sim.AddModelToTask(task_name, one_dof_rotation_profiler) prescribed_rotation_msg_data = messaging.HingedRigidBodyMsgPayload() prescribed_rotation_msg_data.theta = prescribed_theta_ref - prescribed_rotation_msg_data.thetaDot = 0.0 + prescribed_rotation_msg_data.thetaDot = 0.0 # [rad/s] prescribed_rotation_msg = messaging.HingedRigidBodyMsg().write(prescribed_rotation_msg_data) one_dof_rotation_profiler.spinningBodyInMsg.subscribeTo(prescribed_rotation_msg) prescribed_platform.prescribedRotationInMsg.subscribeTo(one_dof_rotation_profiler.prescribedRotationOutMsg) @@ -684,13 +684,13 @@ def run_prescribed_spring_branch(show_plots, hub_omega_name, time_step): one_dof_translation_profiler.setTransHat_M(prescribed_trans_axis_M) one_dof_translation_profiler.setTransAccelMax(prescribed_trans_accel_max) one_dof_translation_profiler.setTransPosInit(prescribed_pos_init) - one_dof_translation_profiler.setCoastOptionBangDuration(1.0) - one_dof_translation_profiler.setSmoothingDuration(1.0) + one_dof_translation_profiler.setCoastOptionBangDuration(1.0) # [s] + one_dof_translation_profiler.setSmoothingDuration(1.0) # [s] sim.AddModelToTask(task_name, one_dof_translation_profiler) prescribed_translation_msg_data = messaging.LinearTranslationRigidBodyMsgPayload() prescribed_translation_msg_data.rho = prescribed_pos_ref - prescribed_translation_msg_data.rhoDot = 0.0 + prescribed_translation_msg_data.rhoDot = 0.0 # [m/s] prescribed_translation_msg = messaging.LinearTranslationRigidBodyMsg().write(prescribed_translation_msg_data) one_dof_translation_profiler.linearTranslationRigidBodyInMsg.subscribeTo(prescribed_translation_msg) prescribed_platform.prescribedTranslationInMsg.subscribeTo(one_dof_translation_profiler.prescribedTranslationOutMsg) @@ -713,7 +713,7 @@ def run_prescribed_spring_branch(show_plots, hub_omega_name, time_step): # Add Earth gravity to the simulation earth_gravity = gravityEffector.GravBodyData() earth_gravity.planetName = "earth_planet_data" - earth_gravity.mu = 0.3986004415E+15 # meters + earth_gravity.mu = 0.3986004415E+15 # [m^3/s^2] earth_gravity.isCentralBody = True sc_object.gravField.gravBodies = spacecraft.GravBodyVector([earth_gravity])