ENH: store the flight solution as per-phase states and derived quantities

Replace the flat 13-state solution table with a Solution container holding
one PhaseSolution per flight phase. A phase declares the states it
integrates through a StateSchema and still reports the full canonical
state, reconstructing or freezing whatever it does not integrate. The
parachute descent now integrates only position and velocity.

Derived quantities (accelerations, aerodynamic forces and moments, net
thrust) are stored per phase alongside the states. A phase's derivative
returns them rather than writing into the flight, and the caller records
them, so the flight no longer has to track which phase is currently
receiving rows. DerivedQuantity says how to label each quantity and what
to report for the phases that do not compute it.

Also in this change:

- pass the descending parachute to the parachute derivative as an argument
  instead of reading it off the flight
- guard the degenerate cases in find_roots_cubic_function, which a cubic
  Hermite fit reaches whenever a quantity changes at a constant rate
- when loading a .rpy, rebuild the solution before the other attributes,
  since the fallbacks read flight outputs computed from it, and stop the
  net_thrust fallback from resampling the motor's own thrust curve

Author: MateusStano

rocketpy/_encoders.py

@@ -161,7 +161,6 @@ def set_minimal_flight_attributes(flight, obj):
         "atol",
         "time_overshoot",
         "name",
-        "solution",
         "out_of_rail_time",
         "apogee_time",
         "apogee",
@@ -189,22 +188,29 @@ def set_minimal_flight_attributes(flight, obj):
         "net_thrust",
     )
 
+    # The solution is stored either as the new phase-based dict or, for older
+    # saved flights, as a flat list of canonical rows. It is restored before the
+    # attributes above because the fallbacks below read flight outputs, which
+    # are all computed from the solution.
+    raw_solution = obj["solution"]
+    if isinstance(raw_solution, dict):
+        flight.solution = Solution.from_dict(raw_solution)
+    else:
+        flight.solution = Solution.from_legacy_list(raw_solution)
+
     for attribute in attributes:
         try:
             setattr(flight, attribute, obj[attribute])
         except KeyError:
             # Manual resolution of new attributes
             if attribute == "net_thrust":
-                flight.net_thrust = obj["rocket"].motor.thrust
-                flight.net_thrust.set_discrete_based_on_model(flight.speed)
-
-    # The solution is stored either as the new segment-based dict or, for older
-    # saved flights, as a flat list of canonical rows.
-    raw_solution = obj["solution"]
-    if isinstance(raw_solution, dict):
-        flight.solution = Solution.from_dict(raw_solution)
-    else:
-        flight.solution = Solution.from_legacy_list(raw_solution)
+                # Resample the thrust curve onto the flight's time grid without
+                # mutating the motor's own thrust Function.
+                flight.net_thrust = obj[
+                    "rocket"
+                ].motor.thrust.set_discrete_based_on_model(
+                    flight.speed, mutate_self=False
+                )
 
     flight.t_initial = flight.initial_solution[0]
 

rocketpy/rocket/parachute.py

@@ -347,14 +347,23 @@ def triggerfunc(**kwargs):
         elif trigger.lower() == "apogee":
 
             def triggerfunc(**kwargs):
-                state_history = kwargs.get("state_history")
-                if not state_history:
-                    return False
-                return state_history[-1][5] > 0 >= kwargs["state"][5]
+                # Deploy when the rocket stops climbing, by comparing this
+                # check against the previous one. Both come from consecutive
+                # checks of this same parachute, so the moment the vertical
+                # velocity turns negative always falls between two of them and
+                # cannot be stepped over. Comparing against the last stored
+                # trajectory point instead would not be safe: those points
+                # advance on their own schedule, and the turn can happen
+                # between two of this parachute's checks without ever showing
+                # up between two stored points.
+                previous = kwargs["previous_check_state"]
+                if previous is None:
+                    return False  # nothing to compare against yet
+                return previous[5] > 0 >= kwargs["state"][5]
 
             self.triggerfunc = triggerfunc
             self._trigger_is_positional = False
-            self._trigger_needs = frozenset({"state_history"})
+            self._trigger_needs = frozenset()
 
         # Case 4: Invalid trigger input
         else:
@@ -392,14 +401,16 @@ def trigger_callback(**kwargs):
             time = kwargs.get("time", None)
             flight = kwargs.get("flight", None)
 
-            flight._active_parachute = self
             flight.parachute_events.append([time, self])
 
-            kwargs["event"].commands.set_derivative(flight.u_dot_parachute)
+            # The descent dynamics carry this parachute, so the phase keeps
+            # using the right one even when another parachute deploys later.
+            kwargs["event"].commands.set_derivative(
+                flight.u_dot_parachute.rebind(parachute=self)
+            )
             kwargs["event"].commands.start_flight_phase(
                 f"{self.name}_parachute_descent",
                 lag=self.lag,
-                parachute=self,
             )
 
         # Resolve effective needs: explicit override > auto-detected default.

rocketpy/simulation/events/commands.py

@@ -19,7 +19,6 @@ def reset(self):
         self.new_flight_phase = None
         self.new_flight_phase_name = None
         self.new_flight_phase_lag = 0
-        self.new_flight_phase_parachute = None
         self._terminate = False
         self.terminate_phase_name = None
 
@@ -39,11 +38,10 @@ def set_derivative(self, derivative):
         self.new_derivative = derivative
         self.new_derivative_set = True
 
-    def start_flight_phase(self, phase_name=None, lag=0, parachute=None):
+    def start_flight_phase(self, phase_name=None, lag=0):
         self.new_flight_phase = True
         self.new_flight_phase_name = phase_name
         self.new_flight_phase_lag = lag
-        self.new_flight_phase_parachute = parachute
 
     def terminate_flight(self):
         self._terminate = True

rocketpy/simulation/events/event.py

@@ -16,7 +16,7 @@
 PRESETS = {
     "apogee": lambda **kwargs: (
         len(kwargs["flight"].solution) >= 2
-        and kwargs["flight"].solution.view(-2)[1]["vz"] > 0 >= kwargs["state"][5]
+        and kwargs["flight"].solution.at_index(-2)["vz"] > 0 >= kwargs["state"][5]
     ),
     "burnout": lambda **kwargs: (
         kwargs.get("time") >= kwargs["rocket"].motor.burn_out_time
@@ -83,7 +83,26 @@ def __init__(  # pylint: disable=too-many-arguments
             ``step_size`` (float, s),
             ``height_agl`` (float, m),
             ``event`` (this :class:`Event` instance),
-            ``sampling_rate`` (float, Hz or ``None``).
+            ``sampling_rate`` (float, Hz or ``None``),
+            ``previous_check_state`` (the ``state`` from the previous time this
+            event was checked, or ``None`` on the first check),
+            ``previous_check_time`` (float, s, or ``None`` on the first check).
+
+            To fire when a value crosses a threshold, compare ``state``
+            against ``previous_check_state``. Those two are always consecutive
+            checks of this same event, so a crossing cannot slip between
+            them::
+
+                def trigger(**kwargs):
+                    previous = kwargs["previous_check_state"]
+                    if previous is None:
+                        return False
+                    return previous[5] > 0 >= kwargs["state"][5]
+
+            Do not use ``state_history[-1]`` for this. That is the last
+            *stored* trajectory point, which advances independently of when
+            this event is checked, so a crossing can fall between the two and
+            never be seen.
             The following keys are only injected when declared via ``needs``:
             ``pressure`` (float, Pa),
             ``state_dot`` (list, time derivative of ``state``),
@@ -215,6 +234,10 @@ def __init__(  # pylint: disable=too-many-arguments
         self.callback_log = []
         self.triggered_times = []
         self._trigger_checked = False
+        # State and time from the last time this event's trigger was checked.
+        # ``None`` until the first check.
+        self._previous_state = None
+        self._previous_time = None
 
         self.is_discrete = self.sampling_rate is not None
         self.sampling_interval = (
@@ -275,6 +298,8 @@ def reset(self):
         self.callback_log.clear()
         self.triggered_times.clear()
         self._trigger_checked = False
+        self._previous_state = None
+        self._previous_time = None
         self.context = deepcopy(self._initial_context)
 
         # Restore enable/disable time history to initial snapshot
@@ -327,14 +352,27 @@ def __call__(self, trigger_only=False, callback_only=False, reset=True, **kwargs
 
         kwargs["event"] = self
         kwargs["sampling_rate"] = self.sampling_rate
+        # What this event saw the last time it was checked. Comparing against
+        # it is the reliable way to notice a value crossing a threshold, since
+        # the two states are always consecutive checks of this same event.
+        kwargs["previous_check_state"] = self._previous_state
+        kwargs["previous_check_time"] = self._previous_time
 
         # --- Trigger Phase ---
         # Skip evaluating triggers if we are only running the callback.
         if callback_only is False:
             if self._call_enable_on(**kwargs) is False:
                 return False
 
-            if self._call_trigger(**kwargs) is False:
+            triggered = self._call_trigger(**kwargs)
+
+            # Remember this check, whether or not it triggered, so the next one
+            # can be compared against it. Recorded before the early return so a
+            # trigger that did not fire is still part of the sequence.
+            self._previous_state = kwargs.get("state")
+            self._previous_time = kwargs.get("time")
+
+            if triggered is False:
                 return False
 
             if trigger_only:
@@ -447,7 +485,7 @@ def _compute_exact_time(self, **kwargs):
         # The exact-time functions read the state by its canonical position, so
         # feed them canonical endpoints and a canonical view of the dense
         # output. The raw dense output is kept so the exact state can be written
-        # back into the (possibly reduced) tail segment of the solution.
+        # back into the (possibly reduced) tail phase of the solution.
         schema = phase.dynamics.schema
         raw_interpolator = phase.solver.dense_output()
         if schema.is_canonical:

rocketpy/simulation/events/event_builders.py

@@ -118,7 +118,7 @@ def apogee_trigger(**kwargs):
     if len(flight.apogee_state) != 1 or len(flight.solution) < 2:
         return False
 
-    previous_vz = flight.solution.view(-2)[1]["vz"]
+    previous_vz = flight.solution.at_index(-2)["vz"]
     current_vz = state[5]
     return previous_vz > 0 >= current_vz