feat: add metadata for segment and sequence solutions

bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Segment.cpp

@@ -22,6 +22,16 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Segment(pybind11::module&
         .def_readonly("dynamics", &Segment::Solution::dynamics)
         .def_readonly("states", &Segment::Solution::states)
         .def_readonly("condition_is_satisfied", &Segment::Solution::conditionIsSatisfied)
+        .def_readonly("segment_type", &Segment::Solution::segmentType)
+
+        .def("access_start_instant", &Segment::Solution::accessStartInstant)
+        .def("access_end_instant", &Segment::Solution::accessEndInstant)
+
+        .def("get_initial_mass", &Segment::Solution::getInitialMass)
+        .def("get_final_mass", &Segment::Solution::getFinalMass)
+
+        .def("compute_delta_mass", &Segment::Solution::computeDeltaMass)
+        .def("compute_delta_v", &Segment::Solution::computeDeltaV, arg("specific_impulse"))
 
         ;
 

bindings/python/src/OpenSpaceToolkitAstrodynamicsPy/Trajectory/Sequence.cpp

@@ -24,7 +24,16 @@ inline void OpenSpaceToolkitAstrodynamicsPy_Trajectory_Sequence(pybind11::module
         .def_readonly("segment_solutions", &Sequence::Solution::segmentSolutions)
         .def_readonly("execution_is_complete", &Sequence::Solution::executionIsComplete)
 
+        .def("access_start_instant", &Sequence::Solution::accessStartInstant)
+        .def("access_end_instant", &Sequence::Solution::accessEndInstant)
+
         .def("get_states", &Sequence::Solution::getStates)
+        .def("get_initial_mass", &Sequence::Solution::getInitialMass)
+        .def("get_final_mass", &Sequence::Solution::getFinalMass)
+        .def("get_propagation_duration", &Sequence::Solution::getPropagationDuration)
+
+        .def("compute_delta_mass", &Sequence::Solution::computeDeltaMass)
+        .def("compute_delta_v", &Sequence::Solution::computeDeltaV, arg("specific_impulse"))
 
         ;
 

bindings/python/test/trajectory/test_segment.py

@@ -7,8 +7,6 @@
 from ostk.physics.time import DateTime
 from ostk.physics.time import Scale
 from ostk.physics.time import Duration
-from ostk.physics.coordinate import Position
-from ostk.physics.coordinate import Velocity
 from ostk.physics.coordinate import Frame
 from ostk.physics.environment.objects.celestial_bodies import Earth
 
@@ -17,11 +15,13 @@
 from ostk.astrodynamics.dynamics import CentralBodyGravity
 from ostk.astrodynamics.dynamics import PositionDerivative
 from ostk.astrodynamics.dynamics.thruster import ConstantThrust
-from ostk.astrodynamics.trajectory import LocalOrbitalFrameDirection
-from ostk.astrodynamics.trajectory import LocalOrbitalFrameFactory
 from ostk.astrodynamics.trajectory import State
 from ostk.astrodynamics.trajectory import Segment
 from ostk.astrodynamics.event_condition import InstantCondition
+from ostk.astrodynamics.trajectory.state import CoordinatesSubset
+from ostk.astrodynamics.trajectory.state import CoordinatesBroker
+from ostk.astrodynamics.trajectory.state.coordinates_subset import CartesianPosition
+from ostk.astrodynamics.trajectory.state.coordinates_subset import CartesianVelocity
 
 
 @pytest.fixture
@@ -31,14 +31,25 @@ def environment() -> Environment:
 
 @pytest.fixture
 def state() -> State:
+    instant: Instant = Instant.date_time(DateTime(2018, 1, 1, 0, 0, 0), Scale.UTC)
+    coordinates: list[float] = [
+        7500000.0,
+        0.0,
+        0.0,
+        0.0,
+        5335.865450622126,
+        5335.865450622126,
+        300.0,
+    ]
     frame: Frame = Frame.GCRF()
-    position: Position = Position.meters([7500000.0, 0.0, 0.0], frame)
-    velocity: Velocity = Velocity.meters_per_second(
-        [0.0, 5335.865450622126, 5335.865450622126], frame
+    coordinates_broker: CoordinatesBroker = CoordinatesBroker(
+        [
+            CartesianPosition.default(),
+            CartesianVelocity.default(),
+            CoordinatesSubset.mass(),
+        ]
     )
-
-    instant: Instant = Instant.date_time(DateTime(2018, 1, 1, 0, 0, 0), Scale.UTC)
-    return State(instant, position, velocity)
+    return State(instant, coordinates, frame, coordinates_broker)
 
 
 @pytest.fixture
@@ -87,21 +98,28 @@ def coast_duration_segment(
     instant_condition: InstantCondition,
     dynamics: list,
     numerical_solver: NumericalSolver,
-):
+) -> Segment:
     return Segment.coast(name, instant_condition, dynamics, numerical_solver)
 
 
 @pytest.fixture
-def thruster_dynamics() -> ConstantThrust:
-    return ConstantThrust(
-        satellite_system=SatelliteSystem.default(),
-        thrust_direction=LocalOrbitalFrameDirection(
-            vector=[1.0, 0.0, 0.0],
-            local_orbital_frame_factory=LocalOrbitalFrameFactory.VNC(Frame.GCRF()),
-        ),
+def maneuver_segment(
+    name: str,
+    instant_condition: InstantCondition,
+    dynamics: list,
+    numerical_solver: NumericalSolver,
+    thruster_dynamics: ConstantThrust,
+) -> Segment:
+    return Segment.maneuver(
+        name, instant_condition, thruster_dynamics, dynamics, numerical_solver
     )
 
 
+@pytest.fixture
+def thruster_dynamics() -> ConstantThrust:
+    return ConstantThrust.intrack(satellite_system=SatelliteSystem.default())
+
+
 class TestSegment:
     def test_get_name(self, coast_duration_segment: Segment, name: str):
         assert coast_duration_segment.get_name() == name
@@ -163,9 +181,9 @@ def test_solve(
         self,
         state: State,
         end_instant: Instant,
-        coast_duration_segment: Segment,
+        maneuver_segment: Segment,
     ):
-        solution = coast_duration_segment.solve(state)
+        solution = maneuver_segment.solve(state)
 
         assert (
             pytest.approx(
@@ -174,5 +192,19 @@ def test_solve(
             )
             == 0.0
         )
+
+        assert solution.name is not None
+        assert solution.dynamics is not None
         assert len(solution.states) > 0
-        assert solution.condition_is_satisfied
+        assert solution.condition_is_satisfied is True
+        assert solution.segment_type is not None
+
+        assert solution.access_start_instant() is not None
+        assert solution.access_end_instant() is not None
+
+        assert solution.get_initial_mass() is not None
+        assert solution.get_final_mass() is not None
+
+        assert solution.get_propagation_duration() is not None
+        assert solution.compute_delta_mass() is not None
+        assert solution.compute_delta_v(1500.0) is not None

bindings/python/test/trajectory/test_sequence.py

@@ -376,5 +376,15 @@ def test_solve(self, state: State, sequence: Sequence, segments: list[Segment]):
 
         assert len(solution.segment_solutions) == len(segments)
 
-        assert solution.get_states() is not None
         assert solution.execution_is_complete
+
+        assert solution.access_start_instant() is not None
+        assert solution.access_end_instant() is not None
+
+        assert solution.get_states() is not None
+        assert solution.get_initial_mass() is not None
+        assert solution.get_final_mass() is not None
+        assert solution.get_propagation_duration() is not None
+
+        assert solution.compute_delta_mass() is not None
+        assert solution.compute_delta_v(1500.0) is not None

include/OpenSpaceToolkit/Astrodynamics/Trajectory/Segment.hpp

@@ -29,6 +29,7 @@ using ostk::core::types::String;
 
 using ostk::physics::time::Instant;
 using ostk::physics::time::Duration;
+using ostk::physics::units::Mass;
 
 using ostk::astro::trajectory::State;
 using ostk::astro::trajectory::state::NumericalSolver;
@@ -54,13 +55,24 @@ class Segment
             const String& aName,
             const Array<Shared<Dynamics>>& aDynamicsArray,
             const Array<State>& aStates,
-            const bool& aConditionIsSatisfied
+            const bool& aConditionIsSatisfied,
+            const Segment::Type& aSegmentType
         );
 
         String name;                       /// Name of the segment.
         Array<Shared<Dynamics>> dynamics;  /// List of dynamics used.
         Array<State> states;               /// Array of states for the segment.
         bool conditionIsSatisfied;         /// True if the event condition is satisfied.
+        Segment::Type segmentType;         /// Type of segment.
+
+        const Instant& accessStartInstant() const;
+        const Instant& accessEndInstant() const;
+
+        Mass getInitialMass() const;
+        Mass getFinalMass() const;
+
+        Real computeDeltaV(const Real& aSpecificImpulse) const;
+        Mass computeDeltaMass() const;
     };
 
     /// @brief                  Output stream operator

include/OpenSpaceToolkit/Astrodynamics/Trajectory/Sequence.hpp

@@ -6,6 +6,7 @@
 #include <OpenSpaceToolkit/Core/Containers/Array.hpp>
 
 #include <OpenSpaceToolkit/Physics/Environment.hpp>
+#include <OpenSpaceToolkit/Physics/Units/Mass.hpp>
 
 #include <OpenSpaceToolkit/Astrodynamics/Flight/System/SatelliteSystem.hpp>
 #include <OpenSpaceToolkit/Astrodynamics/Trajectory/Segment.hpp>
@@ -18,9 +19,11 @@ namespace astro
 namespace trajectory
 {
 
+using ostk::core::types::Real;
 using ostk::core::ctnr::Array;
 
 using ostk::physics::Environment;
+using ostk::physics::units::Mass;
 
 using ostk::astro::trajectory::Segment;
 using ostk::astro::trajectory::State;
@@ -37,7 +40,16 @@ class Sequence
         Array<Segment::Solution> segmentSolutions;
         bool executionIsComplete;
 
+        const Instant& accessStartInstant() const;
+        const Instant& accessEndInstant() const;
+
         Array<State> getStates() const;
+        Mass getInitialMass() const;
+        Mass getFinalMass() const;
+        Duration getPropagationDuration() const;
+
+        Mass computeDeltaMass() const;
+        Real computeDeltaV(const Real& aSpecificImpulse) const;
     };
 
     /// @brief                  Constructor

src/OpenSpaceToolkit/Astrodynamics/Dynamics.cpp

@@ -134,7 +134,7 @@ void Dynamics::applyContribution(
 )
 {
     Index offset = 0;
-    for (Pair<Index, Size> pair : writeInfo)
+    for (const Pair<Index, Size>& pair : writeInfo)
     {
         const Index subsetOffset = pair.first;
         const Size subsetSize = pair.second;

src/OpenSpaceToolkit/Astrodynamics/Trajectory/Segment.cpp

@@ -1,7 +1,10 @@
 /// Apache License 2.0
 
+#include <OpenSpaceToolkit/Physics/Environment/Gravitational/Earth.hpp>
+
 #include <OpenSpaceToolkit/Astrodynamics/Trajectory/Propagator.hpp>
 #include <OpenSpaceToolkit/Astrodynamics/Trajectory/Segment.hpp>
+#include <OpenSpaceToolkit/Astrodynamics/Trajectory/State/CoordinatesSubset.hpp>
 
 namespace ostk
 {
@@ -11,20 +14,86 @@
 {
 
 using ostk::physics::time::Duration;
+using EarthGravitationalModel = ostk::physics::environment::gravitational::Earth;
 
 using ostk::astro::trajectory::Propagator;
+using ostk::astro::trajectory::state::CoordinatesSubset;
 
 Segment::Solution::Solution(
     const String& aName,
     const Array<Shared<Dynamics>>& aDynamicsArray,
     const Array<State>& aStates,
-    const bool& aConditionIsSatisfied
+    const bool& aConditionIsSatisfied,
+    const Segment::Type& aSegmentType
 )
     : name(aName),
       dynamics(aDynamicsArray),
       states(aStates),
-      conditionIsSatisfied(aConditionIsSatisfied)
+      conditionIsSatisfied(aConditionIsSatisfied),
+      segmentType(aSegmentType)
+{
+}
+
+const Instant& Segment::Solution::accessStartInstant() const
+{
+    if (states.isEmpty())
+    {
+        throw ostk::core::error::RuntimeError("No solution available.");
+    }
+
+    return states.accessFirst().accessInstant();
+}
+
+const Instant& Segment::Solution::accessEndInstant() const
+{
+    if (states.isEmpty())
+    {
+        throw ostk::core::error::RuntimeError("No solution available.");
+    }
+
+    return states.accessLast().accessInstant();
+}
+
+Mass Segment::Solution::getInitialMass() const
+{
+    if (states.isEmpty())
+    {
+        throw ostk::core::error::RuntimeError("No solution available.");
+    }
+
+    return Mass::Kilograms(states.accessFirst().extractCoordinates(CoordinatesSubset::Mass())[0]);
+}
+
+Mass Segment::Solution::getFinalMass() const
 {
+    if (states.isEmpty())
+    {
+        throw ostk::core::error::RuntimeError("No solution available.");
+    }
+
+    return Mass::Kilograms(states.accessLast().extractCoordinates(CoordinatesSubset::Mass())[0]);
+}
+
+Real Segment::Solution::computeDeltaV(const Real& aSpecificImpulse) const
+{
+    // TBM: This is only valid for constant thrust, constant Isp
+    if (segmentType != Segment::Type::Maneuver)
+    {
+        return 0.0;
+    }
+
+    return aSpecificImpulse * EarthGravitationalModel::gravityConstant *
+           std::log(getInitialMass().inKilograms() / getFinalMass().inKilograms());
+}
+
+Mass Segment::Solution::computeDeltaMass() const
+{
+    if (segmentType != Segment::Type::Maneuver)
+    {
+        return Mass::Kilograms(0.0);
+    }
+
+    return Mass::Kilograms(getInitialMass().inKilograms() - getFinalMass().inKilograms());
 }
 
 Segment::Segment(
@@ -110,16 +179,16 @@
         dynamics_,
     };
 
-    const Instant startInstant = aState.getInstant();
-
-    const NumericalSolver::ConditionSolution conditionSolution =
-        propagator.calculateStateToCondition(aState, startInstant + maximumPropagationDuration, *eventCondition_);
+    const NumericalSolver::ConditionSolution conditionSolution = propagator.calculateStateToCondition(
+        aState, aState.accessInstant() + maximumPropagationDuration, *eventCondition_
+    );
 
     return {
         name_,
         dynamics_,
         propagator.accessNumericalSolver().accessObservedStates(),
         conditionSolution.conditionIsSatisfied,
+        type_,
     };
 }