From 830ea15445d9f30cf520981388e07557fe2f8da3 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 15 Apr 2024 15:05:31 +0200 Subject: [PATCH 01/28] ENH: refactor retrieve arrays cached properties and controller sim post processing --- rocketpy/simulation/flight.py | 253 +++++++++++++++++++++------------- 1 file changed, 156 insertions(+), 97 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 77632d2f2..ef1b54bbe 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -593,7 +593,6 @@ def __init__( if self.rail_length <= 0: raise ValueError("Rail length must be a positive value.") self.parachutes = self.rocket.parachutes[:] - self._controllers = self.rocket._controllers[:] self.inclination = inclination self.heading = heading self.max_time = max_time @@ -607,6 +606,13 @@ def __init__( self.name = name self.equations_of_motion = equations_of_motion + # Controller initialization + self._controllers = self.rocket._controllers[:] + if self._controllers: + # reset controllable object to initial state (only airbrakes for now) + for air_brakes in self.rocket.air_brakes: + air_brakes.reset() + # Flight initialization self.__init_post_process_variables() self.__init_solution_monitors() @@ -1076,8 +1082,14 @@ def __init__( [self.t, parachute] ) + # If controlled flight, post process must be done on sim time + if self._controllers: + phase.derivative(self.t, self.y_sol, post_processing=True) + self.t_final = self.t self._calculate_pressure_signal() + if self._controllers: + self.__cache_post_process_variables() if verbose: print("Simulation Completed at Time: {:3.4f} s".format(self.t)) @@ -1089,6 +1101,25 @@ def __init_post_process_variables(self): self._bearing = Function(0) self._latitude = Function(0) self._longitude = Function(0) + # Initialize state derivatives, force and atmospheric arrays + self.ax_list = [] + self.ay_list = [] + self.az_list = [] + self.alpha1_list = [] + self.alpha2_list = [] + self.alpha3_list = [] + self.R1_list = [] + self.R2_list = [] + self.R3_list = [] + self.M1_list = [] + self.M2_list = [] + self.M3_list = [] + self.pressure_list = [] + self.density_list = [] + self.dynamic_viscosity_list = [] + self.speed_of_sound_list = [] + self.wind_velocity_x_list = [] + self.wind_velocity_y_list = [] def __init_solution_monitors(self): # Initialize solution monitors @@ -1181,10 +1212,28 @@ def __init_equations_of_motion(self): if self.equations_of_motion == "solid_propulsion": self.u_dot_generalized = self.u_dot - def __init_equations_of_motion(self): - """Initialize equations of motion.""" - if self.equations_of_motion == "solid_propulsion": - self.u_dot_generalized = self.u_dot + def __cache_post_process_variables(self): + """Cache post-process variables for simulations with controllers.""" + self.__retrieve_arrays = [ + self.ax_list, + self.ay_list, + self.az_list, + self.alpha1_list, + self.alpha2_list, + self.alpha3_list, + self.R1_list, + self.R2_list, + self.R3_list, + self.M1_list, + self.M2_list, + self.M3_list, + self.pressure_list, + self.density_list, + self.dynamic_viscosity_list, + self.speed_of_sound_list, + self.wind_velocity_x_list, + self.wind_velocity_y_list, + ] @cached_property def effective_1rl(self): @@ -1261,11 +1310,6 @@ def udot_rail1(self, t, u, post_processing=False): e0dot, e1dot, e2dot, e3dot, alpha1, alpha2, alpha3]. """ - # Check if post processing mode is on - if post_processing: - # Use u_dot post processing code - return self.u_dot_generalized(t, u, True) - # Retrieve integration data x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u @@ -1296,6 +1340,17 @@ def udot_rail1(self, t, u, post_processing=False): else: ax, ay, az = 0, 0, 0 + if post_processing: + # Use u_dot post processing code for forces, moments and env data + self.u_dot_generalized(t, u, post_processing=True) + # Save feasible accelerations + self.ax_list[-1] = [t, ax] + self.ay_list[-1] = [t, ay] + self.az_list[-1] = [t, az] + self.alpha1_list[-1] = [t, 0] + self.alpha2_list[-1] = [t, 0] + self.alpha3_list[-1] = [t, 0] + return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] def udot_rail2(self, t, u, post_processing=False): @@ -1585,6 +1640,13 @@ def u_dot(self, t, u, post_processing=False): ] if post_processing: + # Accelerations + self.ax_list.append([t, ax]) + self.ay_list.append([t, ay]) + self.az_list.append([t, az]) + self.alpha1_list.append([t, alpha1]) + self.alpha2_list.append([t, alpha2]) + self.alpha3_list.append([t, alpha3]) # Dynamics variables self.R1_list.append([t, R1]) self.R2_list.append([t, R2]) @@ -1860,6 +1922,13 @@ def u_dot_generalized(self, t, u, post_processing=False): u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] if post_processing: + # Accelerations + self.ax_list.append([t, v_dot[0]]) + self.ay_list.append([t, v_dot[1]]) + self.az_list.append([t, v_dot[2]]) + self.alpha1_list.append([t, w_dot[0]]) + self.alpha2_list.append([t, w_dot[1]]) + self.alpha3_list.append([t, w_dot[2]]) # Dynamics variables self.R1_list.append([t, R1]) self.R2_list.append([t, R2]) @@ -1944,6 +2013,13 @@ def u_dot_parachute(self, t, u, post_processing=False): az = (Dz - 9.8 * mp) / (mp + ma) if post_processing: + # Accelerations + self.ax_list.append([t, ax]) + self.ay_list.append([t, ay]) + self.az_list.append([t, az]) + self.alpha1_list.append([t, 0]) + self.alpha2_list.append([t, 0]) + self.alpha3_list.append([t, 0]) # Dynamics variables self.R1_list.append([t, Dx]) self.R2_list.append([t, Dy]) @@ -1952,13 +2028,20 @@ def u_dot_parachute(self, t, u, post_processing=False): self.M2_list.append([t, 0]) self.M3_list.append([t, 0]) # Atmospheric Conditions - self.wind_velocity_x_list.append([t, self.env.wind_velocity_x(z)]) - self.wind_velocity_y_list.append([t, self.env.wind_velocity_y(z)]) - self.density_list.append([t, self.env.density(z)]) - self.dynamic_viscosity_list.append([t, self.env.dynamic_viscosity(z)]) - self.pressure_list.append([t, self.env.pressure(z)]) - self.speed_of_sound_list.append([t, self.env.speed_of_sound(z)]) - + self.wind_velocity_x_list.append( + [t, self.env.wind_velocity_x.get_value_opt(z)] + ) + self.wind_velocity_y_list.append( + [t, self.env.wind_velocity_y.get_value_opt(z)] + ) + self.density_list.append([t, self.env.density.get_value_opt(z)]) + self.dynamic_viscosity_list.append( + [t, self.env.dynamic_viscosity.get_value_opt(z)] + ) + self.pressure_list.append([t, self.env.pressure.get_value_opt(z)]) + self.speed_of_sound_list.append( + [t, self.env.speed_of_sound.get_value_opt(z)] + ) return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] @cached_property @@ -2790,13 +2873,64 @@ def longitude(self): return np.column_stack((self.time, longitude)) + @cached_property + def __retrieve_arrays(self): + """post processing function to retrieve arrays from the integration + scheme and store them in lists for further analysis. + + Returns + ------- + temp_values: list + List containing the following arrays: ``ax`` , ``ay`` , ``az`` , + ``alpha1`` , ``alpha2`` , ``alpha3`` , ``R1`` , ``R2`` , ``R3`` , + ``M1`` , ``M2`` , ``M3`` , ``pressure`` , ``density`` , + ``dynamic_viscosity`` , ``speed_of_sound`` , ``wind_velocity_x`` , + ``wind_velocity_y``. + """ + # Go through each time step and calculate forces and atmospheric values + # Get flight phases + for phase_index, phase in self.time_iterator(self.FlightPhases): + init_time = phase.t + final_time = self.FlightPhases[phase_index + 1].t + current_derivative = phase.derivative + # Call callback functions + for callback in phase.callbacks: + callback(self) + # Loop through time steps in flight phase + for step in self.solution: # Can be optimized + if init_time < step[0] <= final_time or ( + init_time == self.t_initial and step[0] == self.t_initial + ): + # Call derivatives in post processing mode + current_derivative(step[0], step[1:], post_processing=True) + + temp_values = [ + self.ax_list, + self.ay_list, + self.az_list, + self.alpha1_list, + self.alpha2_list, + self.alpha3_list, + self.R1_list, + self.R2_list, + self.R3_list, + self.M1_list, + self.M2_list, + self.M3_list, + self.pressure_list, + self.density_list, + self.dynamic_viscosity_list, + self.speed_of_sound_list, + self.wind_velocity_x_list, + self.wind_velocity_y_list, + ] + + return temp_values + @cached_property def retrieve_acceleration_arrays(self): """Retrieve acceleration arrays from the integration scheme - Parameters - ---------- - Returns ------- ax: list @@ -2812,35 +2946,7 @@ def retrieve_acceleration_arrays(self): alpha3: list angular acceleration in z direction """ - # Initialize acceleration arrays - ax, ay, az = [[0, 0]], [[0, 0]], [[0, 0]] - alpha1, alpha2, alpha3 = [[0, 0]], [[0, 0]], [[0, 0]] - # Go through each time step and calculate accelerations - # Get flight phases - for phase_index, phase in self.time_iterator(self.FlightPhases): - init_time = phase.t - final_time = self.FlightPhases[phase_index + 1].t - current_derivative = phase.derivative - # Call callback functions - for callback in phase.callbacks: - callback(self) - # Loop through time steps in flight phase - for step in self.solution: # Can be optimized - if init_time < step[0] <= final_time: - # Get derivatives - u_dot = current_derivative(step[0], step[1:]) - # Get accelerations - ax_value, ay_value, az_value = u_dot[3:6] - alpha1_value, alpha2_value, alpha3_value = u_dot[10:] - # Save accelerations - ax.append([step[0], ax_value]) - ay.append([step[0], ay_value]) - az.append([step[0], az_value]) - alpha1.append([step[0], alpha1_value]) - alpha2.append([step[0], alpha2_value]) - alpha3.append([step[0], alpha3_value]) - - return ax, ay, az, alpha1, alpha2, alpha3 + return self.__retrieve_arrays[:6] @cached_property def retrieve_temporary_values_arrays(self): @@ -2877,54 +2983,7 @@ def retrieve_temporary_values_arrays(self): self.wind_velocity_y_list: list Wind velocity in y direction at each time step. """ - - # Initialize force and atmospheric arrays - self.R1_list = [] - self.R2_list = [] - self.R3_list = [] - self.M1_list = [] - self.M2_list = [] - self.M3_list = [] - self.pressure_list = [] - self.density_list = [] - self.dynamic_viscosity_list = [] - self.speed_of_sound_list = [] - self.wind_velocity_x_list = [] - self.wind_velocity_y_list = [] - - # Go through each time step and calculate forces and atmospheric values - # Get flight phases - for phase_index, phase in self.time_iterator(self.FlightPhases): - init_time = phase.t - final_time = self.FlightPhases[phase_index + 1].t - current_derivative = phase.derivative - # Call callback functions - for callback in phase.callbacks: - callback(self) - # Loop through time steps in flight phase - for step in self.solution: # Can be optimized - if init_time < step[0] <= final_time or ( - init_time == self.t_initial and step[0] == self.t_initial - ): - # Call derivatives in post processing mode - u_dot = current_derivative(step[0], step[1:], post_processing=True) - - temporary_values = [ - self.R1_list, - self.R2_list, - self.R3_list, - self.M1_list, - self.M2_list, - self.M3_list, - self.pressure_list, - self.density_list, - self.dynamic_viscosity_list, - self.speed_of_sound_list, - self.wind_velocity_x_list, - self.wind_velocity_y_list, - ] - - return temporary_values + return self.__retrieve_arrays[6:] def get_controller_observed_variables(self): """Retrieve the observed variables related to air brakes from the From b9a97a454880a3ffbb7194c9807304dfbea0899e Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 15 Apr 2024 15:06:01 +0200 Subject: [PATCH 02/28] ENH: add warning for improper time_overshoot --- rocketpy/simulation/flight.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index ef1b54bbe..178fdf4a0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -609,6 +609,13 @@ def __init__( # Controller initialization self._controllers = self.rocket._controllers[:] if self._controllers: + if self.time_overshoot == True: + warnings.warn( + "time_overshoot is set to True, but controllers are present. " + "Controllers will not work properly. " + "Consider setting time_overshoot=False in the Flight object " + "to use controllers." + ) # reset controllable object to initial state (only airbrakes for now) for air_brakes in self.rocket.air_brakes: air_brakes.reset() From a66fc53323248c74db32874f8dff8ee8e5cd4289 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 15 Apr 2024 15:18:37 +0200 Subject: [PATCH 03/28] ENH: add reset function --- rocketpy/rocket/aero_surface.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py index c5d154f3e..27f87cd5c 100644 --- a/rocketpy/rocket/aero_surface.py +++ b/rocketpy/rocket/aero_surface.py @@ -1978,8 +1978,9 @@ def __init__( brakes drag coefficient will be used for the entire rocket instead. Default is False. deployment_level : float, optional - Current deployment level, ranging from 0 to 1. Deployment level is the - fraction of the total airbrake area that is Deployment. Default is 0. + Initial deployment level, ranging from 0 to 1. Deployment level is + the fraction of the total airbrake area that is Deployment. Default + is 0. name : str, optional Name of the air brakes. Default is "AirBrakes". @@ -1997,6 +1998,7 @@ def __init__( self.reference_area = reference_area self.clamp = clamp self.override_rocket_drag = override_rocket_drag + self.initial_deployment_level = deployment_level self.deployment_level = deployment_level self.prints = _AirBrakesPrints(self) self.plots = _AirBrakesPlots(self) @@ -2023,6 +2025,12 @@ def deployment_level(self, value): ) self._deployment_level = value + def _reset(self): + """Resets the air brakes to their initial state. This is ran at the + beginning of each simulation to ensure the air brakes are in the correct + state.""" + self.deployment_level = self.initial_deployment_level + def evaluate_center_of_pressure(self): """Evaluates the center of pressure of the aerodynamic surface in local coordinates. From 7fe3d0ba3b805f4a134eb399be796f7cb2e20425 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 15 Apr 2024 16:12:39 +0200 Subject: [PATCH 04/28] ENH: optmize post process loop --- 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 178fdf4a0..9a81b81c0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -618,7 +618,7 @@ def __init__( ) # reset controllable object to initial state (only airbrakes for now) for air_brakes in self.rocket.air_brakes: - air_brakes.reset() + air_brakes._reset() # Flight initialization self.__init_post_process_variables() @@ -2903,9 +2903,12 @@ def __retrieve_arrays(self): # Call callback functions for callback in phase.callbacks: callback(self) - # Loop through time steps in flight phase - for step in self.solution: # Can be optimized - if init_time < step[0] <= final_time or ( + # find index of initial and final time of phase in solution array + init_time_index = find_closest(self.time, init_time) + final_time_index = find_closest(self.time, final_time) + 1 + # Loop through time steps solution array + for step in self.solution[init_time_index:final_time_index]: + if init_time != step[0] or ( init_time == self.t_initial and step[0] == self.t_initial ): # Call derivatives in post processing mode From aaa3f7b953d7450eec6056c723f83afa82a67bbb Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 15 Apr 2024 16:52:49 +0200 Subject: [PATCH 05/28] DOCS: change accelerations names in docs --- rocketpy/simulation/flight.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 9a81b81c0..94e7f7fcd 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -2943,17 +2943,17 @@ def retrieve_acceleration_arrays(self): Returns ------- - ax: list + ax_list: list acceleration in x direction - ay: list + ay_list: list acceleration in y direction - az: list + az_list: list acceleration in z direction - alpha1: list + alpha1_list: list angular acceleration in x direction - alpha2: list + alpha2_list: list angular acceleration in y direction - alpha3: list + alpha3_list: list angular acceleration in z direction """ return self.__retrieve_arrays[:6] From a49baf4d0160ca01282b101daa9444b3872b44f7 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 25 Apr 2024 22:52:13 +0200 Subject: [PATCH 06/28] MNT: move controller initialization to private method --- rocketpy/simulation/flight.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 94e7f7fcd..9655b9e01 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -607,18 +607,7 @@ def __init__( self.equations_of_motion = equations_of_motion # Controller initialization - self._controllers = self.rocket._controllers[:] - if self._controllers: - if self.time_overshoot == True: - warnings.warn( - "time_overshoot is set to True, but controllers are present. " - "Controllers will not work properly. " - "Consider setting time_overshoot=False in the Flight object " - "to use controllers." - ) - # reset controllable object to initial state (only airbrakes for now) - for air_brakes in self.rocket.air_brakes: - air_brakes._reset() + self.__init_controllers() # Flight initialization self.__init_post_process_variables() @@ -1219,6 +1208,21 @@ def __init_equations_of_motion(self): if self.equations_of_motion == "solid_propulsion": self.u_dot_generalized = self.u_dot + def __init_controllers(self): + """Initialize controllers""" + self._controllers = self.rocket._controllers[:] + if self._controllers: + if self.time_overshoot == True: + warnings.warn( + "time_overshoot is set to True, but controllers are present. " + "Controllers will not work properly. " + "Consider setting time_overshoot=False in the Flight object " + "to use controllers." + ) + # reset controllable object to initial state (only airbrakes for now) + for air_brakes in self.rocket.air_brakes: + air_brakes._reset() + def __cache_post_process_variables(self): """Cache post-process variables for simulations with controllers.""" self.__retrieve_arrays = [ From 4974e8d7417e0d102088d40007b3c429aed1b896 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 25 Apr 2024 22:52:48 +0200 Subject: [PATCH 07/28] BUG: initialize state derivatives for controllers --- rocketpy/simulation/flight.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 9655b9e01..5fcb7c86e 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1188,6 +1188,11 @@ def __init_flight_state(self): self.out_of_rail_time = self.initial_solution[0] self.out_of_rail_time_index = 0 self.initial_derivative = self.u_dot_generalized + if self._controllers: + # Handle post process during simulation, get initial accel/forces + self.initial_derivative( + self.t_initial, self.initial_solution[1:], post_processing=True + ) def __init_solver_monitors(self): # Initialize solver monitors From 67515f64be6f1702d6497055fbb071f888ac5c09 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 25 Apr 2024 23:02:11 +0200 Subject: [PATCH 08/28] ENH: automatically disable time overshoot for controllers --- rocketpy/simulation/flight.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 3ff87ed52..e90d5d9b8 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1218,11 +1218,9 @@ def __init_controllers(self): self._controllers = self.rocket._controllers[:] if self._controllers: if self.time_overshoot == True: + self.time_overshoot = False warnings.warn( - "time_overshoot is set to True, but controllers are present. " - "Controllers will not work properly. " - "Consider setting time_overshoot=False in the Flight object " - "to use controllers." + "time_overshoot has been set to False due to the presence of controllers. " ) # reset controllable object to initial state (only airbrakes for now) for air_brakes in self.rocket.air_brakes: From 04429a70d22dac3fbed01fdfd53b7d1fbe6f11ad Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:23:51 -0400 Subject: [PATCH 09/28] ENH: add Rocket.evaluate_nozzle_to_center_of_dry_mass_position() --- rocketpy/rocket/rocket.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index c65f1eca1..8d99ee24e 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -734,6 +734,21 @@ def evaluate_inertias(self): self.I_23, ) + def evaluate_nozzle_to_center_of_dry_mass_position(self): + """Evaluates the distance between the nozzle exit and the rocket's + center of dry mass position. + + Returns + ------- + self.nozzle_to_center_of_dry_mass_position : float + Distance between the nozzle exit and the rocket's center of dry + mass position, in meters. + """ + self.nozzle_to_center_of_dry_mass_position = ( + -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys + ) + return self.nozzle_to_center_of_dry_mass_position + def evaluate_nozzle_gyration_tensor(self): pass From 5ef3f04f858c93ccc3b90efa02f9c56d89f924bb Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:24:15 -0400 Subject: [PATCH 10/28] ENH: Add Rocket.evaluate_nozzle_gyration_tensor() --- rocketpy/rocket/rocket.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 8d99ee24e..0afd631c6 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -750,7 +750,28 @@ def evaluate_nozzle_to_center_of_dry_mass_position(self): return self.nozzle_to_center_of_dry_mass_position def evaluate_nozzle_gyration_tensor(self): - pass + """Calculates and returns the nozzle gyration tensor relative to the + rocket's center of dry mass. The gyration tensor is saved and returned + in units of kg*m². + + Returns + ------- + self.nozzle_gyration_tensor : Matrix + Matrix containing the nozzle gyration tensor. + """ + S_noz_33 = 0.5 * self.motor.nozzle_radius**2 + S_noz_11 = S_noz_22 = ( + 0.5 * S_noz_33 + 0.25 * self.nozzle_to_center_of_dry_mass_position**2 + ) + S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 # Due to axis symmetry + self.nozzle_gyration_tensor = Matrix( + [ + [S_noz_11, S_noz_12, S_noz_13], + [S_noz_12, S_noz_22, S_noz_23], + [S_noz_13, S_noz_23, S_noz_33], + ] + ) + return self.nozzle_gyration_tensor def add_motor(self, motor, position): """Adds a motor to the rocket. From cb27b004ce8aa3f1fe3640c979b75dfc73a52cdc Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:24:38 -0400 Subject: [PATCH 11/28] ENH: Adds Rocket.evaluate_z_coordinate_com_to_cdm() --- rocketpy/rocket/rocket.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 0afd631c6..3657c3fb9 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -773,6 +773,36 @@ def evaluate_nozzle_gyration_tensor(self): ) return self.nozzle_gyration_tensor + def evaluate_z_coordinate_com_to_cdm(self): + """Evaluates the z-coordinate of the center of mass (COM) relative to + the center of dry mass (CDM). + + Notes + ----- + 1. The `z_coordinate_com_to_cdm` plus `center_of_mass` should be equal + to `center_of_dry_mass_position` at every time step. + 2. The `z_coordinate_com_to_cdm` is a function of time and will usually + already be discretized. + + Returns + ------- + self.z_coordinate_com_to_cdm : Function + Function of time expressing the z-coordinate of the center of mass + relative to the center of dry mass. + """ + self.z_coordinate_com_to_cdm = ( + -1 + * ( + (self.center_of_propellant_position - self.center_of_dry_mass_position) + * self._csys + ) + * self.motor.propellant_mass + / self.total_mass + ) + self.z_coordinate_com_to_cdm.set_inputs("Time (s)") + self.z_coordinate_com_to_cdm.set_outputs("Z Coordinate COM to CDM (m)") + self.z_coordinate_com_to_cdm.set_title("Z Coordinate COM to CDM") + return self.z_coordinate_com_to_cdm def add_motor(self, motor, position): """Adds a motor to the rocket. From 4da4d19d52fb0697cb96993eb578e22e0c69804d Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:25:35 -0400 Subject: [PATCH 12/28] MNT: updates Rocket.add_motor() to accommodate new evaluations --- rocketpy/rocket/rocket.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 3657c3fb9..57a88c67c 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -844,6 +844,7 @@ def add_motor(self, motor, position): self.evaluate_dry_mass() self.evaluate_total_mass() self.evaluate_center_of_dry_mass() + self.evaluate_nozzle_to_center_of_dry_mass_position() self.evaluate_center_of_mass() self.evaluate_dry_inertias() self.evaluate_inertias() @@ -852,6 +853,8 @@ def add_motor(self, motor, position): self.evaluate_center_of_pressure() self.evaluate_stability_margin() self.evaluate_static_margin() + self.evaluate_z_coordinate_com_to_cdm() + self.evaluate_nozzle_gyration_tensor() def add_surfaces(self, surfaces, positions): """Adds one or more aerodynamic surfaces to the rocket. The aerodynamic From 1a00ca1f971e5157c1e093fb479ba73fd6fe16cd Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:26:19 -0400 Subject: [PATCH 13/28] ENH: Adds methods to get inertia tensors at specific time --- rocketpy/rocket/rocket.py | 61 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 57a88c67c..d3d00612a 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -803,6 +803,67 @@ def evaluate_z_coordinate_com_to_cdm(self): self.z_coordinate_com_to_cdm.set_outputs("Z Coordinate COM to CDM (m)") self.z_coordinate_com_to_cdm.set_title("Z Coordinate COM to CDM") return self.z_coordinate_com_to_cdm + + def get_inertia_tensor_at_time(self, t): + """Returns a Matrix representing the inertia tensor of the rocket with + respect to the rocket's center of mass at a given time. It evaluates + each inertia tensor component at the given time and returns a Matrix + with the computed values. + + Parameters + ---------- + t : float + Time at which the inertia tensor is to be evaluated. + + Returns + ------- + Matrix + Inertia tensor of the rocket at time t. + """ + I_11 = self.I_11.get_value_opt(t) + I_12 = self.I_12.get_value_opt(t) + I_13 = self.I_13.get_value_opt(t) + I_22 = self.I_22.get_value_opt(t) + I_23 = self.I_23.get_value_opt(t) + I_33 = self.I_33.get_value_opt(t) + return Matrix( + [ + [I_11, I_12, I_13], + [I_12, I_22, I_23], + [I_13, I_23, I_33], + ] + ) + + def get_inertia_tensor_derivative_at_time(self, t): + """Returns a Matrix representing the time derivative of the inertia + tensor of the rocket with respect to the rocket's center of mass at a + given time. It evaluates each inertia tensor component's derivative at + the given time and returns a Matrix with the computed values. + + Parameters + ---------- + t : float + Time at which the inertia tensor derivative is to be evaluated. + + Returns + ------- + Matrix + Inertia tensor time derivative of the rocket at time t. + """ + I_11_dot = self.I_11.differentiate_complex_step(t) + I_12_dot = self.I_12.differentiate_complex_step(t) + I_13_dot = self.I_13.differentiate_complex_step(t) + I_22_dot = self.I_22.differentiate_complex_step(t) + I_23_dot = self.I_23.differentiate_complex_step(t) + I_33_dot = self.I_33.differentiate_complex_step(t) + return Matrix( + [ + [I_11_dot, I_12_dot, I_13_dot], + [I_12_dot, I_22_dot, I_23_dot], + [I_13_dot, I_23_dot, I_33_dot], + ] + ) + def add_motor(self, motor, position): """Adds a motor to the rocket. From 74b38b7fc744d8068985e41fe0dcedd037e8e641 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:26:47 -0400 Subject: [PATCH 14/28] MNT: updates docstring of Rocket class --- rocketpy/rocket/rocket.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index d3d00612a..68bd026e8 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -4,6 +4,7 @@ from rocketpy.control.controller import _Controller from rocketpy.mathutils.function import Function +from rocketpy.mathutils.vector_matrix import Matrix from rocketpy.motors.motor import EmptyMotor from rocketpy.plots.rocket_plots import _RocketPlots from rocketpy.prints.rocket_prints import _RocketPrints @@ -72,6 +73,9 @@ class Rocket: for more information regarding the coordinate system. Expressed in meters as a function of time. + Rocket.z_coordinate_com_to_cdm : Function + Function of time expressing the z-coordinate of the center of mass + relative to the center of dry mass. Rocket.reduced_mass : Function Function of time expressing the reduced mass of the rocket, defined as the product of the propellant mass and the mass @@ -146,6 +150,11 @@ class Rocket: defined rocket coordinate system. See :doc:`Positions and Coordinate Systems ` for more information. + Rocket.nozzle_to_center_of_dry_mass_position : float + Distance between the nozzle exit and the rocket's center of dry mass + position, in meters. + Rocket.nozzle_gyration_tensor: Matrix + Matrix representing the nozzle gyration tensor. Rocket.center_of_propellant_position : Function Position of the propellant's center of mass relative to the user defined rocket reference system. See From fac8a023d77705d103a6da6ec4a825cbcddd304f Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:27:20 -0400 Subject: [PATCH 15/28] DOC: Adds Inertia tensors section to Rocket class page --- docs/user/rocket.rst | 51 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/docs/user/rocket.rst b/docs/user/rocket.rst index 28e091e8d..6b4a9f257 100644 --- a/docs/user/rocket.rst +++ b/docs/user/rocket.rst @@ -428,3 +428,54 @@ The lets check all the information available about the rocket: .. jupyter-execute:: calisto.all_info() + +7. Inertia Tensors +------------------ + +The inertia tensor of the rocket at a given time can be obtained using the +``get_inertia_tensor_at_time`` method. This method evaluates each component of +the inertia tensor at the specified time and returns a +:class:`rocketpy.mathutils.Matrix` object. + +The inertia tensor is a matrix that looks like this: + +.. math:: + :label: inertia_tensor + + \mathbf{I} = \begin{bmatrix} + I_{11} & I_{12} & I_{13} \\ + I_{21} & I_{22} & I_{23} \\ + I_{31} & I_{32} & I_{33} + \end{bmatrix} + +For example, to get the inertia tensor of the rocket at time 0.5 seconds, you +can use the following code: + +.. jupyter-execute:: + + calisto.get_inertia_tensor_at_time(0.5) + +Derivative of the Inertia Tensor +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +You can also get the derivative of the inertia tensor at a given time using the +``get_inertia_tensor_derivative_at_time`` method. Here's an example: + +.. jupyter-execute:: + + calisto.get_inertia_tensor_derivative_at_time(0.5) + +Implications from these results +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The inertia tensor reveals important information about the rocket's symmetry +and ease of rotation: + +1. **Axis Symmetry**: If I\ :sub:`11` and I\ :sub:`22` are equal, the rocket is symmetric around the axes perpendicular to the rocket's center axis. In our defined rocket, I\ :sub:`11` and I\ :sub:`22` are indeed equal, indicating that our rocket is axisymmetric. + +2. **Zero Products of Inertia**: The off-diagonal elements of the inertia tensor are zero, which means the products of inertia are zero. This indicates that the rocket is symmetric around its center axis. + +3. **Ease of Rotation**: The I\ :sub:`33` value is significantly lower than the other two. This suggests that the rocket is easier to rotate around its center axis than around the axes perpendicular to the rocket. This is an important factor when considering the rocket's stability and control. + +However, these conclusions are based on the assumption that the inertia tensor is calculated with respect to the rocket's center of mass and aligned with the principal axes of the rocket. If the inertia tensor is calculated with respect to a different point or not aligned with the principal axes, the conclusions may not hold. + From 9a2f842a06b40267142ff7ce9ea3c327e14dbdfd Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:29:53 -0400 Subject: [PATCH 16/28] ENH: updates Flight.u_dot_generalized to use pre-calculated attributes from Rocket class --- rocketpy/simulation/flight.py | 81 ++++++++--------------------------- 1 file changed, 18 insertions(+), 63 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index b762c376f..10c454e2e 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1570,70 +1570,26 @@ def u_dot_generalized(self, t, u, post_processing=False): w = Vector([omega1, omega2, omega3]) # Angular velocity vector # Retrieve necessary quantities - rho = self.env.density.get_value_opt(z) + ## 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(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 = ( - -1 - * ( - ( - self.rocket.center_of_propellant_position - - self.rocket.center_of_dry_mass_position - ) - * self.rocket._csys - ) - * self.rocket.motor.propellant_mass - / total_mass - ) - r_CM = Vector([0, 0, r_CM_z.get_value_opt(t)]) - r_CM_dot = Vector([0, 0, r_CM_z.differentiate(t)]) + r_CM_z = self.rocket.z_coordinate_com_to_cdm + 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_center_of_dry_mass_position]) ## Nozzle gyration tensor - r_NOZ = ( - -(self.rocket.nozzle_position - self.rocket.center_of_dry_mass_position) - * self.rocket._csys - ) - S_noz_33 = 0.5 * self.rocket.motor.nozzle_radius**2 - S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * r_NOZ**2 - S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 - S_nozzle = Matrix( - [ - [S_noz_11, S_noz_12, S_noz_13], - [S_noz_12, S_noz_22, S_noz_23], - [S_noz_13, S_noz_23, S_noz_33], - ] - ) + S_nozzle = self.rocket.nozzle_gyration_tensor ## Inertia tensor - I_11 = self.rocket.I_11.get_value_opt(t) - I_12 = self.rocket.I_12.get_value_opt(t) - I_13 = self.rocket.I_13.get_value_opt(t) - I_22 = self.rocket.I_22.get_value_opt(t) - I_23 = self.rocket.I_23.get_value_opt(t) - I_33 = self.rocket.I_33.get_value_opt(t) - I = Matrix( - [ - [I_11, I_12, I_13], - [I_12, I_22, I_23], - [I_13, I_23, I_33], - ] - ) + I = self.rocket.get_inertia_tensor_at_time(t) ## Inertia tensor time derivative in the body frame - I_11_dot = self.rocket.I_11.differentiate(t) - I_12_dot = self.rocket.I_12.differentiate(t) - I_13_dot = self.rocket.I_13.differentiate(t) - I_22_dot = self.rocket.I_22.differentiate(t) - I_23_dot = self.rocket.I_23.differentiate(t) - I_33_dot = self.rocket.I_33.differentiate(t) - I_dot = Matrix( - [ - [I_11_dot, I_12_dot, I_13_dot], - [I_12_dot, I_22_dot, I_23_dot], - [I_13_dot, I_23_dot, I_33_dot], - ] - ) - ## Inertia tensor relative to CM + 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 = I - H @@ -1723,8 +1679,8 @@ def u_dot_generalized(self, t, u, post_processing=False): R1 += comp_lift_xb R2 += comp_lift_yb # Add to total moment - M1 -= (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_yb - M2 += (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_xb + M1 -= (comp_cpz + r_CM_t) * comp_lift_yb + M2 += (comp_cpz + r_CM_t) * comp_lift_xb # Calculates Roll Moment try: clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters @@ -1749,13 +1705,12 @@ def u_dot_generalized(self, t, u, post_processing=False): pass weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity.get_value_opt(z)]) T00 = total_mass * r_CM - r_NOZ_vector = Vector([0, 0, r_NOZ]) - T03 = 2 * total_mass_dot * (r_NOZ_vector - r_CM) - 2 * total_mass * r_CM_dot + T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot T04 = ( - self.rocket.motor.thrust.get_value_opt(t) * Vector([0, 0, 1]) + Vector([0, 0, self.rocket.motor.thrust.get_value_opt(t)]) - total_mass * r_CM_ddot - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (r_NOZ_vector - r_CM) + + total_mass_ddot * (r_NOZ - r_CM) ) T05 = total_mass_dot * S_nozzle - I_dot From 5c31a0b3cfb070bed528feaaa44ee1a88eed70b0 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:31:38 -0400 Subject: [PATCH 17/28] DEV: Adds #595 to the CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fbe00cfa0..dc7c5a6d6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ### Added +- ENH: Pre-calculate attributes in Rocket class [#595](https://github.com/RocketPy-Team/RocketPy/pull/595) - ENH: Complex step differentiation [#594](https://github.com/RocketPy-Team/RocketPy/pull/594) - ENH: Exponential backoff decorator (fix #449) [#588](https://github.com/RocketPy-Team/RocketPy/pull/588) - ENH: Function Validation Rework & Swap `np.searchsorted` to `bisect_left` [#582](https://github.com/RocketPy-Team/RocketPy/pull/582) From 40a22c4f24dfbd48471eb62368eb3c9c2dda2f67 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Sun, 5 May 2024 04:34:27 +0000 Subject: [PATCH 18/28] Fix code style issues with Black --- rocketpy/rocket/rocket.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 68bd026e8..0bb7f8caf 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -752,12 +752,12 @@ def evaluate_nozzle_to_center_of_dry_mass_position(self): self.nozzle_to_center_of_dry_mass_position : float Distance between the nozzle exit and the rocket's center of dry mass position, in meters. - """ + """ self.nozzle_to_center_of_dry_mass_position = ( -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys ) return self.nozzle_to_center_of_dry_mass_position - + def evaluate_nozzle_gyration_tensor(self): """Calculates and returns the nozzle gyration tensor relative to the rocket's center of dry mass. The gyration tensor is saved and returned From 68d176df04862920a9af702c8252df6243575458 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 01:09:59 -0400 Subject: [PATCH 19/28] TST: add tests for the new Rocket properties --- tests/unit/test_rocket.py | 86 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 7d942ee60..555add5c4 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -413,3 +413,89 @@ def test_evaluate_center_of_mass(calisto): A predefined instance of the calisto Rocket with a motor, used as a base for testing. """ assert isinstance(calisto.evaluate_center_of_mass(), Function) + + +def test_evaluate_nozzle_to_center_of_dry_mass_position(calisto): + expected_distance = 1.255 + atol = 1e-3 # Equivalent to 1mm + assert ( + pytest.approx(expected_distance, atol) + == calisto.nozzle_to_center_of_dry_mass_position + ) + # Test if calling the function returns the same result + res = calisto.evaluate_nozzle_to_center_of_dry_mass_position() + assert pytest.approx(expected_distance, atol) == res + + +def test_evaluate_nozzle_gyration_tensor(calisto): + expected_gyration_tensor = np.array( + [[0.3940207, 0, 0], [0, 0.3940207, 0], [0, 0, 0.0005445]] + ) + atol = 1e-3 * 1e-2 * 1e-2 # Equivalent to 1g * 1cm^2 + assert np.allclose( + expected_gyration_tensor, np.array(calisto.nozzle_gyration_tensor), atol=atol + ) + # Test if calling the function returns the same result + res = calisto.evaluate_nozzle_gyration_tensor() + assert np.allclose(expected_gyration_tensor, np.array(res), atol=atol) + + +def test_evaluate_z_coordinate_com_to_cdm(calisto): + atol = 1e-3 # Equivalent to 1mm + assert np.allclose( + (calisto.center_of_dry_mass_position - calisto.center_of_mass).source, + calisto.z_coordinate_com_to_cdm.source, + atol=atol, + ) + + +def test_get_inertia_tensor_at_time(calisto): + # Expected values (for t = 0) + # TODO: compute these values by hand or using CAD. + Ixx = 10.31379 + Iyy = 10.31379 + Izz = 0.039942 + + # Set tolerance threshold + atol = 1e-5 + + # Get inertia tensor at t = 0 + inertia_tensor = calisto.get_inertia_tensor_at_time(0) + + # Check if the values are close to the expected ones + assert pytest.approx(Ixx, atol) == inertia_tensor.x[0] + assert pytest.approx(Iyy, atol) == inertia_tensor.y[1] + assert pytest.approx(Izz, atol) == inertia_tensor.z[2] + # Check if products of inertia are zero + assert pytest.approx(0, atol) == inertia_tensor.x[1] + assert pytest.approx(0, atol) == inertia_tensor.x[2] + assert pytest.approx(0, atol) == inertia_tensor.y[0] + assert pytest.approx(0, atol) == inertia_tensor.y[2] + assert pytest.approx(0, atol) == inertia_tensor.z[0] + assert pytest.approx(0, atol) == inertia_tensor.z[1] + + +def test_get_inertia_tensor_derivative_at_time(calisto): + # Expected values (for t = 2s) + # TODO: compute these values by hand or using CAD. + Ixx_dot = -0.634805230901143 + Iyy_dot = -0.634805230901143 + Izz_dot = -0.000671493662305 + + # Set tolerance threshold + atol = 1e-3 + + # Get inertia tensor at t = 2s + inertia_tensor = calisto.get_inertia_tensor_derivative_at_time(2) + + # Check if the values are close to the expected ones + assert pytest.approx(Ixx_dot, atol) == inertia_tensor.x[0] + assert pytest.approx(Iyy_dot, atol) == inertia_tensor.y[1] + assert pytest.approx(Izz_dot, atol) == inertia_tensor.z[2] + # Check if products of inertia are zero + assert pytest.approx(0, atol) == inertia_tensor.x[1] + assert pytest.approx(0, atol) == inertia_tensor.x[2] + assert pytest.approx(0, atol) == inertia_tensor.y[0] + assert pytest.approx(0, atol) == inertia_tensor.y[2] + assert pytest.approx(0, atol) == inertia_tensor.z[0] + assert pytest.approx(0, atol) == inertia_tensor.z[1] From 080e9f031ed7c2cc57761dd9e99a1b33f0c58d7e Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 01:17:15 -0400 Subject: [PATCH 20/28] ENH: Update differentiate_complex_step method to use smaller dx value for more accurate differentiation --- 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 7994c355f..b9efdba95 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2427,7 +2427,7 @@ def differentiate(self, x, dx=1e-6, order=1): + self.get_value_opt(x - dx) ) / dx**2 - def differentiate_complex_step(self, x, dx=1e-6): + def differentiate_complex_step(self, x, dx=1e-12): """Differentiate a Function object at a given point using the complex step method. This method can be faster than ``Function.differentiate`` since it requires only one evaluation of the function. However, the From b691ead9e63a256024a378a725be1cab11df4669 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 01:33:09 -0400 Subject: [PATCH 21/28] ENH: use differentiate.complex_time in aero_surface.py --- rocketpy/rocket/aero_surface.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py index d41240ac9..8f0d37f29 100644 --- a/rocketpy/rocket/aero_surface.py +++ b/rocketpy/rocket/aero_surface.py @@ -367,7 +367,7 @@ def evaluate_geometrical_parameters(self): return None def evaluate_nose_shape(self): - """Calculates and saves nose cone's shape as lists and reavulates the + """Calculates and saves nose cone's shape as lists and re-evaluates the nose cone's length for a given bluffness ratio. The shape is saved as two vectors, one for the x coordinates and one for the y coordinates. @@ -381,7 +381,9 @@ def evaluate_nose_shape(self): # Calculate a function to find the tangential intersection point between the circle and nosecone curve. def find_x_intercept(x): - return x + self.y_nosecone(x) * self.y_nosecone.differentiate(x) + return x + self.y_nosecone(x) * self.y_nosecone.differentiate_complex_step( + x + ) # Calculate a function to find the radius of the nosecone curve def find_radius(x): @@ -770,7 +772,9 @@ def evaluate_lift_coefficient(self): ) # Differentiating at alpha = 0 to get cl_alpha - clalpha2D_incompressible = self.airfoil_cl.differentiate(x=1e-3, dx=1e-3) + clalpha2D_incompressible = self.airfoil_cl.differentiate_complex_step( + x=1e-3, dx=1e-3 + ) # Convert to radians if needed if self.airfoil[1] == "degrees": From 3b3549b85b4b4b5aae694c59dd5519a635d597c5 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sun, 5 May 2024 16:33:40 +0200 Subject: [PATCH 22/28] ENH: finalize merge --- rocketpy/simulation/flight.py | 29 +++-------------------------- 1 file changed, 3 insertions(+), 26 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 898b8ddd7..5f282bf81 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1057,6 +1057,7 @@ def __init_solution_monitors(self): self.impact_state = np.array([0]) self.parachute_events = [] self.post_processed = False + self.__post_processed_variables = [] def __init_flight_state(self): """Initialize flight state variables.""" @@ -1148,26 +1149,7 @@ def __init_controllers(self): def __cache_post_process_variables(self): """Cache post-process variables for simulations with controllers.""" - self.__retrieve_arrays = [ - self.ax_list, - self.ay_list, - self.az_list, - self.alpha1_list, - self.alpha2_list, - self.alpha3_list, - self.R1_list, - self.R2_list, - self.R3_list, - self.M1_list, - self.M2_list, - self.M3_list, - self.pressure_list, - self.density_list, - self.dynamic_viscosity_list, - self.speed_of_sound_list, - self.wind_velocity_x_list, - self.wind_velocity_y_list, - ] + self.__evaluate_post_process = np.array(self.__post_processed_variables) @cached_property def effective_1rl(self): @@ -1289,12 +1271,7 @@ def udot_rail1(self, t, u, post_processing=False): # Use u_dot post processing code for forces, moments and env data self.u_dot_generalized(t, u, post_processing=True) # Save feasible accelerations - self.__post_processed_variables[-1][1] = ax - self.__post_processed_variables[-1][2] = ay - self.__post_processed_variables[-1][3] = az - self.__post_processed_variables[-1][4] = 0 - self.__post_processed_variables[-1][5] = 0 - self.__post_processed_variables[-1][6] = 0 + self.__post_processed_variables[-1][1:7] = [ax, ay, az, 0, 0, 0] return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] From 7531ca6a58f034613cf9554edc05031b8be130da Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 15 May 2024 07:16:22 -0400 Subject: [PATCH 23/28] MNT: rename Rocket.nozzle_to_cdm attribute --- rocketpy/rocket/rocket.py | 18 ++++++++---------- rocketpy/simulation/flight.py | 2 +- tests/unit/test_rocket.py | 9 +++------ 3 files changed, 12 insertions(+), 17 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 0bb7f8caf..7d48d8724 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -150,7 +150,7 @@ class Rocket: defined rocket coordinate system. See :doc:`Positions and Coordinate Systems ` for more information. - Rocket.nozzle_to_center_of_dry_mass_position : float + Rocket.nozzle_to_cdm : float Distance between the nozzle exit and the rocket's center of dry mass position, in meters. Rocket.nozzle_gyration_tensor: Matrix @@ -743,20 +743,20 @@ def evaluate_inertias(self): self.I_23, ) - def evaluate_nozzle_to_center_of_dry_mass_position(self): + def evaluate_nozzle_to_cdm(self): """Evaluates the distance between the nozzle exit and the rocket's - center of dry mass position. + center of dry mass. Returns ------- - self.nozzle_to_center_of_dry_mass_position : float + self.nozzle_to_cdm : float Distance between the nozzle exit and the rocket's center of dry mass position, in meters. """ - self.nozzle_to_center_of_dry_mass_position = ( + self.nozzle_to_cdm = ( -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys ) - return self.nozzle_to_center_of_dry_mass_position + return self.nozzle_to_cdm def evaluate_nozzle_gyration_tensor(self): """Calculates and returns the nozzle gyration tensor relative to the @@ -769,9 +769,7 @@ def evaluate_nozzle_gyration_tensor(self): Matrix containing the nozzle gyration tensor. """ S_noz_33 = 0.5 * self.motor.nozzle_radius**2 - S_noz_11 = S_noz_22 = ( - 0.5 * S_noz_33 + 0.25 * self.nozzle_to_center_of_dry_mass_position**2 - ) + S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * self.nozzle_to_cdm**2 S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 # Due to axis symmetry self.nozzle_gyration_tensor = Matrix( [ @@ -914,7 +912,7 @@ def add_motor(self, motor, position): self.evaluate_dry_mass() self.evaluate_total_mass() self.evaluate_center_of_dry_mass() - self.evaluate_nozzle_to_center_of_dry_mass_position() + self.evaluate_nozzle_to_cdm() self.evaluate_center_of_mass() self.evaluate_dry_inertias() self.evaluate_inertias() diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 10c454e2e..373fd16d0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1581,7 +1581,7 @@ def u_dot_generalized(self, t, u, post_processing=False): 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_center_of_dry_mass_position]) + r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm]) ## Nozzle gyration tensor S_nozzle = self.rocket.nozzle_gyration_tensor ## Inertia tensor diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 555add5c4..40ae2aec0 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -415,15 +415,12 @@ def test_evaluate_center_of_mass(calisto): assert isinstance(calisto.evaluate_center_of_mass(), Function) -def test_evaluate_nozzle_to_center_of_dry_mass_position(calisto): +def test_evaluate_nozzle_to_cdm(calisto): expected_distance = 1.255 atol = 1e-3 # Equivalent to 1mm - assert ( - pytest.approx(expected_distance, atol) - == calisto.nozzle_to_center_of_dry_mass_position - ) + assert pytest.approx(expected_distance, atol) == calisto.nozzle_to_cdm # Test if calling the function returns the same result - res = calisto.evaluate_nozzle_to_center_of_dry_mass_position() + res = calisto.evaluate_nozzle_to_cdm() assert pytest.approx(expected_distance, atol) == res From 2eedf980a92f4534869366bd50043f284b4bdab2 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 15 May 2024 07:16:42 -0400 Subject: [PATCH 24/28] DOC: adds Inertia Tensor to the toctree --- docs/user/rocket.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/user/rocket.rst b/docs/user/rocket.rst index 6b4a9f257..fa09f058e 100644 --- a/docs/user/rocket.rst +++ b/docs/user/rocket.rst @@ -12,6 +12,7 @@ Defining a Rocket in RocketPy is simple and requires a few steps: 4. Add, if desired, parachutes; 5. Set, if desired, rail guides; 6. See results. +7. Inertia Tensors. Lets go through each of these steps in detail. From 1b038d5f3ac97354a480e9aedf29ce7bcebcbf60 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 15 May 2024 09:02:27 -0400 Subject: [PATCH 25/28] MNT: use rocket.z_coordinate_com_to_cdm in u_dot --- rocketpy/simulation/flight.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 798a08170..68672d802 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1341,12 +1341,8 @@ def u_dot(self, t, u, post_processing=False): ) * self.rocket._csys ) - # c = -self.rocket.distance_rocket_nozzle - c = ( - -(self.rocket.nozzle_position - self.rocket.center_of_dry_mass_position) - * self.rocket._csys - ) - a = b * Mt / M + c = self.rocket.nozzle_to_cdm + a = self.rocket.z_coordinate_com_to_cdm.get_value_opt(t) rN = self.rocket.motor.nozzle_radius # Prepare transformation matrix a11 = 1 - 2 * (e2**2 + e3**2) From 8b91f0b8333bc884407c858bb7402bb286a51333 Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Thu, 16 May 2024 09:12:54 -0300 Subject: [PATCH 26/28] MNT: remove True check Co-authored-by: Pedro Henrique Marinho Bressan <87212571+phmbressan@users.noreply.github.com> --- rocketpy/simulation/flight.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 64dccae71..1996a4d60 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1138,7 +1138,7 @@ def __init_controllers(self): """Initialize controllers""" self._controllers = self.rocket._controllers[:] if self._controllers: - if self.time_overshoot == True: + if self.time_overshoot: self.time_overshoot = False warnings.warn( "time_overshoot has been set to False due to the presence of controllers. " From 464f2017c46367857582b33b81c7379d4d39cd80 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 16 May 2024 14:40:50 +0200 Subject: [PATCH 27/28] ENH: remove __cache_post_process_variables --- rocketpy/simulation/flight.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 5f282bf81..523402fe1 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1005,7 +1005,8 @@ def __simulate(self, verbose): self.t_final = self.t self.__transform_pressure_signals_lists_to_functions() if self._controllers: - self.__cache_post_process_variables() + # cache post process variables + self.__evaluate_post_process = np.array(self.__post_processed_variables) if verbose: print(f"\n>>> Simulation Completed at Time: {self.t:3.4f} s") @@ -1147,10 +1148,6 @@ def __init_controllers(self): for air_brakes in self.rocket.air_brakes: air_brakes._reset() - def __cache_post_process_variables(self): - """Cache post-process variables for simulations with controllers.""" - self.__evaluate_post_process = np.array(self.__post_processed_variables) - @cached_property def effective_1rl(self): """Original rail length minus the distance measured from nozzle exit From 61ebbc754dc41a7436b9ddc5ccc39310fecc4440 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> Date: Thu, 16 May 2024 20:30:21 +0000 Subject: [PATCH 28/28] MNT: renames z_coordinate_com_to_cdm to com_to_cdm_function --- rocketpy/rocket/rocket.py | 22 +++++++++++----------- rocketpy/simulation/flight.py | 4 ++-- tests/unit/test_rocket.py | 4 ++-- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index c20a35ef2..a0db26eff 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -73,7 +73,7 @@ class Rocket: for more information regarding the coordinate system. Expressed in meters as a function of time. - Rocket.z_coordinate_com_to_cdm : Function + Rocket.com_to_cdm_function : Function Function of time expressing the z-coordinate of the center of mass relative to the center of dry mass. Rocket.reduced_mass : Function @@ -780,24 +780,24 @@ def evaluate_nozzle_gyration_tensor(self): ) return self.nozzle_gyration_tensor - def evaluate_z_coordinate_com_to_cdm(self): + def evaluate_com_to_cdm_function(self): """Evaluates the z-coordinate of the center of mass (COM) relative to the center of dry mass (CDM). Notes ----- - 1. The `z_coordinate_com_to_cdm` plus `center_of_mass` should be equal + 1. The `com_to_cdm_function` plus `center_of_mass` should be equal to `center_of_dry_mass_position` at every time step. - 2. The `z_coordinate_com_to_cdm` is a function of time and will usually + 2. The `com_to_cdm_function` is a function of time and will usually already be discretized. Returns ------- - self.z_coordinate_com_to_cdm : Function + self.com_to_cdm_function : Function Function of time expressing the z-coordinate of the center of mass relative to the center of dry mass. """ - self.z_coordinate_com_to_cdm = ( + self.com_to_cdm_function = ( -1 * ( (self.center_of_propellant_position - self.center_of_dry_mass_position) @@ -806,10 +806,10 @@ def evaluate_z_coordinate_com_to_cdm(self): * self.motor.propellant_mass / self.total_mass ) - self.z_coordinate_com_to_cdm.set_inputs("Time (s)") - self.z_coordinate_com_to_cdm.set_outputs("Z Coordinate COM to CDM (m)") - self.z_coordinate_com_to_cdm.set_title("Z Coordinate COM to CDM") - return self.z_coordinate_com_to_cdm + self.com_to_cdm_function.set_inputs("Time (s)") + self.com_to_cdm_function.set_outputs("Z Coordinate COM to CDM (m)") + self.com_to_cdm_function.set_title("Z Coordinate COM to CDM") + return self.com_to_cdm_function def get_inertia_tensor_at_time(self, t): """Returns a Matrix representing the inertia tensor of the rocket with @@ -921,7 +921,7 @@ def add_motor(self, motor, position): self.evaluate_center_of_pressure() self.evaluate_stability_margin() self.evaluate_static_margin() - self.evaluate_z_coordinate_com_to_cdm() + self.evaluate_com_to_cdm_function() self.evaluate_nozzle_gyration_tensor() def add_surfaces(self, surfaces, positions): diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 68672d802..4cbc752fb 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1342,7 +1342,7 @@ def u_dot(self, t, u, post_processing=False): * self.rocket._csys ) c = self.rocket.nozzle_to_cdm - a = self.rocket.z_coordinate_com_to_cdm.get_value_opt(t) + a = self.rocket.com_to_cdm_function.get_value_opt(t) rN = self.rocket.motor.nozzle_radius # Prepare transformation matrix a11 = 1 - 2 * (e2**2 + e3**2) @@ -1571,7 +1571,7 @@ def u_dot_generalized(self, t, u, post_processing=False): 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.z_coordinate_com_to_cdm + 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)]) diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 40ae2aec0..4d934efef 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -437,11 +437,11 @@ def test_evaluate_nozzle_gyration_tensor(calisto): assert np.allclose(expected_gyration_tensor, np.array(res), atol=atol) -def test_evaluate_z_coordinate_com_to_cdm(calisto): +def test_evaluate_com_to_cdm_function(calisto): atol = 1e-3 # Equivalent to 1mm assert np.allclose( (calisto.center_of_dry_mass_position - calisto.center_of_mass).source, - calisto.z_coordinate_com_to_cdm.source, + calisto.com_to_cdm_function.source, atol=atol, )