From 21331eea9de80f4c3c38c6ef21f83442c3eeb878 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 25 Mar 2024 12:52:54 -0400 Subject: [PATCH 01/71] MNT: use Function.get_value_opt() instead of Function.get_value() - Overall, this already speeds up the flight simulations --- rocketpy/environment/environment.py | 37 +++-- rocketpy/mathutils/function.py | 24 +-- rocketpy/plots/environment_analysis_plots.py | 2 +- rocketpy/rocket/rocket.py | 11 +- rocketpy/simulation/flight.py | 153 ++++++++++++------- 5 files changed, 140 insertions(+), 87 deletions(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 546d446e8..4d9858af3 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -1488,17 +1488,19 @@ def process_custom_atmosphere( # Check maximum height of custom wind input if not callable(self.wind_velocity_x.source): max_expected_height = max(self.wind_velocity_x[-1, 0], max_expected_height) - if not callable(self.wind_velocity_y.source): - max_expected_height = max(self.wind_velocity_y[-1, 0], max_expected_height) - # Compute wind profile direction and heading - wind_heading = ( - lambda h: np.arctan2(self.wind_velocity_x(h), self.wind_velocity_y(h)) - * (180 / np.pi) - % 360 - ) + def wind_heading_func(h): + return ( + np.arctan2( + self.wind_velocity_x.get_value_opt(h), + self.wind_velocity_y.get_value_opt(h), + ) + * (180 / np.pi) + % 360 + ) + self.wind_heading = Function( - wind_heading, + wind_heading_func, inputs="Height Above Sea Level (m)", outputs="Wind Heading (Deg True)", interpolation="linear", @@ -1515,7 +1517,10 @@ def wind_direction(h): ) def wind_speed(h): - return np.sqrt(self.wind_velocity_x(h) ** 2 + self.wind_velocity_y(h) ** 2) + return np.sqrt( + self.wind_velocity_x.get_value_opt(h) ** 2 + + self.wind_velocity_y.get_value_opt(h) ** 2 + ) self.wind_speed = Function( wind_speed, @@ -3142,22 +3147,26 @@ def add_wind_gust(self, wind_gust_x, wind_gust_y): # Reset wind heading and velocity magnitude self.wind_heading = Function( lambda h: (180 / np.pi) - * np.arctan2(self.wind_velocity_x(h), self.wind_velocity_y(h)) + * np.arctan2( + self.wind_velocity_x.get_value_opt(h), + self.wind_velocity_y.get_value_opt(h), + ) % 360, "Height (m)", "Wind Heading (degrees)", extrapolation="constant", ) self.wind_speed = Function( - lambda h: (self.wind_velocity_x(h) ** 2 + self.wind_velocity_y(h) ** 2) + lambda h: ( + self.wind_velocity_x.get_value_opt(h) ** 2 + + self.wind_velocity_y.get_value_opt(h) ** 2 + ) ** 0.5, "Height (m)", "Wind Speed (m/s)", extrapolation="constant", ) - return None - def info(self): """Prints most important data and graphs available about the Environment. diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index cefed044d..7a54ce8e7 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -1988,7 +1988,7 @@ def __add__(self, other): # Create new Function object return Function(source, inputs, outputs, interpolation, extrapolation) else: - return Function(lambda x: (self.get_value(x) + other(x))) + return Function(lambda x: (self.get_value_opt(x) + other(x))) # If other is Float except... except AttributeError: if isinstance(other, NUMERICAL_TYPES): @@ -2009,10 +2009,10 @@ def __add__(self, other): source, inputs, outputs, interpolation, extrapolation ) else: - return Function(lambda x: (self.get_value(x) + other)) + return Function(lambda x: (self.get_value_opt(x) + other)) # Or if it is just a callable elif callable(other): - return Function(lambda x: (self.get_value(x) + other(x))) + return Function(lambda x: (self.get_value_opt(x) + other(x))) def __radd__(self, other): """Sums 'other' and a Function object and returns a new Function @@ -2055,7 +2055,7 @@ def __sub__(self, other): try: return self + (-other) except TypeError: - return Function(lambda x: (self.get_value(x) - other(x))) + return Function(lambda x: (self.get_value_opt(x) - other(x))) def __rsub__(self, other): """Subtracts a Function object from 'other' and returns a new Function @@ -2335,10 +2335,10 @@ def __pow__(self, other): source, inputs, outputs, interpolation, extrapolation ) else: - return Function(lambda x: (self.get_value(x) ** other)) + return Function(lambda x: (self.get_value_opt(x) ** other)) # Or if it is just a callable elif callable(other): - return Function(lambda x: (self.get_value(x) ** other(x))) + return Function(lambda x: (self.get_value_opt(x) ** other(x))) def __rpow__(self, other): """Raises 'other' to the power of a Function object and returns @@ -2371,10 +2371,10 @@ def __rpow__(self, other): # Create new Function object return Function(source, inputs, outputs, interpolation, extrapolation) else: - return Function(lambda x: (other ** self.get_value(x))) + return Function(lambda x: (other ** self.get_value_opt(x))) # Or if it is just a callable elif callable(other): - return Function(lambda x: (other(x) ** self.get_value(x))) + return Function(lambda x: (other(x) ** self.get_value_opt(x))) def __matmul__(self, other): """Operator @ as an alias for composition. Therefore, this @@ -2549,10 +2549,12 @@ def differentiate(self, x, dx=1e-6, order=1): Evaluated derivative. """ if order == 1: - return (self.get_value(x + dx) - self.get_value(x - dx)) / (2 * dx) + return (self.get_value_opt(x + dx) - self.get_value_opt(x - dx)) / (2 * dx) elif order == 2: return ( - self.get_value(x + dx) - 2 * self.get_value(x) + self.get_value(x - dx) + self.get_value_opt(x + dx) + - 2 * self.get_value_opt(x) + + self.get_value_opt(x - dx) ) / dx**2 def identity_function(self): @@ -3262,7 +3264,7 @@ def calc_output(func, inputs): """ output = np.zeros(len(inputs)) for j, value in enumerate(inputs): - output[j] = func.get_value(value) + output[j] = func.get_value_opt(value) return output input_data = [] diff --git a/rocketpy/plots/environment_analysis_plots.py b/rocketpy/plots/environment_analysis_plots.py index 9fbf68d95..26727aba9 100644 --- a/rocketpy/plots/environment_analysis_plots.py +++ b/rocketpy/plots/environment_analysis_plots.py @@ -1454,7 +1454,7 @@ def wind_speed_profile_grid(self, clear_range_limits=False): if self.env_analysis.forecast: forecast = self.env_analysis.forecast y = self.env_analysis.average_wind_speed_profile_by_hour[hour][1] - x = forecast[hour].wind_speed.get_value(y) * convert_units( + x = forecast[hour].wind_speed.get_value_opt(y) * convert_units( 1, "m/s", self.env_analysis.unit_system["wind_speed"] ) ax.plot(x, y, "b--") diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 8154f2967..b07aeeb84 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -542,7 +542,11 @@ def evaluate_stability_margin(self): """ self.stability_margin.set_source( lambda mach, time: ( - (self.center_of_mass(time) - self.cp_position(mach)) / (2 * self.radius) + ( + self.center_of_mass.get_value_opt(time) + - self.cp_position.get_value_opt(mach) + ) + / (2 * self.radius) ) * self._csys ) @@ -560,7 +564,10 @@ def evaluate_static_margin(self): """ # Calculate static margin self.static_margin.set_source( - lambda time: (self.center_of_mass(time) - self.cp_position(0)) + lambda time: ( + self.center_of_mass.get_value_opt(time) + - self.cp_position.get_value_opt(0) + ) / (2 * self.radius) ) # Change sign if coordinate system is upside down diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 77632d2f2..75c31917d 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -714,11 +714,13 @@ def __init__( parachute.noise_signal.append([node.t, noise]) parachute.noisy_pressure_signal.append([node.t, pressure + noise]) # Gets height above ground level considering noise - hAGL = ( - self.env.barometric_height(pressure + noise) + height_above_ground_level = ( + self.env.barometric_height.get_value_opt(pressure + noise) - self.env.elevation ) - if parachute.triggerfunc(pressure + noise, hAGL, self.y_sol): + if parachute.triggerfunc( + pressure + noise, height_above_ground_level, self.y_sol + ): # print('\nEVENT DETECTED') # print('Parachute Triggered') # print('Name: ', parachute.name, ' | Lag: ', parachute.lag) @@ -1020,13 +1022,17 @@ def __init__( [overshootable_node.t, pressure + noise] ) # Gets height above ground level considering noise - hAGL = ( - self.env.barometric_height(pressure + noise) + height_above_ground_level = ( + self.env.barometric_height.get_value_opt( + pressure + noise + ) - self.env.elevation ) if parachute.triggerfunc( - pressure + noise, hAGL, overshootable_node.y + pressure + noise, + height_above_ground_level, + overshootable_node.y, ): # print('\nEVENT DETECTED') # print('Parachute Triggered') @@ -1219,25 +1225,34 @@ def effective_2rl(self): @cached_property def frontal_surface_wind(self): - # Surface wind magnitude in the frontal direction at the rail's elevation - wind_u = self.env.wind_velocity_x(self.env.elevation) - wind_v = self.env.wind_velocity_y(self.env.elevation) + """Frontal wind velocity at the surface level. The frontal wind is + defined as the wind blowing in the direction of the rocket's heading. + + Returns + ------- + float + Wind velocity in the frontal direction at the surface level. + """ + wind_u = self.env.wind_velocity_x.get_value_opt(self.env.elevation) + wind_v = self.env.wind_velocity_y.get_value_opt(self.env.elevation) heading_rad = self.heading * np.pi / 180 - frontal_surface_wind = wind_u * np.sin(heading_rad) + wind_v * np.cos( - heading_rad - ) - return frontal_surface_wind + return wind_u * np.sin(heading_rad) + wind_v * np.cos(heading_rad) @cached_property def lateral_surface_wind(self): - # Surface wind magnitude in the lateral direction at the rail's elevation - wind_u = self.env.wind_velocity_x(self.env.elevation) - wind_v = self.env.wind_velocity_y(self.env.elevation) + """Lateral wind velocity at the surface level. The lateral wind is + defined as the wind blowing perpendicular to the rocket's heading. + + Returns + ------- + float + Wind velocity in the lateral direction at the surface level. + """ + wind_u = self.env.wind_velocity_x.get_value_opt(self.env.elevation) + wind_v = self.env.wind_velocity_y.get_value_opt(self.env.elevation) heading_rad = self.heading * np.pi / 180 - lateral_surface_wind = -wind_u * np.cos(heading_rad) + wind_v * np.sin( - heading_rad - ) - return lateral_surface_wind + + return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad) def udot_rail1(self, t, u, post_processing=False): """Calculates derivative of u state vector with respect to time @@ -1288,7 +1303,9 @@ def udot_rail1(self, t, u, post_processing=False): R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff) # Calculate Linear acceleration - a3 = (R3 + thrust) / M - (e0**2 - e1**2 - e2**2 + e3**2) * self.env.gravity(z) + a3 = (R3 + thrust) / M - ( + e0**2 - e1**2 - e2**2 + e3**2 + ) * self.env.gravity.get_value_opt(z) if a3 > 0: ax = 2 * (e1 * e3 + e0 * e2) * a3 ay = 2 * (e2 * e3 - e0 * e1) * a3 @@ -1389,7 +1406,7 @@ def u_dot(self, t, u, post_processing=False): # b = -self.rocket.distance_rocket_propellant b = ( -( - self.rocket.center_of_propellant_position(0) + self.rocket.center_of_propellant_position.get_value_opt(0) - self.rocket.center_of_dry_mass_position ) * self.rocket._csys @@ -1435,7 +1452,7 @@ def u_dot(self, t, u, post_processing=False): 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( + air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( air_brakes.deployment_level, free_stream_mach ) air_brakes_force = ( @@ -1490,7 +1507,9 @@ def u_dot(self, t, u, post_processing=False): comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed if -1 * comp_stream_vz_bn < 1: comp_attack_angle = np.arccos(-comp_stream_vz_bn) - c_lift = aero_surface.cl(comp_attack_angle, free_stream_mach) + c_lift = aero_surface.cl.get_value_opt( + comp_attack_angle, free_stream_mach + ) # component lift force magnitude comp_lift = ( 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift @@ -1513,14 +1532,14 @@ def u_dot(self, t, u, post_processing=False): * reference_area * 2 * surface_radius - * clf_delta(free_stream_mach) + * clf_delta.get_value_opt(free_stream_mach) * cant_angle_rad ) M3d = ( (1 / 2 * rho * free_stream_speed) * reference_area * (2 * surface_radius) ** 2 - * cld_omega(free_stream_mach) + * cld_omega.get_value_opt(free_stream_mach) * omega3 / 2 ) @@ -1565,7 +1584,7 @@ def u_dot(self, t, u, post_processing=False): (R3 - b * Mt * (alpha2 - omega1 * omega3) + thrust) / M, ] ax, ay, az = np.dot(K, L) - az -= self.env.gravity(z) # Include gravity + az -= self.env.gravity.get_value_opt(z) # Include gravity # Create u_dot u_dot = [ @@ -1734,7 +1753,7 @@ def u_dot_generalized(self, t, u, post_processing=False): 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( + air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( air_brakes.deployment_level, free_stream_mach ) air_brakes_force = ( @@ -1784,7 +1803,9 @@ def u_dot_generalized(self, t, u, post_processing=False): comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed if -1 * comp_stream_vz_bn < 1: comp_attack_angle = np.arccos(-comp_stream_vz_bn) - c_lift = aero_surface.cl(comp_attack_angle, comp_stream_mach) + c_lift = aero_surface.cl.get_value_opt( + comp_attack_angle, comp_stream_mach + ) # Component lift force magnitude comp_lift = ( 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift @@ -1807,28 +1828,28 @@ def u_dot_generalized(self, t, u, post_processing=False): * reference_area * 2 * surface_radius - * clf_delta(comp_stream_mach) + * clf_delta.get_value_opt(comp_stream_mach) * cant_angle_rad ) M3d = ( (1 / 2 * rho * comp_stream_speed) * reference_area * (2 * surface_radius) ** 2 - * cld_omega(comp_stream_mach) + * cld_omega.get_value_opt(comp_stream_mach) * omega3 / 2 ) M3 += M3f - M3d except AttributeError: pass - weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity(z)]) + weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity.get_value_opt(z)]) T00 = total_mass * r_CM T03 = ( 2 * total_mass_dot * (Vector([0, 0, r_NOZ]) - r_CM) - 2 * total_mass * r_CM_dot ) T04 = ( - self.rocket.motor.thrust(t) * Vector([0, 0, 1]) + self.rocket.motor.thrust.get_value_opt(t) * Vector([0, 0, 1]) - total_mass * r_CM_ddot - 2 * total_mass_dot * r_CM_dot + total_mass_ddot * (Vector([0, 0, r_NOZ]) - r_CM) @@ -1952,12 +1973,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] @@ -2197,7 +2226,7 @@ def max_speed_time(self): @cached_property def max_speed(self): """Maximum speed reached by the rocket.""" - return self.speed(self.max_speed_time) + return self.speed.get_value_opt(self.max_speed_time) # Accelerations @funcify_method("Time (s)", "acceleration Magnitude (m/s²)") @@ -2223,7 +2252,7 @@ def max_acceleration_power_on_time(self): @cached_property def max_acceleration_power_on(self): """Maximum acceleration reached by the rocket during motor burn.""" - return self.acceleration(self.max_acceleration_power_on_time) + return self.acceleration.get_value_opt(self.max_acceleration_power_on_time) @cached_property def max_acceleration_power_off_time(self): @@ -2240,7 +2269,7 @@ def max_acceleration_power_off_time(self): @cached_property def max_acceleration_power_off(self): """Maximum acceleration reached by the rocket after motor burn.""" - return self.acceleration(self.max_acceleration_power_off_time) + return self.acceleration.get_value_opt(self.max_acceleration_power_off_time) @cached_property def max_acceleration_time(self): @@ -2388,7 +2417,7 @@ def free_stream_speed(self): @cached_property def apogee_freestream_speed(self): """Freestream speed at apogee in m/s.""" - return self.free_stream_speed(self.apogee_time) + return self.free_stream_speed.get_value_opt(self.apogee_time) # Mach Number @funcify_method("Time (s)", "Mach Number", "spline", "zero") @@ -2405,7 +2434,7 @@ def max_mach_number_time(self): @cached_property def max_mach_number(self): """Maximum Mach number.""" - return self.mach_number(self.max_mach_number_time) + return self.mach_number.get_value_opt(self.max_mach_number_time) # Stability Margin @cached_property @@ -2417,7 +2446,7 @@ def max_stability_margin_time(self): @cached_property def max_stability_margin(self): """Maximum stability margin.""" - return self.stability_margin(self.max_stability_margin_time) + return self.stability_margin.get_value_opt(self.max_stability_margin_time) @cached_property def min_stability_margin_time(self): @@ -2428,7 +2457,7 @@ def min_stability_margin_time(self): @cached_property def min_stability_margin(self): """Minimum stability margin.""" - return self.stability_margin(self.min_stability_margin_time) + return self.stability_margin.get_value_opt(self.min_stability_margin_time) @property def initial_stability_margin(self): @@ -2438,7 +2467,7 @@ def initial_stability_margin(self): ------- float """ - return self.stability_margin(self.time[0]) + return self.stability_margin.get_value_opt(self.time[0]) @property def out_of_rail_stability_margin(self): @@ -2448,7 +2477,7 @@ def out_of_rail_stability_margin(self): ------- float """ - return self.stability_margin(self.out_of_rail_time) + return self.stability_margin.get_value_opt(self.out_of_rail_time) # Reynolds Number @funcify_method("Time (s)", "Reynolds Number", "spline", "zero") @@ -2467,7 +2496,7 @@ def max_reynolds_number_time(self): @cached_property def max_reynolds_number(self): """Maximum Reynolds number.""" - return self.reynolds_number(self.max_reynolds_number_time) + return self.reynolds_number.get_value_opt(self.max_reynolds_number_time) # Dynamic Pressure @funcify_method("Time (s)", "Dynamic Pressure (Pa)", "spline", "zero") @@ -2484,7 +2513,7 @@ def max_dynamic_pressure_time(self): @cached_property def max_dynamic_pressure(self): """Maximum dynamic pressure.""" - return self.dynamic_pressure(self.max_dynamic_pressure_time) + return self.dynamic_pressure.get_value_opt(self.max_dynamic_pressure_time) # Total Pressure @funcify_method("Time (s)", "Total Pressure (Pa)", "spline", "zero") @@ -2500,7 +2529,7 @@ def max_total_pressure_time(self): @cached_property def max_total_pressure(self): """Maximum total pressure.""" - return self.total_pressure(self.max_total_pressure_time) + return self.total_pressure.get_value_opt(self.max_total_pressure_time) # Dynamics functions and variables @@ -2606,7 +2635,7 @@ def angle_of_attack(self): for i in self.time ] # Define freestream speed list - free_stream_speed = [self.free_stream_speed(i) for i in self.time] + free_stream_speed = [self.free_stream_speed.get_value_opt(i) for i in self.time] free_stream_speed = np.nan_to_num(free_stream_speed) # Normalize dot product @@ -3150,12 +3179,12 @@ def export_pressures(self, file_name, time_step): if len(self.rocket.parachutes) == 0: print("No parachutes in the rocket, saving static pressure.") for t in time_points: - file.write(f"{t:f}, {self.pressure(t):.5f}\n") + file.write(f"{t:f}, {self.pressure.get_value_opt(t):.5f}\n") else: for parachute in self.rocket.parachutes: for t in time_points: - p_cl = parachute.clean_pressure_signal_function(t) - p_ns = parachute.noisy_pressure_signal_function(t) + p_cl = parachute.clean_pressure_signal_function.get_value_opt(t) + p_ns = parachute.noisy_pressure_signal_function.get_value_opt(t) file.write(f"{t:f}, {p_cl:.5f}, {p_ns:.5f}\n") # We need to save only 1 parachute data break @@ -3318,9 +3347,9 @@ def export_kml( for t in time_points: coords.append( ( - self.longitude(t), - self.latitude(t), - self.altitude(t), + self.longitude.get_value_opt(t), + self.latitude.get_value_opt(t), + self.altitude.get_value_opt(t), ) ) trajectory.coords = coords @@ -3330,7 +3359,13 @@ def export_kml( # Ensure you use the correct value on self.env.elevation, otherwise # the trajectory path can be offset from ground for t in time_points: - coords.append((self.longitude(t), self.latitude(t), self.z(t))) + coords.append( + ( + self.longitude.get_value_opt(t), + self.latitude.get_value_opt(t), + self.z.get_value_opt(t), + ) + ) trajectory.coords = coords trajectory.altitudemode = simplekml.AltitudeMode.absolute # Modify style of trajectory linestring From 11feb2a3c3d65106bf0344e434fb542e4ab42daf Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 25 Mar 2024 12:56:37 -0400 Subject: [PATCH 02/71] MNT: refactor and improve Function initialization - better input validators - better set_source method - now the function _check_user_input is not called twice! --- rocketpy/mathutils/function.py | 327 ++++++++++++++++----------------- 1 file changed, 154 insertions(+), 173 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 7a54ce8e7..0d88af585 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -110,20 +110,19 @@ def __init__( (II) Fields in CSV files may be enclosed in double quotes. If fields are not quoted, double quotes should not appear inside them. """ - if inputs is None: - inputs = ["Scalar"] - if outputs is None: - outputs = ["Scalar"] - inputs, outputs, interpolation, extrapolation = self._check_user_input( - source, inputs, outputs, interpolation, extrapolation - ) + inputs, outputs = Function._validate_inputs_outputs(inputs, outputs) # initialize variables to avoid errors when being called by other methods self.get_value_opt = None self.__polynomial_coefficients__ = None self.__akima_coefficients__ = None self.__spline_coefficients__ = None + self.EXTRAPOLATION_TYPES = { + "zero": 0, + "natural": 1, + "constant": 2, + } # store variables self.set_inputs(inputs) @@ -212,55 +211,30 @@ def set_source(self, source): self : Function Returns the Function instance. """ - *_, interpolation, extrapolation = self._check_user_input( + source, inputs, outputs, interpolation, extrapolation = self._check_user_input( source, self.__inputs__, self.__outputs__, self.__interpolation__, self.__extrapolation__, ) - # If the source is a Function - if isinstance(source, Function): - source = source.get_source() - # Import CSV if source is a string or Path and convert values to ndarray - if isinstance(source, (str, Path)): - with open(source, "r") as file: - try: - source = np.loadtxt(file, delimiter=",", dtype=float) - except ValueError: - # If an error occurs, headers are present - source = np.loadtxt(source, delimiter=",", dtype=float, skiprows=1) - except Exception as e: - raise ValueError( - "The source file is not a valid csv or txt file." - ) from e - - # Convert to ndarray if source is a list - if isinstance(source, (list, tuple)): - source = np.array(source, dtype=np.float64) - # Convert number source into vectorized lambda function - if isinstance(source, (int, float)): - temp = 1 * source - - def source_function(_): - return temp - - source = source_function + # updates inputs and outputs (could be modified due to csv headers) + self.set_inputs(inputs) + self.set_outputs(outputs) # Handle callable source or number source if callable(source): - # Set source self.source = source - # Set get_value_opt self.get_value_opt = source + self.__interpolation__ = None + self.__extrapolation__ = None + # Set arguments name and domain dimensions parameters = signature(source).parameters self.__dom_dim__ = len(parameters) if self.__inputs__ == ["Scalar"]: self.__inputs__ = list(parameters) - # Set interpolation and extrapolation - self.__interpolation__ = None - self.__extrapolation__ = None + # Handle ndarray source else: # Check to see if dimensions match incoming data set @@ -273,40 +247,23 @@ def source_function(_): self.__dom_dim__ = new_total_dim - 1 self.__inputs__ = self.__dom_dim__ * self.__inputs__ - # Do things if domDim is 1 + # if Function is 1D, sort source by x. If 2D, set z if self.__dom_dim__ == 1: source = source[source[:, 0].argsort()] - - self.x_array = source[:, 0] - self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] - - self.y_array = source[:, 1] - self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] - - # Finally set data source as source - self.source = source - # Do things if function is multivariate - else: - self.x_array = source[:, 0] - self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] - - self.y_array = source[:, 1] - self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] - + elif self.__dom_dim__ == 2: self.z_array = source[:, 2] self.z_initial, self.z_final = self.z_array[0], self.z_array[-1] - # Finally set data source as source - self.source = source - # Update extrapolation method - if self.__extrapolation__ is None: - self.set_extrapolation(extrapolation) - # Set default interpolation for point source if it hasn't - if self.__interpolation__ is None: - self.set_interpolation(interpolation) - else: - # Updates interpolation coefficients - self.set_interpolation(self.__interpolation__) + # Set x and y arrays (common for 1D or multivariate) + self.x_array = source[:, 0] + self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] + self.y_array = source[:, 1] + self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] + + # Finally set source, update extrapolation and interpolation + self.source = source + self.set_extrapolation(extrapolation) + self.set_interpolation(interpolation) return self @cached_property @@ -396,14 +353,12 @@ def set_get_value_opt(self): x_data = self.x_array y_data = self.y_array x_min, x_max = self.x_initial, self.x_final - if self.__extrapolation__ == "zero": - extrapolation = 0 # Extrapolation is zero - elif self.__extrapolation__ == "natural": - extrapolation = 1 # Extrapolation is natural - elif self.__extrapolation__ == "constant": - extrapolation = 2 # Extrapolation is constant - else: - raise ValueError(f"Invalid extrapolation type {self.__extrapolation__}") + try: + extrapolation = self.EXTRAPOLATION_TYPES[self.__extrapolation__] + except KeyError as err: + raise ValueError( + f"Invalid extrapolation type {self.__extrapolation__}" + ) from err # Crete method to interpolate this info for each interpolation type if self.__interpolation__ == "spline": @@ -2998,6 +2953,8 @@ def savetxt( file.write(header_line + newline) np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline) + # Input validators + @staticmethod def _check_user_input( source, @@ -3073,114 +3030,138 @@ def _check_user_input( >>> extrapolation 'zero' """ + if isinstance(source, Function): + source = source.get_source() + + elif isinstance(source, (str, Path)): + # Read csv or txt files and create a numpy array + try: + source = np.loadtxt(source, delimiter=",", dtype=np.float64) + except ValueError: + with open(source, "r") as file: + header, *data = file.read().splitlines() + + header = [label.strip("'").strip('"') for label in header.split(",")] + source = np.loadtxt(data, delimiter=",", dtype=np.float64) + + if len(source[0]) == len(header): + if inputs == ["Scalar"]: + inputs = header[:-1] + if outputs == ["Scalar"]: + outputs = [header[-1]] + except Exception as e: + raise ValueError( + "Could not read the csv or txt file to create Function source." + ) from e + + if isinstance(source, list): + # Triggers an error if source is not a list of numbers + source = np.array(source, dtype=np.float64) + + if isinstance(source, np.ndarray): + inputs, interpolation, extrapolation = ( + Function._validate_interpolation_and_extrapolation( + inputs, interpolation, extrapolation, source + ) + ) + Function._validate_source_dimensions(inputs, outputs, source) + + if isinstance(source, (int, float)): + # Convert number source into vectorized lambda function + temp = 1 * source + + def source_function(_): + return temp + + source = source_function + return source, inputs, outputs, interpolation, extrapolation + + @staticmethod + def _validate_inputs_outputs(inputs, outputs): # None | st | list[str] + if inputs is None: + inputs = ["Scalar"] + if outputs is None: + outputs = ["Scalar"] # check output type and dimensions - if isinstance(outputs, str): - outputs = [outputs] if isinstance(inputs, str): inputs = [inputs] - - if len(outputs) > 1: + if isinstance(outputs, str): + outputs = [outputs] + elif len(outputs) > 1: raise ValueError( "Output must either be a string or have dimension 1, " + f"it currently has dimension ({len(outputs)})." ) + return inputs, outputs - # check source for data type - # if list or ndarray, check for dimensions, interpolation and extrapolation - if isinstance(source, Function): - source = source.get_source() - if isinstance(source, (list, np.ndarray, str, Path)): - # Deal with csv or txt - if isinstance(source, (str, Path)): - # Convert to numpy array - try: - source = np.loadtxt(source, delimiter=",", dtype=float) - except ValueError: - with open(source, "r") as file: - header, *data = file.read().splitlines() - - header = [ - label.strip("'").strip('"') for label in header.split(",") - ] - source = np.loadtxt(data, delimiter=",", dtype=float) - - if len(source[0]) == len(header): - if inputs == ["Scalar"]: - inputs = header[:-1] - if outputs == ["Scalar"]: - outputs = [header[-1]] - except Exception as e: - raise ValueError( - "The source file is not a valid csv or txt file." - ) from e - - else: - # this will also trigger an error if the source is not a list of - # numbers or if the array is not homogeneous - source = np.array(source, dtype=np.float64) - - # check dimensions - source_dim = source.shape[1] - - # check interpolation and extrapolation - - ## single dimension - if source_dim == 2: - # possible interpolation values: linear, polynomial, akima and spline - if interpolation is None: - interpolation = "spline" - elif interpolation.lower() not in [ - "spline", - "linear", - "polynomial", - "akima", - ]: - warnings.warn( - "Interpolation method for single dimensional functions was " - + f"set to 'spline', the {interpolation} method is not supported." - ) - interpolation = "spline" + @staticmethod + def _validate_interpolation_and_extrapolation( + inputs: list[str], interpolation: str, extrapolation: str, source: np.ndarray + ): + source_dim = source.shape[1] + # check interpolation and extrapolation + ## single dimension (1D Functions) + if source_dim == 2: + # possible interpolation values: linear, polynomial, akima and spline + if interpolation is None: + interpolation = "spline" + elif interpolation.lower() not in [ + "spline", + "linear", + "polynomial", + "akima", + ]: + warnings.warn( + "Interpolation method set to 'spline' because the " + f"{interpolation} method is not supported." + ) + interpolation = "spline" # possible extrapolation values: constant, natural, zero - if extrapolation is None: - extrapolation = "constant" - elif extrapolation.lower() not in ["constant", "natural", "zero"]: - warnings.warn( - "Extrapolation method for single dimensional functions was " - + f"set to 'constant', the {extrapolation} method is not supported." - ) - extrapolation = "constant" - - ## multiple dimensions - elif source_dim > 2: - # check for inputs and outputs - if inputs == ["Scalar"]: - inputs = [f"Input {i+1}" for i in range(source_dim - 1)] - - if interpolation not in [None, "shepard"]: - warnings.warn( - ( - "Interpolation method for multidimensional functions was" - "set to 'shepard', other methods are not supported yet." - ), - ) - interpolation = "shepard" - - if extrapolation not in [None, "natural"]: - warnings.warn( - "Extrapolation method for multidimensional functions was set" - "to 'natural', other methods are not supported yet." - ) - extrapolation = "natural" + if extrapolation is None: + extrapolation = "constant" + elif extrapolation.lower() not in ["constant", "natural", "zero"]: + warnings.warn( + "Extrapolation method set to 'constant' because the " + f"{extrapolation} method is not supported." + ) + extrapolation = "constant" + + ## multiple dimensions + elif source_dim > 2: + # check for inputs and outputs + if inputs == ["Scalar"]: + inputs = [f"Input {i+1}" for i in range(source_dim - 1)] + + if interpolation not in [None, "shepard"]: + warnings.warn( + ( + "Interpolation method set to 'shepard'. Other methods " + "are not supported yet." + ), + ) + interpolation = "shepard" - # check input dimensions - in_out_dim = len(inputs) + len(outputs) - if source_dim != in_out_dim: - raise ValueError( - f"Source dimension ({source_dim}) does not match input " - + f"and output dimension ({in_out_dim})." + if extrapolation not in [None, "natural"]: + warnings.warn( + "Extrapolation method set to 'natural'. Other methods " + "are not supported yet." ) - return inputs, outputs, interpolation, extrapolation + extrapolation = "natural" + else: + raise ValueError("Source must have at least 2 dimensions.") + return inputs, interpolation, extrapolation + + @staticmethod + def _validate_source_dimensions(inputs: list, outputs: list, source: np.ndarray): + # check input dimensions + source_dim = source.shape[1] + in_out_dim = len(inputs) + len(outputs) + if source_dim != in_out_dim: + raise ValueError( + f"Source dimension ({source_dim}) does not match input " + f"and output dimension ({in_out_dim})." + ) class PiecewiseFunction(Function): From 5a356cdbdab046ac2895b348688072fce121b2dd Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 25 Mar 2024 12:57:35 -0400 Subject: [PATCH 03/71] MNT: better interpolation and extrapolation for get_value_opt --- rocketpy/mathutils/function.py | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 0d88af585..20fe27a81 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -369,7 +369,7 @@ def get_value_opt(x): # Interval found... interpolate... or extrapolate if x_min <= x <= x_max: # Interpolate - x_interval = x_interval if x_interval != 0 else 1 + x_interval = max(x_interval, 1) a = coeffs[:, x_interval - 1] x = x - x_data[x_interval - 1] y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] @@ -378,8 +378,12 @@ def get_value_opt(x): if extrapolation == 0: # Extrapolation == zero y = 0 elif extrapolation == 1: # Extrapolation == natural - a = coeffs[:, 0] if x < x_min else coeffs[:, -1] - x = x - x_data[0] if x < x_min else x - x_data[-2] + if x < x_min: + a = coeffs[:, 0] + x = x - x_data[0] + else: + a = coeffs[:, -1] + x = x - x_data[-2] y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] else: # Extrapolation is set to constant y = y_data[0] if x < x_min else y_data[-1] @@ -392,22 +396,22 @@ def get_value_opt(x): # Interval found... interpolate... or extrapolate if x_min <= x <= x_max: # Interpolate - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - y = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[ - x_interval - 1 - ] + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + y = (x - x_left) * (dy / dx) + y_left else: # Extrapolate if extrapolation == 0: # Extrapolation == zero y = 0 elif extrapolation == 1: # Extrapolation == natural x_interval = 1 if x < x_min else -1 - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - y = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[ - x_interval - 1 - ] + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + y = (x - x_left) * (dy / dx) + y_left else: # Extrapolation is set to constant y = y_data[0] if x < x_min else y_data[-1] return y From cbaa145d9741fccf297b733b7b85b21b06c515d6 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 25 Mar 2024 12:58:07 -0400 Subject: [PATCH 04/71] MNT: refactors Function.__mul__() --- rocketpy/mathutils/function.py | 87 +++++++++++++++------------------- 1 file changed, 39 insertions(+), 48 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 20fe27a81..fbacd6846 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2054,54 +2054,45 @@ def __mul__(self, other): result : Function A Function object which gives the result of self(x)*other(x). """ - # If other is Function try... - try: - # Check if Function objects source is array or callable - # Check if Function objects have the same domain discretization - if ( - isinstance(other.source, np.ndarray) - and isinstance(self.source, np.ndarray) - and self.__dom_dim__ == other.__dom_dim__ - and np.array_equal(self.x_array, other.x_array) - ): - # Operate on grid values - ys = self.y_array * other.y_array - xs = self.x_array - source = np.concatenate(([xs], [ys])).transpose() - # Retrieve inputs, outputs and interpolation - inputs = self.__inputs__[:] - outputs = self.__outputs__[0] + "*" + other.__outputs__[0] - outputs = "(" + outputs + ")" - interpolation = self.__interpolation__ - extrapolation = self.__extrapolation__ - # Create new Function object - return Function(source, inputs, outputs, interpolation, extrapolation) - else: - return Function(lambda x: (self.get_value(x) * other(x))) - # If other is Float except... - except AttributeError: - if isinstance(other, NUMERICAL_TYPES): - # Check if Function object source is array or callable - if isinstance(self.source, np.ndarray): - # Operate on grid values - ys = self.y_array * other - xs = self.x_array - source = np.concatenate(([xs], [ys])).transpose() - # Retrieve inputs, outputs and interpolation - inputs = self.__inputs__[:] - outputs = self.__outputs__[0] + "*" + str(other) - outputs = "(" + outputs + ")" - interpolation = self.__interpolation__ - extrapolation = self.__extrapolation__ - # Create new Function object - return Function( - source, inputs, outputs, interpolation, extrapolation - ) - else: - return Function(lambda x: (self.get_value(x) * other)) - # Or if it is just a callable - elif callable(other): - return Function(lambda x: (self.get_value(x) * other(x))) + self_source_is_array = isinstance(self.source, np.ndarray) + other_source_is_array = ( + isinstance(other.source, np.ndarray) + if isinstance(other, Function) + else False + ) + inputs = self.__inputs__[:] + interp = self.__interpolation__ + extrap = self.__extrapolation__ + + if ( + self_source_is_array + and other_source_is_array + and self.__dom_dim__ == other.__dom_dim__ + and np.array_equal(self.x_array, other.x_array) + ): + ys = self.y_array * other.y_array + xs = self.x_array + source = np.concatenate(([xs], [ys])).transpose() + outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})" + return Function(source, inputs, outputs, interp, extrap) + elif isinstance(other, NUMERICAL_TYPES): + if not self_source_is_array: + return Function(lambda x: (self.get_value_opt(x) * other), inputs) + ys = np.multiply(self.y_array, other) + xs = self.x_array + source = np.concatenate(([xs], [ys])).transpose() + 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") def __rmul__(self, other): """Multiplies 'other' by a Function object and returns a new Function From 9a637a834e443736e7fa988830d1586aaf1f9201 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 25 Mar 2024 12:58:54 -0400 Subject: [PATCH 05/71] TST: fix tests not passing --- tests/fixtures/environment/environment_fixtures.py | 5 +++-- tests/test_flight.py | 2 +- tests/unit/test_utilities.py | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/fixtures/environment/environment_fixtures.py b/tests/fixtures/environment/environment_fixtures.py index 8949f9973..851be3203 100644 --- a/tests/fixtures/environment/environment_fixtures.py +++ b/tests/fixtures/environment/environment_fixtures.py @@ -1,6 +1,7 @@ from datetime import datetime, timedelta import pytest + from rocketpy import Environment, EnvironmentAnalysis @@ -54,8 +55,8 @@ def env_analysis(): EnvironmentAnalysis """ env_analysis = EnvironmentAnalysis( - start_date=datetime.datetime(2019, 10, 23), - end_date=datetime.datetime(2021, 10, 23), + start_date=datetime(2019, 10, 23), + end_date=datetime(2021, 10, 23), latitude=39.3897, longitude=-8.28896388889, start_hour=6, diff --git a/tests/test_flight.py b/tests/test_flight.py index db882e185..a4b03d101 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -729,7 +729,7 @@ def test_velocities(flight_calisto_custom_wind, flight_time, expected_values): ("t_initial", (1.6542528, 0.65918, -0.067107)), ("out_of_rail_time", (5.05334, 2.01364, -1.7541)), ("apogee_time", (2.35291, -1.8275, -0.87851)), - ("t_final", (0, 0, 159.2212)), + ("t_final", (0, 0, 159.3292416824044)), ], ) def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values): diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py index b07064906..2d5d4c366 100644 --- a/tests/unit/test_utilities.py +++ b/tests/unit/test_utilities.py @@ -153,7 +153,7 @@ def test_fin_flutter_analysis(flight_calisto_custom_wind): assert np.isclose(flutter_mach(np.inf), 1.0048188594647927, atol=5e-3) assert np.isclose(safety_factor(0), 64.78797, atol=5e-3) assert np.isclose(safety_factor(10), 2.1948620401502072, atol=5e-3) - assert np.isclose(safety_factor(np.inf), 61.64222220469224, atol=5e-3) + assert np.isclose(safety_factor(np.inf), 61.62131849619694, atol=5e-3) def test_flutter_prints(flight_calisto_custom_wind): From bb293aa5c07d13c9aa4a1f2346922dbcad5e5ef6 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 25 Mar 2024 13:21:20 -0400 Subject: [PATCH 06/71] DOC: remove forgotten type hints and update docstrings --- rocketpy/mathutils/function.py | 90 ++++++++++++++++++++++++++++------ 1 file changed, 75 insertions(+), 15 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index fbacd6846..c4c079414 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2988,8 +2988,8 @@ def _check_user_input( Returns ------- tuple - A tuple containing the processed inputs, outputs, interpolation, and - extrapolation parameters. + A tuple containing the processed source, inputs, outputs, + interpolation, and extrapolation parameters. Raises ------ @@ -2997,13 +2997,6 @@ def _check_user_input( If the dimensionality of the source does not match the combined dimensions of inputs and outputs. If the outputs list has more than one element. - TypeError - If the source is not a list, np.ndarray, Function object, str or - Path. - Warning - If inputs or outputs do not match for a Function source, or if - defaults are used for inputs, interpolation,and extrapolation for a - multidimensional source. Examples -------- @@ -3072,7 +3065,27 @@ def source_function(_): return source, inputs, outputs, interpolation, extrapolation @staticmethod - def _validate_inputs_outputs(inputs, outputs): # None | st | list[str] + def _validate_inputs_outputs(inputs, outputs): + """Used to validate the inputs and outputs parameters for creating a + Function object. It sets default values if they are not provided. + + Parameters + ---------- + inputs : str, list of str, None + The name(s) of the input variable(s). If None, defaults to "Scalar". + outputs : + The name of the output variables. If None, defaults to "Scalar". + + Returns + ------- + tuple + A tuple containing the validated inputs and outputs parameters. + + Raises + ------ + ValueError + If the output has more than one element. + """ if inputs is None: inputs = ["Scalar"] if outputs is None: @@ -3091,10 +3104,41 @@ def _validate_inputs_outputs(inputs, outputs): # None | st | list[str] @staticmethod def _validate_interpolation_and_extrapolation( - inputs: list[str], interpolation: str, extrapolation: str, source: np.ndarray + inputs, interpolation, extrapolation, source ): + """Used to validate the interpolation and extrapolation methods for + creating a Function object. It sets default values for interpolation + and extrapolation if they are not provided or if they are not supported + for the given source. The inputs and outputs may be modified if the + source is multidimensional. + + Parameters + ---------- + inputs : list of strings + List of inputs, each input is a string. Example: ['x', 'y'] + interpolation : str, None + The type of interpolation to use. The default method is 'spline'. + Currently supported values are 'spline', 'linear', 'polynomial', + 'akima', and 'shepard'. + extrapolation : str, None + The type of extrapolation to use. Currently supported values are + 'constant', 'natural', and 'zero'. The default method is 'constant'. + source : np.ndarray + The source data of the Function object. This has to be a numpy + array. + + Returns + ------- + tuple + A tuple with the validated inputs, interpolation, and extrapolation + parameters (inputs, interpolation, extrapolation). + + Raises + ------ + ValueError + If the source has less than 2 dimensions. + """ source_dim = source.shape[1] - # check interpolation and extrapolation ## single dimension (1D Functions) if source_dim == 2: # possible interpolation values: linear, polynomial, akima and spline @@ -3124,7 +3168,6 @@ def _validate_interpolation_and_extrapolation( ## multiple dimensions elif source_dim > 2: - # check for inputs and outputs if inputs == ["Scalar"]: inputs = [f"Input {i+1}" for i in range(source_dim - 1)] @@ -3148,8 +3191,25 @@ def _validate_interpolation_and_extrapolation( return inputs, interpolation, extrapolation @staticmethod - def _validate_source_dimensions(inputs: list, outputs: list, source: np.ndarray): - # check input dimensions + def _validate_source_dimensions(inputs, outputs, source): + """Used to check whether the source dimensions match the inputs and + outputs. + + Parameters + ---------- + inputs : list of strings + List of inputs, each input is a string. Example: ['x', 'y'] + outputs : list of strings + List of outputs, each output is a string. Example: ['z'] + source : np.ndarray + The source data of the Function object. This has to be a numpy + array. + + Raises + ------ + ValueError + In case the source dimensions do not match the inputs and outputs. + """ source_dim = source.shape[1] in_out_dim = len(inputs) + len(outputs) if source_dim != in_out_dim: From a22d38c8deb8155e387c98579de888d99bfa5eb7 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 02:48:45 -0400 Subject: [PATCH 07/71] MNT: Refactor variable assignment in Function.__mul__ --- rocketpy/mathutils/function.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index c4c079414..b2f4d1dfc 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2067,20 +2067,15 @@ def __mul__(self, other): if ( self_source_is_array and other_source_is_array - and self.__dom_dim__ == other.__dom_dim__ and np.array_equal(self.x_array, other.x_array) ): - ys = self.y_array * other.y_array - xs = self.x_array - source = np.concatenate(([xs], [ys])).transpose() + 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): if not self_source_is_array: return Function(lambda x: (self.get_value_opt(x) * other), inputs) - ys = np.multiply(self.y_array, other) - xs = self.x_array - source = np.concatenate(([xs], [ys])).transpose() + source = np.column_stack((self.x_array, np.multiply(self.y_array, other))) outputs = f"({self.__outputs__[0]}*{other})" return Function( source, From b23e31268e085505b755ddd3e5a414a1b9cc3a34 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 02:51:27 -0400 Subject: [PATCH 08/71] MNT: minor adjustments to Function setters: - set_interpolation - set_extrapolation - set_get_value_opt --- rocketpy/mathutils/function.py | 35 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index b2f4d1dfc..a2dc44e79 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -262,7 +262,7 @@ def set_source(self, source): # Finally set source, update extrapolation and interpolation self.source = source - self.set_extrapolation(extrapolation) + self.__extrapolation__ = extrapolation # to avoid calling set_get_value_opt self.set_interpolation(interpolation) return self @@ -303,7 +303,6 @@ def set_interpolation(self, method="spline"): ------- self : Function """ - # Set interpolation method self.__interpolation__ = method # Spline, akima and polynomial need data processing # Shepard, and linear do not @@ -314,9 +313,7 @@ def set_interpolation(self, method="spline"): elif method == "akima": self.__interpolate_akima__() - # Set get_value_opt self.set_get_value_opt() - return self def set_extrapolation(self, method="constant"): @@ -337,6 +334,7 @@ def set_extrapolation(self, method="constant"): The Function object. """ self.__extrapolation__ = method + self.set_get_value_opt() return self def set_get_value_opt(self): @@ -357,7 +355,7 @@ def set_get_value_opt(self): extrapolation = self.EXTRAPOLATION_TYPES[self.__extrapolation__] except KeyError as err: raise ValueError( - f"Invalid extrapolation type {self.__extrapolation__}" + f"Invalid extrapolation type '{self.__extrapolation__}'" ) from err # Crete method to interpolate this info for each interpolation type @@ -365,10 +363,9 @@ def set_get_value_opt(self): coeffs = self.__spline_coefficients__ def get_value_opt(x): - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate if x_min <= x <= x_max: # Interpolate + x_interval = np.searchsorted(x_data, x) x_interval = max(x_interval, 1) a = coeffs[:, x_interval - 1] x = x - x_data[x_interval - 1] @@ -392,10 +389,9 @@ def get_value_opt(x): elif self.__interpolation__ == "linear": def get_value_opt(x): - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate if x_min <= x <= x_max: # Interpolate + x_interval = np.searchsorted(x_data, x) x_left = x_data[x_interval - 1] y_left = y_data[x_interval - 1] dx = float(x_data[x_interval] - x_left) @@ -420,10 +416,9 @@ def get_value_opt(x): coeffs = np.array(self.__akima_coefficients__) def get_value_opt(x): - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate if x_min <= x <= x_max: # Interpolate + x_interval = np.searchsorted(x_data, x) x_interval = x_interval if x_interval != 0 else 1 a = coeffs[4 * x_interval - 4 : 4 * x_interval] y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] @@ -445,9 +440,7 @@ def get_value_opt(x): # Interpolate... or extrapolate if x_min <= x <= x_max: # Interpolate - y = 0 - for i, coef in enumerate(coeffs): - y += coef * (x**i) + y = np.sum(coeffs * x ** np.arange(len(coeffs))) else: # Extrapolate if extrapolation == 0: # Extrapolation == zero @@ -467,7 +460,13 @@ def get_value_opt_multiple(*args): get_value_opt = get_value_opt_multiple - self.get_value_opt = get_value_opt + try: + self.get_value_opt = get_value_opt + except UnboundLocalError: + warnings.warn( + "Cannot set the get_value_opt method when interpolation is " + f"{self.__interpolation}. Try using the set_interpolation method first." + ) return self def set_discrete( @@ -535,7 +534,7 @@ def set_discrete( ys = func.get_value(xs.tolist()) if one_by_one else func.get_value(xs) func.set_source(np.concatenate(([xs], [ys])).transpose()) func.set_interpolation(interpolation) - func.set_extrapolation(extrapolation) + func.__extrapolation__ = extrapolation # avoid calling set_get_value_opt elif func.__dom_dim__ == 2: lower = 2 * [lower] if isinstance(lower, (int, float)) else lower upper = 2 * [upper] if isinstance(upper, (int, float)) else upper @@ -674,7 +673,7 @@ def set_discrete_based_on_model( ) func.set_interpolation(interp) - func.set_extrapolation(extrap) + func.__extrapolation__ = extrap # avoid calling set_get_value_opt return func @@ -733,7 +732,7 @@ def reset( if interpolation is not None and interpolation != self.__interpolation__: self.set_interpolation(interpolation) if extrapolation is not None and extrapolation != self.__extrapolation__: - self.set_extrapolation(extrapolation) + self.__extrapolation__ = extrapolation self.set_title(title) From c93480f43294640543e6c3f9e4dada54eca632c6 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 02:53:00 -0400 Subject: [PATCH 09/71] MNT: invokes get_value_opt inside get_value - This reduces code lines and improve redability. - performance is similar, but now there's a map() in the code instead of a for loop --- rocketpy/mathutils/function.py | 104 +++------------------------------ 1 file changed, 9 insertions(+), 95 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index a2dc44e79..c887a904a 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -866,103 +866,17 @@ def get_value(self, *args): elif self.__interpolation__ == "shepard": return self.__interpolate_shepard__(args) - # Returns value for polynomial interpolation function type - elif self.__interpolation__ == "polynomial": - if isinstance(args[0], (int, float)): - args = [list(args)] - x = np.array(args[0]) - x_data = self.x_array - y_data = self.y_array - x_min, x_max = self.x_initial, self.x_final - coeffs = self.__polynomial_coefficients__ - matrix = np.zeros((len(args[0]), coeffs.shape[0])) - for i in range(coeffs.shape[0]): - matrix[:, i] = x**i - ans = matrix.dot(coeffs).tolist() - for i, _ in enumerate(x): - if not x_min <= x[i] <= x_max: - if self.__extrapolation__ == "constant": - ans[i] = y_data[0] if x[i] < x_min else y_data[-1] - elif self.__extrapolation__ == "zero": - ans[i] = 0 - return ans if len(ans) > 1 else ans[0] - # Returns value for spline, akima or linear interpolation function type - elif self.__interpolation__ in ["spline", "akima", "linear"]: + # Returns value for other interpolation type + else: # interpolation is "polynomial", "spline", "akima" or "linear" if isinstance(args[0], (int, float, complex, np.integer)): args = [list(args)] - x = list(args[0]) - x_data = self.x_array - y_data = self.y_array - x_intervals = np.searchsorted(x_data, x) - x_min, x_max = self.x_initial, self.x_final - if self.__interpolation__ == "spline": - coeffs = self.__spline_coefficients__ - for i, _ in enumerate(x): - if x[i] == x_min or x[i] == x_max: - x[i] = y_data[x_intervals[i]] - elif x_min < x[i] < x_max or (self.__extrapolation__ == "natural"): - if not x_min < x[i] < x_max: - a = coeffs[:, 0] if x[i] < x_min else coeffs[:, -1] - x[i] = ( - x[i] - x_data[0] if x[i] < x_min else x[i] - x_data[-2] - ) - else: - a = coeffs[:, x_intervals[i] - 1] - x[i] = x[i] - x_data[x_intervals[i] - 1] - x[i] = a[3] * x[i] ** 3 + a[2] * x[i] ** 2 + a[1] * x[i] + a[0] - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x[i] = 0 - else: # Extrapolation is set to constant - x[i] = y_data[0] if x[i] < x_min else y_data[-1] - elif self.__interpolation__ == "linear": - for i, _ in enumerate(x): - # Interval found... interpolate... or extrapolate - inter = x_intervals[i] - if x_min <= x[i] <= x_max: - # Interpolate - dx = float(x_data[inter] - x_data[inter - 1]) - dy = float(y_data[inter] - y_data[inter - 1]) - x[i] = (x[i] - x_data[inter - 1]) * (dy / dx) + y_data[ - inter - 1 - ] - else: - # Extrapolate - if self.__extrapolation__ == "zero": # Extrapolation == zero - x[i] = 0 - elif ( - self.__extrapolation__ == "natural" - ): # Extrapolation == natural - inter = 1 if x[i] < x_min else -1 - dx = float(x_data[inter] - x_data[inter - 1]) - dy = float(y_data[inter] - y_data[inter - 1]) - x[i] = (x[i] - x_data[inter - 1]) * (dy / dx) + y_data[ - inter - 1 - ] - else: # Extrapolation is set to constant - x[i] = y_data[0] if x[i] < x_min else y_data[-1] - else: - coeffs = self.__akima_coefficients__ - for i, _ in enumerate(x): - if x[i] == x_min or x[i] == x_max: - x[i] = y_data[x_intervals[i]] - elif x_min < x[i] < x_max or (self.__extrapolation__ == "natural"): - if not x_min < x[i] < x_max: - a = coeffs[:4] if x[i] < x_min else coeffs[-4:] - else: - a = coeffs[4 * x_intervals[i] - 4 : 4 * x_intervals[i]] - x[i] = a[3] * x[i] ** 3 + a[2] * x[i] ** 2 + a[1] * x[i] + a[0] - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x[i] = 0 - else: # Extrapolation is set to constant - x[i] = y_data[0] if x[i] < x_min else y_data[-1] - if isinstance(args[0], np.ndarray): - return np.array(x) - else: - return x if len(x) > 1 else x[0] + + x = list(args[0]) + x = list(map(self.get_value_opt, x)) + if isinstance(args[0], np.ndarray): + return np.array(x) + else: + return x if len(x) > 1 else x[0] def __getitem__(self, args): """Returns item of the Function source. If the source is not an array, From 99d89482513044a20a26753d7debde143390b7bc Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 02:53:27 -0400 Subject: [PATCH 10/71] BUG: Fix wind direction calculation in Environment class --- 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 4d9858af3..1ecfbb538 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -1507,7 +1507,7 @@ def wind_heading_func(h): ) def wind_direction(h): - return (wind_heading(h) - 180) % 360 + return (wind_heading_func(h) - 180) % 360 self.wind_direction = Function( wind_direction, From 9b47dfe6b660312770a4308f9b29e2034517b384 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 03:02:41 -0400 Subject: [PATCH 11/71] MNT: use local variable instead of class variable in TimeNodes.merge() --- rocketpy/simulation/flight.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 75c31917d..de34076d6 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3664,19 +3664,18 @@ def sort(self): def merge(self): # Initialize temporary list - self.tmp_list = [self.list[0]] - self.copy_list = self.list[1:] + tmp_list = [self.list[0]] # Iterate through all other time nodes - for node in self.copy_list: + for node in self.list[1:]: # If there is already another node with similar time: merge - if abs(node.t - self.tmp_list[-1].t) < 1e-7: - self.tmp_list[-1].parachutes += node.parachutes - self.tmp_list[-1].callbacks += node.callbacks + if abs(node.t - tmp_list[-1].t) < 1e-7: + tmp_list[-1].parachutes += node.parachutes + tmp_list[-1].callbacks += node.callbacks # Add new node to tmp list if there is none with the same time else: - self.tmp_list.append(node) + tmp_list.append(node) # Save tmp list to permanent - self.list = self.tmp_list + self.list = tmp_list def flush_after(self, index): del self.list[index + 1 :] From 6deb23372bf9d6f97536f4cc533bc39ab72b9fc6 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 03:05:08 -0400 Subject: [PATCH 12/71] MNT: speedup Flight._calculate_pressure_signal method - Initializing a Function object is complicated. Instead, let's sum two Functions --- rocketpy/simulation/flight.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index de34076d6..b13f923eb 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3074,17 +3074,13 @@ def _calculate_pressure_signal(self): "Pressure - Without Noise (Pa)", "linear", ) - parachute.noisy_pressure_signal_function = Function( - parachute.noisy_pressure_signal, - "Time (s)", - "Pressure - With Noise (Pa)", - "linear", - ) parachute.noise_signal_function = Function( parachute.noise_signal, "Time (s)", "Pressure Noise (Pa)", "linear" ) - - return None + parachute.noisy_pressure_signal_function = ( + parachute.clean_pressure_signal_function + + parachute.noise_signal_function + ) def post_process(self, interpolation="spline", extrapolation="natural"): """This method is **deprecated** and is only kept here for backwards From 939d294bb470fd9cd71878596c815237fe2acf6c Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 Mar 2024 03:12:48 -0400 Subject: [PATCH 13/71] TST: fix test not passing --- tests/unit/test_utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py index 2d5d4c366..33942b445 100644 --- a/tests/unit/test_utilities.py +++ b/tests/unit/test_utilities.py @@ -153,7 +153,7 @@ def test_fin_flutter_analysis(flight_calisto_custom_wind): assert np.isclose(flutter_mach(np.inf), 1.0048188594647927, atol=5e-3) assert np.isclose(safety_factor(0), 64.78797, atol=5e-3) assert np.isclose(safety_factor(10), 2.1948620401502072, atol=5e-3) - assert np.isclose(safety_factor(np.inf), 61.62131849619694, atol=5e-3) + assert np.isclose(safety_factor(np.inf), 61.64222220697017, atol=5e-3) def test_flutter_prints(flight_calisto_custom_wind): From 26cc6a1b9bb50a3352f1bb0bad62bab6c63c36a5 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 28 Mar 2024 17:46:38 +0100 Subject: [PATCH 14/71] ENH: rework setters and validation --- rocketpy/mathutils/function.py | 632 +++++++++++++++------------------ 1 file changed, 286 insertions(+), 346 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index c887a904a..86135e180 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -20,6 +20,14 @@ from ..tools import cached_property NUMERICAL_TYPES = (float, int, complex, np.ndarray, np.integer, np.floating) +INTERPOLATION_TYPES = { + "linear": 0, + "polynomial": 1, + "akima": 2, + "spline": 3, + "shepard": 4, +} +EXTRAPOLATION_TYPES = {"zero": 0, "natural": 1, "constant": 2} class Function: @@ -110,27 +118,25 @@ def __init__( (II) Fields in CSV files may be enclosed in double quotes. If fields are not quoted, double quotes should not appear inside them. """ - - inputs, outputs = Function._validate_inputs_outputs(inputs, outputs) - # initialize variables to avoid errors when being called by other methods - self.get_value_opt = None self.__polynomial_coefficients__ = None self.__akima_coefficients__ = None self.__spline_coefficients__ = None - self.EXTRAPOLATION_TYPES = { - "zero": 0, - "natural": 1, - "constant": 2, - } - - # store variables - self.set_inputs(inputs) - self.set_outputs(outputs) + + # initialize parameters + self.source = source + self.__inputs__ = inputs + self.__outputs__ = outputs self.__interpolation__ = interpolation self.__extrapolation__ = extrapolation - self.set_source(source) - self.set_title(title) + self.title = title + self.__img_dim__ = 1 # always 1, here for backwards compatibility + + # args must be passed from self. + self.set_source(self.source) + self.set_inputs(self.__inputs__) + self.set_outputs(self.__outputs__) + self.set_title(self.title) # Define all set methods def set_inputs(self, inputs): @@ -145,8 +151,7 @@ def set_inputs(self, inputs): ------- self : Function """ - self.__inputs__ = [inputs] if isinstance(inputs, str) else list(inputs) - self.__dom_dim__ = len(self.__inputs__) + self.__inputs__ = self._validate_inputs(inputs) return self def set_outputs(self, outputs): @@ -161,8 +166,7 @@ def set_outputs(self, outputs): ------- self : Function """ - self.__outputs__ = [outputs] if isinstance(outputs, str) else list(outputs) - self.__img_dim__ = len(self.__outputs__) + self.__outputs__ = self._validate_outputs(outputs) return self def set_source(self, source): @@ -211,20 +215,10 @@ def set_source(self, source): self : Function Returns the Function instance. """ - source, inputs, outputs, interpolation, extrapolation = self._check_user_input( - source, - self.__inputs__, - self.__outputs__, - self.__interpolation__, - self.__extrapolation__, - ) - # updates inputs and outputs (could be modified due to csv headers) - self.set_inputs(inputs) - self.set_outputs(outputs) + source = self._validate_source(source) # Handle callable source or number source if callable(source): - self.source = source self.get_value_opt = source self.__interpolation__ = None self.__extrapolation__ = None @@ -232,13 +226,16 @@ def set_source(self, source): # Set arguments name and domain dimensions parameters = signature(source).parameters self.__dom_dim__ = len(parameters) - if self.__inputs__ == ["Scalar"]: + if self.__inputs__ is None: self.__inputs__ = list(parameters) # Handle ndarray source else: + # Evaluate dimension + self.__dom_dim__ = source.shape[1] - 1 + # Check to see if dimensions match incoming data set - new_total_dim = len(source[0, :]) + new_total_dim = source.shape[1] old_total_dim = self.__dom_dim__ + self.__img_dim__ # If they don't, update default values or throw error @@ -250,20 +247,25 @@ def set_source(self, source): # if Function is 1D, sort source by x. If 2D, set z if self.__dom_dim__ == 1: source = source[source[:, 0].argsort()] + self.x_array = source[:, 0] + self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] + self.y_array = source[:, 1] + self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] + self.get_value_opt = self._get_value_opt_1d elif self.__dom_dim__ == 2: + self.x_array = source[:, 0] + self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] + self.y_array = source[:, 1] + self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] self.z_array = source[:, 2] self.z_initial, self.z_final = self.z_array[0], self.z_array[-1] + self.get_value_opt = self._get_value_opt_nd + else: + self.get_value_opt = self._get_value_opt_nd - # Set x and y arrays (common for 1D or multivariate) - self.x_array = source[:, 0] - self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] - self.y_array = source[:, 1] - self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] - - # Finally set source, update extrapolation and interpolation - self.source = source - self.__extrapolation__ = extrapolation # to avoid calling set_get_value_opt - self.set_interpolation(interpolation) + self.source = source + self.set_interpolation(self.__interpolation__) + self.set_extrapolation(self.__extrapolation__) return self @cached_property @@ -303,18 +305,27 @@ def set_interpolation(self, method="spline"): ------- self : Function """ - self.__interpolation__ = method + if not callable(self.source): + self.__interpolation__ = self._validate_interpolation(method) + self._update_interpolation_coefficients(self.__interpolation__) + self._set_interpolation_function() + return self + + def _update_interpolation_coefficients(self, method): + """Update interpolation coefficients for the given method.""" # Spline, akima and polynomial need data processing # Shepard, and linear do not - if method == "spline": - self.__interpolate_spline__() - elif method == "polynomial": + if method == "polynomial": self.__interpolate_polynomial__() + self._coeffs = self.__polynomial_coefficients__ elif method == "akima": self.__interpolate_akima__() - - self.set_get_value_opt() - return self + self._coeffs = self.__akima_coefficients__ + elif method == "spline" or method is None: + self.__interpolate_spline__() + self._coeffs = self.__spline_coefficients__ + else: + self._coeffs = [] def set_extrapolation(self, method="constant"): """Set extrapolation behavior of data set. @@ -333,142 +344,159 @@ def set_extrapolation(self, method="constant"): self : Function The Function object. """ - self.__extrapolation__ = method - self.set_get_value_opt() + if not callable(self.source): + self.__extrapolation__ = self._validate_extrapolation(method) + self._set_extrapolation_function() return self - def set_get_value_opt(self): - """Crates a method that evaluates interpolations rather quickly - when compared to other options available, such as just calling - the object instance or calling ``Function.get_value`` directly. See - ``Function.get_value_opt`` for documentation. + def _set_interpolation_function(self): + """Return a function that interpolates the Function.""" + interpolation = INTERPOLATION_TYPES[self.__interpolation__] + if interpolation == 0: # linear + + def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = np.searchsorted(x_data, x) + x_left = x_data[x_interval - 1] + y_left = y_data[x_interval - 1] + dx = float(x_data[x_interval] - x_left) + dy = float(y_data[x_interval] - y_left) + y = (x - x_left) * (dy / dx) + y_left + return y - Returns - ------- - self : Function - """ - # Retrieve general info - x_data = self.x_array - y_data = self.y_array - x_min, x_max = self.x_initial, self.x_final - try: - extrapolation = self.EXTRAPOLATION_TYPES[self.__extrapolation__] - except KeyError as err: - raise ValueError( - f"Invalid extrapolation type '{self.__extrapolation__}'" - ) from err + self._interpolation_function = linear_interpolation - # Crete method to interpolate this info for each interpolation type - if self.__interpolation__ == "spline": - coeffs = self.__spline_coefficients__ + elif interpolation == 1: # polynomial - def get_value_opt(x): - if x_min <= x <= x_max: - # Interpolate - x_interval = np.searchsorted(x_data, x) - x_interval = max(x_interval, 1) - a = coeffs[:, x_interval - 1] - x = x - x_data[x_interval - 1] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - if x < x_min: - a = coeffs[:, 0] - x = x - x_data[0] - else: - a = coeffs[:, -1] - x = x - x_data[-2] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] + def polynomial_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + y = np.sum(coeffs * x ** np.arange(len(coeffs))) + return y + + self._interpolation_function = polynomial_interpolation + + elif interpolation == 2: # akima + + def akima_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = np.searchsorted(x_data, x) + x_interval = x_interval if x_interval != 0 else 1 + a = coeffs[4 * x_interval - 4 : 4 * x_interval] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + return y + + self._interpolation_function = akima_interpolation + + elif interpolation == 3: # spline + + def spline_interpolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = np.searchsorted(x_data, x) + x_interval = max(x_interval, 1) + a = coeffs[:, x_interval - 1] + x = x - x_data[x_interval - 1] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] return y - elif self.__interpolation__ == "linear": + self._interpolation_function = spline_interpolation + + elif interpolation == 4: # shepard does not use interpolation function + self._interpolation_function = None + + def _set_extrapolation_function(self): + """Return a function that extrapolates the Function.""" + interpolation = INTERPOLATION_TYPES[self.__interpolation__] + extrapolation = EXTRAPOLATION_TYPES[self.__extrapolation__] + + if interpolation == 4: # shepard does not use extrapolation function + self._extrapolation_function = None + + elif extrapolation == 0: # zero + + def zero_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + return 0 - def get_value_opt(x): - if x_min <= x <= x_max: - # Interpolate - x_interval = np.searchsorted(x_data, x) + self._extrapolation_function = zero_extrapolation + elif extrapolation == 1: # natural + if interpolation == 0: # linear + + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + x_interval = 1 if x < x_min else -1 x_left = x_data[x_interval - 1] y_left = y_data[x_interval - 1] dx = float(x_data[x_interval] - x_left) dy = float(y_data[x_interval] - y_left) y = (x - x_left) * (dy / dx) + y_left - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - x_interval = 1 if x < x_min else -1 - x_left = x_data[x_interval - 1] - y_left = y_data[x_interval - 1] - dx = float(x_data[x_interval] - x_left) - dy = float(y_data[x_interval] - y_left) - y = (x - x_left) * (dy / dx) + y_left - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y + return y + + elif interpolation == 1: # polynomial + + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + y = np.sum(coeffs * x ** np.arange(len(coeffs))) + return y - elif self.__interpolation__ == "akima": - coeffs = np.array(self.__akima_coefficients__) + elif interpolation == 2: # akima - def get_value_opt(x): - if x_min <= x <= x_max: - # Interpolate - x_interval = np.searchsorted(x_data, x) - x_interval = x_interval if x_interval != 0 else 1 - a = coeffs[4 * x_interval - 4 : 4 * x_interval] + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + a = coeffs[:4] if x < x_min else coeffs[-4:] y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - a = coeffs[:4] if x < x_min else coeffs[-4:] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y + return y - elif self.__interpolation__ == "polynomial": - coeffs = self.__polynomial_coefficients__ + elif interpolation == 3: # spline - def get_value_opt(x): - # Interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - y = np.sum(coeffs * x ** np.arange(len(coeffs))) - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - y = 0 - for i, coef in enumerate(coeffs): - y += coef * (x**i) - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y + def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + if x < x_min: + a = coeffs[:, 0] + x = x - x_data[0] + else: + a = coeffs[:, -1] + x = x - x_data[-2] + y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] + return y - elif self.__interpolation__ == "shepard": - # change the function's name to avoid mypy's error - def get_value_opt_multiple(*args): - return self.__interpolate_shepard__(args) + self._extrapolation_function = natural_extrapolation + elif extrapolation == 2: # constant - get_value_opt = get_value_opt_multiple + def constant_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): + return y_data[0] if x < x_min else y_data[-1] - try: - self.get_value_opt = get_value_opt - except UnboundLocalError: - warnings.warn( - "Cannot set the get_value_opt method when interpolation is " - f"{self.__interpolation}. Try using the set_interpolation method first." - ) + self._extrapolation_function = constant_extrapolation + + def set_get_value_opt(self): + """Crates a method that evaluates interpolations rather quickly + when compared to other options available, such as just calling + the object instance or calling ``Function.get_value`` directly. See + ``Function.get_value_opt`` for documentation. + + Returns + ------- + self : Function + """ + if callable(self.source): + self.get_value_opt = self.source + elif self.__dom_dim__ == 1: + self.get_value_opt = self._get_value_opt_1d + elif self.__dom_dim__ > 1: + self.get_value_opt = self._get_value_opt_nd return self + def _get_value_opt_1d(self, x, *args): + # Retrieve general info + x_data = self.x_array + y_data = self.y_array + x_min, x_max = self.x_initial, self.x_final + coeffs = self._coeffs + if x_min <= x <= x_max: + y = self._interpolation_function( + x, x_min, x_max, x_data, y_data, coeffs, *args + ) + else: + y = self._extrapolation_function( + x, x_min, x_max, x_data, y_data, coeffs, *args + ) + return y + + def _get_value_opt_nd(self, *args): + # always use shepard for N-D functions + y = self.__interpolate_shepard__(list(args)) + return y + def set_discrete( self, lower=0, @@ -534,7 +562,7 @@ def set_discrete( ys = func.get_value(xs.tolist()) if one_by_one else func.get_value(xs) func.set_source(np.concatenate(([xs], [ys])).transpose()) func.set_interpolation(interpolation) - func.__extrapolation__ = extrapolation # avoid calling set_get_value_opt + func.set_extrapolation(extrapolation) elif func.__dom_dim__ == 2: lower = 2 * [lower] if isinstance(lower, (int, float)) else lower upper = 2 * [upper] if isinstance(upper, (int, float)) else upper @@ -673,7 +701,7 @@ def set_discrete_based_on_model( ) func.set_interpolation(interp) - func.__extrapolation__ = extrap # avoid calling set_get_value_opt + func.set_extrapolation(extrap) return func @@ -732,7 +760,7 @@ def reset( if interpolation is not None and interpolation != self.__interpolation__: self.set_interpolation(interpolation) if extrapolation is not None and extrapolation != self.__extrapolation__: - self.__extrapolation__ = extrapolation + self.set_extrapolation(extrapolation) self.set_title(title) @@ -862,9 +890,8 @@ def get_value(self, *args): if all(isinstance(arg, Iterable) for arg in args): return [self.source(*arg) for arg in zip(*args)] - # Returns value for shepard interpolation - elif self.__interpolation__ == "shepard": - return self.__interpolate_shepard__(args) + elif self.__dom_dim__ > 1: # deals with nd functions and shepard interp + return self.get_value_opt(*args) # Returns value for other interpolation type else: # interpolation is "polynomial", "spline", "akima" or "linear" @@ -2857,79 +2884,31 @@ def savetxt( np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline) # Input validators - - @staticmethod - def _check_user_input( - source, - inputs=None, - outputs=None, - interpolation=None, - extrapolation=None, - ): - """ - Validates and processes the user input parameters for creating or - modifying a Function object. This function ensures the inputs, outputs, - interpolation, and extrapolation parameters are compatible with the - given source. It converts the source to a numpy array if necessary, sets - default values and raises warnings or errors for incompatible or - ill-defined parameters. + def _validate_source(self, source): + """Used to validate the source parameter for creating a Function object. Parameters ---------- - source : list, np.ndarray, or callable - The source data or Function object. If a list or ndarray, it should - contain numeric data. If a Function, its inputs and outputs are - checked against the provided inputs and outputs. - inputs : list of str or None - The names of the input variables. If None, defaults are generated - based on the dimensionality of the source. - outputs : str or list of str - The name(s) of the output variable(s). If a list is provided, it - must have a single element. - interpolation : str or None - The method of interpolation to be used. For multidimensional sources - it defaults to 'shepard' if not provided. - extrapolation : str or None - The method of extrapolation to be used. For multidimensional sources - it defaults to 'natural' if not provided. + source : np.ndarray, callable, str, Path, Function, list + The source data of the Function object. This can be a numpy array, + a callable function, a string or Path object to a csv or txt file, + a Function object, or a list of numbers. Returns ------- - tuple - A tuple containing the processed source, inputs, outputs, - interpolation, and extrapolation parameters. + np.ndarray, callable + The validated source parameter. Raises ------ ValueError - If the dimensionality of the source does not match the combined - dimensions of inputs and outputs. If the outputs list has more than - one element. - - Examples - -------- - >>> from rocketpy import Function - >>> source = np.array([(1, 1), (2, 4), (3, 9)]) - >>> inputs = "x" - >>> outputs = ["y"] - >>> interpolation = 'linear' - >>> extrapolation = 'zero' - >>> inputs, outputs, interpolation, extrapolation = Function._check_user_input( - ... source, inputs, outputs, interpolation, extrapolation - ... ) - >>> inputs - ['x'] - >>> outputs - ['y'] - >>> interpolation - 'linear' - >>> extrapolation - 'zero' + If the source is not a valid type or if the source is not a 2D array + or a callable function. """ if isinstance(source, Function): - source = source.get_source() + return source.get_source() - elif isinstance(source, (str, Path)): + if isinstance(source, (str, Path)): # Read csv or txt files and create a numpy array try: source = np.loadtxt(source, delimiter=",", dtype=np.float64) @@ -2941,26 +2920,25 @@ def _check_user_input( source = np.loadtxt(data, delimiter=",", dtype=np.float64) if len(source[0]) == len(header): - if inputs == ["Scalar"]: - inputs = header[:-1] - if outputs == ["Scalar"]: - outputs = [header[-1]] + if self.__inputs__ is None: + self.__inputs__ = header[:-1] + if self.__outputs__ is None: + self.__outputs__ = [header[-1]] except Exception as e: raise ValueError( "Could not read the csv or txt file to create Function source." ) from e - if isinstance(source, list): + if isinstance(source, list) or isinstance(source, np.ndarray): # Triggers an error if source is not a list of numbers source = np.array(source, dtype=np.float64) - if isinstance(source, np.ndarray): - inputs, interpolation, extrapolation = ( - Function._validate_interpolation_and_extrapolation( - inputs, interpolation, extrapolation, source + # Checks if 2D array + if len(source.shape) != 2: + raise ValueError( + "Source must be a 2D array in the form [[x1, x2 ..., xn, y], ...]." ) - ) - Function._validate_source_dimensions(inputs, outputs, source) + return source if isinstance(source, (int, float)): # Convert number source into vectorized lambda function @@ -2969,86 +2947,78 @@ def _check_user_input( def source_function(_): return temp - source = source_function - return source, inputs, outputs, interpolation, extrapolation + return source_function - @staticmethod - def _validate_inputs_outputs(inputs, outputs): - """Used to validate the inputs and outputs parameters for creating a - Function object. It sets default values if they are not provided. + # If source is a callable function + return source + + def _validate_inputs(self, inputs): + """Used to validate the inputs parameter for creating a Function object. + It sets a default value if it is not provided. Parameters ---------- - inputs : str, list of str, None + inputs : list of str, None The name(s) of the input variable(s). If None, defaults to "Scalar". - outputs : - The name of the output variables. If None, defaults to "Scalar". Returns ------- - tuple - A tuple containing the validated inputs and outputs parameters. - - Raises - ------ - ValueError - If the output has more than one element. - """ - if inputs is None: - inputs = ["Scalar"] - if outputs is None: - outputs = ["Scalar"] - # check output type and dimensions - if isinstance(inputs, str): - inputs = [inputs] - if isinstance(outputs, str): - outputs = [outputs] - elif len(outputs) > 1: + list + The validated inputs parameter. + """ + if self.__dom_dim__ == 1: + if inputs is None: + return ["Scalar"] + if isinstance(inputs, str): + return [inputs] + if isinstance(inputs, list): + if len(inputs) == 1 and isinstance(inputs[0], str): + return inputs + raise ValueError( + "Inputs must be a string or a list of strings with " + "the length of the domain dimension." + ) + if self.__dom_dim__ > 1: + if inputs is None: + return [f"Input {i+1}" for i in range(self.__dom_dim__)] + if isinstance(inputs, list): + if len(inputs) == self.__dom_dim__ and all( + isinstance(i, str) for i in inputs + ): + return inputs raise ValueError( - "Output must either be a string or have dimension 1, " - + f"it currently has dimension ({len(outputs)})." + "Inputs must be a list of strings with " + "the length of the domain dimension." ) - return inputs, outputs - @staticmethod - def _validate_interpolation_and_extrapolation( - inputs, interpolation, extrapolation, source - ): - """Used to validate the interpolation and extrapolation methods for - creating a Function object. It sets default values for interpolation - and extrapolation if they are not provided or if they are not supported - for the given source. The inputs and outputs may be modified if the - source is multidimensional. + def _validate_outputs(self, outputs): + """Used to validate the outputs parameter for creating a Function object. + It sets a default value if it is not provided. Parameters ---------- - inputs : list of strings - List of inputs, each input is a string. Example: ['x', 'y'] - interpolation : str, None - The type of interpolation to use. The default method is 'spline'. - Currently supported values are 'spline', 'linear', 'polynomial', - 'akima', and 'shepard'. - extrapolation : str, None - The type of extrapolation to use. Currently supported values are - 'constant', 'natural', and 'zero'. The default method is 'constant'. - source : np.ndarray - The source data of the Function object. This has to be a numpy - array. + outputs : str, list of str, None + The name of the output variables. If None, defaults to "Scalar". Returns ------- - tuple - A tuple with the validated inputs, interpolation, and extrapolation - parameters (inputs, interpolation, extrapolation). - - Raises - ------ - ValueError - If the source has less than 2 dimensions. + list + The validated outputs parameter. """ - source_dim = source.shape[1] - ## single dimension (1D Functions) - if source_dim == 2: + if outputs is None: + return ["Scalar"] + if isinstance(outputs, str): + return [outputs] + if isinstance(outputs, list): + if len(outputs) > 1 or not isinstance(outputs[0], str): + raise ValueError( + "Output must either be a string or a list of strings with " + + f"one item. It currently has dimension ({len(outputs)})." + ) + return outputs + + def _validate_interpolation(self, interpolation): + if self.__dom_dim__ == 1: # possible interpolation values: linear, polynomial, akima and spline if interpolation is None: interpolation = "spline" @@ -3063,8 +3033,20 @@ def _validate_interpolation_and_extrapolation( f"{interpolation} method is not supported." ) interpolation = "spline" + ## multiple dimensions + elif self.__dom_dim__ > 1: + if interpolation not in [None, "shepard"]: + warnings.warn( + ( + "Interpolation method set to 'shepard'. Only 'shepard' " + "interpolation is supported for multiple dimensions." + ), + ) + interpolation = "shepard" + return interpolation - # possible extrapolation values: constant, natural, zero + def _validate_extrapolation(self, extrapolation): + if self.__dom_dim__ == 1: if extrapolation is None: extrapolation = "constant" elif extrapolation.lower() not in ["constant", "natural", "zero"]: @@ -3075,56 +3057,14 @@ def _validate_interpolation_and_extrapolation( extrapolation = "constant" ## multiple dimensions - elif source_dim > 2: - if inputs == ["Scalar"]: - inputs = [f"Input {i+1}" for i in range(source_dim - 1)] - - if interpolation not in [None, "shepard"]: - warnings.warn( - ( - "Interpolation method set to 'shepard'. Other methods " - "are not supported yet." - ), - ) - interpolation = "shepard" - + elif self.__dom_dim__ > 1: if extrapolation not in [None, "natural"]: warnings.warn( "Extrapolation method set to 'natural'. Other methods " "are not supported yet." ) extrapolation = "natural" - else: - raise ValueError("Source must have at least 2 dimensions.") - return inputs, interpolation, extrapolation - - @staticmethod - def _validate_source_dimensions(inputs, outputs, source): - """Used to check whether the source dimensions match the inputs and - outputs. - - Parameters - ---------- - inputs : list of strings - List of inputs, each input is a string. Example: ['x', 'y'] - outputs : list of strings - List of outputs, each output is a string. Example: ['z'] - source : np.ndarray - The source data of the Function object. This has to be a numpy - array. - - Raises - ------ - ValueError - In case the source dimensions do not match the inputs and outputs. - """ - source_dim = source.shape[1] - in_out_dim = len(inputs) + len(outputs) - if source_dim != in_out_dim: - raise ValueError( - f"Source dimension ({source_dim}) does not match input " - f"and output dimension ({in_out_dim})." - ) + return extrapolation class PiecewiseFunction(Function): From 387e5326da2bf10f2f7966f8d210fd83efa09f73 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 28 Mar 2024 17:47:08 +0100 Subject: [PATCH 15/71] TST: fix function test for correct behaviour --- tests/test_function.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_function.py b/tests/test_function.py index 2ce94f691..70613671e 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -102,7 +102,8 @@ def test_setters(func_from_csv, func_2d_from_csv): func_2d_from_csv.set_interpolation("shepard") assert func_2d_from_csv.get_interpolation_method() == "shepard" func_2d_from_csv.set_extrapolation("zero") - assert func_2d_from_csv.get_extrapolation_method() == "zero" + # 2d functions do not support zero extrapolation, must change to natural + assert func_2d_from_csv.get_extrapolation_method() == "natural" @patch("matplotlib.pyplot.show") From c175e76d3dad30534ba18a1c7d0e07d635ac8c9d Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 28 Mar 2024 17:47:47 +0100 Subject: [PATCH 16/71] BUG: revert comit https://github.com/RocketPy-Team/RocketPy/pull/581/commits/9a637a834e443736e7fa988830d1586aaf1f9201 --- tests/test_flight.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_flight.py b/tests/test_flight.py index a4b03d101..db882e185 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -729,7 +729,7 @@ def test_velocities(flight_calisto_custom_wind, flight_time, expected_values): ("t_initial", (1.6542528, 0.65918, -0.067107)), ("out_of_rail_time", (5.05334, 2.01364, -1.7541)), ("apogee_time", (2.35291, -1.8275, -0.87851)), - ("t_final", (0, 0, 159.3292416824044)), + ("t_final", (0, 0, 159.2212)), ], ) def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values): From 81a79e3c34b98804645cc10ffd3c6b5895e9cb29 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 28 Mar 2024 18:51:31 +0100 Subject: [PATCH 17/71] ENH: change np.searchsorted to bisect_left --- rocketpy/mathutils/function.py | 43 ++++++++++++++++------------------ 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 86135e180..107c39c61 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -13,6 +13,7 @@ import matplotlib.pyplot as plt import numpy as np from scipy import integrate, linalg, optimize +from bisect import bisect_left try: from functools import cached_property @@ -308,7 +309,7 @@ def set_interpolation(self, method="spline"): if not callable(self.source): self.__interpolation__ = self._validate_interpolation(method) self._update_interpolation_coefficients(self.__interpolation__) - self._set_interpolation_function() + self._set_interpolation_func() return self def _update_interpolation_coefficients(self, method): @@ -346,16 +347,16 @@ def set_extrapolation(self, method="constant"): """ if not callable(self.source): self.__extrapolation__ = self._validate_extrapolation(method) - self._set_extrapolation_function() + self._set_extrapolation_func() return self - def _set_interpolation_function(self): + def _set_interpolation_func(self): """Return a function that interpolates the Function.""" interpolation = INTERPOLATION_TYPES[self.__interpolation__] if interpolation == 0: # linear def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs): - x_interval = np.searchsorted(x_data, x) + x_interval = bisect_left(x_data, x) x_left = x_data[x_interval - 1] y_left = y_data[x_interval - 1] dx = float(x_data[x_interval] - x_left) @@ -363,7 +364,7 @@ def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs): y = (x - x_left) * (dy / dx) + y_left return y - self._interpolation_function = linear_interpolation + self._interpolation_func = linear_interpolation elif interpolation == 1: # polynomial @@ -371,48 +372,48 @@ def polynomial_interpolation(x, x_min, x_max, x_data, y_data, coeffs): y = np.sum(coeffs * x ** np.arange(len(coeffs))) return y - self._interpolation_function = polynomial_interpolation + self._interpolation_func = polynomial_interpolation elif interpolation == 2: # akima def akima_interpolation(x, x_min, x_max, x_data, y_data, coeffs): - x_interval = np.searchsorted(x_data, x) + x_interval = bisect_left(x_data, x) x_interval = x_interval if x_interval != 0 else 1 a = coeffs[4 * x_interval - 4 : 4 * x_interval] y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] return y - self._interpolation_function = akima_interpolation + self._interpolation_func = akima_interpolation elif interpolation == 3: # spline def spline_interpolation(x, x_min, x_max, x_data, y_data, coeffs): - x_interval = np.searchsorted(x_data, x) + x_interval = bisect_left(x_data, x) x_interval = max(x_interval, 1) a = coeffs[:, x_interval - 1] x = x - x_data[x_interval - 1] y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] return y - self._interpolation_function = spline_interpolation + self._interpolation_func = spline_interpolation elif interpolation == 4: # shepard does not use interpolation function - self._interpolation_function = None + self._interpolation_func = None - def _set_extrapolation_function(self): + def _set_extrapolation_func(self): """Return a function that extrapolates the Function.""" interpolation = INTERPOLATION_TYPES[self.__interpolation__] extrapolation = EXTRAPOLATION_TYPES[self.__extrapolation__] if interpolation == 4: # shepard does not use extrapolation function - self._extrapolation_function = None + self._extrapolation_func = None elif extrapolation == 0: # zero def zero_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): return 0 - self._extrapolation_function = zero_extrapolation + self._extrapolation_func = zero_extrapolation elif extrapolation == 1: # natural if interpolation == 0: # linear @@ -450,13 +451,13 @@ def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] return y - self._extrapolation_function = natural_extrapolation + self._extrapolation_func = natural_extrapolation elif extrapolation == 2: # constant def constant_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): return y_data[0] if x < x_min else y_data[-1] - self._extrapolation_function = constant_extrapolation + self._extrapolation_func = constant_extrapolation def set_get_value_opt(self): """Crates a method that evaluates interpolations rather quickly @@ -476,20 +477,16 @@ def set_get_value_opt(self): self.get_value_opt = self._get_value_opt_nd return self - def _get_value_opt_1d(self, x, *args): + def _get_value_opt_1d(self, x): # Retrieve general info x_data = self.x_array y_data = self.y_array x_min, x_max = self.x_initial, self.x_final coeffs = self._coeffs if x_min <= x <= x_max: - y = self._interpolation_function( - x, x_min, x_max, x_data, y_data, coeffs, *args - ) + y = self._interpolation_func(x, x_min, x_max, x_data, y_data, coeffs) else: - y = self._extrapolation_function( - x, x_min, x_max, x_data, y_data, coeffs, *args - ) + y = self._extrapolation_func(x, x_min, x_max, x_data, y_data, coeffs) return y def _get_value_opt_nd(self, *args): From c1d13247133be3402c68216b3009156861f9df20 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 28 Mar 2024 19:08:12 +0100 Subject: [PATCH 18/71] DOC: adds docs to new methods --- rocketpy/mathutils/function.py | 30 ++++++++++++++++++++++++------ 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 107c39c61..ef9b13287 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -351,7 +351,11 @@ def set_extrapolation(self, method="constant"): return self def _set_interpolation_func(self): - """Return a function that interpolates the Function.""" + """Defines interpolation function used by the Function. Each + interpolation method has its own function with exception of shepard, + which has its interpolation/extrapolation function defined in + ``Function.__interpolate_shepard__``. The function is stored in + the attribute _interpolation_func.""" interpolation = INTERPOLATION_TYPES[self.__interpolation__] if interpolation == 0: # linear @@ -401,7 +405,9 @@ def spline_interpolation(x, x_min, x_max, x_data, y_data, coeffs): self._interpolation_func = None def _set_extrapolation_func(self): - """Return a function that extrapolates the Function.""" + """Defines extrapolation function used by the Function. Each + extrapolation method has its own function. The function is stored in + the attribute _extrapolation_func.""" interpolation = INTERPOLATION_TYPES[self.__interpolation__] extrapolation = EXTRAPOLATION_TYPES[self.__extrapolation__] @@ -460,10 +466,7 @@ def constant_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): self._extrapolation_func = constant_extrapolation def set_get_value_opt(self): - """Crates a method that evaluates interpolations rather quickly - when compared to other options available, such as just calling - the object instance or calling ``Function.get_value`` directly. See - ``Function.get_value_opt`` for documentation. + """Defines a method that evaluates interpolations. Returns ------- @@ -478,6 +481,19 @@ def set_get_value_opt(self): return self def _get_value_opt_1d(self, x): + """Evaluate the Function at a single point x. This method is used + when the Function is 1-D. + + Parameters + ---------- + x : scalar + Value where the Function is to be evaluated. + + Returns + ------- + y : scalar + Value of the Function at the specified point. + """ # Retrieve general info x_data = self.x_array y_data = self.y_array @@ -490,6 +506,8 @@ def _get_value_opt_1d(self, x): return y def _get_value_opt_nd(self, *args): + """Evaluate the Function at a single point (x, y, z). This method is + used when the Function is N-D.""" # always use shepard for N-D functions y = self.__interpolate_shepard__(list(args)) return y From d5d9c31e67eff9564f8382d424d5ff0e7ad206b8 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sat, 30 Mar 2024 22:11:47 +0100 Subject: [PATCH 19/71] MNT: run isort --- rocketpy/mathutils/function.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index ef9b13287..0cdb139cc 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -5,6 +5,7 @@ """ import warnings +from bisect import bisect_left from collections.abc import Iterable from copy import deepcopy from inspect import signature From 7cd270b956fa020305a6a6bdc45d52cab7283ce7 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sat, 30 Mar 2024 22:12:02 +0100 Subject: [PATCH 20/71] MNT: remove unecessary varibale intialization --- rocketpy/mathutils/function.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 0cdb139cc..20c9e41d3 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -14,7 +14,6 @@ import matplotlib.pyplot as plt import numpy as np from scipy import integrate, linalg, optimize -from bisect import bisect_left try: from functools import cached_property @@ -120,11 +119,6 @@ def __init__( (II) Fields in CSV files may be enclosed in double quotes. If fields are not quoted, double quotes should not appear inside them. """ - # initialize variables to avoid errors when being called by other methods - self.__polynomial_coefficients__ = None - self.__akima_coefficients__ = None - self.__spline_coefficients__ = None - # initialize parameters self.source = source self.__inputs__ = inputs From 0bd6f6396e8716fb0c76025dd2e29380528cfa0d Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sat, 30 Mar 2024 22:37:48 +0100 Subject: [PATCH 21/71] ENH: x, y and z array for all ndarray functions --- rocketpy/mathutils/function.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 20c9e41d3..731fc7810 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -240,7 +240,7 @@ def set_source(self, source): self.__dom_dim__ = new_total_dim - 1 self.__inputs__ = self.__dom_dim__ * self.__inputs__ - # if Function is 1D, sort source by x. If 2D, set z + # set x and y. If Function is 2D, also set z if self.__dom_dim__ == 1: source = source[source[:, 0].argsort()] self.x_array = source[:, 0] @@ -248,7 +248,7 @@ def set_source(self, source): self.y_array = source[:, 1] self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] self.get_value_opt = self._get_value_opt_1d - elif self.__dom_dim__ == 2: + elif self.__dom_dim__ > 1: self.x_array = source[:, 0] self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] self.y_array = source[:, 1] @@ -256,8 +256,6 @@ def set_source(self, source): self.z_array = source[:, 2] self.z_initial, self.z_final = self.z_array[0], self.z_array[-1] self.get_value_opt = self._get_value_opt_nd - else: - self.get_value_opt = self._get_value_opt_nd self.source = source self.set_interpolation(self.__interpolation__) From ca1ce44c912223c2906341e16077d947027cf7ed Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Sat, 30 Mar 2024 19:28:37 -0300 Subject: [PATCH 22/71] Update rocketpy/mathutils/function.py Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/mathutils/function.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 731fc7810..9dd9ebdf3 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -366,8 +366,7 @@ def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs): elif interpolation == 1: # polynomial def polynomial_interpolation(x, x_min, x_max, x_data, y_data, coeffs): - y = np.sum(coeffs * x ** np.arange(len(coeffs))) - return y + return np.sum(coeffs * x ** np.arange(len(coeffs))) self._interpolation_func = polynomial_interpolation From 8cfaedcae674f317e5b40d5829e5be6e66a1b46b Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Sun, 31 Mar 2024 07:45:26 -0300 Subject: [PATCH 23/71] Update rocketpy/mathutils/function.py Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/mathutils/function.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 9dd9ebdf3..94c9d70ff 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -501,8 +501,7 @@ def _get_value_opt_nd(self, *args): """Evaluate the Function at a single point (x, y, z). This method is used when the Function is N-D.""" # always use shepard for N-D functions - y = self.__interpolate_shepard__(list(args)) - return y + return self.__interpolate_shepard__(list(args)) def set_discrete( self, From 07f9ed50dccf78eda4095042d8868185ecd152dd Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sun, 31 Mar 2024 12:50:25 +0200 Subject: [PATCH 24/71] MNT: remove unecessary list casting --- 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 94c9d70ff..2e3ea0e63 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -501,7 +501,7 @@ def _get_value_opt_nd(self, *args): """Evaluate the Function at a single point (x, y, z). This method is used when the Function is N-D.""" # always use shepard for N-D functions - return self.__interpolate_shepard__(list(args)) + return self.__interpolate_shepard__(args) def set_discrete( self, From fa8a89208fcfb0234df25ed88d7c4d457c307643 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sun, 31 Mar 2024 12:50:44 +0200 Subject: [PATCH 25/71] MNT: return interp and extrap results directly --- rocketpy/mathutils/function.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 2e3ea0e63..02277ce1c 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -358,8 +358,7 @@ def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs): y_left = y_data[x_interval - 1] dx = float(x_data[x_interval] - x_left) dy = float(y_data[x_interval] - y_left) - y = (x - x_left) * (dy / dx) + y_left - return y + return (x - x_left) * (dy / dx) + y_left self._interpolation_func = linear_interpolation @@ -376,8 +375,7 @@ def akima_interpolation(x, x_min, x_max, x_data, y_data, coeffs): x_interval = bisect_left(x_data, x) x_interval = x_interval if x_interval != 0 else 1 a = coeffs[4 * x_interval - 4 : 4 * x_interval] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - return y + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] self._interpolation_func = akima_interpolation @@ -388,8 +386,7 @@ def spline_interpolation(x, x_min, x_max, x_data, y_data, coeffs): x_interval = max(x_interval, 1) a = coeffs[:, x_interval - 1] x = x - x_data[x_interval - 1] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - return y + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] self._interpolation_func = spline_interpolation @@ -421,21 +418,18 @@ def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): y_left = y_data[x_interval - 1] dx = float(x_data[x_interval] - x_left) dy = float(y_data[x_interval] - y_left) - y = (x - x_left) * (dy / dx) + y_left - return y + return (x - x_left) * (dy / dx) + y_left elif interpolation == 1: # polynomial def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): - y = np.sum(coeffs * x ** np.arange(len(coeffs))) - return y + return np.sum(coeffs * x ** np.arange(len(coeffs))) elif interpolation == 2: # akima def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): a = coeffs[:4] if x < x_min else coeffs[-4:] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - return y + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] elif interpolation == 3: # spline @@ -446,8 +440,7 @@ def natural_extrapolation(x, x_min, x_max, x_data, y_data, coeffs): else: a = coeffs[:, -1] x = x - x_data[-2] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - return y + return a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] self._extrapolation_func = natural_extrapolation elif extrapolation == 2: # constant From 29bb5fae988226ad2e119171a78e124e95275f87 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Sun, 31 Mar 2024 13:04:23 +0200 Subject: [PATCH 26/71] ENH: completely private methods --- rocketpy/mathutils/function.py | 44 +++++++++++++++++----------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 02277ce1c..4a7c81154 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -147,7 +147,7 @@ def set_inputs(self, inputs): ------- self : Function """ - self.__inputs__ = self._validate_inputs(inputs) + self.__inputs__ = self.__validate_inputs(inputs) return self def set_outputs(self, outputs): @@ -162,7 +162,7 @@ def set_outputs(self, outputs): ------- self : Function """ - self.__outputs__ = self._validate_outputs(outputs) + self.__outputs__ = self.__validate_outputs(outputs) return self def set_source(self, source): @@ -211,7 +211,7 @@ def set_source(self, source): self : Function Returns the Function instance. """ - source = self._validate_source(source) + source = self.__validate_source(source) # Handle callable source or number source if callable(source): @@ -247,7 +247,7 @@ def set_source(self, source): self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] self.y_array = source[:, 1] self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] - self.get_value_opt = self._get_value_opt_1d + self.get_value_opt = self.__get_value_opt_1d elif self.__dom_dim__ > 1: self.x_array = source[:, 0] self.x_initial, self.x_final = self.x_array[0], self.x_array[-1] @@ -255,7 +255,7 @@ def set_source(self, source): self.y_initial, self.y_final = self.y_array[0], self.y_array[-1] self.z_array = source[:, 2] self.z_initial, self.z_final = self.z_array[0], self.z_array[-1] - self.get_value_opt = self._get_value_opt_nd + self.get_value_opt = self.__get_value_opt_nd self.source = source self.set_interpolation(self.__interpolation__) @@ -300,12 +300,12 @@ def set_interpolation(self, method="spline"): self : Function """ if not callable(self.source): - self.__interpolation__ = self._validate_interpolation(method) - self._update_interpolation_coefficients(self.__interpolation__) - self._set_interpolation_func() + self.__interpolation__ = self.__validate_interpolation(method) + self.__update_interpolation_coefficients(self.__interpolation__) + self.__set_interpolation_func() return self - def _update_interpolation_coefficients(self, method): + def __update_interpolation_coefficients(self, method): """Update interpolation coefficients for the given method.""" # Spline, akima and polynomial need data processing # Shepard, and linear do not @@ -339,11 +339,11 @@ def set_extrapolation(self, method="constant"): The Function object. """ if not callable(self.source): - self.__extrapolation__ = self._validate_extrapolation(method) - self._set_extrapolation_func() + self.__extrapolation__ = self.__validate_extrapolation(method) + self.__set_extrapolation_func() return self - def _set_interpolation_func(self): + def __set_interpolation_func(self): """Defines interpolation function used by the Function. Each interpolation method has its own function with exception of shepard, which has its interpolation/extrapolation function defined in @@ -393,7 +393,7 @@ def spline_interpolation(x, x_min, x_max, x_data, y_data, coeffs): elif interpolation == 4: # shepard does not use interpolation function self._interpolation_func = None - def _set_extrapolation_func(self): + def __set_extrapolation_func(self): """Defines extrapolation function used by the Function. Each extrapolation method has its own function. The function is stored in the attribute _extrapolation_func.""" @@ -460,12 +460,12 @@ def set_get_value_opt(self): if callable(self.source): self.get_value_opt = self.source elif self.__dom_dim__ == 1: - self.get_value_opt = self._get_value_opt_1d + self.get_value_opt = self.__get_value_opt_1d elif self.__dom_dim__ > 1: - self.get_value_opt = self._get_value_opt_nd + self.get_value_opt = self.__get_value_opt_nd return self - def _get_value_opt_1d(self, x): + def __get_value_opt_1d(self, x): """Evaluate the Function at a single point x. This method is used when the Function is 1-D. @@ -490,7 +490,7 @@ def _get_value_opt_1d(self, x): y = self._extrapolation_func(x, x_min, x_max, x_data, y_data, coeffs) return y - def _get_value_opt_nd(self, *args): + def __get_value_opt_nd(self, *args): """Evaluate the Function at a single point (x, y, z). This method is used when the Function is N-D.""" # always use shepard for N-D functions @@ -2883,7 +2883,7 @@ def savetxt( np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline) # Input validators - def _validate_source(self, source): + def __validate_source(self, source): """Used to validate the source parameter for creating a Function object. Parameters @@ -2951,7 +2951,7 @@ def source_function(_): # If source is a callable function return source - def _validate_inputs(self, inputs): + def __validate_inputs(self, inputs): """Used to validate the inputs parameter for creating a Function object. It sets a default value if it is not provided. @@ -2990,7 +2990,7 @@ def _validate_inputs(self, inputs): "the length of the domain dimension." ) - def _validate_outputs(self, outputs): + def __validate_outputs(self, outputs): """Used to validate the outputs parameter for creating a Function object. It sets a default value if it is not provided. @@ -3016,7 +3016,7 @@ def _validate_outputs(self, outputs): ) return outputs - def _validate_interpolation(self, interpolation): + def __validate_interpolation(self, interpolation): if self.__dom_dim__ == 1: # possible interpolation values: linear, polynomial, akima and spline if interpolation is None: @@ -3044,7 +3044,7 @@ def _validate_interpolation(self, interpolation): interpolation = "shepard" return interpolation - def _validate_extrapolation(self, extrapolation): + def __validate_extrapolation(self, extrapolation): if self.__dom_dim__ == 1: if extrapolation is None: extrapolation = "constant" From 3bee49b3c723dcdc33380254513e985ffec0cf97 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 4 Apr 2024 13:17:37 +0200 Subject: [PATCH 27/71] TST: improve tests coverage --- rocketpy/mathutils/function.py | 18 ++++-------------- tests/test_function.py | 27 ++++++++++++++++++++++++++- tests/unit/test_function.py | 28 ++++++++++++++++++++++++++++ 3 files changed, 58 insertions(+), 15 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 4a7c81154..ea98bff78 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -230,16 +230,6 @@ def set_source(self, source): # Evaluate dimension self.__dom_dim__ = source.shape[1] - 1 - # Check to see if dimensions match incoming data set - new_total_dim = source.shape[1] - old_total_dim = self.__dom_dim__ + self.__img_dim__ - - # If they don't, update default values or throw error - if new_total_dim != old_total_dim: - # Update dimensions and inputs - self.__dom_dim__ = new_total_dim - 1 - self.__inputs__ = self.__dom_dim__ * self.__inputs__ - # set x and y. If Function is 2D, also set z if self.__dom_dim__ == 1: source = source[source[:, 0].argsort()] @@ -2970,8 +2960,8 @@ def __validate_inputs(self, inputs): return ["Scalar"] if isinstance(inputs, str): return [inputs] - if isinstance(inputs, list): - if len(inputs) == 1 and isinstance(inputs[0], str): + if isinstance(inputs, (list, tuple)): + if len(inputs) == 1: return inputs raise ValueError( "Inputs must be a string or a list of strings with " @@ -3008,8 +2998,8 @@ def __validate_outputs(self, outputs): return ["Scalar"] if isinstance(outputs, str): return [outputs] - if isinstance(outputs, list): - if len(outputs) > 1 or not isinstance(outputs[0], str): + if isinstance(outputs, (list, tuple)): + if len(outputs) > 1: raise ValueError( "Output must either be a string or a list of strings with " + f"one item. It currently has dimension ({len(outputs)})." diff --git a/tests/test_function.py b/tests/test_function.py index 70613671e..6f4122e47 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -182,7 +182,32 @@ def test_extrapolation_methods(linear_func): assert linear_func.get_extrapolation_method() == "constant" assert np.isclose(linear_func.get_value(-1), 0, atol=1e-6) - # Test natural + # Test natural for linear interpolation + linear_func.set_interpolation("linear") + assert isinstance(linear_func.set_extrapolation("natural"), Function) + linear_func.set_extrapolation("natural") + assert isinstance(linear_func.get_extrapolation_method(), str) + assert linear_func.get_extrapolation_method() == "natural" + assert np.isclose(linear_func.get_value(-1), -1, atol=1e-6) + + # Test natural for spline interpolation + linear_func.set_interpolation("spline") + assert isinstance(linear_func.set_extrapolation("natural"), Function) + linear_func.set_extrapolation("natural") + assert isinstance(linear_func.get_extrapolation_method(), str) + assert linear_func.get_extrapolation_method() == "natural" + assert np.isclose(linear_func.get_value(-1), -1, atol=1e-6) + + # Test natural for akima interpolation + linear_func.set_interpolation("akima") + assert isinstance(linear_func.set_extrapolation("natural"), Function) + linear_func.set_extrapolation("natural") + assert isinstance(linear_func.get_extrapolation_method(), str) + assert linear_func.get_extrapolation_method() == "natural" + assert np.isclose(linear_func.get_value(-1), -1, atol=1e-6) + + # Test natural for polynomial interpolation + linear_func.set_interpolation("polynomial") assert isinstance(linear_func.set_extrapolation("natural"), Function) linear_func.set_extrapolation("natural") assert isinstance(linear_func.get_extrapolation_method(), str) diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index 8bcefb818..17da59498 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -306,3 +306,31 @@ def test_remove_outliers_iqr(x, y, expected_x, expected_y): assert filtered_func.__interpolation__ == func.__interpolation__ assert filtered_func.__extrapolation__ == func.__extrapolation__ assert filtered_func.title == func.title + + +def test_set_get_value_opt(): + """Test the set_value_opt and get_value_opt methods of the Function class.""" + func = Function(lambda x: x**2) + func.source = np.array([[1, 1], [2, 4], [3, 9], [4, 16], [5, 25]]) + func.x_array = np.array([1, 2, 3, 4, 5]) + func.y_array = np.array([1, 4, 9, 16, 25]) + func.x_initial = 1 + func.x_final = 5 + func.set_interpolation("linear") + func.set_get_value_opt() + assert func.get_value_opt(2.5) == 6.5 + + +def test_get_image_dim(linear_func): + """Test the get_img_dim method of the Function class.""" + assert linear_func.get_image_dim() == 1 + + +def test_get_domain_dim(linear_func): + """Test the get_domain_dim method of the Function class.""" + assert linear_func.get_domain_dim() == 1 + + +def test_bool(linear_func): + """Test the __bool__ method of the Function class.""" + assert bool(linear_func) == True From 00a3c1f2971d23933cb86e0ed852dddfe72461c2 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 4 Apr 2024 13:19:14 +0200 Subject: [PATCH 28/71] DEV: changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 90aa086ca..b96efbcc1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ### Added +- ENH: Function Validation Rework & Swap `np.searchsorted` to `bisect_left` [#582](https://github.com/RocketPy-Team/RocketPy/pull/582) - ENH: Add new stability margin properties to Flight class [#572](https://github.com/RocketPy-Team/RocketPy/pull/572) - ENH: adds `Function.remove_outliers` method [#554](https://github.com/RocketPy-Team/RocketPy/pull/554) From 849851255eee91765ea856c53d71c629e63b9289 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 15 Apr 2024 23:45:56 -0400 Subject: [PATCH 29/71] TST: Add unit tests for TimeNodes class and its methods --- tests/unit/test_flight_time_nodes.py | 103 +++++++++++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 tests/unit/test_flight_time_nodes.py diff --git a/tests/unit/test_flight_time_nodes.py b/tests/unit/test_flight_time_nodes.py new file mode 100644 index 000000000..ac22e4490 --- /dev/null +++ b/tests/unit/test_flight_time_nodes.py @@ -0,0 +1,103 @@ +"""Module to test everything related to the TimeNodes class and it's subclass +TimeNode. +""" + +import pytest + +from rocketpy.rocket import Parachute, _Controller + + +def test_time_nodes_init(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + assert len(time_nodes) == 0 + + +def test_time_nodes_getitem(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(1.0, [], []) + assert isinstance(time_nodes[0], flight_calisto.TimeNodes.TimeNode) + assert time_nodes[0].t == 1.0 + + +def test_time_nodes_len(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + assert len(time_nodes) == 0 + + +def test_time_nodes_add(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + example_node = flight_calisto.TimeNodes.TimeNode(1.0, [], []) + time_nodes.add(example_node) + assert len(time_nodes) == 1 + assert isinstance(time_nodes[0], flight_calisto.TimeNodes.TimeNode) + assert time_nodes[0].t == 1.0 + + +def test_time_nodes_add_node(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(2.0, [], []) + assert len(time_nodes) == 1 + assert time_nodes[0].t == 2.0 + assert len(time_nodes[0].parachutes) == 0 + assert len(time_nodes[0].callbacks) == 0 + + +# def test_time_nodes_add_parachutes( +# flight_calisto, calisto_drogue_chute, calisto_main_chute +# ): # TODO: implement this test + + +# def test_time_nodes_add_controllers(flight_calisto): +# TODO: implement this test + + +def test_time_nodes_sort(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add(3.0) + time_nodes.add(1.0) + time_nodes.add(2.0) + time_nodes.sort() + assert len(time_nodes) == 3 + assert time_nodes[0] == 1.0 + assert time_nodes[1] == 2.0 + assert time_nodes[2] == 3.0 + + +def test_time_nodes_merge(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(2.0, [], []) + time_nodes.merge() + assert len(time_nodes) == 2 + assert time_nodes[0].t == 1.0 + assert len(time_nodes[0].parachutes) == 0 + assert len(time_nodes[0].callbacks) == 0 + assert time_nodes[1].t == 2.0 + assert len(time_nodes[1].parachutes) == 0 + assert len(time_nodes[1].callbacks) == 0 + + +def test_time_nodes_flush_after(flight_calisto): + time_nodes = flight_calisto.TimeNodes() + time_nodes.add(1.0) + time_nodes.add(2.0) + time_nodes.add(3.0) + time_nodes.flush_after(1) + assert len(time_nodes) == 2 + assert time_nodes[0] == 1.0 + assert time_nodes[1] == 2.0 + + +def test_time_node_init(flight_calisto): + node = flight_calisto.TimeNodes.TimeNode(1.0, [], []) + assert node.t == 1.0 + assert len(node.parachutes) == 0 + assert len(node.callbacks) == 0 + + +def test_time_node_lt(flight_calisto): + node1 = flight_calisto.TimeNodes.TimeNode(1.0, [], []) + node2 = flight_calisto.TimeNodes.TimeNode(2.0, [], []) + assert node1 < node2 + assert not node2 < node1 From 88eefa28175bb43916b352b6b45b8d6b6b618d40 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 15 Apr 2024 23:47:30 -0400 Subject: [PATCH 30/71] ENH: add Function.calculate_cubic_hermite_coefficients method --- rocketpy/mathutils/function.py | 36 ++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index ea98bff78..652dfd3fd 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2872,6 +2872,42 @@ def savetxt( file.write(header_line + newline) np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline) + # Define auxiliary method + + @staticmethod + def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): + """Calculate the coefficients of a cubic Hermite interpolation function. + The function is defined as ax**3 + bx**2 + cx + d. + + Parameters + ---------- + x0 : float + Position of the first point. + x1 : float + Position of the second point. + y0 : float + Value of the function evaluated at the first point. + yp0 : float + Value of the derivative of the function evaluated at the first + point. + y1 : float + Value of the function evaluated at the second point. + yp1 : float + Value of the derivative of the function evaluated at the second + point. + + Returns + ------- + tuple[float, float, float, float] + The coefficients of the cubic Hermite interpolation function. + """ + dx = x1 - x0 + d = y0 + c = yp0 + b = (3 * y1 - yp1 * dx - 2 * c * dx - 3 * d) / (dx**2) + a = -(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3) + return a, b, c, d + # Input validators def __validate_source(self, source): """Used to validate the source parameter for creating a Function object. From cacd95ffcba88588f92738bfd26cf83e158c0a3c Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 15 Apr 2024 23:48:13 -0400 Subject: [PATCH 31/71] ENH: Adds Function.cardanos_root_finding --- rocketpy/mathutils/function.py | 60 ++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 652dfd3fd..ba90b37f1 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2908,6 +2908,66 @@ def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): a = -(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3) return a, b, c, d + @staticmethod + def cardanos_root_finding(a, b, c, d): + """Calculate the roots of a cubic function using Cardano's method. + + This method applies Cardano's method to find the roots of a cubic + function of the form ax^3 + bx^2 + cx + d. The roots may be complex + numbers. + + Parameters + ---------- + a : float + Coefficient of the cubic term (x^3). + b : float + Coefficient of the quadratic term (x^2). + c : float + Coefficient of the linear term (x). + d : float + Constant term. + + Returns + ------- + tuple[complex, complex, complex] + A tuple containing the real and complex roots of the cubic function. + Note that the roots may be complex numbers. The roots are ordered + in the tuple as x1, x2, x3. + + References + ---------- + - Cardano's method: https://en.wikipedia.org/wiki/Cubic_function#Cardano's_method + + Examples + -------- + >>> from rocketpy import Function + + First we define the coefficients of the function ax**3 + bx**2 + cx + d + >>> a, b, c, d = 1, -3, -1, 3 + >>> x1, x2, x3 = Function.cardanos_root_finding(a, b, c, d) + >>> x1, x2, x3 + ((-1+0j), (3+7.401486830834377e-17j), (1-1.4802973661668753e-16j)) + + To get the real part of the roots, use the real attribute of the complex + number. + >>> x1.real, x2.real, x3.real + (-1.0, 3.0, 1.0) + """ + delta_0 = b**2 - 3 * a * c + delta_1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 + c1 = ((delta_1 + (delta_1**2 - 4 * delta_0**3) ** (0.5)) / 2) ** (1 / 3) + + c2_0 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 0 + x1 = -(1 / (3 * a)) * (b + c2_0 + delta_0 / c2_0) + + c2_1 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 1 + x2 = -(1 / (3 * a)) * (b + c2_1 + delta_0 / c2_1) + + c2_2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 2 + x3 = -(1 / (3 * a)) * (b + c2_2 + delta_0 / c2_2) + + return x1, x2, x3 + # Input validators def __validate_source(self, source): """Used to validate the source parameter for creating a Function object. From c7bbd4537643635877f5c818472be5a7a35dd545 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 15 Apr 2024 23:49:35 -0400 Subject: [PATCH 32/71] TST: Refactor acceptance message to include error threshold --- tests/acceptance/test_bella_lui_rocket.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/tests/acceptance/test_bella_lui_rocket.py b/tests/acceptance/test_bella_lui_rocket.py index a1297e1fb..0041074a8 100644 --- a/tests/acceptance/test_bella_lui_rocket.py +++ b/tests/acceptance/test_bella_lui_rocket.py @@ -228,10 +228,14 @@ def drogue_trigger(p, h, y): apogee_time_measured = time_kalt[np.argmax(altitude_kalt)] apogee_time_simulated = test_flight.apogee_time - assert ( - abs(max(altitude_kalt) - test_flight.apogee + test_flight.env.elevation) - / max(altitude_kalt) - < 0.015 + apogee_error_threshold = 0.015 + apogee_error = abs( + max(altitude_kalt) - test_flight.apogee + test_flight.env.elevation + ) / max(altitude_kalt) + assert apogee_error < apogee_error_threshold, ( + f"Apogee altitude error exceeded the threshold. " + f"Expected the error to be less than {apogee_error_threshold * 100}%, " + f"but got an error of {apogee_error * 100:.1f}%." ) assert abs(max(velocity_rcp) - max(vert_vel_kalt)) / max(vert_vel_kalt) < 0.06 assert ( From 80d271eeb80c87e51c1c67f76bba3f535a4dd424 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 15 Apr 2024 23:51:37 -0400 Subject: [PATCH 33/71] ENH: Allows the comparison of two TimeNode objects based on their initial time. --- rocketpy/simulation/flight.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index b13f923eb..221d9fefb 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3685,9 +3685,22 @@ def __init__(self, t, parachutes, controllers): def __repr__(self): return ( - "{Initial Time: " - + str(self.t) - + " | Parachutes: " - + str(len(self.parachutes)) - + "}" - ) + + def __lt__(self, other): + """Allows the comparison of two TimeNode objects based on their + initial time. This is particularly useful for sorting a list of + TimeNode objects. + + Parameters + ---------- + other : TimeNode + Another TimeNode object to compare with. + + Returns + ------- + bool + True if the initial time of the current TimeNode is less + than the initial time of the other TimeNode, False + otherwise. + """ + return self.t < other.t From ac294c6562d91bc91789465879e003e6f9a305a9 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 15 Apr 2024 23:58:21 -0400 Subject: [PATCH 34/71] ENH: calculates atmospheric attributes from environment --- rocketpy/simulation/flight.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 221d9fefb..5ccc13b74 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -2177,36 +2177,33 @@ def M3(self): @funcify_method("Time (s)", "Pressure (Pa)", "spline", "constant") def pressure(self): """Air pressure felt by the rocket as a Function of time.""" - return self.retrieve_temporary_values_arrays[6] + return [(t, self.env.pressure.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Density (kg/m³)", "spline", "constant") def density(self): """Air density felt by the rocket as a Function of time.""" - return self.retrieve_temporary_values_arrays[7] + return [(t, self.env.density.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Dynamic Viscosity (Pa s)", "spline", "constant") def dynamic_viscosity(self): """Air dynamic viscosity felt by the rocket as a Function of time.""" - return self.retrieve_temporary_values_arrays[8] + return [(t, self.env.dynamic_viscosity.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Speed of Sound (m/s)", "spline", "constant") def speed_of_sound(self): - """Speed of sound in the air felt by the rocket as a Function - of time.""" - return self.retrieve_temporary_values_arrays[9] + """Speed of sound in the air felt by the rocket as a Function of time.""" + return [(t, self.env.speed_of_sound.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Wind Velocity X (East) (m/s)", "spline", "constant") def wind_velocity_x(self): - """Wind velocity in the X direction (east) as a Function of - time.""" - return self.retrieve_temporary_values_arrays[10] + """Wind velocity in the X direction (east) as a Function of time.""" + return [(t, self.env.wind_velocity_x.get_value_opt(z)) for t, z in self.z] @funcify_method("Time (s)", "Wind Velocity Y (North) (m/s)", "spline", "constant") def wind_velocity_y(self): - """Wind velocity in the y direction (north) as a Function of - time.""" - return self.retrieve_temporary_values_arrays[11] + """Wind velocity in the y direction (north) as a Function of time.""" + return [(t, self.env.wind_velocity_y.get_value_opt(z)) for t, z in self.z] # Process fourth type of output - values calculated from previous outputs From cadcf10000325019731419c01379f18abd8d7d21 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 16 Apr 2024 00:45:18 -0400 Subject: [PATCH 35/71] ENH: Add property to calculate differences in function evaluations per time step --- rocketpy/simulation/flight.py | 45 ++++++++++++++++++++++------------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 5ccc13b74..271b7ff2b 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -759,29 +759,22 @@ def __init__( # Step through simulation while phase.solver.status == "running": - # Step + # Execute solver step, log solution and function evaluations phase.solver.step() - # Save step result self.solution += [[phase.solver.t, *phase.solver.y]] - # Step step metrics - self.function_evaluations_per_time_step.append( - phase.solver.nfev - self.function_evaluations[-1] - ) self.function_evaluations.append(phase.solver.nfev) - self.time_steps.append(phase.solver.step_size) + # Update time and state self.t = phase.solver.t self.y_sol = phase.solver.y - if verbose: - print( - "Current Simulation Time: {:3.4f} s".format(self.t), - end="\r", - ) - # print('\n\t\t\tCurrent Step Details') - # print('\t\t\tIState: ', phase.solver._lsoda_solver._integrator.istate) - # print('\t\t\tTime: ', phase.solver.t) - # print('\t\t\tAltitude: ', phase.solver.y[2]) - # print('\t\t\tEvals: ', self.function_evaluations_per_time_step[-1]) + if verbose: # TODO: change verbose to a logging. + print(f"Current Simulation Time: {self.t:3.4f} s", end="\r") + # print("\n\t\t\tCurrent Step Details:") + # print( + # "\t\t\tIState: ", phase.solver._lsoda_solver._integrator.istate + # ) + # print("\t\t\tTime: ", phase.solver.t) + # print("\t\t\tAltitude (ASL): ", phase.solver.y[2]) # Check for first out of rail event if len(self.out_of_rail_state) == 1 and ( @@ -1995,11 +1988,29 @@ def solution_array(self): """Returns solution array of the rocket flight.""" return np.array(self.solution) + @property + def function_evaluations_per_time_step(self): + """Get the number of function evaluations per time step. This method + calculates the difference between consecutive function evaluations + during numerical integration and returns it as a list. + + Returns + ------- + list + The list of differences in function evaluations per time step. + """ + return np.diff(self.function_evaluations).tolist() + @cached_property def time(self): """Returns time array from solution.""" return self.solution_array[:, 0] + @cached_property + def time_steps(self): + """Returns time step array.""" + return np.diff(self.time) + def get_solution_at_time(self, t, atol=1e-3): """Returns the solution state vector at a given time. If the time is not found in the solution, the closest time is used and a warning is From 77be755532ae2dda1ecb14f909e3af837b055178 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 16 Apr 2024 00:55:58 -0400 Subject: [PATCH 36/71] ENH: small fixes to the TimeNode() class --- rocketpy/simulation/flight.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 271b7ff2b..a94b8744e 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3685,7 +3685,26 @@ def flush_after(self, index): del self.list[index + 1 :] class TimeNode: + """TimeNode is a class that represents a time node in the time + nodes list. It stores the time, the parachutes and the controllers + that are active at that time. This class is supposed to work + exclusively within the TimeNodes class. + """ + def __init__(self, t, parachutes, controllers): + """Create a TimeNode object. + + Parameters + ---------- + t : float + Initial time of the time node. + parachutes : list[Parachute] + List containing all the parachutes that should be evaluated + at this time node. + controllers : list[_Controller] + List containing all the controllers that should be evaluated + at this time node. + """ self.t = t self.parachutes = parachutes self.callbacks = [] @@ -3693,6 +3712,11 @@ def __init__(self, t, parachutes, controllers): def __repr__(self): return ( + f"" + ) def __lt__(self, other): """Allows the comparison of two TimeNode objects based on their From dc6ff4aa939b4730b17a8a77cad1a64a2a585f67 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 16 Apr 2024 01:09:53 -0400 Subject: [PATCH 37/71] ENH: Improve Flight.__init_xxxx methods --- rocketpy/simulation/flight.py | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a94b8744e..5e127aea4 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -608,9 +608,10 @@ def __init__( self.equations_of_motion = equations_of_motion # Flight initialization - self.__init_post_process_variables() self.__init_solution_monitors() self.__init_equations_of_motion() + self.__init_solver_monitors() + # Initialize prints and plots objects self.prints = _FlightPrints(self) @@ -1094,12 +1095,8 @@ def __init_solution_monitors(self): self.out_of_rail_time = 0 self.out_of_rail_time_index = 0 self.out_of_rail_state = np.array([0]) - self.out_of_rail_velocity = 0 self.apogee_state = np.array([0]) - self.apogee_time = 0 - self.apogee_x = 0 - self.apogee_y = 0 - self.apogee = 0 + # self.apogee_time = 0 self.x_impact = 0 self.y_impact = 0 self.impact_velocity = 0 @@ -1107,8 +1104,6 @@ def __init_solution_monitors(self): self.parachute_events = [] self.post_processed = False - return None - def __init_flight_state(self): """Initialize flight state variables.""" if self.initial_solution is None: @@ -1164,8 +1159,6 @@ def __init_flight_state(self): def __init_solver_monitors(self): # Initialize solver monitors self.function_evaluations = [] - self.function_evaluations_per_time_step = [] - self.time_steps = [] # Initialize solution state self.solution = [] self.__init_flight_state() @@ -1180,11 +1173,6 @@ 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 - @cached_property def effective_1rl(self): """Original rail length minus the distance measured from nozzle exit From 826331f694f63bf5c2ebb3f2f57f8b452a25f79a Mon Sep 17 00:00:00 2001 From: Lint Action Date: Tue, 16 Apr 2024 05:10:56 +0000 Subject: [PATCH 38/71] Fix code style issues with Black --- rocketpy/simulation/flight.py | 1 - 1 file changed, 1 deletion(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 5e127aea4..92b0a5432 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -612,7 +612,6 @@ def __init__( self.__init_equations_of_motion() self.__init_solver_monitors() - # Initialize prints and plots objects self.prints = _FlightPrints(self) self.plots = _FlightPlots(self) From 25093f6744221d3752f2726a1365241f602280d4 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 16:42:04 -0400 Subject: [PATCH 39/71] BUG: restore tests, roll back errors --- rocketpy/simulation/flight.py | 42 +++++++++++++++------------- tests/unit/test_flight_time_nodes.py | 22 +++++++-------- 2 files changed, 33 insertions(+), 31 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 51c5b82e5..6e25b479c 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -616,9 +616,6 @@ def __init__( self.prints = _FlightPrints(self) self.plots = _FlightPlots(self) - # Initialize solver monitors - self.__init_solver_monitors() - # Create known flight phases self.FlightPhases = FlightPhases() self.FlightPhases.add_phase( @@ -767,7 +764,7 @@ def __init__( # Update time and state self.t = phase.solver.t self.y_sol = phase.solver.y - if verbose: # TODO: change verbose to a logging. + if verbose: print(f"Current Simulation Time: {self.t:3.4f} s", end="\r") # print("\n\t\t\tCurrent Step Details:") # print( @@ -852,9 +849,6 @@ def __init__( self.out_of_rail_time = self.t self.out_of_rail_time_index = len(self.solution) - 1 self.out_of_rail_state = self.y_sol - self.out_of_rail_velocity = ( - self.y_sol[3] ** 2 + self.y_sol[4] ** 2 + self.y_sol[5] ** 2 - ) ** (0.5) # Create new flight phase self.FlightPhases.add_phase( self.t, @@ -882,9 +876,6 @@ def __init__( self.apogee_state = interpolator(t_root) # Store apogee data self.apogee_time = t_root - self.apogee_x = self.apogee_state[0] - self.apogee_y = self.apogee_state[1] - self.apogee = self.apogee_state[2] if self.terminate_on_apogee: # print('Terminate on Apogee Activated!') self.t = t_root @@ -1080,22 +1071,13 @@ def __init__( if verbose: print("Simulation Completed at Time: {:3.4f} s".format(self.t)) - def __init_post_process_variables(self): - """Initialize post-process variables.""" - # Initialize all variables calculated after initialization. - # Important to do so that MATLAB® can access them - self._drift = Function(0) - self._bearing = Function(0) - self._latitude = Function(0) - self._longitude = Function(0) - def __init_solution_monitors(self): # Initialize solution monitors self.out_of_rail_time = 0 self.out_of_rail_time_index = 0 self.out_of_rail_state = np.array([0]) self.apogee_state = np.array([0]) - # self.apogee_time = 0 + self.apogee_time = 0 self.x_impact = 0 self.y_impact = 0 self.impact_velocity = 0 @@ -2046,6 +2028,21 @@ def z(self): """Rocket z position as a Function of time.""" return self.solution_array[:, [0, 3]] + @property + def apogee_x(self): + """Rocket x position at apogee.""" + return self.x.get_value_opt(self.apogee_time) + + @property + def apogee_y(self): + """Rocket y position at apogee.""" + return self.y.get_value_opt(self.apogee_time) + + @property + def apogee(self): + """Rocket z position at apogee.""" + return self.z.get_value_opt(self.apogee_time) + @funcify_method("Time (s)", "Altitude AGL (m)", "spline", "constant") def altitude(self): """Rocket altitude above ground level as a Function of time.""" @@ -2212,6 +2209,11 @@ def speed(self): """Rocket speed, or velocity magnitude, as a Function of time.""" return (self.vx**2 + self.vy**2 + self.vz**2) ** 0.5 + @property + def out_of_rail_velocity(self): + """Velocity at which the rocket leaves the launch rail.""" + return self.speed.get_value_opt(self.out_of_rail_time) + @cached_property def max_speed_time(self): """Time at which the rocket reaches its maximum speed.""" diff --git a/tests/unit/test_flight_time_nodes.py b/tests/unit/test_flight_time_nodes.py index ac22e4490..10f6b6c30 100644 --- a/tests/unit/test_flight_time_nodes.py +++ b/tests/unit/test_flight_time_nodes.py @@ -53,14 +53,14 @@ def test_time_nodes_add_node(flight_calisto): def test_time_nodes_sort(flight_calisto): time_nodes = flight_calisto.TimeNodes() - time_nodes.add(3.0) - time_nodes.add(1.0) - time_nodes.add(2.0) + time_nodes.add_node(3.0, [], []) + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(2.0, [], []) time_nodes.sort() assert len(time_nodes) == 3 - assert time_nodes[0] == 1.0 - assert time_nodes[1] == 2.0 - assert time_nodes[2] == 3.0 + assert time_nodes[0].t == 1.0 + assert time_nodes[1].t == 2.0 + assert time_nodes[2].t == 3.0 def test_time_nodes_merge(flight_calisto): @@ -80,13 +80,13 @@ def test_time_nodes_merge(flight_calisto): def test_time_nodes_flush_after(flight_calisto): time_nodes = flight_calisto.TimeNodes() - time_nodes.add(1.0) - time_nodes.add(2.0) - time_nodes.add(3.0) + time_nodes.add_node(1.0, [], []) + time_nodes.add_node(2.0, [], []) + time_nodes.add_node(3.0, [], []) time_nodes.flush_after(1) assert len(time_nodes) == 2 - assert time_nodes[0] == 1.0 - assert time_nodes[1] == 2.0 + assert time_nodes[0].t == 1.0 + assert time_nodes[1].t == 2.0 def test_time_node_init(flight_calisto): From 13b3d571eeedc326e736677c30ee5593f5843ef7 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 16:48:39 -0400 Subject: [PATCH 40/71] MNT: renames FlightPhases to flight_phases --- rocketpy/simulation/flight.py | 41 +++++++++++++++++------------------ 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 6e25b479c..db09831e1 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -188,7 +188,7 @@ class Flight: Flight.time_steps : array List of time steps taking during numerical integration in seconds. - Flight.FlightPhases : Flight.FlightPhases + Flight.flight_phases : Flight.FlightPhases Stores and manages flight phases. Flight.wind_velocity_x : Function Wind velocity X (East) experienced by the rocket as a @@ -581,7 +581,6 @@ def __init__( None """ # Fetch helper classes and functions - FlightPhases = self.FlightPhases TimeNodes = self.TimeNodes time_iterator = self.time_iterator @@ -617,20 +616,20 @@ def __init__( self.plots = _FlightPlots(self) # Create known flight phases - self.FlightPhases = FlightPhases() - self.FlightPhases.add_phase( + self.flight_phases = self.FlightPhases() + self.flight_phases.add_phase( self.t_initial, self.initial_derivative, clear=False ) - self.FlightPhases.add_phase(self.max_time) + self.flight_phases.add_phase(self.max_time) # Simulate flight - for phase_index, phase in time_iterator(self.FlightPhases): + for phase_index, phase in time_iterator(self.flight_phases): # print('\nCurrent Flight Phase List') - # print(self.FlightPhases) + # print(self.flight_phases) # print('\n\tCurrent Flight Phase') # print('\tIndex: ', phase_index, ' | Phase: ', phase) # Determine maximum time for this flight phase - phase.time_bound = self.FlightPhases[phase_index + 1].t + phase.time_bound = self.flight_phases[phase_index + 1].t # Evaluate callbacks for callback in phase.callbacks: @@ -727,7 +726,7 @@ def __init__( # Must only be created if parachute has any lag i = 1 if parachute.lag != 0: - self.FlightPhases.add_phase( + self.flight_phases.add_phase( node.t, phase.derivative, clear=True, @@ -740,7 +739,7 @@ def __init__( self, "parachute_cd_s", parachute_cd_s ) ] - self.FlightPhases.add_phase( + self.flight_phases.add_phase( node.t + parachute.lag, self.u_dot_parachute, callbacks, @@ -850,7 +849,7 @@ def __init__( self.out_of_rail_time_index = len(self.solution) - 1 self.out_of_rail_state = self.y_sol # Create new flight phase - self.FlightPhases.add_phase( + self.flight_phases.add_phase( self.t, self.u_dot_generalized, index=phase_index + 1, @@ -884,8 +883,8 @@ def __init__( # Roll back solution self.solution[-1] = [self.t, *self.state] # Set last flight phase - self.FlightPhases.flush_after(phase_index) - self.FlightPhases.add_phase(self.t) + self.flight_phases.flush_after(phase_index) + self.flight_phases.add_phase(self.t) # Prepare to leave loops and start new flight phase phase.TimeNodes.flush_after(node_index) phase.TimeNodes.add_node(self.t, [], []) @@ -947,8 +946,8 @@ def __init__( self.impact_velocity = self.impact_state[5] self.t_final = self.t # Set last flight phase - self.FlightPhases.flush_after(phase_index) - self.FlightPhases.add_phase(self.t) + self.flight_phases.flush_after(phase_index) + self.flight_phases.add_phase(self.t) # Prepare to leave loops and start new flight phase phase.TimeNodes.flush_after(node_index) phase.TimeNodes.add_node(self.t, [], []) @@ -1027,7 +1026,7 @@ def __init__( # Must only be created if parachute has any lag i = 1 if parachute.lag != 0: - self.FlightPhases.add_phase( + self.flight_phases.add_phase( overshootable_node.t, phase.derivative, clear=True, @@ -1040,7 +1039,7 @@ def __init__( self, "parachute_cd_s", parachute_cd_s ) ] - self.FlightPhases.add_phase( + self.flight_phases.add_phase( overshootable_node.t + parachute.lag, self.u_dot_parachute, callbacks, @@ -2843,9 +2842,9 @@ def retrieve_acceleration_arrays(self): 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): + for phase_index, phase in self.time_iterator(self.flight_phases): init_time = phase.t - final_time = self.FlightPhases[phase_index + 1].t + final_time = self.flight_phases[phase_index + 1].t current_derivative = phase.derivative # Call callback functions for callback in phase.callbacks: @@ -2920,9 +2919,9 @@ def retrieve_temporary_values_arrays(self): # Go through each time step and calculate forces and atmospheric values # Get flight phases - for phase_index, phase in self.time_iterator(self.FlightPhases): + for phase_index, phase in self.time_iterator(self.flight_phases): init_time = phase.t - final_time = self.FlightPhases[phase_index + 1].t + final_time = self.flight_phases[phase_index + 1].t current_derivative = phase.derivative # Call callback functions for callback in phase.callbacks: From 9c6c7dbba30dd1d2c68bdfe4a78d4bf821c3f910 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 16:51:17 -0400 Subject: [PATCH 41/71] MNT: renames TimeNodes to time_nodes --- rocketpy/simulation/flight.py | 48 +++++++++++++++++------------------ 1 file changed, 23 insertions(+), 25 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index db09831e1..fe04b5df3 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -580,8 +580,6 @@ def __init__( ------- None """ - # Fetch helper classes and functions - TimeNodes = self.TimeNodes time_iterator = self.time_iterator # Save rocket, parachutes, environment, maximum simulation time @@ -656,37 +654,37 @@ def __init__( # print('\tTolerances: ', self.rtol, self.atol) # Initialize phase time nodes - phase.TimeNodes = TimeNodes() + phase.time_nodes = self.TimeNodes() # Add first time node to permanent list - phase.TimeNodes.add_node(phase.t, [], []) + phase.time_nodes.add_node(phase.t, [], []) # Add non-overshootable parachute time nodes if self.time_overshoot is False: - phase.TimeNodes.add_parachutes( + phase.time_nodes.add_parachutes( self.parachutes, phase.t, phase.time_bound ) - phase.TimeNodes.add_controllers( + phase.time_nodes.add_controllers( self._controllers, phase.t, phase.time_bound ) # Add lst time node to permanent list - phase.TimeNodes.add_node(phase.time_bound, [], []) + phase.time_nodes.add_node(phase.time_bound, [], []) # Sort time nodes - phase.TimeNodes.sort() + phase.time_nodes.sort() # Merge equal time nodes - phase.TimeNodes.merge() + phase.time_nodes.merge() # Clear triggers from first time node if necessary if phase.clear: - phase.TimeNodes[0].parachutes = [] - phase.TimeNodes[0].callbacks = [] + phase.time_nodes[0].parachutes = [] + phase.time_nodes[0].callbacks = [] # print('\n\tPhase Time Nodes') - # print('\tTime Nodes Length: ', str(len(phase.TimeNodes)), ' | Time Nodes Preview: ', phase.TimeNodes[0:3]) + # print('\tTime Nodes Length: ', str(len(phase.time_nodes)), ' | Time Nodes Preview: ', phase.time_nodes[0:3]) # Iterate through time nodes - for node_index, node in time_iterator(phase.TimeNodes): + for node_index, node in time_iterator(phase.time_nodes): # print('\n\t\tCurrent Time Node') # print('\t\tIndex: ', node_index, ' | Time Node: ', node) # Determine time bound for this time node - node.time_bound = phase.TimeNodes[node_index + 1].t + node.time_bound = phase.time_nodes[node_index + 1].t phase.solver.t_bound = node.time_bound phase.solver._lsoda_solver._integrator.rwork[0] = phase.solver.t_bound phase.solver._lsoda_solver._integrator.call_args[4] = ( @@ -747,8 +745,8 @@ def __init__( index=phase_index + i, ) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Save parachute event self.parachute_events.append([self.t, parachute]) @@ -855,8 +853,8 @@ def __init__( index=phase_index + 1, ) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Check for apogee event @@ -886,8 +884,8 @@ def __init__( self.flight_phases.flush_after(phase_index) self.flight_phases.add_phase(self.t) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Check for impact event if self.y_sol[2] < self.env.elevation: @@ -949,14 +947,14 @@ def __init__( self.flight_phases.flush_after(phase_index) self.flight_phases.add_phase(self.t) # Prepare to leave loops and start new flight phase - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # List and feed overshootable time nodes if self.time_overshoot: # Initialize phase overshootable time nodes - overshootable_nodes = TimeNodes() + overshootable_nodes = self.TimeNodes() # Add overshootable parachute time nodes overshootable_nodes.add_parachutes( self.parachutes, self.solution[-2][0], self.t @@ -1057,8 +1055,8 @@ def __init__( overshootable_nodes.flush_after( overshootable_index ) - phase.TimeNodes.flush_after(node_index) - phase.TimeNodes.add_node(self.t, [], []) + phase.time_nodes.flush_after(node_index) + phase.time_nodes.add_node(self.t, [], []) phase.solver.status = "finished" # Save parachute event self.parachute_events.append( From 41cbcc3c5494449cd2babbf3dd3dbb4cfdb2ff90 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 16:55:21 -0400 Subject: [PATCH 42/71] MNT: don't re-assign the time_iterator function --- rocketpy/simulation/flight.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index fe04b5df3..13aa40659 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -580,10 +580,7 @@ def __init__( ------- None """ - time_iterator = self.time_iterator - - # Save rocket, parachutes, environment, maximum simulation time - # and termination events + # Save arguments self.env = environment self.rocket = rocket self.rail_length = rail_length @@ -621,7 +618,10 @@ def __init__( self.flight_phases.add_phase(self.max_time) # Simulate flight - for phase_index, phase in time_iterator(self.flight_phases): + self.__simulate__(verbose) + + def __simulate__(self, verbose): + for phase_index, phase in self.time_iterator(self.flight_phases): # print('\nCurrent Flight Phase List') # print(self.flight_phases) # print('\n\tCurrent Flight Phase') @@ -680,7 +680,7 @@ def __init__( # print('\tTime Nodes Length: ', str(len(phase.time_nodes)), ' | Time Nodes Preview: ', phase.time_nodes[0:3]) # Iterate through time nodes - for node_index, node in time_iterator(phase.time_nodes): + for node_index, node in self.time_iterator(phase.time_nodes): # print('\n\t\tCurrent Time Node') # print('\t\tIndex: ', node_index, ' | Time Node: ', node) # Determine time bound for this time node @@ -978,7 +978,7 @@ def __init__( for ( overshootable_index, overshootable_node, - ) in time_iterator(overshootable_nodes): + ) in self.time_iterator(overshootable_nodes): # print('\n\t\t\t\tCurrent Overshootable Node') # print('\t\t\t\tIndex: ', overshootable_index, ' | Overshootable Node: ', overshootable_node) # Calculate state at node time From cabdcd2aefb64ea412cef70a2361e0027bf197c2 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 18:35:14 -0400 Subject: [PATCH 43/71] MNT: Refactor Flight class root finding algorithm for rail exit and impact time calculations --- rocketpy/simulation/flight.py | 34 ++++++++++++---------------------- 1 file changed, 12 insertions(+), 22 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 13aa40659..a3ab6fda3 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -818,18 +818,13 @@ def __simulate__(self, verbose): b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + 1e-5 # Find roots - d0 = b**2 - 3 * a * c - d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 - c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3) - t_roots = [] - for k in [0, 1, 2]: - c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k - t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2)) + t_roots = Function.cardanos_root_finding(a, b, c, d) # Find correct root - valid_t_root = [] - for t_root in t_roots: - if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001: - valid_t_root.append(t_root.real) + valid_t_root = [ + t_root.real + for t_root in t_roots + if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001 + ] if len(valid_t_root) > 1: raise ValueError( "Multiple roots found when solving for rail exit time." @@ -914,18 +909,13 @@ def __simulate__(self, verbose): b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) # Find roots - d0 = b**2 - 3 * a * c - d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 - c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3) - t_roots = [] - for k in [0, 1, 2]: - c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k - t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2)) + t_roots = Function.cardanos_root_finding(a, b, c, d) # Find correct root - valid_t_root = [] - for t_root in t_roots: - if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001: - valid_t_root.append(t_root.real) + valid_t_root = [ + t_root.real + for t_root in t_roots + if abs(t_root.imag) < 0.001 and 0 < t_root.real < t1 + ] if len(valid_t_root) > 1: raise ValueError( "Multiple roots found when solving for impact time." From ce32cb14f5eb5587a72ef3cdad290e3a3a77f522 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 19:20:02 -0400 Subject: [PATCH 44/71] MNT: Refactor simulation loop to use Function.calculate_cubic_hermite_coefficients --- rocketpy/simulation/flight.py | 54 ++++++++++++++--------------------- 1 file changed, 22 insertions(+), 32 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a3ab6fda3..a360fab19 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -812,11 +812,15 @@ def __simulate__(self, verbose): self.solution[-2][3] += self.env.elevation self.solution[-1][3] += self.env.elevation # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) - D = float(phase.solver.step_size) - d = float(y0) - c = float(yp0) - b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) - a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + 1e-5 + a, b, c, d = Function.calculate_cubic_hermite_coefficients( + 0, + float(phase.solver.step_size), + y0, + yp0, + y1, + yp1, + ) + a += 1e-5 # TODO: why?? # Find roots t_roots = Function.cardanos_root_finding(a, b, c, d) # Find correct root @@ -884,33 +888,21 @@ def __simulate__(self, verbose): phase.solver.status = "finished" # Check for impact event if self.y_sol[2] < self.env.elevation: - # print('\nPASSIVE EVENT DETECTED') - # print('Rocket Has Reached Ground!') - # Impact reported - # Check exactly when it went out using root finding - # States before and after - # t0 -> 0 - # Disconsider elevation - self.solution[-2][3] -= self.env.elevation - self.solution[-1][3] -= self.env.elevation - # Get points - y0 = self.solution[-2][3] - yp0 = self.solution[-2][6] - t1 = self.solution[-1][0] - self.solution[-2][0] - y1 = self.solution[-1][3] - yp1 = self.solution[-1][6] - # Put elevation back - self.solution[-2][3] += self.env.elevation - self.solution[-1][3] += self.env.elevation + # print('\n>>>PASSIVE EVENT DETECTED: Touchdown!') + # Check exactly when it happened using root finding # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) - D = float(phase.solver.step_size) - d = float(y0) - c = float(yp0) - b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2)) - a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + a, b, c, d = Function.calculate_cubic_hermite_coefficients( + x0=0, # t0 + x1=float(phase.solver.step_size), # t1 - t0 + y0=float(self.solution[-2][3] - self.env.elevation), # z0 + yp0=float(self.solution[-2][6]), # vz0 + y1=float(self.solution[-1][3] - self.env.elevation), # z1 + yp1=float(self.solution[-1][6]), # vz1 + ) # Find roots t_roots = Function.cardanos_root_finding(a, b, c, d) # Find correct root + t1 = self.solution[-1][0] - self.solution[-2][0] valid_t_root = [ t_root.real for t_root in t_roots @@ -921,18 +913,16 @@ def __simulate__(self, verbose): "Multiple roots found when solving for impact time." ) # Determine impact state at t_root - self.t = valid_t_root[0] + self.solution[-2][0] + self.t = self.t_final = valid_t_root[0] + self.solution[-2][0] interpolator = phase.solver.dense_output() - self.y_sol = interpolator(self.t) + self.y_sol = self.impact_state = interpolator(self.t) # Roll back solution self.solution[-1] = [self.t, *self.y_sol] # Save impact state - self.impact_state = self.y_sol self.x_impact = self.impact_state[0] self.y_impact = self.impact_state[1] self.z_impact = self.impact_state[2] self.impact_velocity = self.impact_state[5] - self.t_final = self.t # Set last flight phase self.flight_phases.flush_after(phase_index) self.flight_phases.add_phase(self.t) From 48ac2f18bc7ad2ebda5c395339a80b397f6aa095 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 19:21:15 -0400 Subject: [PATCH 45/71] TST: adds tests for new Function calculations --- tests/unit/test_function.py | 40 +++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index 17da59498..633da2cf6 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -334,3 +334,43 @@ def test_get_domain_dim(linear_func): def test_bool(linear_func): """Test the __bool__ method of the Function class.""" assert bool(linear_func) == True + + +def test_calculate_cubic_hermite_coefficients(): + """Test the calculate_cubic_hermite_coefficients method of the Function class.""" + # Function: f(x) = x**3 + 2x**2 -1 ; derivative: f'(x) = 3x**2 + 4x + x = np.array([-3, -2, -1, 0, 1]) + y = np.array([-10, -1, 0, -1, 2]) + + # Selects two points as x0 and x1 + x0, x1 = 0, 1 + y0, y1 = -1, 2 + yp0, yp1 = 0, 7 + + a, b, c, d = Function.calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1) + + assert np.isclose(a, 1) + assert np.isclose(b, 2) + assert np.isclose(c, 0) + assert np.isclose(d, -1) + assert np.allclose( + a * x**3 + b * x**2 + c * x + d, + y, + ) + + +def test_cardanos_root_finding(): + """Tests the cardanos_root_finding method of the Function class.""" + # Function: f(x) = x**3 + 2x**2 -1 + # roots: (-1 - 5**0.5) / 2; -1; (-1 + 5**0.5) / 2 + + roots = list(Function.cardanos_root_finding(a=1, b=2, c=0, d=-1)) + roots.sort(key=lambda x: x.real) + + assert np.isclose(roots[0].real, (-1 - 5**0.5) / 2) + assert np.isclose(roots[1].real, -1) + assert np.isclose(roots[2].real, (-1 + 5**0.5) / 2) + + assert np.isclose(roots[0].imag, 0) + assert np.isclose(roots[1].imag, 0) + assert np.isclose(roots[2].imag, 0) From af31b94d1e1389137a524b0d16fe1def299d958a Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 19:21:41 -0400 Subject: [PATCH 46/71] ENH: cast float type in the Hermite interpolation function in Function class --- rocketpy/mathutils/function.py | 58 +++++++++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 4 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index ba90b37f1..d2cf6c0c5 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2902,10 +2902,10 @@ def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): The coefficients of the cubic Hermite interpolation function. """ dx = x1 - x0 - d = y0 - c = yp0 - b = (3 * y1 - yp1 * dx - 2 * c * dx - 3 * d) / (dx**2) - a = -(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3) + d = float(y0) + c = float(yp0) + b = float((3 * y1 - yp1 * dx - 2 * c * dx - 3 * d) / (dx**2)) + a = float(-(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3)) return a, b, c, d @staticmethod @@ -2966,8 +2966,58 @@ def cardanos_root_finding(a, b, c, d): c2_2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 2 x3 = -(1 / (3 * a)) * (b + c2_2 + delta_0 / c2_2) + # Alternative method (same results) + + # Q = (3 * a * c - b**2) / (9 * a**2) + # R = (9 * a * b * c - 27 * a**2 * d - 2 * b**3) / (54 * a**3) + # S = (R + (Q**3 + R**2)**0.5)**(1 / 3) + # T = (R - (Q**3 + R**2)**0.5)**(1 / 3) + + # x1 = S + T - b / (3 * a) + # x2 = -(S + T) / 2 - b / (3 * a) + 1j * (S - T) * (3**0.5) / 2 + # x3 = -(S + T) / 2 - b / (3 * a) - 1j * (S - T) * (3**0.5) / 2 + return x1, x2, x3 + @staticmethod + def find_root_linear_interpolation(x0, x1, y0, y1, y): + """Calculate the root of a linear interpolation function. + + This method calculates the root of a linear interpolation function + given two points (x0, y0) and (x1, y1) and a value y. The function + is defined as y = m*x + c. + + Parameters + ---------- + x0 : float + Position of the first point. + x1 : float + Position of the second point. + y0 : float + Value of the function evaluated at the first point. + y1 : float + Value of the function evaluated at the second point. + y : float + Value of the function at the desired point. + + Returns + ------- + float + The root of the linear interpolation function. This represents the + value of x at which the function evaluates to y. + + Examples + -------- + >>> from rocketpy import Function + >>> x0, x1, y0, y1, y = 0, 1, 0, 1, 0.5 + >>> x = Function.linear_interpolation_root_finding(x0, x1, y0, y1, y) + >>> x + 0.5 + """ + m = (y1 - y0) / (x1 - x0) + c = y0 - m * x0 + return (y - c) / m + # Input validators def __validate_source(self, source): """Used to validate the source parameter for creating a Function object. From 2602e4d1454ee10d81559b28f1ea4b774bb89ff3 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 19:43:12 -0400 Subject: [PATCH 47/71] ENH: improve Flight.TimeNodes.merge method --- rocketpy/simulation/flight.py | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a360fab19..419ab75fe 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3633,19 +3633,22 @@ def sort(self): self.list.sort(key=(lambda node: node.t)) def merge(self): - # Initialize temporary list - tmp_list = [self.list[0]] - # Iterate through all other time nodes - for node in self.list[1:]: - # If there is already another node with similar time: merge - if abs(node.t - tmp_list[-1].t) < 1e-7: - tmp_list[-1].parachutes += node.parachutes - tmp_list[-1].callbacks += node.callbacks - # Add new node to tmp list if there is none with the same time - else: - tmp_list.append(node) - # Save tmp list to permanent - self.list = tmp_list + """Merge all the time nodes that have the same time. This is made to + avoid multiple evaluations of the same time node. This method does + not guarantee the order of the nodes in the list, so it is + recommended to sort the list before or after using this method. + """ + tmp_dict = {} + for node in self.list: + time = round(node.t, 7) + try: + # Try to access the node and merge if it exists + tmp_dict[time].parachutes += node.parachutes + tmp_dict[time].callbacks += node.callbacks + except KeyError: + # If the node does not exist, add it to the dictionary + tmp_dict[time] = node + self.list = list(tmp_dict.values()) def flush_after(self, index): del self.list[index + 1 :] From b55236e78a243b94f036b19d77a92d85342e47e2 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 19:44:54 -0400 Subject: [PATCH 48/71] MNT: rename cardanos_root_finding to find_roots_cubic_function in Function class --- rocketpy/mathutils/function.py | 4 ++-- rocketpy/simulation/flight.py | 4 ++-- tests/unit/test_function.py | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index d2cf6c0c5..8f96d49c5 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2909,7 +2909,7 @@ def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): return a, b, c, d @staticmethod - def cardanos_root_finding(a, b, c, d): + def find_roots_cubic_function(a, b, c, d): """Calculate the roots of a cubic function using Cardano's method. This method applies Cardano's method to find the roots of a cubic @@ -2944,7 +2944,7 @@ def cardanos_root_finding(a, b, c, d): First we define the coefficients of the function ax**3 + bx**2 + cx + d >>> a, b, c, d = 1, -3, -1, 3 - >>> x1, x2, x3 = Function.cardanos_root_finding(a, b, c, d) + >>> x1, x2, x3 = Function.find_roots_cubic_function(a, b, c, d) >>> x1, x2, x3 ((-1+0j), (3+7.401486830834377e-17j), (1-1.4802973661668753e-16j)) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 419ab75fe..2c24e70e5 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -822,7 +822,7 @@ def __simulate__(self, verbose): ) a += 1e-5 # TODO: why?? # Find roots - t_roots = Function.cardanos_root_finding(a, b, c, d) + t_roots = Function.find_roots_cubic_function(a, b, c, d) # Find correct root valid_t_root = [ t_root.real @@ -900,7 +900,7 @@ def __simulate__(self, verbose): yp1=float(self.solution[-1][6]), # vz1 ) # Find roots - t_roots = Function.cardanos_root_finding(a, b, c, d) + t_roots = Function.find_roots_cubic_function(a, b, c, d) # Find correct root t1 = self.solution[-1][0] - self.solution[-2][0] valid_t_root = [ diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index 633da2cf6..35a5b21e9 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -360,11 +360,11 @@ def test_calculate_cubic_hermite_coefficients(): def test_cardanos_root_finding(): - """Tests the cardanos_root_finding method of the Function class.""" + """Tests the find_roots_cubic_function method of the Function class.""" # Function: f(x) = x**3 + 2x**2 -1 # roots: (-1 - 5**0.5) / 2; -1; (-1 + 5**0.5) / 2 - roots = list(Function.cardanos_root_finding(a=1, b=2, c=0, d=-1)) + roots = list(Function.find_roots_cubic_function(a=1, b=2, c=0, d=-1)) roots.sort(key=lambda x: x.real) assert np.isclose(roots[0].real, (-1 - 5**0.5) / 2) From 28d0c62cc85edf4333bc232aea3be8623e57c68c Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 20:55:47 -0400 Subject: [PATCH 49/71] DOC: representation methods, TODOs and docstring --- rocketpy/simulation/flight.py | 188 +++++++++++++++++----------------- 1 file changed, 92 insertions(+), 96 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 2c24e70e5..f929a8244 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -182,7 +182,7 @@ class Flight: Flight.function_evaluations : array List that stores number of derivative function evaluations during numerical integration in cumulative manner. - Flight.function_evaluations_per_time_step : array + Flight.function_evaluations_per_time_step : list List that stores number of derivative function evaluations per time step during numerical integration. Flight.time_steps : array @@ -606,10 +606,6 @@ def __init__( self.__init_equations_of_motion() self.__init_solver_monitors() - # Initialize prints and plots objects - self.prints = _FlightPrints(self) - self.plots = _FlightPlots(self) - # Create known flight phases self.flight_phases = self.FlightPhases() self.flight_phases.add_phase( @@ -618,14 +614,29 @@ def __init__( self.flight_phases.add_phase(self.max_time) # Simulate flight - self.__simulate__(verbose) + self.__simulate(verbose) + + # Initialize prints and plots objects + self.prints = _FlightPrints(self) + self.plots = _FlightPlots(self) - def __simulate__(self, verbose): + def __repr__(self): + return ( + f"" + ) + + def __simulate(self, verbose): + """Simulate the flight trajectory.""" for phase_index, phase in self.time_iterator(self.flight_phases): - # print('\nCurrent Flight Phase List') + # print("\nCurrent Flight Phase List:") # print(self.flight_phases) - # print('\n\tCurrent Flight Phase') - # print('\tIndex: ', phase_index, ' | Phase: ', phase) + # print("\n\tCurrent Flight Phase:") + # print("\t\tIndex: ", phase_index, " | Phase: ", phase) # Determine maximum time for this flight phase phase.time_bound = self.flight_phases[phase_index + 1].t @@ -633,7 +644,7 @@ def __simulate__(self, verbose): for callback in phase.callbacks: callback(self) - # Create solver for this flight phase + # Create solver for this flight phase # TODO: allow different integrators self.function_evaluations.append(0) phase.solver = integrate.LSODA( phase.derivative, @@ -645,46 +656,45 @@ def __simulate__(self, verbose): rtol=self.rtol, atol=self.atol, ) - # print('\n\tSolver Initialization Details') - # print('\t_initial Time: ', phase.t) - # print('\t_initial State: ', self.y_sol) - # print('\tTime Bound: ', phase.time_bound) - # print('\tMin Step: ', self.min_time_step) - # print('\tMax Step: ', self.max_time_step) - # print('\tTolerances: ', self.rtol, self.atol) + # print(f"\n\tSolver Initialization Details:") + # print(f"\t\tInitial Time (t): {phase.t}") + # print(f"\t\tInitial State (y_sol): {self.y_sol}") + # print(f"\t\tTime Bound: {phase.time_bound}") + # print(f"\t\tMin Step: {self.min_time_step}") + # print(f"\t\tMax Step: {self.max_time_step}") + # print(f"\t\tTolerances: rtol = {self.rtol}, atol = {self.atol}") # Initialize phase time nodes phase.time_nodes = self.TimeNodes() - # Add first time node to permanent list + # Add first time node to the time_nodes list phase.time_nodes.add_node(phase.t, [], []) # Add non-overshootable parachute time nodes if self.time_overshoot is False: + # TODO: move parachutes to controllers phase.time_nodes.add_parachutes( self.parachutes, phase.t, phase.time_bound ) phase.time_nodes.add_controllers( self._controllers, phase.t, phase.time_bound ) - # Add lst time node to permanent list + # Add last time node to the time_nodes list phase.time_nodes.add_node(phase.time_bound, [], []) - # Sort time nodes + # Organize time nodes with sort() and merge() phase.time_nodes.sort() - # Merge equal time nodes phase.time_nodes.merge() # Clear triggers from first time node if necessary if phase.clear: phase.time_nodes[0].parachutes = [] phase.time_nodes[0].callbacks = [] - # print('\n\tPhase Time Nodes') - # print('\tTime Nodes Length: ', str(len(phase.time_nodes)), ' | Time Nodes Preview: ', phase.time_nodes[0:3]) - # Iterate through time nodes for node_index, node in self.time_iterator(phase.time_nodes): - # print('\n\t\tCurrent Time Node') - # print('\t\tIndex: ', node_index, ' | Time Node: ', node) + # print("\n\t\tCurrent Time Node") + # print("\t\tIndex: ", node_index, " | Time Node: ", node) # Determine time bound for this time node node.time_bound = phase.time_nodes[node_index + 1].t + # NOTE: Setting the time bound and status for the phase solver, + # and updating its internal state for the next integration step. phase.solver.t_bound = node.time_bound phase.solver._lsoda_solver._integrator.rwork[0] = phase.solver.t_bound phase.solver._lsoda_solver._integrator.call_args[4] = ( @@ -693,6 +703,7 @@ def __simulate__(self, verbose): phase.solver.status = "running" # Feed required parachute and discrete controller triggers + # TODO: parachutes should be moved to controllers for callback in node.callbacks: callback(self) @@ -715,9 +726,8 @@ def __simulate__(self, verbose): if parachute.triggerfunc( pressure + noise, height_above_ground_level, self.y_sol ): - # print('\nEVENT DETECTED') - # print('Parachute Triggered') - # print('Name: ', parachute.name, ' | Lag: ', parachute.lag) + # print("\nEVENT DETECTED: Parachute Triggered") + # print("Name: ", parachute.name, " | Lag: ", parachute.lag) # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create flight phase for time after detection and before inflation @@ -777,12 +787,8 @@ def __simulate__(self, verbose): + (self.y_sol[2] - self.env.elevation) ** 2 >= self.effective_1rl**2 ): - # Rocket is out of rail + # print("\n>>> EVENT DETECTED: Rocket is Out of Rail!") # Check exactly when it went out using root finding - # States before and after - # t0 -> 0 - # print('\nEVENT DETECTED') - # print('Rocket is Out of Rail!') # Disconsider elevation self.solution[-2][3] -= self.env.elevation self.solution[-1][3] -= self.env.elevation @@ -858,9 +864,7 @@ def __simulate__(self, verbose): # Check for apogee event if len(self.apogee_state) == 1 and self.y_sol[5] < 0: - # print('\nPASSIVE EVENT DETECTED') - # print('Rocket Has Reached Apogee!') - # Apogee reported + # print("\n>>> EVENT DETECTED: Rocket has reached apogee") # Assume linear vz(t) to detect when vz = 0 vz0 = self.solution[-2][6] t0 = self.solution[-2][0] @@ -873,10 +877,7 @@ def __simulate__(self, verbose): # Store apogee data self.apogee_time = t_root if self.terminate_on_apogee: - # print('Terminate on Apogee Activated!') - self.t = t_root - self.t_final = self.t - self.state = self.apogee_state + # print(">>> Terminate on apogee activated!") # Roll back solution self.solution[-1] = [self.t, *self.state] # Set last flight phase @@ -942,25 +943,34 @@ def __simulate__(self, verbose): # Add last time node (always skipped) overshootable_nodes.add_node(self.t, [], []) if len(overshootable_nodes) > 1: - # Sort overshootable time nodes + # Sort and merge equal overshootable time nodes overshootable_nodes.sort() - # Merge equal overshootable time nodes overshootable_nodes.merge() # Clear if necessary if overshootable_nodes[0].t == phase.t and phase.clear: overshootable_nodes[0].parachutes = [] overshootable_nodes[0].callbacks = [] - # print('\n\t\t\t\tOvershootable Time Nodes') - # print('\t\t\t\tInterval: ', self.solution[-2][0], self.t) - # print('\t\t\t\tOvershootable Nodes Length: ', str(len(overshootable_nodes)), ' | Overshootable Nodes: ', overshootable_nodes) + # print("\n\t\t\t\tOvershootable Time Nodes") + # print("\t\t\t\tInterval: ", self.solution[-2][0], self.t) + # print( + # "\t\t\t\tOvershootable Nodes Length: ", + # str(len(overshootable_nodes)), + # " | Overshootable Nodes: ", + # overshootable_nodes, + # ) # Feed overshootable time nodes trigger interpolator = phase.solver.dense_output() for ( overshootable_index, overshootable_node, ) in self.time_iterator(overshootable_nodes): - # print('\n\t\t\t\tCurrent Overshootable Node') - # print('\t\t\t\tIndex: ', overshootable_index, ' | Overshootable Node: ', overshootable_node) + # print("\n\t\t\t\tCurrent Overshootable Node:") + # print( + # "\t\t\t\tIndex: ", + # overshootable_index, + # " | Overshootable Node: ", + # overshootable_node, + # ) # Calculate state at node time overshootable_node.y = interpolator( overshootable_node.t @@ -995,9 +1005,12 @@ def __simulate__(self, verbose): height_above_ground_level, overshootable_node.y, ): - # print('\nEVENT DETECTED') - # print('Parachute Triggered') - # print('Name: ', parachute.name, ' | Lag: ', parachute.lag) + # print( + # "\n>>> EVENT DETECTED: Parachute Triggered!" + # ) + # print( + # f"\tParachute Name: {parachute.name} | Lag: {parachute.lag}" + # ) # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create flight phase for time after detection and before inflation @@ -1325,7 +1338,6 @@ def u_dot(self, t, u, post_processing=False): M2 -= self.rocket.thrust_eccentricity_y * thrust else: # Motor stopped - # Retrieve important motor quantities # Inertias Tz, Ti, Tzdot, Tidot = 0, 0, 0, 0 # Mass @@ -1405,7 +1417,6 @@ def u_dot(self, t, u, post_processing=False): R3 = air_brakes_force # Substitutes rocket drag coefficient else: R3 += air_brakes_force - # R3 += self.__computeDragForce(z, Vector(vx, vy, vz)) # Off center moment M1 += self.rocket.cp_eccentricity_y * R3 M2 -= self.rocket.cp_eccentricity_x * R3 @@ -2303,6 +2314,7 @@ def attitude_angle(self): @funcify_method("Time (s)", "Lateral Attitude Angle (°)") def lateral_attitude_angle(self): """Rocket lateral attitude angle as a Function of time.""" + # TODO: complex method, it should be defined elsewhere. lateral_vector_angle = (np.pi / 180) * (self.heading - 90) lateral_vector_x = np.sin(lateral_vector_angle) lateral_vector_y = np.cos(lateral_vector_angle) @@ -2560,16 +2572,15 @@ def potential_energy(self): """Potential energy as a Function of time in relation to sea level.""" # Constants - GM = 3.986004418e14 + GM = 3.986004418e14 # TODO: this constant should come from Environment. # Redefine total_mass time grid to allow for efficient Function algebra total_mass = deepcopy(self.rocket.total_mass) total_mass.set_discrete_based_on_model(self.z) - potential_energy = ( + return ( GM * total_mass * (1 / (self.z + self.env.earth_radius) - 1 / self.env.earth_radius) ) - return potential_energy # Total Mechanical Energy @funcify_method("Time (s)", "Mechanical Energy (J)", "spline", "constant") @@ -2771,6 +2782,7 @@ def latitude(self): ) return np.column_stack((self.time, latitude)) + # TODO: haversine should be defined in tools.py so we just invoke it in here. @funcify_method("Time (s)", "Longitude (°)", "linear", "constant") def longitude(self): """Rocket longitude coordinate, in degrees, as a Function of @@ -2944,7 +2956,7 @@ def get_controller_observed_variables(self): ) @cached_property - def __calculate_rail_button_forces(self): + def __calculate_rail_button_forces(self): # TODO: complex method. """Calculate the forces applied to the rail buttons while rocket is still on the launch rail. It will return 0 if no rail buttons are defined. @@ -3070,12 +3082,10 @@ def post_process(self, interpolation="spline", extrapolation="natural"): ------- None """ - # Register post processing + # TODO: add a deprecation warning maybe? self.post_processed = True - return None - - def calculate_stall_wind_velocity(self, stall_angle): + def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities """Function to calculate the maximum wind velocity before the angle of attack exceeds a desired angle, at the instant of departing rail launch. Can be helpful if you know the exact stall angle of all aerodynamics @@ -3086,6 +3096,7 @@ def calculate_stall_wind_velocity(self, stall_angle): stall_angle : float Angle, in degrees, for which you would like to know the maximum wind speed before the angle of attack exceeds it + Return ------ None @@ -3240,13 +3251,13 @@ class attributes which are instances of the Function class. Usage try: obj = getattr(self.__class__, variable) variable_function = obj.__get__(self, self.__class__) - except AttributeError: + except AttributeError as exc: raise AttributeError( - "Variable '{}' not found in Flight class".format(variable) - ) + f"Variable '{variable}' not found in Flight class" + ) from exc variable_points = variable_function(time_points) exported_matrix += [variable_points] - exported_header += ", " + variable_function.__outputs__[0] + exported_header += f", {variable_function.__outputs__[0]}" exported_matrix = np.array(exported_matrix).T # Fix matrix orientation @@ -3259,9 +3270,7 @@ class attributes which are instances of the Function class. Usage encoding="utf-8", ) - return - - def export_kml( + def export_kml( # TODO: should be moved out of this class. self, file_name="trajectory.kml", time_step=None, @@ -3347,29 +3356,13 @@ def export_kml( kml.save(file_name) print("File ", file_name, " saved with success!") - return None - def info(self): - """Prints out a summary of the data available about the Flight. - - Returns - ------- - None - """ + """Prints out a summary of the data available about the Flight.""" self.prints.all() - return None def all_info(self): - """Prints out all data and graphs available about the Flight. - - Returns - ------- - None - """ - - # Print a summary of data about the flight + """Prints out all data and graphs available about the Flight.""" self.info() - self.plots.all() return None @@ -3409,7 +3402,7 @@ def display_warning(self, *messages): """A simple function to print a warning message.""" print("WARNING:", *messages) - def add(self, flight_phase, index=None): + def add(self, flight_phase, index=None): # TODO: quite complex method """Add a flight phase to the list. It will be inserted in the correct position, according to its initial time. If no index is provided, it will be appended to the end of the list. If by any @@ -3566,6 +3559,8 @@ class FlightPhase: A flag indicating whether to clear the solution after the phase. """ + # TODO: add a "name" optional argument to the FlightPhase. Really helps. + def __init__(self, t, derivative=None, callbacks=None, clear=True): self.t = t self.derivative = derivative @@ -3573,17 +3568,19 @@ def __init__(self, t, derivative=None, callbacks=None, clear=True): self.clear = clear def __repr__(self): - if self.derivative is None: - return "{Initial Time: " + str(self.t) + " | Derivative: None}" + name = "None" if self.derivative is None else self.derivative.__name__ return ( - "{Initial Time: " - + str(self.t) - + " | Derivative: " - + self.derivative.__name__ - + "}" + f"" ) class TimeNodes: + """TimeNodes is a class that stores all the time nodes of a simulation. + It is meant to work like a python list, but it has some additional + methods that are useful for the simulation. Them items stored in are + TimeNodes object are instances of the TimeNode class. + """ + def __init__(self, init_list=[]): self.list = init_list[:] @@ -3616,7 +3613,6 @@ def add_parachutes(self, parachutes, t_init, t_end): self.list += parachute_node_list def add_controllers(self, controllers, t_init, t_end): - # Iterate over controllers for controller in controllers: # Calculate start of sampling time nodes controller_time_step = 1 / controller.sampling_rate From 1b37951e203dda4b2c1a188b84a8bd2a63375cfe Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 20:56:18 -0400 Subject: [PATCH 50/71] ENH: providing lambda to sort() is no longer needed --- 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 f929a8244..c519a437f 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3626,7 +3626,7 @@ def add_controllers(self, controllers, t_init, t_end): self.list += controller_node_list def sort(self): - self.list.sort(key=(lambda node: node.t)) + self.list.sort() def merge(self): """Merge all the time nodes that have the same time. This is made to From ce681506c4fcb1c6146370db7e92de8ce9696f7f Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 20:58:19 -0400 Subject: [PATCH 51/71] ENH: new Flight.__calculate_and_save_pressure_signals method --- rocketpy/simulation/flight.py | 91 ++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 38 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index c519a437f..570796aaf 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -712,19 +712,13 @@ def __simulate(self, verbose): for parachute in node.parachutes: # Calculate and save pressure signal - pressure = self.env.pressure.get_value_opt(self.y_sol[2]) - parachute.clean_pressure_signal.append([node.t, pressure]) - # Calculate and save noise - noise = parachute.noise_function() - parachute.noise_signal.append([node.t, noise]) - parachute.noisy_pressure_signal.append([node.t, pressure + noise]) - # Gets height above ground level considering noise - height_above_ground_level = ( - self.env.barometric_height.get_value_opt(pressure + noise) - - self.env.elevation + noisy_pressure, height_above_ground_level = ( + self.__calculate_and_save_pressure_signals( + parachute, node.t, self.y_sol[2] + ) ) if parachute.triggerfunc( - pressure + noise, height_above_ground_level, self.y_sol + noisy_pressure, height_above_ground_level, self.y_sol ): # print("\nEVENT DETECTED: Parachute Triggered") # print("Name: ", parachute.name, " | Lag: ", parachute.lag) @@ -972,38 +966,24 @@ def __simulate(self, verbose): # overshootable_node, # ) # Calculate state at node time - overshootable_node.y = interpolator( + overshootable_node.y_sol = interpolator( overshootable_node.t ) - # Calculate and save pressure signal - pressure = self.env.pressure.get_value_opt( - overshootable_node.y[2] - ) for parachute in overshootable_node.parachutes: - # Save pressure signal - parachute.clean_pressure_signal.append( - [overshootable_node.t, pressure] - ) - # Calculate and save noise - noise = parachute.noise_function() - parachute.noise_signal.append( - [overshootable_node.t, noise] - ) - parachute.noisy_pressure_signal.append( - [overshootable_node.t, pressure + noise] - ) - # Gets height above ground level considering noise - height_above_ground_level = ( - self.env.barometric_height.get_value_opt( - pressure + noise + # Calculate and save pressure signal + noisy_pressure, height_above_ground_level = ( + self.__calculate_and_save_pressure_signals( + parachute, + overshootable_node.t, + overshootable_node.y_sol[2], ) - - self.env.elevation ) + # Check for parachute trigger if parachute.triggerfunc( - pressure + noise, + noisy_pressure, height_above_ground_level, - overshootable_node.y, + overshootable_node.y_sol, ): # print( # "\n>>> EVENT DETECTED: Parachute Triggered!" @@ -1039,10 +1019,10 @@ def __simulate(self, verbose): ) # Rollback history self.t = overshootable_node.t - self.y_sol = overshootable_node.y + self.y_sol = overshootable_node.y_sol self.solution[-1] = [ overshootable_node.t, - *overshootable_node.y, + *overshootable_node.y_sol, ] # Prepare to leave loops and start new flight phase overshootable_nodes.flush_after( @@ -1059,7 +1039,42 @@ def __simulate(self, verbose): self.t_final = self.t self._calculate_pressure_signal() if verbose: - print("Simulation Completed at Time: {:3.4f} s".format(self.t)) + print(f"\n>>> Simulation Completed at Time: {self.t:3.4f} s") + + def __calculate_and_save_pressure_signals(self, parachute, t, z): + """Gets noise and pressure signals and saves them in the parachute + object given the current time and altitude. + + Parameters + ---------- + parachute : Parachute + The parachute object to calculate signals for. + t : float + The current time in seconds. + z : float + The altitude above sea level in meters. + + Returns + ------- + tuple[float, float] + The noisy pressure and height above ground level. + """ + # Calculate pressure and noise + pressure = self.env.pressure.get_value_opt(z) + noise = parachute.noise_function() + noisy_pressure = pressure + noise + + # Stores in the parachute object + parachute.clean_pressure_signal.append([t, pressure]) + parachute.noise_signal.append([t, noise]) + + # Gets height above ground level considering noise + height_above_ground_level = ( + self.env.barometric_height.get_value_opt(noisy_pressure) + - self.env.elevation + ) + + return noisy_pressure, height_above_ground_level def __init_solution_monitors(self): # Initialize solution monitors From 8f6efdf6b03120fc16076ecf5cf723915a33f4df Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 20:58:48 -0400 Subject: [PATCH 52/71] ENH: use Function.find_root_linear_interpolation in the simulation loop --- rocketpy/simulation/flight.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 570796aaf..cf479b22d 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -860,11 +860,11 @@ def __simulate(self, verbose): if len(self.apogee_state) == 1 and self.y_sol[5] < 0: # print("\n>>> EVENT DETECTED: Rocket has reached apogee") # Assume linear vz(t) to detect when vz = 0 - vz0 = self.solution[-2][6] - t0 = self.solution[-2][0] - vz1 = self.solution[-1][6] - t1 = self.solution[-1][0] - t_root = -(t1 - t0) * vz0 / (vz1 - vz0) + t0 + t0, vz0 = self.solution[-2][0], self.solution[-2][6] + t1, vz1 = self.solution[-1][0], self.solution[-1][6] + t_root = Function.find_root_linear_interpolation( + t0, t1, vz0, vz1, 0 + ) # Fetch state at t_root interpolator = phase.solver.dense_output() self.apogee_state = interpolator(t_root) @@ -872,8 +872,9 @@ def __simulate(self, verbose): self.apogee_time = t_root if self.terminate_on_apogee: # print(">>> Terminate on apogee activated!") + self.t = self.t_final = t_root # Roll back solution - self.solution[-1] = [self.t, *self.state] + self.solution[-1] = [self.t, *self.apogee_state] # Set last flight phase self.flight_phases.flush_after(phase_index) self.flight_phases.add_phase(self.t) From 39003cd5552e4d2629827f24300a75234f7efbe9 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 20:59:21 -0400 Subject: [PATCH 53/71] MNT: rename __transform_pressure_signals_lists_to_functions method --- 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 cf479b22d..97cc885fb 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1038,7 +1038,7 @@ def __simulate(self, verbose): ) self.t_final = self.t - self._calculate_pressure_signal() + self.__transform_pressure_signals_lists_to_functions() if verbose: print(f"\n>>> Simulation Completed at Time: {self.t:3.4f} s") @@ -3059,7 +3059,7 @@ def __calculate_rail_button_forces(self): # TODO: complex method. rail_button2_shear_force, ) - def _calculate_pressure_signal(self): + def __transform_pressure_signals_lists_to_functions(self): """Calculate the pressure signal from the pressure sensor. It creates a signal_function attribute in the parachute object. Parachute works as a subclass of Rocket class. From b8660dbe84a67df1ed8532f007468cb5cb25e62d Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 21:00:50 -0400 Subject: [PATCH 54/71] ENH: optimize post_process mode --- rocketpy/simulation/flight.py | 106 +++++----------------------------- 1 file changed, 13 insertions(+), 93 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 97cc885fb..061533695 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1570,28 +1570,7 @@ def u_dot(self, t, u, post_processing=False): ] if post_processing: - # Dynamics variables - self.R1_list.append([t, R1]) - self.R2_list.append([t, R2]) - self.R3_list.append([t, R3]) - self.M1_list.append([t, M1]) - self.M2_list.append([t, M2]) - self.M3_list.append([t, M3]) - # Atmospheric Conditions - 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)] - ) + u_dot.extend([R1, R2, R3, M1, M2, M3]) return u_dot @@ -1847,28 +1826,7 @@ def u_dot_generalized(self, t, u, post_processing=False): u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] if post_processing: - # Dynamics variables - self.R1_list.append([t, R1]) - self.R2_list.append([t, R2]) - self.R3_list.append([t, R3]) - self.M1_list.append([t, M1]) - self.M2_list.append([t, M2]) - self.M3_list.append([t, M3]) - # Atmospheric Conditions - 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)] - ) + u_dot.extend([R1, R2, R3, M1, M2, M3]) return u_dot @@ -1931,28 +1889,7 @@ def u_dot_parachute(self, t, u, post_processing=False): az = (Dz - 9.8 * mp) / (mp + ma) if post_processing: - # Dynamics variables - self.R1_list.append([t, Dx]) - self.R2_list.append([t, Dy]) - self.R3_list.append([t, Dz]) - self.M1_list.append([t, 0]) - 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.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, Dx, Dy, Dz, 0, 0, 0] return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] @@ -2910,18 +2847,14 @@ def retrieve_temporary_values_arrays(self): """ # 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 = [] + temporary_values = [ + [], # R1_list + [], # R2_list + [], # R3_list + [], # M1_list + [], # M2_list + [], # M3_list + ] # Go through each time step and calculate forces and atmospheric values # Get flight phases @@ -2939,21 +2872,8 @@ def retrieve_temporary_values_arrays(self): ): # 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, - ] + for i, value in enumerate(temporary_values): + value.append(u_dot[i + 13]) return temporary_values From f84a03fadf17be8e4e6c4901a34fb5f8d187faa9 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 21:01:22 -0400 Subject: [PATCH 55/71] MNT: small fixes to the Flight class --- rocketpy/simulation/flight.py | 114 ++++++++++++++-------------------- 1 file changed, 47 insertions(+), 67 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 061533695..43ab16ef5 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1333,8 +1333,7 @@ def u_dot(self, t, u, post_processing=False): # Retrieve integration data x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Determine lift force and moment - R1, R2 = 0, 0 - M1, M2, M3 = 0, 0, 0 + R1, R2, M1, M2, M3 = 0, 0, 0, 0, 0 # Determine current behavior if t < self.rocket.motor.burn_out_time: # Motor burning @@ -1397,8 +1396,6 @@ def u_dot(self, t, u, post_processing=False): a33 = 1 - 2 * (e1**2 + e2**2) # Transformation matrix: (123) -> (XYZ) K = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]] - # Transformation matrix: (XYZ) -> (123) or K transpose - Kt = [[a11, a21, a31], [a12, a22, a32], [a13, a23, a33]] # Calculate Forces and Moments # Get freestream speed @@ -1633,11 +1630,8 @@ def u_dot_generalized(self, t, u, post_processing=False): * self.rocket._csys ) S_noz_33 = 0.5 * self.rocket.motor.nozzle_radius**2 - S_noz_11 = 0.5 * S_noz_33 + 0.25 * r_NOZ**2 - S_noz_22 = S_noz_11 - S_noz_12 = 0 - S_noz_13 = 0 - S_noz_23 = 0 + 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], @@ -1690,7 +1684,9 @@ def u_dot_generalized(self, t, u, post_processing=False): 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))) - free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z) + speed_of_sound = self.env.speed_of_sound.get_value_opt(z) + free_stream_mach = free_stream_speed / speed_of_sound + if t < self.rocket.motor.burn_out_time: drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach) else: @@ -1737,13 +1733,11 @@ def u_dot_generalized(self, t, u, post_processing=False): comp_stream_velocity = comp_wind_vb - comp_vb comp_stream_vx_b, comp_stream_vy_b, comp_stream_vz_b = comp_stream_velocity comp_stream_speed = abs(comp_stream_velocity) - comp_stream_mach = ( - comp_stream_speed / self.env.speed_of_sound.get_value_opt(z) - ) + comp_stream_mach = comp_stream_speed / speed_of_sound # Component attack angle and lift force comp_attack_angle = 0 comp_lift, comp_lift_xb, comp_lift_yb = 0, 0, 0 - if comp_stream_vx_b**2 + comp_stream_vy_b**2 != 0: + if comp_stream_vx_b**2 + comp_stream_vy_b**2 != 0: # TODO: maybe try/except # Normalize component stream velocity in body frame comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed if -1 * comp_stream_vz_bn < 1: @@ -1789,15 +1783,13 @@ 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 - T03 = ( - 2 * total_mass_dot * (Vector([0, 0, r_NOZ]) - r_CM) - - 2 * total_mass * r_CM_dot - ) + r_NOZ_vector = Vector([0, 0, r_NOZ]) + T03 = 2 * total_mass_dot * (r_NOZ_vector - r_CM) - 2 * total_mass * r_CM_dot T04 = ( self.rocket.motor.thrust.get_value_opt(t) * Vector([0, 0, 1]) - total_mass * r_CM_ddot - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (Vector([0, 0, r_NOZ]) - r_CM) + + total_mass_ddot * (r_NOZ_vector - r_CM) ) T05 = total_mass_dot * S_nozzle - I_dot @@ -2075,19 +2067,19 @@ def alpha3(self): def R1(self): """Aerodynamic force along the first axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[0] + return np.column_stack([self.time, self.retrieve_temporary_values_arrays[0]]) @funcify_method("Time (s)", "R2 (N)", "spline", "zero") def R2(self): """Aerodynamic force along the second axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[1] + return np.column_stack([self.time, self.retrieve_temporary_values_arrays[1]]) @funcify_method("Time (s)", "R3 (N)", "spline", "zero") def R3(self): """Aerodynamic force along the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[2] + return np.column_stack([self.time, self.retrieve_temporary_values_arrays[2]]) @funcify_method("Time (s)", "M1 (Nm)", "spline", "zero") def M1(self): @@ -2095,20 +2087,20 @@ def M1(self): perpendicular to the rocket's axis of symmetry as a Function of time. """ - return self.retrieve_temporary_values_arrays[3] + return np.column_stack([self.time, self.retrieve_temporary_values_arrays[3]]) @funcify_method("Time (s)", "M2 (Nm)", "spline", "zero") def M2(self): """Aerodynamic bending moment in the same direction as the axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[4] + return np.column_stack([self.time, self.retrieve_temporary_values_arrays[4]]) @funcify_method("Time (s)", "M3 (Nm)", "spline", "zero") def M3(self): """Aerodynamic bending moment in the same direction as the rocket's axis of symmetry as a Function of time.""" - return self.retrieve_temporary_values_arrays[5] + return np.column_stack([self.time, self.retrieve_temporary_values_arrays[5]]) @funcify_method("Time (s)", "Pressure (Pa)", "spline", "constant") def pressure(self): @@ -2260,8 +2252,7 @@ def attitude_angle(self): attitude_angle = (180 / np.pi) * np.arctan2( self.attitude_vector_z[:, 1], horizontal_attitude_proj[:, 1] ) - attitude_angle = np.column_stack([self.time, attitude_angle]) - return attitude_angle + return np.column_stack([self.time, attitude_angle]) # Lateral Attitude Angle @funcify_method("Time (s)", "Lateral Attitude Angle (°)") @@ -2323,24 +2314,17 @@ def theta(self): @funcify_method("Time (s)", "Freestream Velocity X (m/s)", "spline", "constant") def stream_velocity_x(self): """Freestream velocity X component as a Function of time.""" - stream_velocity_x = np.column_stack( - (self.time, self.wind_velocity_x[:, 1] - self.vx[:, 1]) - ) - return stream_velocity_x + return np.column_stack((self.time, self.wind_velocity_x[:, 1] - self.vx[:, 1])) @funcify_method("Time (s)", "Freestream Velocity Y (m/s)", "spline", "constant") def stream_velocity_y(self): """Freestream velocity Y component as a Function of time.""" - stream_velocity_y = np.column_stack( - (self.time, self.wind_velocity_y[:, 1] - self.vy[:, 1]) - ) - return stream_velocity_y + return np.column_stack((self.time, self.wind_velocity_y[:, 1] - self.vy[:, 1])) @funcify_method("Time (s)", "Freestream Velocity Z (m/s)", "spline", "constant") def stream_velocity_z(self): """Freestream velocity Z component as a Function of time.""" - stream_velocity_z = np.column_stack((self.time, -self.vz[:, 1])) - return stream_velocity_z + return np.column_stack((self.time, -self.vz[:, 1])) @funcify_method("Time (s)", "Freestream Speed (m/s)", "spline", "constant") def free_stream_speed(self): @@ -2511,7 +2495,7 @@ def translational_energy(self): # Redefine total_mass time grid to allow for efficient Function algebra total_mass = deepcopy(self.rocket.total_mass) total_mass.set_discrete_based_on_model(self.vz) - translational_energy = 0.5 * total_mass * (self.vx**2 + self.vy**2 + self.vz**2) + translational_energy = 0.5 * total_mass * (self.speed**2) return translational_energy @funcify_method("Time (s)", "Kinetic Energy (J)", "spline", "zero") @@ -2585,9 +2569,8 @@ def angle_of_attack(self): # Calculate angle of attack and convert to degrees angle_of_attack = np.rad2deg(np.arccos(dot_product_normalized)) - angle_of_attack = np.column_stack([self.time, angle_of_attack]) - return angle_of_attack + return np.column_stack([self.time, angle_of_attack]) # Frequency response and stability variables @funcify_method("Frequency (Hz)", "ω1 Fourier Amplitude", "spline", "zero") @@ -3039,9 +3022,8 @@ def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities """ v_f = self.out_of_rail_velocity - # Convert angle to radians - theta = self.inclination * 3.14159265359 / 180 - stall_angle = stall_angle * 3.14159265359 / 180 + theta = np.radians(self.inclination) + stall_angle = np.radians(stall_angle) c = (math.cos(stall_angle) ** 2 - math.cos(theta) ** 2) / math.sin( stall_angle @@ -3055,16 +3037,13 @@ def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities ** 0.5 ) / 2 - # Convert stall_angle to degrees - stall_angle = stall_angle * 180 / np.pi + stall_angle = np.degrees(stall_angle) print( "Maximum wind velocity at Rail Departure time before angle" + f" of attack exceeds {stall_angle:.3f}°: {w_v:.3f} m/s" ) - return None - - def export_pressures(self, file_name, time_step): + def export_pressures(self, file_name, time_step): # TODO: move out """Exports the pressure experienced by the rocket during the flight to an external file, the '.csv' format is recommended, as the columns will be separated by commas. It can handle flights with or without @@ -3129,6 +3108,7 @@ class attributes which are instances of the Function class. Usage will be exported. Otherwise, linear interpolation is carried out to calculate values at the desired time steps. Example: 0.001. """ + # TODO: we should move this method to outside of class. # Fast evaluation for the most basic scenario if time_step is None and len(variables) == 0: @@ -3254,33 +3234,33 @@ def export_kml( # TODO: should be moved out of this class. # Open kml file with simplekml library kml = simplekml.Kml(open=1) trajectory = kml.newlinestring(name="Rocket Trajectory - Powered by RocketPy") - coords = [] + if altitude_mode == "relativetoground": # In this mode the elevation data will be the Above Ground Level # elevation. Only works properly if the ground level is similar to # a plane, i.e. it might not work well if the terrain has mountains - for t in time_points: - coords.append( - ( - self.longitude.get_value_opt(t), - self.latitude.get_value_opt(t), - self.altitude.get_value_opt(t), - ) + coords = [ + ( + self.longitude.get_value_opt(t), + self.latitude.get_value_opt(t), + self.altitude.get_value_opt(t), ) + for t in time_points + ] trajectory.coords = coords trajectory.altitudemode = simplekml.AltitudeMode.relativetoground else: # altitude_mode == 'absolute' # In this case the elevation data will be the Above Sea Level elevation # Ensure you use the correct value on self.env.elevation, otherwise # the trajectory path can be offset from ground - for t in time_points: - coords.append( - ( - self.longitude.get_value_opt(t), - self.latitude.get_value_opt(t), - self.z.get_value_opt(t), - ) + coords = [ + ( + self.longitude.get_value_opt(t), + self.latitude.get_value_opt(t), + self.z.get_value_opt(t), ) + for t in time_points + ] trajectory.coords = coords trajectory.altitudemode = simplekml.AltitudeMode.absolute # Modify style of trajectory linestring @@ -3536,14 +3516,14 @@ def add_node(self, t, parachutes, callbacks): self.list.append(self.TimeNode(t, parachutes, callbacks)) def add_parachutes(self, parachutes, t_init, t_end): - # Iterate over parachutes for parachute in parachutes: # Calculate start of sampling time nodes - pcDt = 1 / parachute.sampling_rate + sampling_interval = 1 / parachute.sampling_rate parachute_node_list = [ - self.TimeNode(i * pcDt, [parachute], []) + self.TimeNode(i * sampling_interval, [parachute], []) for i in range( - math.ceil(t_init / pcDt), math.floor(t_end / pcDt) + 1 + math.ceil(t_init / sampling_interval), + math.floor(t_end / sampling_interval) + 1, ) ] self.list += parachute_node_list From ac01264285d2fbb8f792f5a01172eb1154bcd5cf Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 21:06:38 -0400 Subject: [PATCH 56/71] DEV: Adds PR 581 to CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 70b466134..bf02f98e2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -38,6 +38,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ### Changed +- ENH: Flight simulation speed up [#581] (https://github.com/RocketPy-Team/RocketPy/pull/581) - DOC: Improvements of Environment docstring phrasing [#565](https://github.com/RocketPy-Team/RocketPy/pull/565) - MNT: Refactor flight prints module [#579](https://github.com/RocketPy-Team/RocketPy/pull/579) - DOC: Convert CompareFlights example notebooks to .rst files [#576](https://github.com/RocketPy-Team/RocketPy/pull/576) From 5f2a4dc23a39ac42e5418a5c89ec9d5648aa8ac7 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 18 Apr 2024 21:15:40 -0400 Subject: [PATCH 57/71] BUG: fix linear interpolation root finding method name in Function class --- 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 8f96d49c5..3207b3142 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -3010,7 +3010,7 @@ def find_root_linear_interpolation(x0, x1, y0, y1, y): -------- >>> from rocketpy import Function >>> x0, x1, y0, y1, y = 0, 1, 0, 1, 0.5 - >>> x = Function.linear_interpolation_root_finding(x0, x1, y0, y1, y) + >>> x = Function.find_root_linear_interpolation(x0, x1, y0, y1, y) >>> x 0.5 """ From 8d0986490fff76726ce09aa64715427db225a0ce Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Fri, 26 Apr 2024 15:30:50 -0400 Subject: [PATCH 58/71] MNT: move auxiliary functions from Function to tools.py --- rocketpy/mathutils/function.py | 146 ------------------------------- rocketpy/simulation/flight.py | 15 ++-- rocketpy/tools.py | 152 +++++++++++++++++++++++++++++++++ tests/unit/test_function.py | 40 --------- tests/unit/test_tools.py | 46 ++++++++++ 5 files changed, 206 insertions(+), 193 deletions(-) create mode 100644 tests/unit/test_tools.py diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index c8688aae5..6dc1c764b 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2868,152 +2868,6 @@ def savetxt( file.write(header_line + newline) np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline) - # Define auxiliary method - - @staticmethod - def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): - """Calculate the coefficients of a cubic Hermite interpolation function. - The function is defined as ax**3 + bx**2 + cx + d. - - Parameters - ---------- - x0 : float - Position of the first point. - x1 : float - Position of the second point. - y0 : float - Value of the function evaluated at the first point. - yp0 : float - Value of the derivative of the function evaluated at the first - point. - y1 : float - Value of the function evaluated at the second point. - yp1 : float - Value of the derivative of the function evaluated at the second - point. - - Returns - ------- - tuple[float, float, float, float] - The coefficients of the cubic Hermite interpolation function. - """ - dx = x1 - x0 - d = float(y0) - c = float(yp0) - b = float((3 * y1 - yp1 * dx - 2 * c * dx - 3 * d) / (dx**2)) - a = float(-(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3)) - return a, b, c, d - - @staticmethod - def find_roots_cubic_function(a, b, c, d): - """Calculate the roots of a cubic function using Cardano's method. - - This method applies Cardano's method to find the roots of a cubic - function of the form ax^3 + bx^2 + cx + d. The roots may be complex - numbers. - - Parameters - ---------- - a : float - Coefficient of the cubic term (x^3). - b : float - Coefficient of the quadratic term (x^2). - c : float - Coefficient of the linear term (x). - d : float - Constant term. - - Returns - ------- - tuple[complex, complex, complex] - A tuple containing the real and complex roots of the cubic function. - Note that the roots may be complex numbers. The roots are ordered - in the tuple as x1, x2, x3. - - References - ---------- - - Cardano's method: https://en.wikipedia.org/wiki/Cubic_function#Cardano's_method - - Examples - -------- - >>> from rocketpy import Function - - First we define the coefficients of the function ax**3 + bx**2 + cx + d - >>> a, b, c, d = 1, -3, -1, 3 - >>> x1, x2, x3 = Function.find_roots_cubic_function(a, b, c, d) - >>> x1, x2, x3 - ((-1+0j), (3+7.401486830834377e-17j), (1-1.4802973661668753e-16j)) - - To get the real part of the roots, use the real attribute of the complex - number. - >>> x1.real, x2.real, x3.real - (-1.0, 3.0, 1.0) - """ - delta_0 = b**2 - 3 * a * c - delta_1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 - c1 = ((delta_1 + (delta_1**2 - 4 * delta_0**3) ** (0.5)) / 2) ** (1 / 3) - - c2_0 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 0 - x1 = -(1 / (3 * a)) * (b + c2_0 + delta_0 / c2_0) - - c2_1 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 1 - x2 = -(1 / (3 * a)) * (b + c2_1 + delta_0 / c2_1) - - c2_2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 2 - x3 = -(1 / (3 * a)) * (b + c2_2 + delta_0 / c2_2) - - # Alternative method (same results) - - # Q = (3 * a * c - b**2) / (9 * a**2) - # R = (9 * a * b * c - 27 * a**2 * d - 2 * b**3) / (54 * a**3) - # S = (R + (Q**3 + R**2)**0.5)**(1 / 3) - # T = (R - (Q**3 + R**2)**0.5)**(1 / 3) - - # x1 = S + T - b / (3 * a) - # x2 = -(S + T) / 2 - b / (3 * a) + 1j * (S - T) * (3**0.5) / 2 - # x3 = -(S + T) / 2 - b / (3 * a) - 1j * (S - T) * (3**0.5) / 2 - - return x1, x2, x3 - - @staticmethod - def find_root_linear_interpolation(x0, x1, y0, y1, y): - """Calculate the root of a linear interpolation function. - - This method calculates the root of a linear interpolation function - given two points (x0, y0) and (x1, y1) and a value y. The function - is defined as y = m*x + c. - - Parameters - ---------- - x0 : float - Position of the first point. - x1 : float - Position of the second point. - y0 : float - Value of the function evaluated at the first point. - y1 : float - Value of the function evaluated at the second point. - y : float - Value of the function at the desired point. - - Returns - ------- - float - The root of the linear interpolation function. This represents the - value of x at which the function evaluates to y. - - Examples - -------- - >>> from rocketpy import Function - >>> x0, x1, y0, y1, y = 0, 1, 0, 1, 0.5 - >>> x = Function.find_root_linear_interpolation(x0, x1, y0, y1, y) - >>> x - 0.5 - """ - m = (y1 - y0) / (x1 - x0) - c = y0 - m * x0 - return (y - c) / m - # Input validators def __validate_source(self, source): """Used to validate the source parameter for creating a Function object. diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 43ab16ef5..96deb165b 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -12,7 +12,10 @@ from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints from ..tools import ( + calculate_cubic_hermite_coefficients, find_closest, + find_root_linear_interpolation, + find_roots_cubic_function, quaternions_to_nutation, quaternions_to_precession, quaternions_to_spin, @@ -812,7 +815,7 @@ def __simulate(self, verbose): self.solution[-2][3] += self.env.elevation self.solution[-1][3] += self.env.elevation # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) - a, b, c, d = Function.calculate_cubic_hermite_coefficients( + a, b, c, d = calculate_cubic_hermite_coefficients( 0, float(phase.solver.step_size), y0, @@ -822,7 +825,7 @@ def __simulate(self, verbose): ) a += 1e-5 # TODO: why?? # Find roots - t_roots = Function.find_roots_cubic_function(a, b, c, d) + t_roots = find_roots_cubic_function(a, b, c, d) # Find correct root valid_t_root = [ t_root.real @@ -862,9 +865,7 @@ def __simulate(self, verbose): # Assume linear vz(t) to detect when vz = 0 t0, vz0 = self.solution[-2][0], self.solution[-2][6] t1, vz1 = self.solution[-1][0], self.solution[-1][6] - t_root = Function.find_root_linear_interpolation( - t0, t1, vz0, vz1, 0 - ) + t_root = find_root_linear_interpolation(t0, t1, vz0, vz1, 0) # Fetch state at t_root interpolator = phase.solver.dense_output() self.apogee_state = interpolator(t_root) @@ -887,7 +888,7 @@ def __simulate(self, verbose): # print('\n>>>PASSIVE EVENT DETECTED: Touchdown!') # Check exactly when it happened using root finding # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) - a, b, c, d = Function.calculate_cubic_hermite_coefficients( + a, b, c, d = calculate_cubic_hermite_coefficients( x0=0, # t0 x1=float(phase.solver.step_size), # t1 - t0 y0=float(self.solution[-2][3] - self.env.elevation), # z0 @@ -896,7 +897,7 @@ def __simulate(self, verbose): yp1=float(self.solution[-1][6]), # vz1 ) # Find roots - t_roots = Function.find_roots_cubic_function(a, b, c, d) + t_roots = find_roots_cubic_function(a, b, c, d) # Find correct root t1 = self.solution[-1][0] - self.solution[-2][0] valid_t_root = [ diff --git a/rocketpy/tools.py b/rocketpy/tools.py index ccecffd4f..716f3a0ab 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -1,3 +1,11 @@ +"""The module rocketpy.tools contains a set of functions that are used +throughout the rocketpy package. These functions are not specific to any +particular class or module, and are used to perform general tasks that are +required by multiple classes or modules. These functions can be modified or +expanded to suit the needs of other modules and may present breaking changes +between minor versions if necessary, although this will be always avoided. +""" + import functools import importlib import importlib.metadata @@ -42,6 +50,150 @@ def tuple_handler(value): raise ValueError("value must be a list or tuple of length 1 or 2.") +def calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1): + """Calculate the coefficients of a cubic Hermite interpolation function. + The function is defined as ax**3 + bx**2 + cx + d. + + Parameters + ---------- + x0 : float + Position of the first point. + x1 : float + Position of the second point. + y0 : float + Value of the function evaluated at the first point. + yp0 : float + Value of the derivative of the function evaluated at the first + point. + y1 : float + Value of the function evaluated at the second point. + yp1 : float + Value of the derivative of the function evaluated at the second + point. + + Returns + ------- + tuple[float, float, float, float] + The coefficients of the cubic Hermite interpolation function. + """ + dx = x1 - x0 + d = float(y0) + c = float(yp0) + b = float((3 * y1 - yp1 * dx - 2 * c * dx - 3 * d) / (dx**2)) + a = float(-(2 * y1 - yp1 * dx - c * dx - 2 * d) / (dx**3)) + return a, b, c, d + + +def find_roots_cubic_function(a, b, c, d): + """Calculate the roots of a cubic function using Cardano's method. + + This method applies Cardano's method to find the roots of a cubic + function of the form ax^3 + bx^2 + cx + d. The roots may be complex + numbers. + + Parameters + ---------- + a : float + Coefficient of the cubic term (x^3). + b : float + Coefficient of the quadratic term (x^2). + c : float + Coefficient of the linear term (x). + d : float + Constant term. + + Returns + ------- + tuple[complex, complex, complex] + A tuple containing the real and complex roots of the cubic function. + Note that the roots may be complex numbers. The roots are ordered + in the tuple as x1, x2, x3. + + References + ---------- + - Cardano's method: https://en.wikipedia.org/wiki/Cubic_function#Cardano's_method + + Examples + -------- + >>> from rocketpy import Function + + First we define the coefficients of the function ax**3 + bx**2 + cx + d + >>> a, b, c, d = 1, -3, -1, 3 + >>> x1, x2, x3 = Function.find_roots_cubic_function(a, b, c, d) + >>> x1, x2, x3 + ((-1+0j), (3+7.401486830834377e-17j), (1-1.4802973661668753e-16j)) + + To get the real part of the roots, use the real attribute of the complex + number. + >>> x1.real, x2.real, x3.real + (-1.0, 3.0, 1.0) + """ + delta_0 = b**2 - 3 * a * c + delta_1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 + c1 = ((delta_1 + (delta_1**2 - 4 * delta_0**3) ** (0.5)) / 2) ** (1 / 3) + + c2_0 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 0 + x1 = -(1 / (3 * a)) * (b + c2_0 + delta_0 / c2_0) + + c2_1 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 1 + x2 = -(1 / (3 * a)) * (b + c2_1 + delta_0 / c2_1) + + c2_2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 2 + x3 = -(1 / (3 * a)) * (b + c2_2 + delta_0 / c2_2) + + # Alternative method (same results) + + # Q = (3 * a * c - b**2) / (9 * a**2) + # R = (9 * a * b * c - 27 * a**2 * d - 2 * b**3) / (54 * a**3) + # S = (R + (Q**3 + R**2)**0.5)**(1 / 3) + # T = (R - (Q**3 + R**2)**0.5)**(1 / 3) + + # x1 = S + T - b / (3 * a) + # x2 = -(S + T) / 2 - b / (3 * a) + 1j * (S - T) * (3**0.5) / 2 + # x3 = -(S + T) / 2 - b / (3 * a) - 1j * (S - T) * (3**0.5) / 2 + + return x1, x2, x3 + + +def find_root_linear_interpolation(x0, x1, y0, y1, y): + """Calculate the root of a linear interpolation function. + + This method calculates the root of a linear interpolation function + given two points (x0, y0) and (x1, y1) and a value y. The function + is defined as y = m*x + c. + + Parameters + ---------- + x0 : float + Position of the first point. + x1 : float + Position of the second point. + y0 : float + Value of the function evaluated at the first point. + y1 : float + Value of the function evaluated at the second point. + y : float + Value of the function at the desired point. + + Returns + ------- + float + The root of the linear interpolation function. This represents the + value of x at which the function evaluates to y. + + Examples + -------- + >>> from rocketpy import Function + >>> x0, x1, y0, y1, y = 0, 1, 0, 1, 0.5 + >>> x = Function.find_root_linear_interpolation(x0, x1, y0, y1, y) + >>> x + 0.5 + """ + m = (y1 - y0) / (x1 - x0) + c = y0 - m * x0 + return (y - c) / m + + def bilinear_interpolation(x, y, x1, x2, y1, y2, z11, z12, z21, z22): """Bilinear interpolation. It considers the values of the four points around the point to be interpolated and returns the interpolated value. diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index 35a5b21e9..17da59498 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -334,43 +334,3 @@ def test_get_domain_dim(linear_func): def test_bool(linear_func): """Test the __bool__ method of the Function class.""" assert bool(linear_func) == True - - -def test_calculate_cubic_hermite_coefficients(): - """Test the calculate_cubic_hermite_coefficients method of the Function class.""" - # Function: f(x) = x**3 + 2x**2 -1 ; derivative: f'(x) = 3x**2 + 4x - x = np.array([-3, -2, -1, 0, 1]) - y = np.array([-10, -1, 0, -1, 2]) - - # Selects two points as x0 and x1 - x0, x1 = 0, 1 - y0, y1 = -1, 2 - yp0, yp1 = 0, 7 - - a, b, c, d = Function.calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1) - - assert np.isclose(a, 1) - assert np.isclose(b, 2) - assert np.isclose(c, 0) - assert np.isclose(d, -1) - assert np.allclose( - a * x**3 + b * x**2 + c * x + d, - y, - ) - - -def test_cardanos_root_finding(): - """Tests the find_roots_cubic_function method of the Function class.""" - # Function: f(x) = x**3 + 2x**2 -1 - # roots: (-1 - 5**0.5) / 2; -1; (-1 + 5**0.5) / 2 - - roots = list(Function.find_roots_cubic_function(a=1, b=2, c=0, d=-1)) - roots.sort(key=lambda x: x.real) - - assert np.isclose(roots[0].real, (-1 - 5**0.5) / 2) - assert np.isclose(roots[1].real, -1) - assert np.isclose(roots[2].real, (-1 + 5**0.5) / 2) - - assert np.isclose(roots[0].imag, 0) - assert np.isclose(roots[1].imag, 0) - assert np.isclose(roots[2].imag, 0) diff --git a/tests/unit/test_tools.py b/tests/unit/test_tools.py new file mode 100644 index 000000000..75a526aac --- /dev/null +++ b/tests/unit/test_tools.py @@ -0,0 +1,46 @@ +import numpy as np + +from rocketpy.tools import ( + calculate_cubic_hermite_coefficients, + find_roots_cubic_function, +) + + +def test_calculate_cubic_hermite_coefficients(): + """Test the calculate_cubic_hermite_coefficients method of the Function class.""" + # Function: f(x) = x**3 + 2x**2 -1 ; derivative: f'(x) = 3x**2 + 4x + x = np.array([-3, -2, -1, 0, 1]) + y = np.array([-10, -1, 0, -1, 2]) + + # Selects two points as x0 and x1 + x0, x1 = 0, 1 + y0, y1 = -1, 2 + yp0, yp1 = 0, 7 + + a, b, c, d = calculate_cubic_hermite_coefficients(x0, x1, y0, yp0, y1, yp1) + + assert np.isclose(a, 1) + assert np.isclose(b, 2) + assert np.isclose(c, 0) + assert np.isclose(d, -1) + assert np.allclose( + a * x**3 + b * x**2 + c * x + d, + y, + ) + + +def test_cardanos_root_finding(): + """Tests the find_roots_cubic_function method of the Function class.""" + # Function: f(x) = x**3 + 2x**2 -1 + # roots: (-1 - 5**0.5) / 2; -1; (-1 + 5**0.5) / 2 + + roots = list(find_roots_cubic_function(a=1, b=2, c=0, d=-1)) + roots.sort(key=lambda x: x.real) + + assert np.isclose(roots[0].real, (-1 - 5**0.5) / 2) + assert np.isclose(roots[1].real, -1) + assert np.isclose(roots[2].real, (-1 + 5**0.5) / 2) + + assert np.isclose(roots[0].imag, 0) + assert np.isclose(roots[1].imag, 0) + assert np.isclose(roots[2].imag, 0) From 505275a134f06e5a471b331aa021bcda74024a02 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Fri, 26 Apr 2024 15:31:24 -0400 Subject: [PATCH 59/71] MNT: Fix apogee position calculation in Flight class --- rocketpy/simulation/flight.py | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 96deb165b..5341d07b2 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -871,8 +871,11 @@ def __simulate(self, verbose): self.apogee_state = interpolator(t_root) # Store apogee data self.apogee_time = t_root + self.apogee_x = self.apogee_state[0] + self.apogee_y = self.apogee_state[1] + self.apogee = self.apogee_state[2] + if self.terminate_on_apogee: - # print(">>> Terminate on apogee activated!") self.t = self.t_final = t_root # Roll back solution self.solution[-1] = [self.t, *self.apogee_state] @@ -1962,21 +1965,6 @@ def z(self): """Rocket z position as a Function of time.""" return self.solution_array[:, [0, 3]] - @property - def apogee_x(self): - """Rocket x position at apogee.""" - return self.x.get_value_opt(self.apogee_time) - - @property - def apogee_y(self): - """Rocket y position at apogee.""" - return self.y.get_value_opt(self.apogee_time) - - @property - def apogee(self): - """Rocket z position at apogee.""" - return self.z.get_value_opt(self.apogee_time) - @funcify_method("Time (s)", "Altitude AGL (m)", "spline", "constant") def altitude(self): """Rocket altitude above ground level as a Function of time.""" From 2690fb293290addd647744b5d01a941029446a8a Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Fri, 26 Apr 2024 15:31:38 -0400 Subject: [PATCH 60/71] MNT: delete unused comments in the Flight class --- rocketpy/simulation/flight.py | 45 ----------------------------------- 1 file changed, 45 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 5341d07b2..1d7329bee 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -636,10 +636,6 @@ def __repr__(self): def __simulate(self, verbose): """Simulate the flight trajectory.""" for phase_index, phase in self.time_iterator(self.flight_phases): - # print("\nCurrent Flight Phase List:") - # print(self.flight_phases) - # print("\n\tCurrent Flight Phase:") - # print("\t\tIndex: ", phase_index, " | Phase: ", phase) # Determine maximum time for this flight phase phase.time_bound = self.flight_phases[phase_index + 1].t @@ -659,13 +655,6 @@ def __simulate(self, verbose): rtol=self.rtol, atol=self.atol, ) - # print(f"\n\tSolver Initialization Details:") - # print(f"\t\tInitial Time (t): {phase.t}") - # print(f"\t\tInitial State (y_sol): {self.y_sol}") - # print(f"\t\tTime Bound: {phase.time_bound}") - # print(f"\t\tMin Step: {self.min_time_step}") - # print(f"\t\tMax Step: {self.max_time_step}") - # print(f"\t\tTolerances: rtol = {self.rtol}, atol = {self.atol}") # Initialize phase time nodes phase.time_nodes = self.TimeNodes() @@ -692,8 +681,6 @@ def __simulate(self, verbose): # Iterate through time nodes for node_index, node in self.time_iterator(phase.time_nodes): - # print("\n\t\tCurrent Time Node") - # print("\t\tIndex: ", node_index, " | Time Node: ", node) # Determine time bound for this time node node.time_bound = phase.time_nodes[node_index + 1].t # NOTE: Setting the time bound and status for the phase solver, @@ -723,8 +710,6 @@ def __simulate(self, verbose): if parachute.triggerfunc( noisy_pressure, height_above_ground_level, self.y_sol ): - # print("\nEVENT DETECTED: Parachute Triggered") - # print("Name: ", parachute.name, " | Lag: ", parachute.lag) # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create flight phase for time after detection and before inflation @@ -770,12 +755,6 @@ def __simulate(self, verbose): self.y_sol = phase.solver.y if verbose: print(f"Current Simulation Time: {self.t:3.4f} s", end="\r") - # print("\n\t\t\tCurrent Step Details:") - # print( - # "\t\t\tIState: ", phase.solver._lsoda_solver._integrator.istate - # ) - # print("\t\t\tTime: ", phase.solver.t) - # print("\t\t\tAltitude (ASL): ", phase.solver.y[2]) # Check for first out of rail event if len(self.out_of_rail_state) == 1 and ( @@ -784,7 +763,6 @@ def __simulate(self, verbose): + (self.y_sol[2] - self.env.elevation) ** 2 >= self.effective_1rl**2 ): - # print("\n>>> EVENT DETECTED: Rocket is Out of Rail!") # Check exactly when it went out using root finding # Disconsider elevation self.solution[-2][3] -= self.env.elevation @@ -861,7 +839,6 @@ def __simulate(self, verbose): # Check for apogee event if len(self.apogee_state) == 1 and self.y_sol[5] < 0: - # print("\n>>> EVENT DETECTED: Rocket has reached apogee") # Assume linear vz(t) to detect when vz = 0 t0, vz0 = self.solution[-2][0], self.solution[-2][6] t1, vz1 = self.solution[-1][0], self.solution[-1][6] @@ -888,7 +865,6 @@ def __simulate(self, verbose): phase.solver.status = "finished" # Check for impact event if self.y_sol[2] < self.env.elevation: - # print('\n>>>PASSIVE EVENT DETECTED: Touchdown!') # Check exactly when it happened using root finding # Cubic Hermite interpolation (ax**3 + bx**2 + cx + d) a, b, c, d = calculate_cubic_hermite_coefficients( @@ -949,27 +925,12 @@ def __simulate(self, verbose): if overshootable_nodes[0].t == phase.t and phase.clear: overshootable_nodes[0].parachutes = [] overshootable_nodes[0].callbacks = [] - # print("\n\t\t\t\tOvershootable Time Nodes") - # print("\t\t\t\tInterval: ", self.solution[-2][0], self.t) - # print( - # "\t\t\t\tOvershootable Nodes Length: ", - # str(len(overshootable_nodes)), - # " | Overshootable Nodes: ", - # overshootable_nodes, - # ) # Feed overshootable time nodes trigger interpolator = phase.solver.dense_output() for ( overshootable_index, overshootable_node, ) in self.time_iterator(overshootable_nodes): - # print("\n\t\t\t\tCurrent Overshootable Node:") - # print( - # "\t\t\t\tIndex: ", - # overshootable_index, - # " | Overshootable Node: ", - # overshootable_node, - # ) # Calculate state at node time overshootable_node.y_sol = interpolator( overshootable_node.t @@ -990,12 +951,6 @@ def __simulate(self, verbose): height_above_ground_level, overshootable_node.y_sol, ): - # print( - # "\n>>> EVENT DETECTED: Parachute Triggered!" - # ) - # print( - # f"\tParachute Name: {parachute.name} | Lag: {parachute.lag}" - # ) # Remove parachute from flight parachutes self.parachutes.remove(parachute) # Create flight phase for time after detection and before inflation From 0df1a4685e35e05500a597c7e6a19ad8593366dd Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 29 Apr 2024 14:35:27 -0400 Subject: [PATCH 61/71] BUG: fix Flight.__repr__ method --- 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 1d7329bee..27eb7fbef 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -625,8 +625,8 @@ def __init__( def __repr__(self): return ( - f" Date: Mon, 29 Apr 2024 14:39:09 -0400 Subject: [PATCH 62/71] MNT: simplify post process evaluation with __evaluate_post_process method - This introduces small breaking changes since some public variables were deleted, however they are not crucial for the class to work. --- rocketpy/simulation/flight.py | 182 +++++++++------------------------- 1 file changed, 48 insertions(+), 134 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 27eb7fbef..29d579e41 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1526,7 +1526,9 @@ def u_dot(self, t, u, post_processing=False): ] if post_processing: - u_dot.extend([R1, R2, R3, M1, M2, M3]) + self.__post_processed_variables.append( + [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3] + ) return u_dot @@ -1777,7 +1779,9 @@ def u_dot_generalized(self, t, u, post_processing=False): u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] if post_processing: - u_dot.extend([R1, R2, R3, M1, M2, M3]) + self.__post_processed_variables.append( + [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3] + ) return u_dot @@ -1840,7 +1844,9 @@ def u_dot_parachute(self, t, u, post_processing=False): az = (Dz - 9.8 * mp) / (mp + ma) if post_processing: - return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0] + return self.__post_processed_variables.append( + [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0] + ) return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] @@ -1979,51 +1985,51 @@ def w3(self): @funcify_method("Time (s)", "Ax (m/s²)", "spline", "zero") def ax(self): """Rocket x acceleration as a Function of time.""" - return self.retrieve_acceleration_arrays[0] + return self.__evaluate_post_process[:, [0, 1]] @funcify_method("Time (s)", "Ay (m/s²)", "spline", "zero") def ay(self): """Rocket y acceleration as a Function of time.""" - return self.retrieve_acceleration_arrays[1] + return self.__evaluate_post_process[:, [0, 2]] @funcify_method("Time (s)", "Az (m/s²)", "spline", "zero") def az(self): """Rocket z acceleration as a Function of time.""" - return self.retrieve_acceleration_arrays[2] + return self.__evaluate_post_process[:, [0, 3]] @funcify_method("Time (s)", "α1 (rad/s²)", "spline", "zero") def alpha1(self): """Rocket angular acceleration α1 as a Function of time.""" - return self.retrieve_acceleration_arrays[3] + return self.__evaluate_post_process[:, [0, 4]] @funcify_method("Time (s)", "α2 (rad/s²)", "spline", "zero") def alpha2(self): """Rocket angular acceleration α2 as a Function of time.""" - return self.retrieve_acceleration_arrays[4] + return self.__evaluate_post_process[:, [0, 5]] @funcify_method("Time (s)", "α3 (rad/s²)", "spline", "zero") def alpha3(self): """Rocket angular acceleration α3 as a Function of time.""" - return self.retrieve_acceleration_arrays[5] + return self.__evaluate_post_process[:, [0, 6]] # Process third type of outputs - Temporary values @funcify_method("Time (s)", "R1 (N)", "spline", "zero") def R1(self): """Aerodynamic force along the first axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return np.column_stack([self.time, self.retrieve_temporary_values_arrays[0]]) + return self.__evaluate_post_process[:, [0, 7]] @funcify_method("Time (s)", "R2 (N)", "spline", "zero") def R2(self): """Aerodynamic force along the second axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return np.column_stack([self.time, self.retrieve_temporary_values_arrays[1]]) + return self.__evaluate_post_process[:, [0, 8]] @funcify_method("Time (s)", "R3 (N)", "spline", "zero") def R3(self): """Aerodynamic force along the rocket's axis of symmetry as a Function of time.""" - return np.column_stack([self.time, self.retrieve_temporary_values_arrays[2]]) + return self.__evaluate_post_process[:, [0, 9]] @funcify_method("Time (s)", "M1 (Nm)", "spline", "zero") def M1(self): @@ -2031,20 +2037,20 @@ def M1(self): perpendicular to the rocket's axis of symmetry as a Function of time. """ - return np.column_stack([self.time, self.retrieve_temporary_values_arrays[3]]) + return self.__evaluate_post_process[:, [0, 10]] @funcify_method("Time (s)", "M2 (Nm)", "spline", "zero") def M2(self): """Aerodynamic bending moment in the same direction as the axis that is perpendicular to the rocket's axis of symmetry as a Function of time.""" - return np.column_stack([self.time, self.retrieve_temporary_values_arrays[4]]) + return self.__evaluate_post_process[:, [0, 11]] @funcify_method("Time (s)", "M3 (Nm)", "spline", "zero") def M3(self): """Aerodynamic bending moment in the same direction as the rocket's axis of symmetry as a Function of time.""" - return np.column_stack([self.time, self.retrieve_temporary_values_arrays[5]]) + return self.__evaluate_post_process[:, [0, 12]] @funcify_method("Time (s)", "Pressure (Pa)", "spline", "constant") def pressure(self): @@ -2685,125 +2691,6 @@ def longitude(self): return np.column_stack((self.time, longitude)) - @cached_property - def retrieve_acceleration_arrays(self): - """Retrieve acceleration arrays from the integration scheme - - Parameters - ---------- - - Returns - ------- - ax: list - acceleration in x direction - ay: list - acceleration in y direction - az: list - acceleration in z direction - alpha1: list - angular acceleration in x direction - alpha2: list - angular acceleration in y direction - 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.flight_phases): - init_time = phase.t - final_time = self.flight_phases[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 - - @cached_property - def retrieve_temporary_values_arrays(self): - """Retrieve temporary values arrays from the integration scheme. - Currently, the following temporary values are retrieved: ``R1`` , ``R2`` - ``R3`` , ``M1`` , ``M2`` , ``M3`` , ``pressure`` , ``density`` , - ``dynamic_viscosity`` , ``speed_of_sound`` . - - Returns - ------- - self.R1_list: list - R1 values. - self.R2_list: list - R2 values. - self.R3_list: list - R3 values are the aerodynamic force values in the rocket's axis - direction. - self.M1_list: list - M1 values. - self.M2_list: list - Aerodynamic bending moment in e2 direction at each time step. - self.M3_list: list - Aerodynamic bending moment in e3 direction at each time step. - self.pressure_list: list - Air pressure at each time step. - self.density_list: list - Air density at each time step. - self.dynamic_viscosity_list: list - Dynamic viscosity at each time step. - self.speed_of_sound_list: list - Speed of sound at each time step. - self.wind_velocity_x_list: list - Wind velocity in x direction at each time step. - self.wind_velocity_y_list: list - Wind velocity in y direction at each time step. - """ - - # Initialize force and atmospheric arrays - temporary_values = [ - [], # R1_list - [], # R2_list - [], # R3_list - [], # M1_list - [], # M2_list - [], # M3_list - ] - - # Go through each time step and calculate forces and atmospheric values - # Get flight phases - for phase_index, phase in self.time_iterator(self.flight_phases): - init_time = phase.t - final_time = self.flight_phases[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) - for i, value in enumerate(temporary_values): - value.append(u_dot[i + 13]) - - return temporary_values - def get_controller_observed_variables(self): """Retrieve the observed variables related to air brakes from the controllers. If there is only one set of observed variables, it is @@ -2931,6 +2818,33 @@ def __transform_pressure_signals_lists_to_functions(self): + parachute.noise_signal_function ) + @cached_property + def __evaluate_post_process(self): + """Evaluate all post-processing variables by running the simulation + again but with the post-processing flag set to True. + + Returns + ------- + np.array + An array containing all post-processed variables evaluated at each + time step. Each element of the array is a list containing: + [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3] + """ + self.__post_processed_variables = [] + for phase_index, phase in self.time_iterator(self.flight_phases): + init_time = phase.t + final_time = self.flight_phases[phase_index + 1].t + current_derivative = phase.derivative + for callback in phase.callbacks: + callback(self) + for step in self.solution: + if init_time < step[0] <= final_time or ( + init_time == self.t_initial and step[0] == self.t_initial + ): + current_derivative(step[0], step[1:], post_processing=True) + + return np.array(self.__post_processed_variables) + def post_process(self, interpolation="spline", extrapolation="natural"): """This method is **deprecated** and is only kept here for backwards compatibility. All attributes that need to be post processed are From 79be71439d1f43f9d6b2a5e20a9a8eef701807d9 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 29 Apr 2024 14:51:19 -0400 Subject: [PATCH 63/71] TST: update test values for accelerations calculations --- tests/test_flight.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_flight.py b/tests/test_flight.py index 4fb4036eb..2379489d5 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -619,8 +619,8 @@ def test_max_values(flight_calisto_robust): """ test = flight_calisto_robust atol = 1e-2 - assert pytest.approx(105.2774, abs=atol) == test.max_acceleration_power_on - assert pytest.approx(105.2774, abs=atol) == test.max_acceleration + assert pytest.approx(105.1599, abs=atol) == test.max_acceleration_power_on + assert pytest.approx(105.1599, abs=atol) == test.max_acceleration assert pytest.approx(0.85999, abs=atol) == test.max_mach_number assert pytest.approx(285.94948, abs=atol) == test.max_speed @@ -649,8 +649,8 @@ def test_rail_buttons_forces(flight_calisto_custom_wind): @pytest.mark.parametrize( "flight_time, expected_values", [ - ("t_initial", (0, 0, 0)), - ("out_of_rail_time", (0, 7.8068, 89.2325)), + ("t_initial", (0.07757, 0.07106, -9.33378)), + ("out_of_rail_time", (0.24249, 8.75150, 89.13007)), ("apogee_time", (0.07534, -0.058127, -9.614386)), ("t_final", (0, 0, 0.0017346294117130806)), ], From 4b3d884e287485f50f2ed86e82e78ce6136aa8ff Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 29 Apr 2024 14:56:31 -0400 Subject: [PATCH 64/71] MNT: Refactor find_roots_cubic_function method in tools.py --- rocketpy/tools.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 716f3a0ab..e3696af84 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -132,7 +132,7 @@ def find_roots_cubic_function(a, b, c, d): delta_1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2 c1 = ((delta_1 + (delta_1**2 - 4 * delta_0**3) ** (0.5)) / 2) ** (1 / 3) - c2_0 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 0 + c2_0 = c1 x1 = -(1 / (3 * a)) * (b + c2_0 + delta_0 / c2_0) c2_1 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 1 @@ -141,17 +141,6 @@ def find_roots_cubic_function(a, b, c, d): c2_2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** 2 x3 = -(1 / (3 * a)) * (b + c2_2 + delta_0 / c2_2) - # Alternative method (same results) - - # Q = (3 * a * c - b**2) / (9 * a**2) - # R = (9 * a * b * c - 27 * a**2 * d - 2 * b**3) / (54 * a**3) - # S = (R + (Q**3 + R**2)**0.5)**(1 / 3) - # T = (R - (Q**3 + R**2)**0.5)**(1 / 3) - - # x1 = S + T - b / (3 * a) - # x2 = -(S + T) / 2 - b / (3 * a) + 1j * (S - T) * (3**0.5) / 2 - # x3 = -(S + T) / 2 - b / (3 * a) - 1j * (S - T) * (3**0.5) / 2 - return x1, x2, x3 From a1737c89bddc8a8fd23570fe97089f2f8fbc0490 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 29 Apr 2024 16:37:00 -0400 Subject: [PATCH 65/71] ENH: defaults to u_dot in case the rocket has a solid propulsion --- rocketpy/simulation/flight.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 29d579e41..7ad091e84 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -9,6 +9,7 @@ from ..mathutils.function import Function, funcify_method from ..mathutils.vector_matrix import Matrix, Vector +from ..motors.solid_motor import SolidMotor from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints from ..tools import ( @@ -1116,7 +1117,11 @@ def __init_solver_monitors(self): def __init_equations_of_motion(self): """Initialize equations of motion.""" - if self.equations_of_motion == "solid_propulsion": + if ( + isinstance(self.rocket.motor, SolidMotor) + or 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 @cached_property From 327d9b3e016cc59867dccc37f3c2859fae1c47a7 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 1 May 2024 22:37:54 -0400 Subject: [PATCH 66/71] REV: reverts commit a1737c8 --- rocketpy/simulation/flight.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 7ad091e84..e5079445b 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1117,10 +1117,7 @@ def __init_solver_monitors(self): def __init_equations_of_motion(self): """Initialize equations of motion.""" - if ( - isinstance(self.rocket.motor, SolidMotor) - or self.equations_of_motion == "solid_propulsion" - ): + 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 From 0eb4ae8b9c0779d232b126de58f8394771ff3435 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 2 May 2024 01:28:12 -0400 Subject: [PATCH 67/71] TST: test_max_values now receives relative tolerance --- tests/test_flight.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_flight.py b/tests/test_flight.py index 2379489d5..c0b995555 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -618,11 +618,11 @@ def test_max_values(flight_calisto_robust): regarding this pytest fixture. """ test = flight_calisto_robust - atol = 1e-2 - assert pytest.approx(105.1599, abs=atol) == test.max_acceleration_power_on - assert pytest.approx(105.1599, abs=atol) == test.max_acceleration - assert pytest.approx(0.85999, abs=atol) == test.max_mach_number - assert pytest.approx(285.94948, abs=atol) == test.max_speed + rtol = 1e-2 + assert pytest.approx(105.1599, rel=rtol) == test.max_acceleration_power_on + assert pytest.approx(105.1599, rel=rtol) == test.max_acceleration + assert pytest.approx(0.85999, rel=rtol) == test.max_mach_number + assert pytest.approx(285.94948, rel=rtol) == test.max_speed def test_rail_buttons_forces(flight_calisto_custom_wind): From 14eca0979b49357d677a0c05f0869c1a3216523e Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 2 May 2024 01:33:29 -0400 Subject: [PATCH 68/71] TST: fix the documentation tests of tools.py --- rocketpy/tools.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rocketpy/tools.py b/rocketpy/tools.py index e3696af84..86ad7f17e 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -115,11 +115,11 @@ def find_roots_cubic_function(a, b, c, d): Examples -------- - >>> from rocketpy import Function + >>> from rocketpy.tools import find_roots_cubic_function First we define the coefficients of the function ax**3 + bx**2 + cx + d >>> a, b, c, d = 1, -3, -1, 3 - >>> x1, x2, x3 = Function.find_roots_cubic_function(a, b, c, d) + >>> x1, x2, x3 = find_roots_cubic_function(a, b, c, d) >>> x1, x2, x3 ((-1+0j), (3+7.401486830834377e-17j), (1-1.4802973661668753e-16j)) @@ -172,9 +172,9 @@ def find_root_linear_interpolation(x0, x1, y0, y1, y): Examples -------- - >>> from rocketpy import Function + >>> from rocketpy.tools import find_root_linear_interpolation >>> x0, x1, y0, y1, y = 0, 1, 0, 1, 0.5 - >>> x = Function.find_root_linear_interpolation(x0, x1, y0, y1, y) + >>> x = find_root_linear_interpolation(x0, x1, y0, y1, y) >>> x 0.5 """ From dd9bc159383ce104bcd341d906ad1e870ff6b8b0 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> Date: Thu, 2 May 2024 19:59:24 -0400 Subject: [PATCH 69/71] BUG: fix append function in rocketpy/simulation/flight.py Co-authored-by: MateusStano <69485049+MateusStano@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 e5079445b..c5c74c7f1 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1846,7 +1846,7 @@ def u_dot_parachute(self, t, u, post_processing=False): az = (Dz - 9.8 * mp) / (mp + ma) if post_processing: - return self.__post_processed_variables.append( + self.__post_processed_variables.append( [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0] ) From 92e93a7f616f9c7a02b281fba09418e96f387b61 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Thu, 2 May 2024 20:20:46 -0400 Subject: [PATCH 70/71] ENH: modify the post_processing version of udot_rail1 --- rocketpy/simulation/flight.py | 15 +++++++++------ tests/test_flight.py | 6 +++--- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index c5c74c7f1..b762c376f 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -9,7 +9,6 @@ from ..mathutils.function import Function, funcify_method from ..mathutils.vector_matrix import Matrix, Vector -from ..motors.solid_motor import SolidMotor from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints from ..tools import ( @@ -1205,11 +1204,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 @@ -1242,6 +1236,15 @@ def udot_rail1(self, t, u, post_processing=False): else: ax, ay, az = 0, 0, 0 + if post_processing: + self.u_dot_generalized(t, u, post_processing=True) + 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 + return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0] def udot_rail2(self, t, u, post_processing=False): diff --git a/tests/test_flight.py b/tests/test_flight.py index c0b995555..8b9d338f9 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -618,7 +618,7 @@ def test_max_values(flight_calisto_robust): regarding this pytest fixture. """ test = flight_calisto_robust - rtol = 1e-2 + rtol = 5e-3 assert pytest.approx(105.1599, rel=rtol) == test.max_acceleration_power_on assert pytest.approx(105.1599, rel=rtol) == test.max_acceleration assert pytest.approx(0.85999, rel=rtol) == test.max_mach_number @@ -649,8 +649,8 @@ def test_rail_buttons_forces(flight_calisto_custom_wind): @pytest.mark.parametrize( "flight_time, expected_values", [ - ("t_initial", (0.07757, 0.07106, -9.33378)), - ("out_of_rail_time", (0.24249, 8.75150, 89.13007)), + ("t_initial", (0, 0, 0)), + ("out_of_rail_time", (0, 7.8068, 89.2325)), ("apogee_time", (0.07534, -0.058127, -9.614386)), ("t_final", (0, 0, 0.0017346294117130806)), ], From 80d85bb0c9d1d1136fb3b74b14b8f5e3402c3829 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Fri, 3 May 2024 00:21:16 +0000 Subject: [PATCH 71/71] Fix code style issues with Black --- tests/test_flight.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_flight.py b/tests/test_flight.py index 8b9d338f9..d4ed1e114 100644 --- a/tests/test_flight.py +++ b/tests/test_flight.py @@ -650,7 +650,7 @@ def test_rail_buttons_forces(flight_calisto_custom_wind): "flight_time, expected_values", [ ("t_initial", (0, 0, 0)), - ("out_of_rail_time", (0, 7.8068, 89.2325)), + ("out_of_rail_time", (0, 7.8068, 89.2325)), ("apogee_time", (0.07534, -0.058127, -9.614386)), ("t_final", (0, 0, 0.0017346294117130806)), ],