From 9fe3a37c93069156ac45a6f241f5935832ddb21d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Fri, 14 Mar 2025 18:46:09 +0100 Subject: [PATCH 01/37] wind factor bug corrected the wind factor wasn't applied to the env.wind_velocity properties --- rocketpy/stochastic/stochastic_environment.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rocketpy/stochastic/stochastic_environment.py b/rocketpy/stochastic/stochastic_environment.py index 58afe0fed..e0fc33eec 100644 --- a/rocketpy/stochastic/stochastic_environment.py +++ b/rocketpy/stochastic/stochastic_environment.py @@ -184,5 +184,6 @@ def create_object(self): # get original attribute value and multiply by factor attribute_name = f"_{key.replace('_factor', '')}" value = getattr(self, attribute_name) * value + key = f"{key.replace('_factor', '')}" setattr(self.obj, key, value) return self.obj From f6efa816442a740279ebaf880458049d28799b3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sat, 15 Mar 2025 11:46:06 +0100 Subject: [PATCH 02/37] BUG: StochasticModel visualize attributes of a uniform distribution It showed the nominal and the standard deviation values and it doesn't make sense in a uniform distribution. In a np.random.uniform the 'nominal value' is the lower bound of the distribution, and the 'standard deviation' value is the upper bound. Now, a new condition has been added for the uniform distributions where the mean and semi range are calculated and showed. This way the visualize_attribute function will show the whole range where the random values are uniformly taken in --- rocketpy/stochastic/stochastic_model.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index c82232fcb..e44e2613d 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -486,11 +486,20 @@ def format_attribute(attr, value): ) elif isinstance(value, tuple): nominal_value, std_dev, dist_func = value - return ( - f"\t{attr.ljust(max_str_length)} " - f"{nominal_value:.5f} ± " - f"{std_dev:.5f} ({dist_func.__name__})" - ) + if dist_func == np.random.uniform: + mean = (std_dev + nominal_value) / 2 + semi_range = (std_dev - nominal_value) / 2 + return ( + f"\t{attr.ljust(max_str_length)} " + f"{mean:.5f} ± " + f"{semi_range:.5f} ({dist_func.__name__})" + ) + else: + return ( + f"\t{attr.ljust(max_str_length)} " + f"{nominal_value:.5f} ± " + f"{std_dev:.5f} ({dist_func.__name__})" + ) return None attributes = {k: v for k, v in self.__dict__.items() if not k.startswith("_")} From 04337ab79982f28196d4b3950d36407cac6752b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sat, 15 Mar 2025 16:11:27 +0100 Subject: [PATCH 03/37] variable names corrections --- rocketpy/stochastic/stochastic_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index e44e2613d..e833e01bf 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -488,11 +488,11 @@ def format_attribute(attr, value): nominal_value, std_dev, dist_func = value if dist_func == np.random.uniform: mean = (std_dev + nominal_value) / 2 - semi_range = (std_dev - nominal_value) / 2 + half_range = (std_dev - nominal_value) / 2 return ( f"\t{attr.ljust(max_str_length)} " f"{mean:.5f} ± " - f"{semi_range:.5f} ({dist_func.__name__})" + f"{half_range:.5f} ({dist_func.__name__})" ) else: return ( From 18f8323f39798ae90caada573660dcf63b113769 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sat, 15 Mar 2025 17:46:04 +0100 Subject: [PATCH 04/37] Corrections requested by the pylint test --- rocketpy/stochastic/stochastic_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index e833e01bf..39fc478b7 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -486,7 +486,7 @@ def format_attribute(attr, value): ) elif isinstance(value, tuple): nominal_value, std_dev, dist_func = value - if dist_func == np.random.uniform: + if callable(dist_func) and dist_func.__name__ == "uniform": mean = (std_dev + nominal_value) / 2 half_range = (std_dev - nominal_value) / 2 return ( From a46de46eaa87bd8de797e030738facdbe4492587 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sun, 16 Mar 2025 13:21:19 +0100 Subject: [PATCH 05/37] ENH: add multiplication for 2D functions in rocketpy.function Added the ability to multiply functions with 2D domains in the __mul__ function --- rocketpy/mathutils/function.py | 89 ++++++++++++++++++++++------------ 1 file changed, 59 insertions(+), 30 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 143912357..869f5a7b9 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2222,7 +2222,7 @@ def __rsub__(self, other): def __mul__(self, other): """Multiplies a Function object and returns a new Function object which gives the result of the multiplication. Only implemented for 1D - domains. + and 2D domains. Parameters ---------- @@ -2238,7 +2238,7 @@ def __mul__(self, other): Returns ------- result : Function - A Function object which gives the result of self(x)*other(x). + A Function object which gives the result of self(x,y)*other(x,y). """ self_source_is_array = isinstance(self.source, np.ndarray) other_source_is_array = ( @@ -2250,37 +2250,66 @@ def __mul__(self, other): interp = self.__interpolation__ extrap = self.__extrapolation__ - if ( - self_source_is_array - and other_source_is_array - and np.array_equal(self.x_array, other.x_array) - ): - source = np.column_stack((self.x_array, self.y_array * other.y_array)) - outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})" - return Function(source, inputs, outputs, interp, extrap) - elif isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array( - other - ): - if not self_source_is_array: - return Function(lambda x: (self.get_value_opt(x) * other), inputs) - source = np.column_stack((self.x_array, np.multiply(self.y_array, other))) - outputs = f"({self.__outputs__[0]}*{other})" - return Function( - source, - inputs, - outputs, - interp, - extrap, - ) - elif callable(other): - return Function(lambda x: (self.get_value_opt(x) * other(x)), inputs) - else: - raise TypeError("Unsupported type for multiplication") + if self.__dom_dim__ == 1: + if ( + self_source_is_array + and other_source_is_array + and np.array_equal(self.x_array, other.x_array) + ): + source = np.column_stack((self.x_array, self.y_array * other.y_array)) + outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})" + return Function(source, inputs, outputs, interp, extrap) + elif isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array( + other + ): + if not self_source_is_array: + return Function(lambda x: (self.get_value_opt(x) * other), inputs) + source = np.column_stack((self.x_array, np.multiply(self.y_array, other))) + outputs = f"({self.__outputs__[0]}*{other})" + return Function( + source, + inputs, + outputs, + interp, + extrap, + ) + elif callable(other): + return Function(lambda x: (self.get_value_opt(x) * other(x)), inputs) + else: + raise TypeError("Unsupported type for multiplication") + elif self.__dom_dim__ == 2: + if ( + self_source_is_array + and other_source_is_array + and np.array_equal(self.x_array, other.x_array) + and np.array_equal(self.y_array, other.y_array) + ): + source = np.column_stack((self.x_array, self.y_array, self.z_array * other.z_array)) + outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})" + return Function(source, inputs, outputs, interp, extrap) + elif isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array( + other + ): + if not self_source_is_array: + return Function(lambda x,y: (self.get_value_opt(x,y) * other), inputs) + source = np.column_stack((self.x_array, self.y_array, np.multiply(self.z_array, other))) + outputs = f"({self.__outputs__[0]}*{other})" + return Function( + source, + inputs, + outputs, + interp, + extrap, + ) + elif callable(other): + return Function(lambda x,y: (self.get_value_opt(x,y) * other(x)), inputs) + else: + raise TypeError("Unsupported type for multiplication") def __rmul__(self, other): """Multiplies 'other' by a Function object and returns a new Function object which gives the result of the multiplication. Only implemented for - 1D domains. + 1D and 2D domains. Parameters ---------- @@ -2290,7 +2319,7 @@ def __rmul__(self, other): Returns ------- result : Function - A Function object which gives the result of other(x)*self(x). + A Function object which gives the result of other(x,y)*self(x,y). """ return self * other From 718136045b7678a2c3175aee7a1821e8c5387f6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sun, 16 Mar 2025 13:30:46 +0100 Subject: [PATCH 06/37] ENH: StochasticAirBrakes class created The StochasticAirBrakes class has been created. The __init__.py files in the stochastic and rocketpy folders have also been modified accordingly to incorporate this new class --- rocketpy/__init__.py | 1 + rocketpy/stochastic/__init__.py | 1 + .../stochastic/stochastic_aero_surfaces.py | 90 +++++++++++++++++++ 3 files changed, 92 insertions(+) diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index 539b8b2cb..e0436bf02 100644 --- a/rocketpy/__init__.py +++ b/rocketpy/__init__.py @@ -44,6 +44,7 @@ from .sensors import Accelerometer, Barometer, GnssReceiver, Gyroscope from .simulation import Flight, MonteCarlo from .stochastic import ( + StochasticAirBrakes, StochasticEllipticalFins, StochasticEnvironment, StochasticFlight, diff --git a/rocketpy/stochastic/__init__.py b/rocketpy/stochastic/__init__.py index 692eca85f..b1e146246 100644 --- a/rocketpy/stochastic/__init__.py +++ b/rocketpy/stochastic/__init__.py @@ -6,6 +6,7 @@ """ from .stochastic_aero_surfaces import ( + StochasticAirBrakes, StochasticEllipticalFins, StochasticNoseCone, StochasticRailButtons, diff --git a/rocketpy/stochastic/stochastic_aero_surfaces.py b/rocketpy/stochastic/stochastic_aero_surfaces.py index 31c1e9efc..17f22337b 100644 --- a/rocketpy/stochastic/stochastic_aero_surfaces.py +++ b/rocketpy/stochastic/stochastic_aero_surfaces.py @@ -9,6 +9,7 @@ RailButtons, Tail, TrapezoidalFins, + AirBrakes, ) from .stochastic_model import StochasticModel @@ -432,3 +433,92 @@ def create_object(self): """ generated_dict = next(self.dict_generator()) return RailButtons(**generated_dict) + + +class StochasticAirBrakes(StochasticModel): + """A Stochastic Air Brakes class that inherits from StochasticModel. + + See Also + -------- + :ref:`stochastic_model` and + :class:`AirBrakes ` + + Attributes + ---------- + object : AirBrakes + AirBrakes object to be used for validation. + drag_coefficient_curve : list + The drag coefficient curve of the air brakes can account for + either the air brakes' drag alone or the combined drag of both + the rocket and the air brakes. + drag_coefficient_curve_factor : tuple, list, int, float + The drag curve factor of the air brakes. + reference_area : tuple, list, int, float + Reference area used to non-dimensionalize the drag coefficients. + deployment_level : tuple, list, int, float + Initial deployment level, ranging from 0 to 1. + name : list[str] + List with the air brakes object name. This attribute can't be randomized. + """ + + def __init__( + self, + air_brakes, + drag_coefficient_curve=None, + drag_coefficient_curve_factor=(1, 0), + reference_area=None, + clamp=None, + override_rocket_drag=None, + deployment_level=(0,0), + ): + """Initializes the Stochastic AirBrakes class. + + See Also + -------- + :ref:`stochastic_model` + + Parameters + ---------- + air_brakes : AirBrakes + AirBrakes object to be used for validation. + drag_coefficient_curve : list, optional + The drag coefficient curve of the air brakes can account for + either the air brakes' drag alone or the combined drag of both + the rocket and the air brakes. + drag_coefficient_curve_factor : tuple, list, int, float, optional + The drag curve factor of the air brakes. + reference_area : tuple, list, int, float, optional + Reference area used to non-dimensionalize the drag coefficients. + deployment_level : tuple, list, int, float, optional + Initial deployment level, ranging from 0 to 1. + """ + super().__init__( + air_brakes, + drag_coefficient_curve=drag_coefficient_curve, + drag_coefficient_curve_factor=drag_coefficient_curve_factor, + reference_area=reference_area, + clamp=clamp, + override_rocket_drag=override_rocket_drag, + deployment_level=deployment_level, + name=None, + ) + + def create_object(self): + """Creates and returns an AirBrakes object from the randomly generated + input arguments. + + Returns + ------- + air_brake : AirBrakes + AirBrakes object with the randomly generated input arguments. + """ + generated_dict = next(self.dict_generator()) + air_brakes = AirBrakes( + drag_coefficient_curve=generated_dict['drag_coefficient_curve'], + reference_area=generated_dict['reference_area'], + clamp=generated_dict['clamp'], + override_rocket_drag=generated_dict['override_rocket_drag'], + deployment_level=generated_dict['deployment_level'], + ) + air_brakes.drag_coefficient *= generated_dict['drag_coefficient_curve_factor'] + return air_brakes From e1549e0a2ba12302b617939544ea48f8c895e219 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 10:33:59 +0100 Subject: [PATCH 07/37] ENH: set_air_brakes function created This functions appends an airbrake and controller objects previuosly created to the rocket --- rocketpy/rocket/rocket.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 609e1b76c..e40b71a9c 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1534,6 +1534,15 @@ def add_sensor(self, sensor, position): except KeyError: sensor._attached_rockets[self] = 1 + def set_air_brakes( + self, + air_brakes, + controller, + ): + + self.air_brakes.append(air_brakes) + self._add_controllers(controller) + def add_air_brakes( self, drag_coefficient_curve, From 2d2a0e870d5e1741460de132ed010e621164a90c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 11:08:00 +0100 Subject: [PATCH 08/37] ENH: add StochasticAirBrake to rocketpy.stochastic_rocket Some functions has been modified and other has been created in order to include the new StochasticAirBrakes feature into the StochasticRocket class. A new function named 'add_air_brakes' has been created to append a StochasticAirBrakes and Controller objects to the StochasticRocket object. A new function '_create_air_brake' has been introduced to create a sample of an AirBrake object through a StochasticAirBrake object. Enventually, the 'create_object' function has been modified to add the sampled AirBrakes to the sampled Rocket --- rocketpy/stochastic/stochastic_rocket.py | 47 ++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 224262ff1..c8f91b342 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -3,11 +3,13 @@ import warnings from random import choice +from rocketpy.control import _Controller from rocketpy.mathutils.vector_matrix import Vector from rocketpy.motors.empty_motor import EmptyMotor from rocketpy.motors.motor import GenericMotor, Motor from rocketpy.motors.solid_motor import SolidMotor from rocketpy.rocket.aero_surface import ( + AirBrakes, EllipticalFins, NoseCone, RailButtons, @@ -21,6 +23,7 @@ from rocketpy.stochastic.stochastic_motor_model import StochasticMotorModel from .stochastic_aero_surfaces import ( + StochasticAirBrakes, StochasticEllipticalFins, StochasticNoseCone, StochasticRailButtons, @@ -148,6 +151,7 @@ def __init__( self.motors = Components() self.aerodynamic_surfaces = Components() self.rail_buttons = Components() + self.air_brakes = [] self.parachutes = [] self.__components_map = {} super().__init__( @@ -389,6 +393,36 @@ def set_rail_buttons( rail_buttons, self._validate_position(rail_buttons, lower_button_position) ) + def add_air_brakes( + self, + air_brakes, + controller_function, + sampling_rate, + initial_observed_variables=None, + return_controller=False, + controller_name="AirBrakes Controller", + ): + + # checks if input is a StochasticAirbrakes type + if not isinstance(air_brakes, (AirBrakes, StochasticAirBrakes)): + raise TypeError( + "`air_brake` must be of AirBrakes or StochasticAirBrakes type" + ) + if isinstance(air_brakes, AirBrakes): + air_brakes = StochasticAirBrakes(air_brakes=air_brakes) + + self.air_brakes.append(air_brakes) + _controller = _Controller( + interactive_objects=air_brakes, + controller_function=controller_function, + sampling_rate=sampling_rate, + initial_observed_variables=initial_observed_variables, + name=controller_name, + ) + self.air_brake_controller = _controller + if return_controller: + return _controller + def _validate_position(self, validated_object, position): """Validate the position argument. @@ -531,6 +565,7 @@ def dict_generator(self): generated_dict["motors"] = [] generated_dict["aerodynamic_surfaces"] = [] generated_dict["rail_buttons"] = [] + generated_dict["air_brakes"] = [] generated_dict["parachutes"] = [] self.last_rnd_dict = generated_dict yield generated_dict @@ -572,6 +607,11 @@ def _create_rail_buttons(self, component_stochastic_rail_buttons): upper_button_position_rnd ) return rail_buttons, lower_button_position_rnd, upper_button_position_rnd + + def _create_air_brake(self, stochastic_air_brake): + air_brake = stochastic_air_brake.create_object() + self.last_rnd_dict["air_brakes"].append(stochastic_air_brake.last_rnd_dict) + return air_brake def _create_parachute(self, stochastic_parachute): parachute = stochastic_parachute.create_object() @@ -617,6 +657,13 @@ def create_object(self): surface, position_rnd = self._create_surface(component_surface) rocket.add_surfaces(surface, position_rnd) + for air_brake in self.air_brakes: + air_brake = self._create_air_brake(air_brake) + rocket.set_air_brakes( + air_brakes=air_brake, + controller=self.air_brake_controller, + ) + for component_rail_buttons in self.rail_buttons: ( rail_buttons, From b2066478521064c0d3ffd83552b4610c596e1488 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 12:36:33 +0100 Subject: [PATCH 09/37] BUG: StochasticAirBrake object input in _Controller When defining the _Controller object a StochasticAirBrake was input. This is already corrected and a AirBrake object is now introduced --- rocketpy/stochastic/stochastic_rocket.py | 25 +++++++++--------------- 1 file changed, 9 insertions(+), 16 deletions(-) diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index c8f91b342..399ea568d 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -396,11 +396,7 @@ def set_rail_buttons( def add_air_brakes( self, air_brakes, - controller_function, - sampling_rate, - initial_observed_variables=None, - return_controller=False, - controller_name="AirBrakes Controller", + controller ): # checks if input is a StochasticAirbrakes type @@ -412,16 +408,7 @@ def add_air_brakes( air_brakes = StochasticAirBrakes(air_brakes=air_brakes) self.air_brakes.append(air_brakes) - _controller = _Controller( - interactive_objects=air_brakes, - controller_function=controller_function, - sampling_rate=sampling_rate, - initial_observed_variables=initial_observed_variables, - name=controller_name, - ) - self.air_brake_controller = _controller - if return_controller: - return _controller + self.air_brake_controller = controller def _validate_position(self, validated_object, position): """Validate the position argument. @@ -659,9 +646,15 @@ def create_object(self): for air_brake in self.air_brakes: air_brake = self._create_air_brake(air_brake) + _controller = _Controller( + interactive_objects=air_brake, + controller_function=self.air_brake_controller.base_controller_function, + sampling_rate=self.air_brake_controller.sampling_rate, + initial_observed_variables=self.air_brake_controller.initial_observed_variables, + ) rocket.set_air_brakes( air_brakes=air_brake, - controller=self.air_brake_controller, + controller=_controller, ) for component_rail_buttons in self.rail_buttons: From e45ae76bf1b8fa4db78cd2c8c818ee1f4fd2d9f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 13:47:43 +0100 Subject: [PATCH 10/37] ENH: add time_overshoot option to rocketpy.stochastic_flight Since the new StochasticAirBrake class is defined, we need the 'time_overshoot' option in the Flight class to ensure that the time step defined in the simulation is the controller sampling rate. The MonteCarlo class has had to be modified as well to include this option. --- rocketpy/simulation/monte_carlo.py | 1 + rocketpy/stochastic/stochastic_flight.py | 13 +++++++++++++ 2 files changed, 14 insertions(+) diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index 3914f1e54..ff86db7b3 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -458,6 +458,7 @@ def __run_single_simulation(self): heading=self.flight._randomize_heading(), initial_solution=self.flight.initial_solution, terminate_on_apogee=self.flight.terminate_on_apogee, + time_overshoot=self.flight.time_overshoot, ) def __evaluate_flight_inputs(self, sim_idx): diff --git a/rocketpy/stochastic/stochastic_flight.py b/rocketpy/stochastic/stochastic_flight.py index 4ed778eb5..ed9037088 100644 --- a/rocketpy/stochastic/stochastic_flight.py +++ b/rocketpy/stochastic/stochastic_flight.py @@ -29,6 +29,9 @@ class StochasticFlight(StochasticModel): terminate_on_apogee : bool Whether or not the flight should terminate on apogee. This attribute can not be randomized. + time_overshoot : bool + If False, the simulation will run at the time step defined by the controller + sampling rate. Be aware that this will make the simulation run much slower. """ def __init__( @@ -39,6 +42,7 @@ def __init__( heading=None, initial_solution=None, terminate_on_apogee=None, + time_overshoot=None, ): """Initializes the Stochastic Flight class. @@ -63,11 +67,18 @@ def __init__( terminate_on_apogee : bool, optional Whether or not the flight should terminate on apogee. This attribute can not be randomized. + time_overshoot : bool + If False, the simulation will run at the time step defined by the controller + sampling rate. Be aware that this will make the simulation run much slower. """ if terminate_on_apogee is not None: assert isinstance(terminate_on_apogee, bool), ( "`terminate_on_apogee` must be a boolean" ) + if time_overshoot is not None: + assert isinstance(time_overshoot, bool), ( + "`time_overshoot` must be a boolean" + ) super().__init__( flight, rail_length=rail_length, @@ -77,6 +88,7 @@ def __init__( self.initial_solution = initial_solution self.terminate_on_apogee = terminate_on_apogee + self.time_overshoot = time_overshoot def _validate_initial_solution(self, initial_solution): if initial_solution is not None: @@ -128,4 +140,5 @@ def create_object(self): heading=generated_dict["heading"], initial_solution=self.initial_solution, terminate_on_apogee=self.terminate_on_apogee, + time_overshoot=self.time_overshoot, ) From 7f082d2cbcba9a8342b4e5324b99f86b999bd68e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 17:35:52 +0100 Subject: [PATCH 11/37] DOC: StochasticAirBrakes related documentation added Documentation related to the StochasticAirBrakes implementation has been added in StochasticAirBrakes, StochasticRocket and Rocket classes. --- rocketpy/rocket/rocket.py | 16 +++++++++++-- .../stochastic/stochastic_aero_surfaces.py | 24 +++++++++++++++++-- rocketpy/stochastic/stochastic_rocket.py | 10 +++++++- 3 files changed, 45 insertions(+), 5 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index e40b71a9c..520747f77 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1539,9 +1539,21 @@ def set_air_brakes( air_brakes, controller, ): - + """Adds an air brake with a controller to the rocket. + + Parameters + ---------- + air_brakes : AirBrakes + Air brake to be added to the rocket. + controller : _Controller + Air brake controller + + Returns + ------- + None + """ self.air_brakes.append(air_brakes) - self._add_controllers(controller) + self._add_controllers(controller) def add_air_brakes( self, diff --git a/rocketpy/stochastic/stochastic_aero_surfaces.py b/rocketpy/stochastic/stochastic_aero_surfaces.py index 17f22337b..4241130c8 100644 --- a/rocketpy/stochastic/stochastic_aero_surfaces.py +++ b/rocketpy/stochastic/stochastic_aero_surfaces.py @@ -447,7 +447,7 @@ class StochasticAirBrakes(StochasticModel): ---------- object : AirBrakes AirBrakes object to be used for validation. - drag_coefficient_curve : list + drag_coefficient_curve : list, str The drag coefficient curve of the air brakes can account for either the air brakes' drag alone or the combined drag of both the rocket and the air brakes. @@ -455,6 +455,16 @@ class StochasticAirBrakes(StochasticModel): The drag curve factor of the air brakes. reference_area : tuple, list, int, float Reference area used to non-dimensionalize the drag coefficients. + clamp : bool + If True, the simulation will clamp the deployment level to 0 or 1 if + the deployment level is out of bounds. If False, the simulation will + not clamp the deployment level and will instead raise a warning if + the deployment level is out of bounds. + override_rocket_drag : bool + If False, the air brakes drag coefficient will be added to the + rocket's power off drag coefficient curve. If True, during the + simulation, the rocket's power off drag will be ignored and the air + brakes drag coefficient will be used for the entire rocket instead. deployment_level : tuple, list, int, float Initial deployment level, ranging from 0 to 1. name : list[str] @@ -481,7 +491,7 @@ def __init__( ---------- air_brakes : AirBrakes AirBrakes object to be used for validation. - drag_coefficient_curve : list, optional + drag_coefficient_curve : list, str, optional The drag coefficient curve of the air brakes can account for either the air brakes' drag alone or the combined drag of both the rocket and the air brakes. @@ -489,6 +499,16 @@ def __init__( The drag curve factor of the air brakes. reference_area : tuple, list, int, float, optional Reference area used to non-dimensionalize the drag coefficients. + clamp : bool, optional + If True, the simulation will clamp the deployment level to 0 or 1 if + the deployment level is out of bounds. If False, the simulation will + not clamp the deployment level and will instead raise a warning if + the deployment level is out of bounds. + override_rocket_drag : bool, optional + If False, the air brakes drag coefficient will be added to the + rocket's power off drag coefficient curve. If True, during the + simulation, the rocket's power off drag will be ignored and the air + brakes drag coefficient will be used for the entire rocket instead. deployment_level : tuple, list, int, float, optional Initial deployment level, ranging from 0 to 1. """ diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 399ea568d..3e29e5963 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -398,7 +398,15 @@ def add_air_brakes( air_brakes, controller ): - + """Adds an air brake to the stochastic rocket. + + Parameters + ---------- + air_brakes : StochasticAirBrakes or Airbrakes + The air brake to be added to the stochastic rocket. + controller : _Controller + Deterministic air brake controller. + """ # checks if input is a StochasticAirbrakes type if not isinstance(air_brakes, (AirBrakes, StochasticAirBrakes)): raise TypeError( From 1c4c79122958e9fc994cae904fa10efc28c8cfea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 23:00:00 +0100 Subject: [PATCH 12/37] ENH: pylint recommendations done --- rocketpy/mathutils/function.py | 2 +- rocketpy/stochastic/stochastic_aero_surfaces.py | 4 ++-- rocketpy/stochastic/stochastic_rocket.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 869f5a7b9..209544b10 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2219,7 +2219,7 @@ def __rsub__(self, other): """ return other + (-self) - def __mul__(self, other): + def __mul__(self, other): # pylint: disable=too-many-statements """Multiplies a Function object and returns a new Function object which gives the result of the multiplication. Only implemented for 1D and 2D domains. diff --git a/rocketpy/stochastic/stochastic_aero_surfaces.py b/rocketpy/stochastic/stochastic_aero_surfaces.py index 4241130c8..e9745ac01 100644 --- a/rocketpy/stochastic/stochastic_aero_surfaces.py +++ b/rocketpy/stochastic/stochastic_aero_surfaces.py @@ -433,7 +433,7 @@ def create_object(self): """ generated_dict = next(self.dict_generator()) return RailButtons(**generated_dict) - + class StochasticAirBrakes(StochasticModel): """A Stochastic Air Brakes class that inherits from StochasticModel. @@ -541,4 +541,4 @@ def create_object(self): deployment_level=generated_dict['deployment_level'], ) air_brakes.drag_coefficient *= generated_dict['drag_coefficient_curve_factor'] - return air_brakes + return air_brakes diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 3e29e5963..6dcc7963a 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -602,7 +602,7 @@ def _create_rail_buttons(self, component_stochastic_rail_buttons): upper_button_position_rnd ) return rail_buttons, lower_button_position_rnd, upper_button_position_rnd - + def _create_air_brake(self, stochastic_air_brake): air_brake = stochastic_air_brake.create_object() self.last_rnd_dict["air_brakes"].append(stochastic_air_brake.last_rnd_dict) From 85e82e639931e45b6bac1d35484a93e64368d6dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 18 Mar 2025 23:11:33 +0100 Subject: [PATCH 13/37] ENH: Reformatted files to pass Ruff linting checks --- rocketpy/mathutils/function.py | 20 ++++++++++++++----- .../stochastic/stochastic_aero_surfaces.py | 18 ++++++++--------- rocketpy/stochastic/stochastic_rocket.py | 6 +----- 3 files changed, 25 insertions(+), 19 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 209544b10..aa717b2c3 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2264,7 +2264,9 @@ def __mul__(self, other): # pylint: disable=too-many-statements ): if not self_source_is_array: return Function(lambda x: (self.get_value_opt(x) * other), inputs) - source = np.column_stack((self.x_array, np.multiply(self.y_array, other))) + source = np.column_stack( + (self.x_array, np.multiply(self.y_array, other)) + ) outputs = f"({self.__outputs__[0]}*{other})" return Function( source, @@ -2284,15 +2286,21 @@ def __mul__(self, other): # pylint: disable=too-many-statements and np.array_equal(self.x_array, other.x_array) and np.array_equal(self.y_array, other.y_array) ): - source = np.column_stack((self.x_array, self.y_array, self.z_array * other.z_array)) + source = np.column_stack( + (self.x_array, self.y_array, self.z_array * other.z_array) + ) outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})" return Function(source, inputs, outputs, interp, extrap) elif isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array( other ): if not self_source_is_array: - return Function(lambda x,y: (self.get_value_opt(x,y) * other), inputs) - source = np.column_stack((self.x_array, self.y_array, np.multiply(self.z_array, other))) + return Function( + lambda x, y: (self.get_value_opt(x, y) * other), inputs + ) + source = np.column_stack( + (self.x_array, self.y_array, np.multiply(self.z_array, other)) + ) outputs = f"({self.__outputs__[0]}*{other})" return Function( source, @@ -2302,7 +2310,9 @@ def __mul__(self, other): # pylint: disable=too-many-statements extrap, ) elif callable(other): - return Function(lambda x,y: (self.get_value_opt(x,y) * other(x)), inputs) + return Function( + lambda x, y: (self.get_value_opt(x, y) * other(x)), inputs + ) else: raise TypeError("Unsupported type for multiplication") diff --git a/rocketpy/stochastic/stochastic_aero_surfaces.py b/rocketpy/stochastic/stochastic_aero_surfaces.py index e9745ac01..863c97d64 100644 --- a/rocketpy/stochastic/stochastic_aero_surfaces.py +++ b/rocketpy/stochastic/stochastic_aero_surfaces.py @@ -448,7 +448,7 @@ class StochasticAirBrakes(StochasticModel): object : AirBrakes AirBrakes object to be used for validation. drag_coefficient_curve : list, str - The drag coefficient curve of the air brakes can account for + The drag coefficient curve of the air brakes can account for either the air brakes' drag alone or the combined drag of both the rocket and the air brakes. drag_coefficient_curve_factor : tuple, list, int, float @@ -479,7 +479,7 @@ def __init__( reference_area=None, clamp=None, override_rocket_drag=None, - deployment_level=(0,0), + deployment_level=(0, 0), ): """Initializes the Stochastic AirBrakes class. @@ -492,7 +492,7 @@ def __init__( air_brakes : AirBrakes AirBrakes object to be used for validation. drag_coefficient_curve : list, str, optional - The drag coefficient curve of the air brakes can account for + The drag coefficient curve of the air brakes can account for either the air brakes' drag alone or the combined drag of both the rocket and the air brakes. drag_coefficient_curve_factor : tuple, list, int, float, optional @@ -534,11 +534,11 @@ def create_object(self): """ generated_dict = next(self.dict_generator()) air_brakes = AirBrakes( - drag_coefficient_curve=generated_dict['drag_coefficient_curve'], - reference_area=generated_dict['reference_area'], - clamp=generated_dict['clamp'], - override_rocket_drag=generated_dict['override_rocket_drag'], - deployment_level=generated_dict['deployment_level'], + drag_coefficient_curve=generated_dict["drag_coefficient_curve"], + reference_area=generated_dict["reference_area"], + clamp=generated_dict["clamp"], + override_rocket_drag=generated_dict["override_rocket_drag"], + deployment_level=generated_dict["deployment_level"], ) - air_brakes.drag_coefficient *= generated_dict['drag_coefficient_curve_factor'] + air_brakes.drag_coefficient *= generated_dict["drag_coefficient_curve_factor"] return air_brakes diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 6dcc7963a..3a298dd61 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -393,11 +393,7 @@ def set_rail_buttons( rail_buttons, self._validate_position(rail_buttons, lower_button_position) ) - def add_air_brakes( - self, - air_brakes, - controller - ): + def add_air_brakes(self, air_brakes, controller): """Adds an air brake to the stochastic rocket. Parameters From a52ad3e3f3dc43ce3f87e69d6f75e21ce98b9ace Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Wed, 19 Mar 2025 21:52:56 +0100 Subject: [PATCH 14/37] ENH: Update rocketpy/stochastic/stochastic_rocket.py Unnecessary comment Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/stochastic/stochastic_rocket.py | 1 - 1 file changed, 1 deletion(-) diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 3a298dd61..67dd4cc8d 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -403,7 +403,6 @@ def add_air_brakes(self, air_brakes, controller): controller : _Controller Deterministic air brake controller. """ - # checks if input is a StochasticAirbrakes type if not isinstance(air_brakes, (AirBrakes, StochasticAirBrakes)): raise TypeError( "`air_brake` must be of AirBrakes or StochasticAirBrakes type" From 7c68241adf11fe0d8b08b8047eb3c9f511dd242d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sat, 22 Mar 2025 10:06:13 +0100 Subject: [PATCH 15/37] ENH: more intuitive uniform distribution display in StochasticModel Co-authored-by: MateusStano <69485049+MateusStano@users.noreply.github.com> --- rocketpy/stochastic/stochastic_model.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index b5be5293e..fe47a252a 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -515,12 +515,11 @@ def format_attribute(attr, value): elif isinstance(value, tuple): nominal_value, std_dev, dist_func = value if callable(dist_func) and dist_func.__name__ == "uniform": - mean = (std_dev + nominal_value) / 2 - half_range = (std_dev - nominal_value) / 2 + lower_bound = nominal_value + upper_bound = std_dev return ( f"\t{attr.ljust(max_str_length)} " - f"{mean:.5f} ± " - f"{half_range:.5f} ({dist_func.__name__})" + f"{lower_bound:.5f}, {upper_bound:.5f} ({dist_func.__name__})" ) else: return ( From 91f83c89fed3d90b7b5993f5ac8cb124fcb8bfc9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sun, 23 Mar 2025 11:38:14 +0100 Subject: [PATCH 16/37] DOC: improve drag curve factor definition in StochasticAirBrakes --- rocketpy/stochastic/stochastic_aero_surfaces.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/rocketpy/stochastic/stochastic_aero_surfaces.py b/rocketpy/stochastic/stochastic_aero_surfaces.py index 863c97d64..5dda716bb 100644 --- a/rocketpy/stochastic/stochastic_aero_surfaces.py +++ b/rocketpy/stochastic/stochastic_aero_surfaces.py @@ -452,7 +452,8 @@ class StochasticAirBrakes(StochasticModel): either the air brakes' drag alone or the combined drag of both the rocket and the air brakes. drag_coefficient_curve_factor : tuple, list, int, float - The drag curve factor of the air brakes. + The drag curve factor of the air brakes. This value scales the + drag coefficient curve to introduce stochastic variability. reference_area : tuple, list, int, float Reference area used to non-dimensionalize the drag coefficients. clamp : bool @@ -496,7 +497,8 @@ def __init__( either the air brakes' drag alone or the combined drag of both the rocket and the air brakes. drag_coefficient_curve_factor : tuple, list, int, float, optional - The drag curve factor of the air brakes. + The drag curve factor of the air brakes. This value scales the + drag coefficient curve to introduce stochastic variability. reference_area : tuple, list, int, float, optional Reference area used to non-dimensionalize the drag coefficients. clamp : bool, optional From 95a69faa7996dc9006fd73fd2d4542910647b252 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sun, 23 Mar 2025 11:52:11 +0100 Subject: [PATCH 17/37] ENH: Change assert statement to if Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/stochastic/stochastic_flight.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/rocketpy/stochastic/stochastic_flight.py b/rocketpy/stochastic/stochastic_flight.py index ed9037088..729313736 100644 --- a/rocketpy/stochastic/stochastic_flight.py +++ b/rocketpy/stochastic/stochastic_flight.py @@ -76,9 +76,8 @@ def __init__( "`terminate_on_apogee` must be a boolean" ) if time_overshoot is not None: - assert isinstance(time_overshoot, bool), ( - "`time_overshoot` must be a boolean" - ) + if not isinstance(time_overshoot, bool): + raise TypeError("`time_overshoot` must be a boolean") super().__init__( flight, rail_length=rail_length, From ac958d637db35acf6a30c5167a3f2f73ed67ae8f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sun, 23 Mar 2025 19:38:52 +0100 Subject: [PATCH 18/37] DOC: better explanation of __mul__ function Co-authored-by: MateusStano <69485049+MateusStano@users.noreply.github.com> --- rocketpy/mathutils/function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index aa717b2c3..c7d11910b 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2238,7 +2238,7 @@ def __mul__(self, other): # pylint: disable=too-many-statements Returns ------- result : Function - A Function object which gives the result of self(x,y)*other(x,y). + A Function object which gives the result of self(x)*other(x) or self(x,y)*other(x,y). """ self_source_is_array = isinstance(self.source, np.ndarray) other_source_is_array = ( From d51f15b3cb163eb79881d2a4fbaedf7711edc05e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Sun, 23 Mar 2025 20:12:30 +0100 Subject: [PATCH 19/37] ENH: delete set_air_brakes function for simplicity --- rocketpy/rocket/rocket.py | 21 --------------------- rocketpy/stochastic/stochastic_rocket.py | 6 ++---- 2 files changed, 2 insertions(+), 25 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 520747f77..609e1b76c 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1534,27 +1534,6 @@ def add_sensor(self, sensor, position): except KeyError: sensor._attached_rockets[self] = 1 - def set_air_brakes( - self, - air_brakes, - controller, - ): - """Adds an air brake with a controller to the rocket. - - Parameters - ---------- - air_brakes : AirBrakes - Air brake to be added to the rocket. - controller : _Controller - Air brake controller - - Returns - ------- - None - """ - self.air_brakes.append(air_brakes) - self._add_controllers(controller) - def add_air_brakes( self, drag_coefficient_curve, diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 67dd4cc8d..718176e36 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/rocketpy/stochastic/stochastic_rocket.py @@ -655,10 +655,8 @@ def create_object(self): sampling_rate=self.air_brake_controller.sampling_rate, initial_observed_variables=self.air_brake_controller.initial_observed_variables, ) - rocket.set_air_brakes( - air_brakes=air_brake, - controller=_controller, - ) + rocket.air_brakes.append(air_brake) + rocket._add_controllers(_controller) for component_rail_buttons in self.rail_buttons: ( From 2f42a50fe14a706c9f7af5bb7e0e8eb176e57c34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 25 Mar 2025 08:58:50 +0100 Subject: [PATCH 20/37] ENH: inertial foreces added to u_dot_generalized --- rocketpy/simulation/flight.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 9d53ae2d6..292b5d1c0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1740,10 +1740,10 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. """ # Retrieve integration data - _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u + x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Create necessary vectors - # r = Vector([x, y, z]) # CDM position vector + r = Vector([x, y, z]) # CDM position vector v = Vector([vx, vy, vz]) # CDM velocity vector e = [e0, e1, e2, e3] # Euler parameters/quaternions w = Vector([omega1, omega2, omega3]) # Angular velocity vector @@ -1896,8 +1896,11 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too # Angular velocity derivative w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) - # Velocity vector derivative - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) + # Velocity vector derivative + Inertial forces + wE = Kt @ self.env.earth_rotation_vector + v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - ( + 2*(wE ^ v) + wE ^ (wE ^ r) + ) # Euler parameters derivative e_dot = [ From 299decc9fd12e18a91f76cffdc2bb6c741fcb31b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 25 Mar 2025 10:02:16 +0100 Subject: [PATCH 21/37] ENH: define Earth's angular velocity vector in Environment --- rocketpy/environment/environment.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 1a0c044ed..48da835b3 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -353,6 +353,7 @@ def __init__( self.set_location(latitude, longitude) self.__initialize_earth_geometry(datum) self.__initialize_utm_coordinates() + self.__set_earth_rotation_vector() # Set the gravity model self.gravity = self.set_gravity_model(gravity) @@ -583,6 +584,23 @@ def __reset_wind_direction_function(self): self.wind_direction.set_inputs("Height Above Sea Level (m)") self.wind_direction.set_outputs("Wind Direction (Deg True)") self.wind_direction.set_title("Wind Direction Profile") + + def __set_earth_rotation_vector(self): + """Calculates and stores the Earth's angular velocity vector in the local + reference frame, which is essential for evaluating inertial forces. + """ + # Sidereal day + T = 86164.1 # seconds + + # Earth's angular velocity magnitude + wE = 2*np.pi/T + + # Vector in the local reference frame + lat = np.radians(self.latitude) + wL = [0, wE * np.cos(lat), wE * np.sin(lat)] + + # Store the attribute + self.earth_rotation_vector = wL # Validators (used to verify an attribute is being set correctly.) From 54aeab8a97adfb8e6ef8d444d7341e4f48d86fb4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 25 Mar 2025 10:33:18 +0100 Subject: [PATCH 22/37] ENH: some corrections to the Flight class --- rocketpy/simulation/flight.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 292b5d1c0..b6d6e2ce0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1743,7 +1743,7 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Create necessary vectors - r = Vector([x, y, z]) # CDM position vector + r = Vector([x, y, z]) # CDM position vector v = Vector([vx, vy, vz]) # CDM velocity vector e = [e0, e1, e2, e3] # Euler parameters/quaternions w = Vector([omega1, omega2, omega3]) # Angular velocity vector @@ -1897,7 +1897,7 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) # Velocity vector derivative + Inertial forces - wE = Kt @ self.env.earth_rotation_vector + wE = Kt @ Vector(self.env.earth_rotation_vector) v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - ( 2*(wE ^ v) + wE ^ (wE ^ r) ) From 653f0a65d155eacaba875960c096608132df2da0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 25 Mar 2025 12:24:00 +0100 Subject: [PATCH 23/37] ENH: modifications in the Flight class --- rocketpy/simulation/flight.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index b6d6e2ce0..e61a56504 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1740,10 +1740,10 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. """ # Retrieve integration data - x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u + _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Create necessary vectors - r = Vector([x, y, z]) # CDM position vector + # r = Vector([x, y, z]) # CDM position vector v = Vector([vx, vy, vz]) # CDM velocity vector e = [e0, e1, e2, e3] # Euler parameters/quaternions w = Vector([omega1, omega2, omega3]) # Angular velocity vector @@ -1896,11 +1896,9 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too # Angular velocity derivative w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) - # Velocity vector derivative + Inertial forces + # Velocity vector derivative + Coriolis acceleration wE = Kt @ Vector(self.env.earth_rotation_vector) - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - ( - 2*(wE ^ v) + wE ^ (wE ^ r) - ) + v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2*(wE ^ v) # Euler parameters derivative e_dot = [ From a1387465b8d85a07a6a70913c0a523a7537d272d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 25 Mar 2025 14:32:06 +0100 Subject: [PATCH 24/37] DOC: improving Environment documentation --- rocketpy/environment/environment.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 48da835b3..35ee2d35e 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -586,8 +586,8 @@ def __reset_wind_direction_function(self): self.wind_direction.set_title("Wind Direction Profile") def __set_earth_rotation_vector(self): - """Calculates and stores the Earth's angular velocity vector in the local - reference frame, which is essential for evaluating inertial forces. + """Calculates and stores the Earth's angular velocity vector in the Flight + Coordinate System, which is essential for evaluating inertial forces. """ # Sidereal day T = 86164.1 # seconds @@ -595,7 +595,7 @@ def __set_earth_rotation_vector(self): # Earth's angular velocity magnitude wE = 2*np.pi/T - # Vector in the local reference frame + # Vector in the Flight Coordinate System lat = np.radians(self.latitude) wL = [0, wE * np.cos(lat), wE * np.sin(lat)] From 41977ba5d5fa5e70d1a6ef41b6d6d4cbfb154cb4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Tue, 25 Mar 2025 14:36:16 +0100 Subject: [PATCH 25/37] DOC: more improvements in the Environment class --- rocketpy/environment/environment.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 35ee2d35e..8e2cb777a 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -248,6 +248,8 @@ class Environment: Number of ensemble members. Only defined when using Ensembles. Environment.ensemble_member : int Current selected ensemble member. Only defined when using Ensembles. + Environment.earth_rotation_vector : list + Earth's angular velocity vector in the Flight Coordinate System. Notes ----- From 38cad32da028b1e7870e66ef3c9d228acb21ea0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Mon, 31 Mar 2025 08:12:18 +0200 Subject: [PATCH 26/37] ENH: format changes done --- rocketpy/environment/environment.py | 10 +++++----- rocketpy/simulation/flight.py | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 8e2cb777a..4b1e4e182 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -586,23 +586,23 @@ def __reset_wind_direction_function(self): self.wind_direction.set_inputs("Height Above Sea Level (m)") self.wind_direction.set_outputs("Wind Direction (Deg True)") self.wind_direction.set_title("Wind Direction Profile") - + def __set_earth_rotation_vector(self): """Calculates and stores the Earth's angular velocity vector in the Flight Coordinate System, which is essential for evaluating inertial forces. """ # Sidereal day - T = 86164.1 # seconds + T = 86164.1 # seconds # Earth's angular velocity magnitude - wE = 2*np.pi/T + w_earth = 2 * np.pi / T # Vector in the Flight Coordinate System lat = np.radians(self.latitude) - wL = [0, wE * np.cos(lat), wE * np.sin(lat)] + w_local = [0, w_earth * np.cos(lat), w_earth * np.sin(lat)] # Store the attribute - self.earth_rotation_vector = wL + self.earth_rotation_vector = w_local # Validators (used to verify an attribute is being set correctly.) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index e61a56504..717dfa350 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1897,8 +1897,8 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) # Velocity vector derivative + Coriolis acceleration - wE = Kt @ Vector(self.env.earth_rotation_vector) - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2*(wE ^ v) + w_earth = Kt @ Vector(self.env.earth_rotation_vector) + v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v) # Euler parameters derivative e_dot = [ From 3d482213594f0d042ee7369f6203069e5dd6b3e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Fri, 4 Apr 2025 20:22:58 +0200 Subject: [PATCH 27/37] DOC: env.earth_rotation_vector improved Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/environment/environment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 4b1e4e182..d3fbfb60a 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -248,7 +248,7 @@ class Environment: Number of ensemble members. Only defined when using Ensembles. Environment.ensemble_member : int Current selected ensemble member. Only defined when using Ensembles. - Environment.earth_rotation_vector : list + Environment.earth_rotation_vector : list[float] Earth's angular velocity vector in the Flight Coordinate System. Notes From 057e7fe9b5ac4e812095b6b1ffbcbab228552cab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Mon, 7 Apr 2025 12:29:10 +0200 Subject: [PATCH 28/37] ENH: Coriolis acceleration added to u_dot --- rocketpy/simulation/flight.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 717dfa350..67a08155d 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1692,6 +1692,13 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals ax, ay, az = K @ Vector(L) az -= self.env.gravity.get_value_opt(z) # Include gravity + # Coriolis acceleration + print('u_dot') + _, w_earth_y, w_earth_z = self.env.earth_rotation_vector + ax -= 2 * (-a23*vy*w_earth_y + a22*vz*w_earth_y - a33*vy*w_earth_z + a32*vz*w_earth_z) + ay -= 2 * (a23*vx*w_earth_y - a21*vz*w_earth_y + a33*vx*w_earth_z - a31*vz*w_earth_z) + az -= 2 * (-a22*vx*w_earth_y + a21*vy*w_earth_y - a32*vx*w_earth_z + a31*vy*w_earth_z) + # Create u_dot u_dot = [ vx, From 6a6cc77b4c88a124aaff0f53856182e781772d6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Mon, 7 Apr 2025 12:30:41 +0200 Subject: [PATCH 29/37] BUG: print left --- rocketpy/simulation/flight.py | 1 - 1 file changed, 1 deletion(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 67a08155d..1e98a8def 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1693,7 +1693,6 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals az -= self.env.gravity.get_value_opt(z) # Include gravity # Coriolis acceleration - print('u_dot') _, w_earth_y, w_earth_z = self.env.earth_rotation_vector ax -= 2 * (-a23*vy*w_earth_y + a22*vz*w_earth_y - a33*vy*w_earth_z + a32*vz*w_earth_z) ay -= 2 * (a23*vx*w_earth_y - a21*vz*w_earth_y + a33*vx*w_earth_z - a31*vz*w_earth_z) From 79ed52f37aeb153d7ebc5b6c0b8d74f5c75465a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Mon, 7 Apr 2025 12:32:29 +0200 Subject: [PATCH 30/37] ENH: ruff changes --- rocketpy/simulation/flight.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 1e98a8def..11f1f7455 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1694,9 +1694,24 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Coriolis acceleration _, w_earth_y, w_earth_z = self.env.earth_rotation_vector - ax -= 2 * (-a23*vy*w_earth_y + a22*vz*w_earth_y - a33*vy*w_earth_z + a32*vz*w_earth_z) - ay -= 2 * (a23*vx*w_earth_y - a21*vz*w_earth_y + a33*vx*w_earth_z - a31*vz*w_earth_z) - az -= 2 * (-a22*vx*w_earth_y + a21*vy*w_earth_y - a32*vx*w_earth_z + a31*vy*w_earth_z) + ax -= 2 * ( + -a23 * vy * w_earth_y + + a22 * vz * w_earth_y + - a33 * vy * w_earth_z + + a32 * vz * w_earth_z + ) + ay -= 2 * ( + a23 * vx * w_earth_y + - a21 * vz * w_earth_y + + a33 * vx * w_earth_z + - a31 * vz * w_earth_z + ) + az -= 2 * ( + -a22 * vx * w_earth_y + + a21 * vy * w_earth_y + - a32 * vx * w_earth_z + + a31 * vy * w_earth_z + ) # Create u_dot u_dot = [ From 01d981b4d26c92ec1308fefa62a22b092bf10d31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kevin=20Alca=C3=B1iz?= Date: Mon, 7 Apr 2025 12:41:27 +0200 Subject: [PATCH 31/37] ENH: CHANGELOG updated --- CHANGELOG.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a0d664c8..7a7d1fc5f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,19 @@ Attention: The newest changes should be on top --> ### Added +- ENH: Support for the RSE file format has been added to the library [#798](https://github.com/RocketPy-Team/RocketPy/pull/798) +- ENH: Add the Coriolis Force to the Flight class [#799](https://github.com/RocketPy-Team/RocketPy/pull/799) + +### Changed + + +### Fixed + + +## v1.9.0 - 2025-03-24 + +### Added + - ENH: Parallel mode for monte-carlo simulations 2 [#768](https://github.com/RocketPy-Team/RocketPy/pull/768) - DOC: ASTRA Flight Example [#770](https://github.com/RocketPy-Team/RocketPy/pull/770) - ENH: Add Eccentricity to Stochastic Simulations [#792](https://github.com/RocketPy-Team/RocketPy/pull/792) @@ -39,6 +52,7 @@ Attention: The newest changes should be on top --> ### Changed +- DEP: Remove Pending Deprecations and Add Warnings Where Needed [#794](https://github.com/RocketPy-Team/RocketPy/pull/794) - DOCS: reshape docs (closes #659) [#781](https://github.com/RocketPy-Team/RocketPy/pull/781) - MNT: EmptyMotor class inherits from Motor(ABC) [#779](https://github.com/RocketPy-Team/RocketPy/pull/779) From 950265cd2cc7b9b4cd9d5674f09d8941b0489904 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 13 Jun 2025 18:01:24 +0200 Subject: [PATCH 32/37] ENH: remove unecessary frame rotation --- rocketpy/simulation/flight.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 203505b32..808151438 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1925,10 +1925,6 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too # Angular velocity derivative w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) - # Velocity vector derivative + Coriolis acceleration - w_earth = Kt @ Vector(self.env.earth_rotation_vector) - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v) - # Euler parameters derivative e_dot = [ 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3), @@ -1937,6 +1933,10 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2), ] + # Velocity vector derivative + Coriolis acceleration + w_earth = Vector(self.env.earth_rotation_vector) + v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v) + # Position vector derivative r_dot = [vx, vy, vz] From 07d1a6f6a5514a00d03dde76f71dbcf12566735d Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 13 Jun 2025 19:24:08 +0200 Subject: [PATCH 33/37] ENH: remove rotation from solid prop udot --- rocketpy/simulation/flight.py | 21 +++------------------ 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 808151438..d3256c297 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1681,24 +1681,9 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Coriolis acceleration _, w_earth_y, w_earth_z = self.env.earth_rotation_vector - ax -= 2 * ( - -a23 * vy * w_earth_y - + a22 * vz * w_earth_y - - a33 * vy * w_earth_z - + a32 * vz * w_earth_z - ) - ay -= 2 * ( - a23 * vx * w_earth_y - - a21 * vz * w_earth_y - + a33 * vx * w_earth_z - - a31 * vz * w_earth_z - ) - az -= 2 * ( - -a22 * vx * w_earth_y - + a21 * vy * w_earth_y - - a32 * vx * w_earth_z - + a31 * vy * w_earth_z - ) + ax -= 2 * (vz * w_earth_y - vy * w_earth_z) + ay -= 2 * (vx * w_earth_z) + az -= 2 * (-vx * w_earth_y) # Create u_dot u_dot = [ From d985eabd1f9837d105680eb5f047a01865f89af0 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 13 Jun 2025 19:24:43 +0200 Subject: [PATCH 34/37] ENH: add coriolis to parachute --- rocketpy/simulation/flight.py | 220 ++++++++++++++++++++++++++++++++++ 1 file changed, 220 insertions(+) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index d3256c297..454531e34 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1197,6 +1197,10 @@ def __init_equations_of_motion(self): if self.equations_of_motion == "solid_propulsion": # NOTE: The u_dot is faster, but only works for solid propulsion self.u_dot_generalized = self.u_dot + if self.equations_of_motion == "2": + self.u_dot_generalized = self.u_dot_generalized_std + if self.equations_of_motion == "3": + self.u_dot_generalized = self.u_dot_generalized_moving_omega def __init_controllers(self): """Initialize controllers and sensors""" @@ -1935,6 +1939,216 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too return u_dot + def u_dot_generalized_std(self, t, u, post_processing=False): # pylint: disable=too-many-locals,too-many-statements + """Calculates derivative of u state vector with respect to time when the + rocket is flying in 6 DOF motion in space and significant mass variation + effects exist. Typical flight phases include powered ascent after launch + rail. + + Parameters + ---------- + t : float + Time in seconds + u : list + State vector defined by u = [x, y, z, vx, vy, vz, q0, q1, + q2, q3, omega1, omega2, omega3]. + post_processing : bool, optional + If True, adds flight data information directly to self variables + such as self.angle_of_attack, by default False. + + Returns + ------- + u_dot : list + State vector defined by u_dot = [vx, vy, vz, ax, ay, az, + e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. + """ + # Retrieve integration data + _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u + + # Create necessary vectors + # r = Vector([x, y, z]) # CDM position vector + v = Vector([vx, vy, vz]) # CDM velocity vector + e = [e0, e1, e2, e3] # Euler parameters/quaternions + w = Vector([omega1, omega2, omega3]) # Angular velocity vector + + # Retrieve necessary quantities + ## Rocket mass + total_mass = self.rocket.total_mass.get_value_opt(t) + total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) + total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t) + ## CM position vector and time derivatives relative to CDM in body frame + r_CM_z = self.rocket.com_to_cdm_function + r_CM_t = r_CM_z.get_value_opt(t) + r_CM = Vector([0, 0, r_CM_t]) + r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) + r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)]) + ## Nozzle position vector + r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm]) + ## Nozzle gyration tensor + S_nozzle = self.rocket.nozzle_gyration_tensor + ## Inertia tensor + inertia_tensor = self.rocket.get_inertia_tensor_at_time(t) + ## Inertia tensor time derivative in the body frame + I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) + + # Calculate the Inertia tensor relative to CM + H = (r_CM.cross_matrix @ -r_CM.cross_matrix) * total_mass + I_CM = inertia_tensor - H + + # Prepare transformation matrices + K = Matrix.transformation(e) + Kt = K.transpose + + # Compute aerodynamic forces and moments + R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0 + + ## Drag force + rho = self.env.density.get_value_opt(z) + wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) + wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) + wind_velocity = Vector([wind_velocity_x, wind_velocity_y, 0]) + free_stream_speed = abs((wind_velocity - Vector(v))) + speed_of_sound = self.env.speed_of_sound.get_value_opt(z) + free_stream_mach = free_stream_speed / speed_of_sound + + if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time: + pressure = self.env.pressure.get_value_opt(z) + net_thrust = max( + self.rocket.motor.thrust.get_value_opt(t) + + self.rocket.motor.pressure_thrust(pressure), + 0, + ) + drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach) + else: + net_thrust = 0 + drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach) + R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff + for air_brakes in self.rocket.air_brakes: + if air_brakes.deployment_level > 0: + air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( + air_brakes.deployment_level, free_stream_mach + ) + air_brakes_force = ( + -0.5 + * rho + * (free_stream_speed**2) + * air_brakes.reference_area + * air_brakes_cd + ) + if air_brakes.override_rocket_drag: + R3 = air_brakes_force # Substitutes rocket drag coefficient + else: + R3 += air_brakes_force + # Get rocket velocity in body frame + velocity_in_body_frame = Kt @ v + # Calculate lift and moment for each component of the rocket + for aero_surface, _ in self.rocket.aerodynamic_surfaces: + # Component cp relative to CDM in body frame + comp_cp = self.rocket.surfaces_cp_to_cdm[aero_surface] + # Component absolute velocity in body frame + comp_vb = velocity_in_body_frame + (w ^ comp_cp) + # Wind velocity at component altitude + comp_z = z + (K @ comp_cp).z + comp_wind_vx = self.env.wind_velocity_x.get_value_opt(comp_z) + comp_wind_vy = self.env.wind_velocity_y.get_value_opt(comp_z) + # Component freestream velocity in body frame + comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0]) + comp_stream_velocity = comp_wind_vb - comp_vb + comp_stream_speed = abs(comp_stream_velocity) + comp_stream_mach = comp_stream_speed / speed_of_sound + # Reynolds at component altitude + # TODO: Reynolds is only used in generic surfaces. This calculation + # should be moved to the surface class for efficiency + comp_reynolds = ( + self.env.density.get_value_opt(comp_z) + * comp_stream_speed + * aero_surface.reference_length + / self.env.dynamic_viscosity.get_value_opt(comp_z) + ) + # Forces and moments + X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( + comp_stream_velocity, + comp_stream_speed, + comp_stream_mach, + rho, + comp_cp, + w, + comp_reynolds, + ) + R1 += X + R2 += Y + R3 += Z + M1 += M + M2 += N + M3 += L + + # Off center moment + M1 += ( + self.rocket.cp_eccentricity_y * R3 + + self.rocket.thrust_eccentricity_y * net_thrust + ) + M2 -= ( + self.rocket.cp_eccentricity_x * R3 + + self.rocket.thrust_eccentricity_x * net_thrust + ) + M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1 + + weight_in_body_frame = Kt @ Vector( + [0, 0, -total_mass * self.env.gravity.get_value_opt(z)] + ) + + T00 = total_mass * r_CM + T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot + T04 = ( + Vector([0, 0, net_thrust]) + - total_mass * r_CM_ddot + - 2 * total_mass_dot * r_CM_dot + + total_mass_ddot * (r_NOZ - r_CM) + ) + T05 = total_mass_dot * S_nozzle - I_dot + + T20 = ( + ((w ^ T00) ^ w) + + (w ^ T03) + + T04 + + weight_in_body_frame + + Vector([R1, R2, R3]) + ) + + T21 = ( + ((inertia_tensor @ w) ^ w) + + T05 @ w + - (weight_in_body_frame ^ r_CM) + + Vector([M1, M2, M3]) + ) + + # Angular velocity derivative + w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) + + # Velocity vector derivative + Coriolis acceleration + v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) + + # Euler parameters derivative + e_dot = [ + 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3), + 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3), + 0.5 * (omega2 * e0 - omega3 * e1 + omega1 * e3), + 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2), + ] + + # Position vector derivative + r_dot = [vx, vy, vz] + + # Create u_dot + u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] + + if post_processing: + self.__post_processed_variables.append( + [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3, net_thrust] + ) + + return u_dot + def u_dot_parachute(self, t, u, post_processing=False): """Calculates derivative of u state vector with respect to time when rocket is flying under parachute. A 3 DOF approximation is @@ -2001,6 +2215,12 @@ def u_dot_parachute(self, t, u, post_processing=False): ay = Dy / (mp + ma) az = (Dz - 9.8 * mp) / (mp + ma) + # Add coriolis acceleration + _, w_earth_y, w_earth_z = self.env.earth_rotation_vector + ax -= 2 * (vz * w_earth_y - vy * w_earth_z) + ay -= 2 * (vx * w_earth_z) + az -= 2 * (-vx * w_earth_y) + if post_processing: self.__post_processed_variables.append( [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0, 0] From ef38d40af13b72902afc756d86f8116cfa9f2f97 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 13 Jun 2025 20:24:38 +0200 Subject: [PATCH 35/37] TST: fix tests values --- tests/integration/test_flight.py | 18 +++++++++++++----- tests/unit/test_flight.py | 28 ++++++++++++++-------------- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/tests/integration/test_flight.py b/tests/integration/test_flight.py index 8fddb6486..c47b9b124 100644 --- a/tests/integration/test_flight.py +++ b/tests/integration/test_flight.py @@ -491,8 +491,9 @@ def test_empty_motor_flight(mock_show, example_plain_env, calisto_motorless): # def test_freestream_speed_at_apogee(example_plain_env, calisto): """ - Asserts that a rocket at apogee has a free stream speed of 0.0 m/s in all - directions given that the environment doesn't have any wind. + Asserts that a rocket at apogee has a free stream speed of near 0.0 m/s + in all directions given that the environment doesn't have any wind. Any + speed values comes from coriolis effect. """ # NOTE: this rocket doesn't move in x or z direction. There's no wind. hard_atol = 1e-12 @@ -508,7 +509,9 @@ def test_freestream_speed_at_apogee(example_plain_env, calisto): ) assert np.isclose( - test_flight.stream_velocity_x(test_flight.apogee_time), 0.0, atol=hard_atol + test_flight.stream_velocity_x(test_flight.apogee_time), + 0.4641492104717301, + atol=hard_atol, ) assert np.isclose( test_flight.stream_velocity_y(test_flight.apogee_time), 0.0, atol=hard_atol @@ -518,9 +521,13 @@ def test_freestream_speed_at_apogee(example_plain_env, calisto): test_flight.stream_velocity_z(test_flight.apogee_time), 0.0, atol=soft_atol ) assert np.isclose( - test_flight.free_stream_speed(test_flight.apogee_time), 0.0, atol=soft_atol + test_flight.free_stream_speed(test_flight.apogee_time), + 0.4641492104717798, + atol=hard_atol, + ) + assert np.isclose( + test_flight.apogee_freestream_speed, 0.4641492104717798, atol=hard_atol ) - assert np.isclose(test_flight.apogee_freestream_speed, 0.0, atol=soft_atol) def test_rocket_csys_equivalence( @@ -546,6 +553,7 @@ def test_rocket_csys_equivalence( assert np.isclose( flight_calisto_robust.x_impact, flight_calisto_nose_to_tail_robust.x_impact, + atol=1e-3, ) assert np.isclose( flight_calisto_robust.y_impact, diff --git a/tests/unit/test_flight.py b/tests/unit/test_flight.py index 260a2b138..04bd97c40 100644 --- a/tests/unit/test_flight.py +++ b/tests/unit/test_flight.py @@ -109,15 +109,15 @@ def test_get_solution_at_time(flight_calisto): flight_calisto.get_solution_at_time(flight_calisto.t_final), np.array( [ - 48.4313533, - 0.0, - 985.7665845, - -0.00000229951048, - 0.0, - 11.2223284, - -341.028803, - 0.999048222, - -0.0436193874, + 48.43719482805657, + -14.836008075478597, + 985.9858934483618, + -3.4415459237894554e-05, + 0.0007572309307800201, + 11.21695000766671, + -341.1460775169661, + 0.9990482215818578, + -0.043619387365336, 0.0, 0.0, 0.0, @@ -212,7 +212,7 @@ def test_export_sensor_data(flight_calisto_with_sensors): [ ("t_initial", (0.25886, -0.649623, 0)), ("out_of_rail_time", (0.792028, -1.987634, 0)), - ("apogee_time", (-0.522875, -0.741825, 0)), + ("apogee_time", (-0.509420, -0.732933, -2.089120e-14)), ("t_final", (0, 0, 0)), ], ) @@ -251,8 +251,8 @@ def test_aerodynamic_moments(flight_calisto_custom_wind, flight_time, expected_v [ ("t_initial", (1.654150, 0.659142, -0.067103)), ("out_of_rail_time", (5.052628, 2.013361, -1.75370)), - ("apogee_time", (2.339424, -1.648934, -0.938867)), - ("t_final", (0, 0, 159.2210)), + ("apogee_time", (2.321838, -1.613641, -0.962108)), + ("t_final", (-0.025792, 0.012030, 159.202481)), ], ) def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values): @@ -292,7 +292,7 @@ def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_va ("out_of_rail_time", (0, 2.248540, 25.700928)), ( "apogee_time", - (-14.488364, 15.638049, -0.000191), + (-14.826350, 15.670022, -0.000264), ), ("t_final", (5, 2, -5.660155)), ], @@ -540,7 +540,7 @@ def test_lat_lon_conversion_from_origin(mock_show, example_plain_env, calisto_ro heading=0, ) - assert abs(test_flight.longitude(test_flight.t_final) - 0) < 1e-12 + assert abs(test_flight.longitude(test_flight.t_final)) < 1e-4 assert test_flight.latitude(test_flight.t_final) > 0 From 5ca314dc2615e85052a4518c349be413372f4826 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 13 Jun 2025 20:30:18 +0200 Subject: [PATCH 36/37] MNT: remove debug functions --- rocketpy/simulation/flight.py | 214 ---------------------------------- 1 file changed, 214 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 454531e34..c6c75ad5d 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1197,10 +1197,6 @@ def __init_equations_of_motion(self): if self.equations_of_motion == "solid_propulsion": # NOTE: The u_dot is faster, but only works for solid propulsion self.u_dot_generalized = self.u_dot - if self.equations_of_motion == "2": - self.u_dot_generalized = self.u_dot_generalized_std - if self.equations_of_motion == "3": - self.u_dot_generalized = self.u_dot_generalized_moving_omega def __init_controllers(self): """Initialize controllers and sensors""" @@ -1939,216 +1935,6 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too return u_dot - def u_dot_generalized_std(self, t, u, post_processing=False): # pylint: disable=too-many-locals,too-many-statements - """Calculates derivative of u state vector with respect to time when the - rocket is flying in 6 DOF motion in space and significant mass variation - effects exist. Typical flight phases include powered ascent after launch - rail. - - Parameters - ---------- - t : float - Time in seconds - u : list - State vector defined by u = [x, y, z, vx, vy, vz, q0, q1, - q2, q3, omega1, omega2, omega3]. - post_processing : bool, optional - If True, adds flight data information directly to self variables - such as self.angle_of_attack, by default False. - - Returns - ------- - u_dot : list - State vector defined by u_dot = [vx, vy, vz, ax, ay, az, - e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. - """ - # Retrieve integration data - _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u - - # Create necessary vectors - # r = Vector([x, y, z]) # CDM position vector - v = Vector([vx, vy, vz]) # CDM velocity vector - e = [e0, e1, e2, e3] # Euler parameters/quaternions - w = Vector([omega1, omega2, omega3]) # Angular velocity vector - - # Retrieve necessary quantities - ## Rocket mass - total_mass = self.rocket.total_mass.get_value_opt(t) - total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) - total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t) - ## CM position vector and time derivatives relative to CDM in body frame - r_CM_z = self.rocket.com_to_cdm_function - r_CM_t = r_CM_z.get_value_opt(t) - r_CM = Vector([0, 0, r_CM_t]) - r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) - r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)]) - ## Nozzle position vector - r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm]) - ## Nozzle gyration tensor - S_nozzle = self.rocket.nozzle_gyration_tensor - ## Inertia tensor - inertia_tensor = self.rocket.get_inertia_tensor_at_time(t) - ## Inertia tensor time derivative in the body frame - I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) - - # Calculate the Inertia tensor relative to CM - H = (r_CM.cross_matrix @ -r_CM.cross_matrix) * total_mass - I_CM = inertia_tensor - H - - # Prepare transformation matrices - K = Matrix.transformation(e) - Kt = K.transpose - - # Compute aerodynamic forces and moments - R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0 - - ## Drag force - rho = self.env.density.get_value_opt(z) - wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) - wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) - wind_velocity = Vector([wind_velocity_x, wind_velocity_y, 0]) - free_stream_speed = abs((wind_velocity - Vector(v))) - speed_of_sound = self.env.speed_of_sound.get_value_opt(z) - free_stream_mach = free_stream_speed / speed_of_sound - - if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time: - pressure = self.env.pressure.get_value_opt(z) - net_thrust = max( - self.rocket.motor.thrust.get_value_opt(t) - + self.rocket.motor.pressure_thrust(pressure), - 0, - ) - drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach) - else: - net_thrust = 0 - drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach) - R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff - for air_brakes in self.rocket.air_brakes: - if air_brakes.deployment_level > 0: - air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( - air_brakes.deployment_level, free_stream_mach - ) - air_brakes_force = ( - -0.5 - * rho - * (free_stream_speed**2) - * air_brakes.reference_area - * air_brakes_cd - ) - if air_brakes.override_rocket_drag: - R3 = air_brakes_force # Substitutes rocket drag coefficient - else: - R3 += air_brakes_force - # Get rocket velocity in body frame - velocity_in_body_frame = Kt @ v - # Calculate lift and moment for each component of the rocket - for aero_surface, _ in self.rocket.aerodynamic_surfaces: - # Component cp relative to CDM in body frame - comp_cp = self.rocket.surfaces_cp_to_cdm[aero_surface] - # Component absolute velocity in body frame - comp_vb = velocity_in_body_frame + (w ^ comp_cp) - # Wind velocity at component altitude - comp_z = z + (K @ comp_cp).z - comp_wind_vx = self.env.wind_velocity_x.get_value_opt(comp_z) - comp_wind_vy = self.env.wind_velocity_y.get_value_opt(comp_z) - # Component freestream velocity in body frame - comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0]) - comp_stream_velocity = comp_wind_vb - comp_vb - comp_stream_speed = abs(comp_stream_velocity) - comp_stream_mach = comp_stream_speed / speed_of_sound - # Reynolds at component altitude - # TODO: Reynolds is only used in generic surfaces. This calculation - # should be moved to the surface class for efficiency - comp_reynolds = ( - self.env.density.get_value_opt(comp_z) - * comp_stream_speed - * aero_surface.reference_length - / self.env.dynamic_viscosity.get_value_opt(comp_z) - ) - # Forces and moments - X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( - comp_stream_velocity, - comp_stream_speed, - comp_stream_mach, - rho, - comp_cp, - w, - comp_reynolds, - ) - R1 += X - R2 += Y - R3 += Z - M1 += M - M2 += N - M3 += L - - # Off center moment - M1 += ( - self.rocket.cp_eccentricity_y * R3 - + self.rocket.thrust_eccentricity_y * net_thrust - ) - M2 -= ( - self.rocket.cp_eccentricity_x * R3 - + self.rocket.thrust_eccentricity_x * net_thrust - ) - M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1 - - weight_in_body_frame = Kt @ Vector( - [0, 0, -total_mass * self.env.gravity.get_value_opt(z)] - ) - - T00 = total_mass * r_CM - T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot - T04 = ( - Vector([0, 0, net_thrust]) - - total_mass * r_CM_ddot - - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (r_NOZ - r_CM) - ) - T05 = total_mass_dot * S_nozzle - I_dot - - T20 = ( - ((w ^ T00) ^ w) - + (w ^ T03) - + T04 - + weight_in_body_frame - + Vector([R1, R2, R3]) - ) - - T21 = ( - ((inertia_tensor @ w) ^ w) - + T05 @ w - - (weight_in_body_frame ^ r_CM) - + Vector([M1, M2, M3]) - ) - - # Angular velocity derivative - w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) - - # Velocity vector derivative + Coriolis acceleration - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - - # Euler parameters derivative - e_dot = [ - 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3), - 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3), - 0.5 * (omega2 * e0 - omega3 * e1 + omega1 * e3), - 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2), - ] - - # Position vector derivative - r_dot = [vx, vy, vz] - - # Create u_dot - u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] - - if post_processing: - self.__post_processed_variables.append( - [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3, net_thrust] - ) - - return u_dot - def u_dot_parachute(self, t, u, post_processing=False): """Calculates derivative of u state vector with respect to time when rocket is flying under parachute. A 3 DOF approximation is From ff82200f659872d057c71cec7dc3b85c342a666e Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 4 Jul 2025 11:56:28 +0200 Subject: [PATCH 37/37] DEV: changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index dce368bc8..f78b34edb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,11 +32,11 @@ Attention: The newest changes should be on top --> ### Added +- ENH: Add the Coriolis Force to the Flight class [#799](https://github.com/RocketPy-Team/RocketPy/pull/799) ### Changed - ENH: _MotorPrints inheritance - issue #460 [#828](https://github.com/RocketPy-Team/RocketPy/pull/828) - - MNT: fix deprecations and warnings [#829](https://github.com/RocketPy-Team/RocketPy/pull/829) ### Fixed