diff --git a/rocketpy/motors/ClusterMotor b/rocketpy/motors/ClusterMotor new file mode 100644 index 000000000..0c1624524 --- /dev/null +++ b/rocketpy/motors/ClusterMotor @@ -0,0 +1,405 @@ +import numpy as np +from rocketpy.mathutils.function import Function +from rocketpy.mathutils.vector_matrix import Vector +import matplotlib.pyplot as plt +import matplotlib.patches as patches +from rocketpy.tools import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I22, + parallel_axis_theorem_I33, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I23, +) + + +class ClusterMotor: + """ + Manages a cluster of motors. + + This class behaves like a single motor by aggregating the properties + of several individual motors and implementing the full interface + expected by the Rocket class. It handles the calculation of + aggregated thrust, mass, center of mass, and inertia tensor as + functions of time, considering the 3D position and orientation + of each motor. + """ + + def __init__(self, motors, positions, orientations=None): + """Initializes the ClusterMotor, aggregating multiple motors. + + Parameters + ---------- + motors : list[Motor] + A list of instantiated RocketPy Motor objects to be part of the cluster. + positions : list[tuple] or list[list] or list[np.array] + A list of 3D position vectors [x, y, z] for each motor. + The position is relative to the rocket's coordinate system origin, + which is also the cluster's origin. The coordinate system + orientation is assumed to be 'tail_to_nose'. + orientations : list[tuple] or list[list] or list[np.array], optional + A list of 3D unit vectors [x, y, z] specifying the thrust + direction for each motor in the rocket's coordinate system. + If None, all motors are assumed to thrust along the rocket's + positive Z-axis (e.g., [0, 0, 1]). Defaults to None. + """ + self._validate_inputs(motors, positions, orientations) + + self.coordinate_system_orientation = "tail_to_nose" + self._csys = 1 + + self._initialize_basic_properties() + self._initialize_thrust_and_mass() + self._initialize_center_of_mass() + self._initialize_inertia_properties() + + def _validate_inputs(self, motors, positions, orientations): + """Validates and stores the primary inputs for the cluster.""" + if not motors: + raise ValueError("The list of motors cannot be empty.") + + self.motors = motors + self.positions = [np.array(pos) for pos in positions] + + if orientations is None: + self.orientations = [np.array([0, 0, 1.0]) for _ in motors] + else: + self.orientations = [ + np.array(ori) / np.linalg.norm(ori) for ori in orientations + ] + + if not len(self.motors) == len(self.positions) == len(self.orientations): + raise ValueError( + "The 'motors', 'positions', and 'orientations' lists must have the same length." + ) + + def _initialize_basic_properties(self): + """Calculates simple aggregated scalar properties.""" + self.propellant_initial_mass = sum( + m.propellant_initial_mass for m in self.motors + ) + self.dry_mass = sum(m.dry_mass for m in self.motors) + self.total_impulse = sum(m.total_impulse for m in self.motors) + + self.burn_start_time = min(motor.burn_start_time for motor in self.motors) + self.burn_out_time = max(motor.burn_out_time for motor in self.motors) + self.burn_time = (self.burn_start_time, self.burn_out_time) + + self.nozzle_radius = np.sqrt(sum(m.nozzle_radius**2 for m in self.motors)) + self.throat_radius = np.sqrt(sum(m.throat_radius**2 for m in self.motors)) + + self.nozzle_position = np.mean([ + pos[2] + motor.nozzle_position * motor._csys + for motor, pos in zip(self.motors, self.positions) + ]) + + def _initialize_thrust_and_mass(self): + """Initializes thrust and mass-related Function objects.""" + self.thrust = Function( + self._calculate_total_thrust_scalar, inputs="t", outputs="Thrust (N)" + ) + if self.burn_time[1] > self.burn_time[0]: + self.thrust.set_discrete(self.burn_start_time, self.burn_out_time, 100) + + if self.propellant_initial_mass > 1e-9: + average_exhaust_velocity = self.total_impulse / self.propellant_initial_mass + self.total_mass_flow_rate = self.thrust / -average_exhaust_velocity + else: + self.total_mass_flow_rate = Function(0) + + self.mass_flow_rate = self.total_mass_flow_rate + + self.propellant_mass = Function( + self._calculate_propellant_mass, inputs="t", outputs="Propellant Mass (kg)" + ) + self.total_mass = Function( + self._calculate_total_mass, inputs="t", outputs="Total Mass (kg)" + ) + + def _initialize_center_of_mass(self): + """Initializes center of mass Function objects.""" + self.center_of_mass = Function( + self._calculate_center_of_mass, inputs="t", outputs="Cluster CoM Vector (m)" + ) + self.center_of_propellant_mass = Function( + self._calculate_center_of_propellant_mass, + inputs="t", + outputs="Cluster CoPM Vector (m)", + ) + + com_time = self.burn_out_time if self.burn_out_time > self.burn_start_time else self.burn_start_time + self.center_of_dry_mass_position = self._calculate_center_of_mass( + com_time + ) + + def _initialize_inertia_properties(self): + """Initializes dry and propellant inertia properties.""" + ( + self.dry_I_11, + self.dry_I_22, + self.dry_I_33, + self.dry_I_12, + self.dry_I_13, + self.dry_I_23, + ) = self._calculate_total_dry_inertia() + + def propellant_inertia_func(t): + return self._calculate_total_propellant_inertia(t) + + self.propellant_I_11 = Function( + lambda t: propellant_inertia_func(t)[0], inputs="t", outputs="I_11 (kg*m^2)" + ) + self.propellant_I_22 = Function( + lambda t: propellant_inertia_func(t)[1], inputs="t", outputs="I_22 (kg*m^2)" + ) + self.propellant_I_33 = Function( + lambda t: propellant_inertia_func(t)[2], inputs="t", outputs="I_33 (kg*m^2)" + ) + self.propellant_I_12 = Function( + lambda t: propellant_inertia_func(t)[3], inputs="t", outputs="I_12 (kg*m^2)" + ) + self.propellant_I_13 = Function( + lambda t: propellant_inertia_func(t)[4], inputs="t", outputs="I_13 (kg*m^2)" + ) + self.propellant_I_23 = Function( + lambda t: propellant_inertia_func(t)[5], inputs="t", outputs="I_23 (kg*m^2)" + ) + + def _calculate_total_mass(self, t): + """Calculates total cluster mass at time t.""" + return self.dry_mass + self._calculate_propellant_mass(t) + + def _calculate_propellant_mass(self, t): + """ + Calculates the total propellant mass at time t by integrating the + total mass flow rate. This ensures consistency between mass and thrust. + """ + return self.propellant_initial_mass + self.total_mass_flow_rate.integral( + self.burn_start_time, t + ) + + def get_total_thrust_vector(self, t): + """ + Calculates the total thrust vector of the cluster at a given time t. + This vector is the sum of all individual motor thrust vectors, + considering their orientation. + """ + total_thrust = np.zeros(3) + for motor, orientation in zip(self.motors, self.orientations): + scalar_thrust = motor.thrust.get_value_opt(t) + total_thrust += scalar_thrust * orientation + return Vector(total_thrust) + + def _calculate_total_thrust_scalar(self, t): + """ + Calculates the magnitude of the total thrust vector. + This is what is wrapped by the `self.thrust` Function. + """ + return abs(self.get_total_thrust_vector(t)) + + def get_total_moment(self, t, ref_point): + """ + Calculates the total moment (torque) generated by the cluster + about a given reference point (e.g., the rocket's CoM). + """ + total_moment = np.zeros(3, dtype=np.float64) + ref_point_arr = np.array(ref_point, dtype=np.float64) + + for motor, pos, orientation in zip( + self.motors, self.positions, self.orientations + ): + force_magnitude = motor.thrust.get_value_opt(t) + force = force_magnitude * orientation + arm = pos - ref_point_arr + total_moment += np.cross(arm, force) + + if hasattr(motor, "thrust_eccentricity_y") and hasattr( + motor, "thrust_eccentricity_x" + ): + total_moment[0] += motor.thrust_eccentricity_y * force_magnitude + total_moment[1] -= motor.thrust_eccentricity_x * force_magnitude + + return Vector(total_moment) + + def pressure_thrust(self, pressure): + """Calculates the total pressure thrust correction for the cluster.""" + return sum(motor.pressure_thrust(pressure) for motor in self.motors) + + def _calculate_center_of_mass(self, t): + """Calculates the aggregated center of mass of the cluster at time t.""" + total_mass_val = self._calculate_total_mass(t) + if total_mass_val < 1e-9: + return Vector( + np.mean(self.positions, axis=0) if self.positions else np.zeros(3) + ) + + weighted_sum = np.zeros(3, dtype=np.float64) + for motor, pos in zip(self.motors, self.positions): + motor_com_local_z = motor.center_of_mass.get_value_opt(t) + motor_com_global = pos + np.array([0, 0, motor_com_local_z * motor._csys]) + weighted_sum += motor.total_mass.get_value_opt(t) * motor_com_global + + return Vector(weighted_sum / total_mass_val) + + def _calculate_center_of_propellant_mass(self, t): + """ + Calculates the aggregated center of mass of the cluster's propellant. + This calculation is based on the individual motor properties. + """ + total_prop_mass = 0.0 + weighted_sum = np.zeros(3, dtype=np.float64) + + for motor in self.motors: + total_prop_mass += motor.propellant_mass.get_value_opt(t) + + if total_prop_mass < 1e-9: + return self.center_of_dry_mass_position + + for motor, pos in zip(self.motors, self.positions): + prop_mass_t = motor.propellant_mass.get_value_opt(t) + if prop_mass_t > 1e-9: + prop_com_local_z = motor.center_of_propellant_mass.get_value_opt(t) + prop_com_global = pos + np.array([0, 0, prop_com_local_z * motor._csys]) + weighted_sum += prop_mass_t * prop_com_global + + return Vector(weighted_sum / total_prop_mass) + + def _calculate_total_dry_inertia(self): + """ + Calculates the cluster's total dry inertia tensor relative to the + cluster's center of dry mass. + """ + I_11, I_22, I_33 = 0.0, 0.0, 0.0 + I_12, I_13, I_23 = 0.0, 0.0, 0.0 + + ref_point = self.center_of_dry_mass_position + + for motor, pos in zip(self.motors, self.positions): + m = motor.dry_mass + motor_cdm_local_z = motor.center_of_dry_mass_position + motor_cdm_global = pos + np.array([0, 0, motor_cdm_local_z * motor._csys]) + r_vec = Vector(motor_cdm_global) - ref_point + + I_11 += parallel_axis_theorem_I11(motor.dry_I_11, m, r_vec) + I_22 += parallel_axis_theorem_I22(motor.dry_I_22, m, r_vec) + I_33 += parallel_axis_theorem_I33(motor.dry_I_33, m, r_vec) + I_12 += parallel_axis_theorem_I12(motor.dry_I_12, m, r_vec) + I_13 += parallel_axis_theorem_I13(motor.dry_I_13, m, r_vec) + I_23 += parallel_axis_theorem_I23(motor.dry_I_23, m, r_vec) + + return I_11, I_22, I_33, I_12, I_13, I_23 + + def _calculate_total_propellant_inertia(self, t): + """ + Calculates the cluster's total propellant inertia tensor relative to the + cluster's center of propellant mass at time t. + """ + I_11, I_22, I_33 = 0.0, 0.0, 0.0 + I_12, I_13, I_23 = 0.0, 0.0, 0.0 + + ref_point = self._calculate_center_of_propellant_mass(t) + + for motor, pos in zip(self.motors, self.positions): + m = motor.propellant_mass.get_value_opt(t) + if m < 1e-9: + continue + + prop_com_local_z = motor.center_of_propellant_mass.get_value_opt(t) + prop_com_global = pos + np.array([0, 0, prop_com_local_z * motor._csys]) + r_vec = Vector(prop_com_global) - ref_point + + I_11 += parallel_axis_theorem_I11( + motor.propellant_I_11.get_value_opt(t), m, r_vec + ) + I_22 += parallel_axis_theorem_I22( + motor.propellant_I_22.get_value_opt(t), m, r_vec + ) + I_33 += parallel_axis_theorem_I33( + motor.propellant_I_33.get_value_opt(t), m, r_vec + ) + I_12 += parallel_axis_theorem_I12( + motor.propellant_I_12.get_value_opt(t), m, r_vec + ) + I_13 += parallel_axis_theorem_I13( + motor.propellant_I_13.get_value_opt(t), m, r_vec + ) + I_23 += parallel_axis_theorem_I23( + motor.propellant_I_23.get_value_opt(t), m, r_vec + ) + + return I_11, I_22, I_33, I_12, I_13, I_23 + + def draw_rear_view(self, rocket_radius, tail_radius=None, filename=None): + """ + Plots a 2D rear view of the motor cluster, showing the main + rocket body and optionally the tail cone diameter. + + Args: + rocket_radius (float): The main radius of the rocket body. + tail_radius (float, optional): The rocket's radius at the tail (boattail). + If provided, a second circle will be drawn. + filename (str, optional): If provided, saves the plot to a file. + Otherwise, shows the plot. + """ + _, ax = plt.subplots(figsize=(8.5, 8.5)) + + rocket_body = patches.Circle( + (0, 0), + radius=rocket_radius, + facecolor="lightgrey", + edgecolor="black", + linewidth=2, + label=f"Main Body ({rocket_radius * 200:.0f} cm)", + ) + ax.add_patch(rocket_body) + + if tail_radius is not None: + tail_body = patches.Circle( + (0, 0), + radius=tail_radius, + facecolor="silver", + edgecolor="black", + linestyle="--", + linewidth=1.5, + label=f"Shrunk ({tail_radius * 200:.0f} cm)", + ) + ax.add_patch(tail_body) + + for i, (motor, pos) in enumerate(zip(self.motors, self.positions)): + motor_body = patches.Circle( + (pos[0], pos[1]), + radius=getattr(motor, "grain_outer_radius", motor.nozzle_radius/2), + facecolor="dimgrey", + edgecolor="black", + linewidth=1.5, + ) + ax.add_patch(motor_body) + ax.text( + pos[0], + pos[1], + f"M{i + 1}", + ha="center", + va="center", + color="white", + fontsize=12, + fontweight="bold", + ) + + ax.set_aspect("equal", adjustable="box") + plot_limit = rocket_radius * 1.2 + ax.set_xlim(-plot_limit, plot_limit) + ax.set_ylim(-plot_limit, plot_limit) + ax.set_title("Detailed Rear View of the Cluster", fontsize=16) + ax.set_xlabel("Axis X (m)") + ax.set_ylabel("Axis Y (m)") + ax.grid(True, linestyle="--", alpha=0.6) + ax.axhline(0, color="black", linewidth=0.5) + ax.axvline(0, color="black", linewidth=0.5) + + plt.legend() + + if filename: + plt.savefig(filename) + print(f"Rear view plot saved to: {filename}") + else: + plt.show() diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py index 7cb28670c..a46112347 100644 --- a/rocketpy/motors/hybrid_motor.py +++ b/rocketpy/motors/hybrid_motor.py @@ -1,13 +1,17 @@ from functools import cached_property - -from rocketpy.tools import parallel_axis_theorem_from_com - from ..mathutils.function import Function, funcify_method, reset_funcified_methods from ..plots.hybrid_motor_plots import _HybridMotorPlots -from ..prints.hybrid_motor_prints import _HybridMotorPrints -from .liquid_motor import LiquidMotor from .motor import Motor from .solid_motor import SolidMotor +from ..tools import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I22, + parallel_axis_theorem_I33, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I23, +) +from ..mathutils.vector_matrix import Vector class HybridMotor(Motor): @@ -195,192 +199,421 @@ class HybridMotor(Motor): It will allow to obtain the net thrust in the Flight class. """ - def __init__( # pylint: disable=too-many-arguments + # pylint: disable=too-many-locals, too-many-statements + def __init__( self, thrust_source, dry_mass, dry_inertia, nozzle_radius, + burn_time, + center_of_dry_mass_position, grain_number, + grain_separation, grain_density, grain_outer_radius, grain_initial_inner_radius, grain_initial_height, - grain_separation, grains_center_of_mass_position, - center_of_dry_mass_position, + tanks_mass, + oxidizer_initial_mass, + oxidizer_mass_flow_rate_curve, + oxidizer_density, + oxidizer_tanks_geometries, + oxidizer_tanks_positions, + oxidizer_initial_inertia=(0, 0, 0), nozzle_position=0, - burn_time=None, - throat_radius=0.01, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, ): - """Initialize Motor class, process thrust curve and geometrical + """Initializes HybridMotor class, process thrust curve and geometrical parameters and store results. Parameters ---------- - thrust_source : int, float, callable, string, array, Function - Motor's thrust curve. Can be given as an int or float, in which - case the thrust will be considered constant in time. It can - also be given as a callable function, whose argument is time in - seconds and returns the thrust supplied by the motor in the - instant. If a string is given, it must point to a .csv or .eng file. - The .csv file can contain a single line header and the first column - must specify time in seconds, while the second column specifies - thrust. Arrays may also be specified, following rules set by the - class Function. Thrust units are Newtons. - - .. seealso:: :doc:`Thrust Source Details ` + thrust_source : int, float, callable, string, numpy.ndarray, list + Motor's thrust curve. Can be given as Thrust-Time pairs array or + file path (csv or eng format). Is passed to the Function class, see + help(Function) for more information. Thrust units are Newtons, time + units are seconds. dry_mass : int, float - Same as in Motor class. See the :class:`Motor ` docs + Motor's dry mass in kg. This is the mass of the motor without + propellant. dry_inertia : tuple, list Tuple or list containing the motor's dry mass inertia tensor components, in kg*m^2. This inertia is defined with respect to the - the `center_of_dry_mass_position` position. - Assuming e_3 is the rocket's axis of symmetry, e_1 and e_2 are - orthogonal and form a plane perpendicular to e_3, the dry mass - inertia tensor components must be given in the following order: - (I_11, I_22, I_33, I_12, I_13, I_23), where I_ij is the - component of the inertia tensor in the direction of e_i x e_j. - Alternatively, the inertia tensor can be given as - (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. + motor's center of dry mass position. Assuming e_3 is the motor's axis + of symmetry, e_1 and e_2 are orthogonal and form a plane + perpendicular to e_3, the dry mass inertia tensor components must be + given in the following order: (I_11, I_22, I_33, I_12, I_13, I_23). + Alternatively, the inertia tensor can be given as (I_11, I_22, I_33), + where I_12 = I_13 = I_23 = 0. nozzle_radius : int, float - Motor's nozzle outlet radius in meters. + Motor's nozzle radius in meters. + burn_time: int, float, tuple of int, float + Motor's burn time. + If a tuple is passed, the first value is the ignition time and the + second value is the end of burn time. If a single number is passed, + the ignition time is assumed to be 0 and the end of burn time is + the number passed. + center_of_dry_mass_position : int, float + Position of the motor's center of dry mass (i.e. center of mass + without propellant) relative to the motor's coordinate system + origin, in meters. See the ``coordinate_system_orientation`` + parameter for details on the coordinate system. grain_number : int - Number of solid grains + Number of solid grains. + grain_separation : int, float + Distance between grains, in meters. grain_density : int, float - Solid grain density in kg/m3. + Solid grain density in kg/m³. grain_outer_radius : int, float Solid grain outer radius in meters. grain_initial_inner_radius : int, float Solid grain initial inner radius in meters. grain_initial_height : int, float Solid grain initial height in meters. - grain_separation : int, float - Distance between grains, in meters. - grains_center_of_mass_position : float - Position of the center of mass of the grains in meters. More - specifically, the coordinate of the center of mass specified in the - motor's coordinate system. - See :doc:`Positions and Coordinate Systems ` for - more information. - center_of_dry_mass_position : int, float - The position, in meters, of the motor's center of mass with respect - to the motor's coordinate system when it is devoid of propellant. - See :doc:`Positions and Coordinate Systems `. + grains_center_of_mass_position : int, float + Position of the center of mass of the grains relative to the motor's + coordinate system origin in meters. Generally equal to + ``center_of_dry_mass_position``. + grain_initial_inertia : tuple, list, optional + Tuple or list containing the initial inertia tensor components of a + single grain, in kg*m^2. This inertia is defined with respect to the + the grains_center_of_mass_position position. If not specified, the + grain is assumed to be a hollow cylinder with the initial dimensions. + Assuming e_3 is the grain's axis of symmetry, e_1 and e_2 are + orthogonal and form a plane perpendicular to e_3, the initial inertia + tensor components must be given in the following order: + (I_11, I_22, I_33, I_12, I_13, I_23). Alternatively, the inertia + tensor can be given as (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. + Default is (0, 0, 0). + tanks_mass : float + Total mass of the oxidizer tanks structures in kg. Includes the mass + of the tanks themselves, valves, pipes, etc. It is assumed constant + over time. + oxidizer_initial_mass : float + Initial mass of the oxidizer, including liquid and gas phases, in kg. + oxidizer_mass_flow_rate_curve : int, float, callable, string, numpy.ndarray, list + Oxidizer mass flow rate curve. Can be given as MassFlowRate-Time + pairs array or file path (csv format). It is used to calculate the + oxidizer mass and center of mass position as a function of time. + If int or float is given, it is assumed constant. Mass flow rate + units are kg/s, time units are seconds. Passed to the Function + class, see help(Function) for more information. + oxidizer_density : float + Density of the oxidizer in kg/m³. It is used to calculate the volume + and height of the oxidizer in the tanks. It is assumed constant over + time. + oxidizer_tanks_geometries : list + List of tuples, where each tuple represents the geometry of an + oxidizer tank. Accepted geometries are: + ('cylinder', (top_radius, bottom_radius, height)) + ('sphere', radius) + ('ullage', volume) + Dimensions should be in meters and volume in cubic meters. + The list must contain at least one tank geometry. Ullage tanks can only be + placed at the top or bottom of the tanks stack. + Example: [('ullage', 0.01), ('cylinder', (0.1, 0.1, 0.5)), ('cylinder', (0.1, 0.05, 0.2))] + oxidizer_tanks_positions : list + List of floats, representing the position of the centroid of each + oxidizer tank's geometry with respect to the motor's coordinate system + origin, in meters. The list must have the same length as + ``oxidizer_tanks_geometries``. + See the ``coordinate_system_orientation`` parameter for details on the coordinate system. + oxidizer_tanks_initial_liquid_level : float, optional + Initial liquid level in the tanks, measured in meters from the bottom + of the tanks stack. If specified, this parameter overrides the initial + oxidizer mass calculation based on ``oxidizer_initial_mass``, allowing + precise control over the starting volume of the liquid oxidizer. If + None, the initial liquid level is derived from ``oxidizer_initial_mass``. + Default is None. + oxidizer_tanks_initial_ullage_mass : float, optional + Initial mass of the ullage gas in kg. If not specified, it is assumed + to be 0. Default is 0. + oxidizer_tanks_initial_ullage_volume : float, optional + Initial volume of the ullage gas in cubic meters. If not specified, it + is automatically calculated based on the tanks geometries and the + initial liquid level. Default is None. + oxidizer_initial_inertia : tuple, list, optional + Tuple or list containing the initial inertia tensor components of the + oxidizer (liquid + gas), in kg*m^2. This inertia is defined with + respect to the initial oxidizer center of mass position. If not + specified, the oxidizer is assumed to be a point mass. Default is (0, 0, 0). + Assuming e_3 is the motor's axis of symmetry, e_1 and e_2 are + orthogonal and form a plane perpendicular to e_3, the initial inertia + tensor components must be given in the following order: + (I_11, I_22, I_33, I_12, I_13, I_23). Alternatively, the inertia + tensor can be given as (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. nozzle_position : int, float, optional - Motor's nozzle outlet position in meters, in the motor's coordinate - system. See :doc:`Positions and Coordinate Systems ` - for details. Default is 0, in which case the origin of the - coordinate system is placed at the motor's nozzle outlet. - burn_time: float, tuple of float, optional - Motor's burn time. - If a float is given, the burn time is assumed to be between 0 and - the given float, in seconds. - If a tuple of float is given, the burn time is assumed to be between - the first and second elements of the tuple, in seconds. - If not specified, automatically sourced as the range between the - first and last-time step of the motor's thrust curve. This can only - be used if the motor's thrust is defined by a list of points, such - as a .csv file, a .eng file or a Function instance whose source is - a list. - throat_radius : int, float, optional - Motor's nozzle throat radius in meters. Used to calculate Kn curve. - Optional if the Kn curve is not interesting. Its value does not - impact trajectory simulation. + Motor's nozzle outlet position in meters, specified in the motor's + coordinate system. Default is 0, which corresponds to the motor's + origin. See the ``coordinate_system_orientation`` parameter for + details on the coordinate system. reshape_thrust_curve : boolean, tuple, optional - If False, the original thrust curve supplied is not altered. If a - tuple is given, whose first parameter is a new burn out time and - whose second parameter is a new total impulse in Ns, the thrust - curve is reshaped to match the new specifications. May be useful - for motors whose thrust curve shape is expected to remain similar - in case the impulse and burn time varies slightly. Default is - False. + If False, the original thrust curve supplied is used. If a tuple is + given, the thrust curve is reshaped to match the new grain mass + flow rate and burn time. The tuple should contain the initial grain + mass and the final grain mass, in kg. Default is False. interpolation_method : string, optional Method of interpolation to be used in case thrust curve is given - by data set in .csv or .eng, or as an array. Options are 'spline' - 'akima' and 'linear'. Default is "linear". + by data points. Options are 'spline', 'akima' and 'linear'. + Default is "linear". coordinate_system_orientation : string, optional Orientation of the motor's coordinate system. The coordinate system - is defined by the motor's axis of symmetry. The origin of the - coordinate system may be placed anywhere along such axis, such as - at the nozzle area, and must be kept the same for all other - positions specified. Options are "nozzle_to_combustion_chamber" and - "combustion_chamber_to_nozzle". Default is - "nozzle_to_combustion_chamber". + has its origin at the motor's center of mass and is oriented + according to the following options: + "nozzle_to_combustion_chamber" : the coordinate system is oriented + with the z-axis pointing from the nozzle towards the combustion + chamber. + "combustion_chamber_to_nozzle" : the coordinate system is oriented + with the z-axis pointing from the combustion chamber towards the + nozzle. + Default is "nozzle_to_combustion_chamber". reference_pressure : int, float, optional - Atmospheric pressure in Pa at which the thrust data was recorded. + Reference pressure in Pa used to calculate vacuum thrust and + pressure thrust. This corresponds to the atmospheric pressure + measured during the static test of the motor. Default is None, + which means no pressure correction is applied. Returns ------- None """ - super().__init__( + # Call SolidMotor init to initialize grain parameters + # Note: dry_mass and dry_inertia are temporarily set to 0, they will be + # calculated later considering the tanks mass. + SolidMotor.__init__( + self, thrust_source=thrust_source, - dry_inertia=dry_inertia, + dry_mass=0, + dry_inertia=(0, 0, 0), nozzle_radius=nozzle_radius, + burn_time=burn_time, center_of_dry_mass_position=center_of_dry_mass_position, - dry_mass=dry_mass, + grain_number=grain_number, + grain_separation=grain_separation, + grain_density=grain_density, + grain_outer_radius=grain_outer_radius, + grain_initial_inner_radius=grain_initial_inner_radius, + grain_initial_height=grain_initial_height, + grains_center_of_mass_position=grains_center_of_mass_position, nozzle_position=nozzle_position, - burn_time=burn_time, reshape_thrust_curve=reshape_thrust_curve, interpolation_method=interpolation_method, coordinate_system_orientation=coordinate_system_orientation, reference_pressure=reference_pressure, ) - self.liquid = LiquidMotor( - thrust_source, - dry_mass, - dry_inertia, - nozzle_radius, - center_of_dry_mass_position, - nozzle_position, - burn_time, - reshape_thrust_curve, - interpolation_method, - coordinate_system_orientation, - reference_pressure, + + # Oxidizer parameters initialization + self.tanks_mass = tanks_mass + self.oxidizer_initial_mass = oxidizer_initial_mass + self.oxidizer_density = oxidizer_density + self.oxidizer_tanks_geometries = oxidizer_tanks_geometries + self.oxidizer_tanks_positions = oxidizer_tanks_positions + self.oxidizer_initial_inertia = ( + (*oxidizer_initial_inertia, 0, 0, 0) + if len(oxidizer_initial_inertia) == 3 + else oxidizer_initial_inertia ) - self.solid = SolidMotor( - thrust_source, - dry_mass, - dry_inertia, - nozzle_radius, - grain_number, - grain_density, - grain_outer_radius, - grain_initial_inner_radius, - grain_initial_height, - grain_separation, - grains_center_of_mass_position, - center_of_dry_mass_position, - nozzle_position, - burn_time, - throat_radius, - reshape_thrust_curve, + + # Oxidizer mass flow rate definition and processing + self.oxidizer_mass_flow_rate = Function( + oxidizer_mass_flow_rate_curve, + "Time (s)", + "Oxidizer Mass Flow Rate (kg/s)", interpolation_method, - coordinate_system_orientation, - reference_pressure, + extrapolation="zero", + ) + + # Correct dry mass and dry inertia to include tanks mass + self.dry_mass = dry_mass + tanks_mass + dry_inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia + self.dry_I_11 = dry_inertia[0] + self.dry_I_22 = dry_inertia[1] + self.dry_I_33 = dry_inertia[2] + self.dry_I_12 = dry_inertia[3] + self.dry_I_13 = dry_inertia[4] + self.dry_I_23 = dry_inertia[5] + # TODO: Calculate tanks inertia tensor based on their geometry and mass + """ + # Initialize Tanks object + self.tanks = Tank( + geometries=oxidizer_tanks_geometries, + positions=oxidizer_tanks_positions, + fluid_mass=self.oxidizer_mass, + fluid_density=self.oxidizer_density, + initial_liquid_level=oxidizer_tanks_initial_liquid_level, + initial_ullage_mass=oxidizer_tanks_initial_ullage_mass, + initial_ullage_volume=oxidizer_tanks_initial_ullage_volume, + ) + """ + # Store important functions + self.liquid_propellant_mass = self.tanks.liquid_mass + self.gas_propellant_mass = self.tanks.gas_mass + self.center_of_liquid_propellant_mass = self.tanks.liquid_center_of_mass + self.center_of_gas_propellant_mass = self.tanks.gas_center_of_mass + self.liquid_propellant_I_11 = self.tanks.liquid_I_11 + self.liquid_propellant_I_22 = self.tanks.liquid_I_22 + self.liquid_propellant_I_33 = self.tanks.liquid_I_33 + self.gas_propellant_I_11 = self.tanks.gas_I_11 + self.gas_propellant_I_22 = self.tanks.gas_I_22 + self.gas_propellant_I_33 = self.tanks.gas_I_33 + + # Rename grain attributes for clarity + self.grain_propellant_mass = self.grain_mass + self.center_of_grain_propellant_mass = self.grains_center_of_mass_position + self.grain_propellant_I_11 = self.grains_I_11 + self.grain_propellant_I_22 = self.grains_I_22 + self.grain_propellant_I_33 = self.grains_I_33 + + # Overall propellant inertia tensor components relative to propellant CoM + # We need to recalculate the total propellant CoM function first + # (Assuming self.liquid_propellant_mass and self.grain_propellant_mass exist and are Functions) + # (Assuming self.center_of_liquid_propellant_mass and self.center_of_grain_propellant_mass exist and are Functions returning scalars) + self._propellant_mass = self.liquid_propellant_mass + self.grain_propellant_mass + self._center_of_propellant_mass = ( + self.center_of_liquid_propellant_mass * self.liquid_propellant_mass + + self.center_of_grain_propellant_mass * self.grain_propellant_mass + ) / self._propellant_mass + # Ensure division by zero is handled if needed, although propellant mass shouldn't be zero initially + + # Create Functions returning distance vectors relative to the overall propellant CoM + liquid_com_to_prop_com = ( + self.center_of_liquid_propellant_mass - self._center_of_propellant_mass + ) + grain_com_to_prop_com = ( + self.center_of_grain_propellant_mass - self._center_of_propellant_mass + ) + + # Convert scalar distances to 3D vectors for PAT functions + # The distance is along the Z-axis in the motor's coordinate system + liquid_dist_vec_func = Function( + lambda t: Vector([0, 0, liquid_com_to_prop_com(t)]), inputs="t" + ) + grain_dist_vec_func = Function( + lambda t: Vector([0, 0, grain_com_to_prop_com(t)]), inputs="t" + ) + + # Apply PAT using the new specific functions + # Inertias relative to component CoMs are needed (e.g., self.liquid_propellant_I_11_from_liquid_CM) + # Assuming these exist, otherwise adjust the first argument of the PAT functions + + # --- I_11 --- + # Assuming self.liquid_propellant_I_11 refers to inertia relative to liquid CoM + liquid_I_11_prop_com = parallel_axis_theorem_I11( + self.liquid_propellant_I_11, # Inertia relative to liquid's own CoM + self.liquid_propellant_mass, + liquid_dist_vec_func, # Distance from total prop CoM to liquid CoM + ) + # Assuming self.grain_propellant_I_11 refers to inertia relative to grain CoM + grain_I_11_prop_com = parallel_axis_theorem_I11( + self.grain_propellant_I_11, # Inertia relative to grain's own CoM + self.grain_propellant_mass, + grain_dist_vec_func, # Distance from total prop CoM to grain CoM + ) + self.propellant_I_11_from_propellant_CM = ( + liquid_I_11_prop_com + grain_I_11_prop_com + ) + + # --- I_22 --- + liquid_I_22_prop_com = parallel_axis_theorem_I22( + self.liquid_propellant_I_22, # Inertia relative to liquid's own CoM + self.liquid_propellant_mass, + liquid_dist_vec_func, + ) + grain_I_22_prop_com = parallel_axis_theorem_I22( + self.grain_propellant_I_22, # Inertia relative to grain's own CoM + self.grain_propellant_mass, + grain_dist_vec_func, + ) + self.propellant_I_22_from_propellant_CM = ( + liquid_I_22_prop_com + grain_I_22_prop_com + ) + + # --- I_33 --- + liquid_I_33_prop_com = parallel_axis_theorem_I33( + self.liquid_propellant_I_33, # Inertia relative to liquid's own CoM + self.liquid_propellant_mass, + liquid_dist_vec_func, + ) + grain_I_33_prop_com = parallel_axis_theorem_I33( + self.grain_propellant_I_33, # Inertia relative to grain's own CoM + self.grain_propellant_mass, + grain_dist_vec_func, + ) + self.propellant_I_33_from_propellant_CM = ( + liquid_I_33_prop_com + grain_I_33_prop_com + ) + + # --- Products of Inertia (I_12, I_13, I_23) --- + # Assume components PoI are 0 relative to their own CoM due to axisymmetry + # PAT calculation will correctly handle the axisymmetry (result should be 0) + + # I_12 + liquid_I_12_prop_com = parallel_axis_theorem_I12( + Function(0), self.liquid_propellant_mass, liquid_dist_vec_func + ) + grain_I_12_prop_com = parallel_axis_theorem_I12( + Function(0), self.grain_propellant_mass, grain_dist_vec_func + ) + # Store intermediate result if needed by Motor.__init__ later, prefix with '_' if not part of public API + self._propellant_I_12_from_propellant_CM = ( + liquid_I_12_prop_com + grain_I_12_prop_com + ) + + # I_13 + liquid_I_13_prop_com = parallel_axis_theorem_I13( + Function(0), self.liquid_propellant_mass, liquid_dist_vec_func + ) + grain_I_13_prop_com = parallel_axis_theorem_I13( + Function(0), self.grain_propellant_mass, grain_dist_vec_func + ) + self._propellant_I_13_from_propellant_CM = ( + liquid_I_13_prop_com + grain_I_13_prop_com + ) + + # I_23 + liquid_I_23_prop_com = parallel_axis_theorem_I23( + Function(0), self.liquid_propellant_mass, liquid_dist_vec_func + ) + grain_I_23_prop_com = parallel_axis_theorem_I23( + Function(0), self.grain_propellant_mass, grain_dist_vec_func + ) + self._propellant_I_23_from_propellant_CM = ( + liquid_I_23_prop_com + grain_I_23_prop_com ) - self.positioned_tanks = self.liquid.positioned_tanks - self.grain_number = grain_number - self.grain_density = grain_density - self.grain_outer_radius = grain_outer_radius - self.grain_initial_inner_radius = grain_initial_inner_radius - self.grain_initial_height = grain_initial_height - self.grain_separation = grain_separation - self.grains_center_of_mass_position = grains_center_of_mass_position - self.throat_radius = throat_radius - - # Initialize plots and prints object - self.prints = _HybridMotorPrints(self) + # IMPORTANT: Call the parent __init__ AFTER calculating component inertias + # because the parent __init__ uses these calculated values. + super().__init__( + thrust_source=thrust_source, + dry_mass=self.dry_mass, # Use the corrected dry mass + dry_inertia=( + self.dry_I_11, + self.dry_I_22, + self.dry_I_33, + self.dry_I_12, + self.dry_I_13, + self.dry_I_23, + ), # Use corrected dry inertia + nozzle_radius=nozzle_radius, + burn_time=burn_time, + center_of_dry_mass_position=center_of_dry_mass_position, + nozzle_position=nozzle_position, + reshape_thrust_curve=reshape_thrust_curve, + interpolation_method=interpolation_method, + coordinate_system_orientation=coordinate_system_orientation, + reference_pressure=reference_pressure, + ) + # The parent __init__ will now correctly use the calculated + # self.propellant_I_xx_from_propellant_CM values. + + # Initialize plots object specific to HybridMotor self.plots = _HybridMotorPlots(self) - @funcify_method("Time (s)", "Exhaust velocity (m/s)") def exhaust_velocity(self): """Exhaust velocity by assuming it as a constant. The formula used is total impulse/propellant initial mass. @@ -458,145 +691,51 @@ def center_of_propellant_mass(self): ) return mass_balance / self.propellant_mass + @funcify_method("Time (s)", "Inertia I_11 (kg m²)") + @property @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def propellant_I_11(self): """Inertia tensor 11 component of the propellant, the inertia is relative to the e_1 axis, centered at the instantaneous propellant center of mass. - - Returns - ------- - Function - Propellant inertia tensor 11 component at time t. - - Notes - ----- - The e_1 direction is assumed to be the direction perpendicular to the - motor body axis. - - References - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ + # Returns the value calculated in __init__ + return self.propellant_I_11_from_propellant_CM - solid_mass = self.solid.propellant_mass - liquid_mass = self.liquid.propellant_mass - - cm = self.center_of_propellant_mass - solid_cm_to_cm = self.solid.center_of_propellant_mass - cm - liquid_cm_to_cm = self.liquid.center_of_propellant_mass - cm - - solid_prop_inertia = self.solid.propellant_I_11 - liquid_prop_inertia = self.liquid.propellant_I_11 - - I_11 = parallel_axis_theorem_from_com( - solid_prop_inertia, solid_mass, solid_cm_to_cm - ) + parallel_axis_theorem_from_com( - liquid_prop_inertia, liquid_mass, liquid_cm_to_cm - ) - - return I_11 - + @property @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): - """Inertia tensor 22 component of the propellant, the inertia is - relative to the e_2 axis, centered at the instantaneous propellant - center of mass. - - Returns - ------- - Function - Propellant inertia tensor 22 component at time t. - - Notes - ----- - The e_2 direction is assumed to be the direction perpendicular to the - motor body axis, and perpendicular to e_1. + """Inertia tensor 22 component of the propellant... (Identical to I_11)""" - References - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor - """ - return self.propellant_I_11 + return self.propellant_I_22_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_33 (kg m²)") def propellant_I_33(self): - """Inertia tensor 33 component of the propellant, the inertia is - relative to the e_3 axis, centered at the instantaneous propellant - center of mass. + """Inertia tensor 33 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 33 component at time t. - - Notes - ----- - The e_3 direction is assumed to be the axial direction of the rocket - motor. - - References - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor - """ - return self.solid.propellant_I_33 + self.liquid.propellant_I_33 + return self.propellant_I_33_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): - """Inertia tensor 12 component of the propellant, the inertia is - relative to the e_1 and e_2 axes, centered at the instantaneous - propellant center of mass. + """Inertia tensor 12 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 12 component at time t. - - Notes - ----- - This is assumed to be zero due to axial symmetry of the motor. This - could be improved in the future to account for the fact that the - motor is not perfectly symmetric. - """ - return 0 + return self._propellant_I_12_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): - """Inertia tensor 13 component of the propellant, the inertia is - relative to the e_1 and e_3 axes, centered at the instantaneous - propellant center of mass. + """Inertia tensor 13 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 13 component at time t. - - Notes - ----- - This is assumed to be zero due to axial symmetry of the motor. This - could be improved in the future to account for the fact that the - motor is not perfectly symmetric. - """ - return 0 + return self._propellant_I_13_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): - """Inertia tensor 23 component of the propellant, the inertia is - relative to the e_2 and e_3 axes, centered at the instantaneous - propellant center of mass. + """Inertia tensor 23 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 23 component at time t. - - Notes - ----- - This is assumed to be zero due to axial symmetry of the motor. This - could be improved in the future to account for the fact that the - motor is not perfectly symmetric. - """ - return 0 + return self._propellant_I_23_from_propellant_CM def add_tank(self, tank, position): """Adds a tank to the motor. @@ -641,8 +780,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "grain_number": self.grain_number, @@ -660,18 +799,13 @@ def to_dict(self, **kwargs): } ) - if kwargs.get("include_outputs", False): - burn_rate = self.solid.burn_rate - if kwargs.get("discretize", False): - burn_rate = burn_rate.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) + if include_outputs: data.update( { "grain_inner_radius": self.solid.grain_inner_radius, "grain_height": self.solid.grain_height, "burn_area": self.solid.burn_area, - "burn_rate": burn_rate, + "burn_rate": self.solid.burn_rate, } ) @@ -703,11 +837,12 @@ def from_dict(cls, data): grain_separation=data["grain_separation"], grains_center_of_mass_position=data["grains_center_of_mass_position"], nozzle_position=data["nozzle_position"], - throat_radius=data["throat_radius"], + tanks_mass=data["tanks_mass"], + oxidizer_initial_mass=data["oxidizer_initial_mass"], + oxidizer_mass_flow_rate_curve=data["oxidizer_mass_flow_rate_curve"], + oxidizer_density=data["oxidizer_density"], + oxidizer_tanks_geometries=data["oxidizer_tanks_geometries"], + oxidizer_tanks_positions=data["oxidizer_tanks_positions"], reference_pressure=data.get("reference_pressure"), ) - - for tank in data["positioned_tanks"]: - motor.add_tank(tank["tank"], tank["position"]) - return motor diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py index 1cd76a52f..88c695e8d 100644 --- a/rocketpy/motors/liquid_motor.py +++ b/rocketpy/motors/liquid_motor.py @@ -1,23 +1,16 @@ from functools import cached_property - import numpy as np - from rocketpy.mathutils.function import funcify_method, reset_funcified_methods -from rocketpy.tools import parallel_axis_theorem_from_com - from ..plots.liquid_motor_plots import _LiquidMotorPlots -from ..prints.liquid_motor_prints import _LiquidMotorPrints from .motor import Motor class LiquidMotor(Motor): """Class to specify characteristics and useful operations for Liquid motors. This class inherits from the Motor class. - See Also -------- Motor - Attributes ---------- LiquidMotor.coordinate_system_orientation : str @@ -157,113 +150,44 @@ class LiquidMotor(Motor): It will allow to obtain the net thrust in the Flight class. """ + # pylint: disable=too-many-locals, too-many-statements, too-many-arguments def __init__( self, thrust_source, dry_mass, dry_inertia, nozzle_radius, + burn_time, center_of_dry_mass_position, nozzle_position=0, - burn_time=None, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, + tanks_mass=0, ): - """Initialize LiquidMotor class, process thrust curve and geometrical - parameters and store results. - Parameters - ---------- - thrust_source : int, float, callable, string, array, Function - Motor's thrust curve. Can be given as an int or float, in which - case the thrust will be considered constant in time. It can - also be given as a callable function, whose argument is time in - seconds and returns the thrust supplied by the motor in the - instant. If a string is given, it must point to a .csv or .eng file. - The .csv file can contain a single line header and the first column - must specify time in seconds, while the second column specifies - thrust. Arrays may also be specified, following rules set by the - class Function. Thrust units are Newtons. - - .. seealso:: :doc:`Thrust Source Details ` - dry_mass : int, float - Same as in Motor class. See the :class:`Motor ` docs. - dry_inertia : tuple, list - Tuple or list containing the motor's dry mass inertia tensor - components, in kg*m^2. This inertia is defined with respect to the - the ``center_of_dry_mass_position`` position. - Assuming e_3 is the rocket's axis of symmetry, e_1 and e_2 are - orthogonal and form a plane perpendicular to e_3, the dry mass - inertia tensor components must be given in the following order: - (I_11, I_22, I_33, I_12, I_13, I_23), where I_ij is the - component of the inertia tensor in the direction of e_i x e_j. - Alternatively, the inertia tensor can be given as - (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. - nozzle_radius : int, float - Motor's nozzle outlet radius in meters. - center_of_dry_mass_position : int, float - The position, in meters, of the motor's center of mass with respect - to the motor's coordinate system when it is devoid of propellant. - See :doc:`Positions and Coordinate Systems ` - nozzle_position : float - Motor's nozzle outlet position in meters, specified in the motor's - coordinate system. See - :doc:`Positions and Coordinate Systems ` for - more information. - burn_time: float, tuple of float, optional - Motor's burn time. - If a float is given, the burn time is assumed to be between 0 and - the given float, in seconds. - If a tuple of float is given, the burn time is assumed to be between - the first and second elements of the tuple, in seconds. - If not specified, automatically sourced as the range between the - first and last-time step of the motor's thrust curve. This can only - be used if the motor's thrust is defined by a list of points, such - as a .csv file, a .eng file or a Function instance whose source is - a list. - reshape_thrust_curve : boolean, tuple, optional - If False, the original thrust curve supplied is not altered. If a - tuple is given, whose first parameter is a new burn out time and - whose second parameter is a new total impulse in Ns, the thrust - curve is reshaped to match the new specifications. May be useful - for motors whose thrust curve shape is expected to remain similar - in case the impulse and burn time varies slightly. Default is - False. - interpolation_method : string, optional - Method of interpolation to be used in case thrust curve is given - by data set in .csv or .eng, or as an array. Options are 'spline' - 'akima' and 'linear'. Default is "linear". - coordinate_system_orientation : string, optional - Orientation of the motor's coordinate system. The coordinate system - is defined by the motor's axis of symmetry. The origin of the - coordinate system may be placed anywhere along such axis, such as - at the nozzle area, and must be kept the same for all other - positions specified. Options are "nozzle_to_combustion_chamber" - and "combustion_chamber_to_nozzle". Default is - "nozzle_to_combustion_chamber". - reference_pressure : int, float, optional - Atmospheric pressure in Pa at which the thrust data was recorded. - """ + # Initialize the list of tanks + self.positioned_tanks = [] + + # Correct the dry mass to include the mass of the tanks + dry_mass = dry_mass + tanks_mass + dry_inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia + super().__init__( thrust_source=thrust_source, + dry_mass=dry_mass, dry_inertia=dry_inertia, nozzle_radius=nozzle_radius, + burn_time=burn_time, center_of_dry_mass_position=center_of_dry_mass_position, - dry_mass=dry_mass, nozzle_position=nozzle_position, - burn_time=burn_time, reshape_thrust_curve=reshape_thrust_curve, interpolation_method=interpolation_method, coordinate_system_orientation=coordinate_system_orientation, reference_pressure=reference_pressure, ) - self.positioned_tanks = [] - - # Initialize plots and prints object - self.prints = _LiquidMotorPrints(self) self.plots = _LiquidMotorPlots(self) @funcify_method("Time (s)", "Exhaust Velocity (m/s)") @@ -374,36 +298,16 @@ def center_of_propellant_mass(self): @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def propellant_I_11(self): - """Inertia tensor 11 component of the propellant, the inertia is - relative to the e_1 axis, centered at the instantaneous propellant - center of mass. + """Inertia tensor 11 component of the total propellant, the inertia is + relative to the e_1 axis, centered at the instantaneous total propellant + center of mass. Recalculated here relative to the instantaneous CoM. Returns ------- Function - Propellant inertia tensor 11 component at time t. - - Notes - ----- - The e_1 direction is assumed to be the direction perpendicular to the - motor body axis. - - References - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor + Total propellant inertia tensor 11 component at time t relative to total propellant CoM. """ - I_11 = 0 - center_of_mass = self.center_of_propellant_mass - - for positioned_tank in self.positioned_tanks: - tank = positioned_tank.get("tank") - tank_position = positioned_tank.get("position") - distance = tank_position + tank.center_of_mass - center_of_mass - I_11 += parallel_axis_theorem_from_com( - tank.inertia, tank.fluid_mass, distance - ) - - return I_11 + return self.propellant_I_11_from_propellant_CM @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): @@ -497,8 +401,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "positioned_tanks": [ diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py index 7930ed52b..89e33b47d 100644 --- a/rocketpy/motors/motor.py +++ b/rocketpy/motors/motor.py @@ -4,13 +4,20 @@ from abc import ABC, abstractmethod from functools import cached_property from os import path - import numpy as np - from ..mathutils.function import Function, funcify_method +from ..mathutils.vector_matrix import Vector from ..plots.motor_plots import _MotorPlots from ..prints.motor_prints import _MotorPrints -from ..tools import parallel_axis_theorem_from_com, tuple_handler +from ..tools import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I22, + parallel_axis_theorem_I33, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I23, + tuple_handler, +) # pylint: disable=too-many-public-methods @@ -165,177 +172,153 @@ class Motor(ABC): def __init__( self, thrust_source, + dry_mass, dry_inertia, nozzle_radius, + burn_time, center_of_dry_mass_position, - dry_mass=None, nozzle_position=0, - burn_time=None, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, ): - """Initialize Motor class, process thrust curve and geometrical + """Initializes Motor class, process thrust curve and geometrical parameters and store results. - Parameters - ---------- - thrust_source : int, float, callable, string, array, Function - Motor's thrust curve. Can be given as an int or float, in which - case the thrust will be considered constant in time. It can - also be given as a callable function, whose argument is time in - seconds and returns the thrust supplied by the motor in the - instant. If a string is given, it must point to a .csv or .eng file. - The .csv file can contain a single line header and the first column - must specify time in seconds, while the second column specifies - thrust. Arrays may also be specified, following rules set by the - class Function. Thrust units are Newtons. - - .. seealso:: :doc:`Thrust Source Details ` - - dry_mass : int, float - Same as in Motor class. See the :class:`Motor ` docs - center_of_dry_mass_position : int, float - The position, in meters, of the motor's center of mass with respect - to the motor's coordinate system when it is devoid of propellant. - See :doc:`Positions and Coordinate Systems ` - dry_inertia : tuple, list - Tuple or list containing the motor's dry mass inertia tensor - components, in kg*m^2. This inertia is defined with respect to the - the `center_of_dry_mass_position` position. - Assuming e_3 is the rocket's axis of symmetry, e_1 and e_2 are - orthogonal and form a plane perpendicular to e_3, the dry mass - inertia tensor components must be given in the following order: - (I_11, I_22, I_33, I_12, I_13, I_23), where I_ij is the - component of the inertia tensor in the direction of e_i x e_j. - Alternatively, the inertia tensor can be given as - (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. - nozzle_radius : int, float, optional - Motor's nozzle outlet radius in meters. - burn_time: float, tuple of float, optional - Motor's burn time. - If a float is given, the burn time is assumed to be between 0 and - the given float, in seconds. - If a tuple of float is given, the burn time is assumed to be between - the first and second elements of the tuple, in seconds. - If not specified, automatically sourced as the range between the - first and last-time step of the motor's thrust curve. This can only - be used if the motor's thrust is defined by a list of points, such - as a .csv file, a .eng file or a Function instance whose source is a - list. - nozzle_position : int, float, optional - Motor's nozzle outlet position in meters, in the motor's coordinate - system. See :doc:`Positions and Coordinate Systems ` - for details. Default is 0, in which case the origin of the - coordinate system is placed at the motor's nozzle outlet. - reshape_thrust_curve : boolean, tuple, optional - If False, the original thrust curve supplied is not altered. If a - tuple is given, whose first parameter is a new burn out time and - whose second parameter is a new total impulse in Ns, the thrust - curve is reshaped to match the new specifications. May be useful - for motors whose thrust curve shape is expected to remain similar - in case the impulse and burn time varies slightly. Default is - False. Note that the Motor burn_time parameter must include the new - reshaped burn time. - interpolation_method : string, optional - Method of interpolation to be used in case thrust curve is given - by data set in .csv or .eng, or as an array. Options are 'spline' - 'akima' and 'linear'. Default is "linear". - coordinate_system_orientation : string, optional - Orientation of the motor's coordinate system. The coordinate system - is defined by the motor's axis of symmetry. The origin of the - coordinate system may be placed anywhere along such axis, such as - at the nozzle area, and must be kept the same for all other - positions specified. Options are "nozzle_to_combustion_chamber" and - "combustion_chamber_to_nozzle". Default is - "nozzle_to_combustion_chamber". - reference_pressure : int, float, optional - Atmospheric pressure in Pa at which the thrust data was recorded. - - Returns - ------- - None + [... Docstring inchangée ...] """ + # Motor attributes + self.nozzle_radius = nozzle_radius + self.nozzle_position = nozzle_position + self.dry_mass = dry_mass + self.center_of_dry_mass_position = center_of_dry_mass_position + dry_inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia + self.dry_I_11 = dry_inertia[0] + self.dry_I_22 = dry_inertia[1] + self.dry_I_33 = dry_inertia[2] + self.dry_I_12 = dry_inertia[3] + self.dry_I_13 = dry_inertia[4] + self.dry_I_23 = dry_inertia[5] + # Define coordinate system orientation - self.coordinate_system_orientation = coordinate_system_orientation if coordinate_system_orientation == "nozzle_to_combustion_chamber": self._csys = 1 elif coordinate_system_orientation == "combustion_chamber_to_nozzle": self._csys = -1 - else: # pragma: no cover - raise ValueError( - "Invalid coordinate system orientation. Options are " - "'nozzle_to_combustion_chamber' and 'combustion_chamber_to_nozzle'." + else: + raise TypeError( + "Invalid coordinate system orientation. Please choose between " + + '"nozzle_to_combustion_chamber" and ' + + '"combustion_chamber_to_nozzle".' ) + self.coordinate_system_orientation = coordinate_system_orientation - # Motor parameters - self.interpolate = interpolation_method - self.nozzle_position = nozzle_position - self.nozzle_radius = nozzle_radius - self.nozzle_area = np.pi * nozzle_radius**2 - self.center_of_dry_mass_position = center_of_dry_mass_position self.reference_pressure = reference_pressure - # Inertia tensor setup - inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia - self.dry_I_11 = inertia[0] - self.dry_I_22 = inertia[1] - self.dry_I_33 = inertia[2] - self.dry_I_12 = inertia[3] - self.dry_I_13 = inertia[4] - self.dry_I_23 = inertia[5] - - # Handle .eng or .rse file inputs - self.description_eng_file = None - self.rse_motor_data = None - if isinstance(thrust_source, str): - if ( - path.exists(thrust_source) - and path.splitext(path.basename(thrust_source))[1] == ".eng" - ): - _, self.description_eng_file, points = Motor.import_eng(thrust_source) - thrust_source = points - elif ( - path.exists(thrust_source) - and path.splitext(path.basename(thrust_source))[1] == ".rse" - ): - self.rse_motor_data, points = Motor.import_rse(thrust_source) - thrust_source = points - # Evaluate raw thrust source - self.thrust_source = thrust_source + # Thrust curve definition and processing self.thrust = Function( - thrust_source, "Time (s)", "Thrust (N)", self.interpolate, "zero" + thrust_source, + "Time (s)", + "Thrust (N)", + interpolation_method, + extrapolation="zero", ) - # Handle dry_mass input - self.dry_mass = dry_mass + self.interpolate = interpolation_method + # Burn time definition + self.burn_time = tuple_handler(burn_time) + self.burn_start_time = self.burn_time[0] + self.burn_out_time = self.burn_time[1] - # Handle burn_time input - self.burn_time = burn_time + # Reshape thrust curve if needed + self.reshape_thrust_curve = reshape_thrust_curve + if reshape_thrust_curve: + self._reshape_thrust_curve(*reshape_thrust_curve) - if callable(self.thrust.source): - self.thrust.set_discrete(*self.burn_time, 50, self.interpolate, "zero") + # Basic calculations and attributes + self.burn_duration = self.burn_out_time - self.burn_start_time + self.total_impulse = self.thrust.integral( + self.burn_start_time, self.burn_out_time + ) + self.max_thrust = self.thrust.max - # Reshape thrust_source if needed - if reshape_thrust_curve: - # Overwrites burn_time and thrust - self.thrust = Motor.reshape_thrust_curve(self.thrust, *reshape_thrust_curve) - self.burn_time = (self.thrust.x_array[0], self.thrust.x_array[-1]) + self.average_thrust = self.total_impulse / self.burn_duration - # Post process thrust - self.thrust = Motor.clip_thrust(self.thrust, self.burn_time) + # Abstract methods - must be implemented by subclasses + self._propellant_initial_mass = 0 + self.exhaust_velocity = Function(0) + self.propellant_mass = Function(0) + self.total_mass = Function(0) + self.propellant_I_11_from_propellant_CM = Function(0) + self.propellant_I_22_from_propellant_CM = Function(0) + self.propellant_I_33_from_propellant_CM = Function(0) + self.center_of_propellant_mass = Function(0) + self.center_of_mass = Function(0) + + # Inertia tensor for propellant referenced to the MOTOR's origin + # Note: distance_vec_3d is the vector from the motor origin to the propellant CoM + propellant_com_func = self.center_of_propellant_mass + + # Create a Function that returns the distance vector [0, 0, center_of_propellant_mass(t)] + propellant_com_vector_func = Function( + lambda t: Vector([0, 0, propellant_com_func(t)]), + inputs="t", + outputs="Vector (m)", + ) - # Auxiliary quantities - self.burn_start_time = self.burn_time[0] - self.burn_out_time = self.burn_time[1] - self.burn_duration = self.burn_time[1] - self.burn_time[0] + # Use the new specific PAT functions + self.propellant_I_11 = parallel_axis_theorem_I11( + self.propellant_I_11_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_11.set_outputs("Propellant I_11 (kg*m^2)") - # Compute thrust metrics - self.max_thrust = np.amax(self.thrust.y_array) - max_thrust_index = np.argmax(self.thrust.y_array) - self.max_thrust_time = self.thrust.source[max_thrust_index, 0] - self.average_thrust = self.total_impulse / self.burn_duration + self.propellant_I_22 = parallel_axis_theorem_I22( + self.propellant_I_22_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_22.set_outputs("Propellant I_22 (kg*m^2)") + + self.propellant_I_33 = parallel_axis_theorem_I33( + self.propellant_I_33_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_33.set_outputs("Propellant I_33 (kg*m^2)") + + self.propellant_I_12 = parallel_axis_theorem_I12( + Function(0), self.propellant_mass, propellant_com_vector_func + ) + self.propellant_I_12.set_outputs("Propellant I_12 (kg*m^2)") + + self.propellant_I_13 = parallel_axis_theorem_I13( + Function(0), self.propellant_mass, propellant_com_vector_func + ) + self.propellant_I_13.set_outputs("Propellant I_13 (kg*m^2)") + + self.propellant_I_23 = parallel_axis_theorem_I23( + Function(0), self.propellant_mass, propellant_com_vector_func + ) + self.propellant_I_23.set_outputs("Propellant I_23 (kg*m^2)") + + # Calculate total motor inertia relative to motor's origin + self.I_11 = Function(lambda t: self.dry_I_11) + self.propellant_I_11 + self.I_11.set_outputs("Motor I_11 (kg*m^2)") + self.I_22 = Function(lambda t: self.dry_I_22) + self.propellant_I_22 + self.I_22.set_outputs("Motor I_22 (kg*m^2)") + self.I_33 = Function(lambda t: self.dry_I_33) + self.propellant_I_33 + self.I_33.set_outputs("Motor I_33 (kg*m^2)") + + # Calculate total products of inertia relative to motor's origin + self.I_12 = Function(lambda t: self.dry_I_12) + self.propellant_I_12 + self.I_12.set_outputs("Motor I_12 (kg*m^2)") + self.I_13 = Function(lambda t: self.dry_I_13) + self.propellant_I_13 + self.I_13.set_outputs("Motor I_13 (kg*m^2)") + self.I_23 = Function(lambda t: self.dry_I_23) + self.propellant_I_23 + self.I_23.set_outputs("Motor I_23 (kg*m^2)") # Initialize plots and prints object self.prints = _MotorPrints(self) @@ -575,36 +558,30 @@ def center_of_propellant_mass(self): @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def I_11(self): - """Inertia tensor 11 component, which corresponds to the inertia - relative to the e_1 axis, centered at the instantaneous center of mass. - - Returns - ------- - Function - Propellant inertia tensor 11 component at time t. - - Notes - ----- - The e_1 direction is assumed to be the direction perpendicular to the - motor body axis. Also, due to symmetry, I_11 = I_22. - - See Also - -------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor + """Builds a Function for the total I_11 inertia component relative + to the instantaneous center of mass of the motor. """ + # Inertias relative to the motor origin (calculated in __init__) + prop_I_11_origin = self.propellant_I_11 + dry_I_11_origin = Function(lambda t: self.dry_I_11) # Dry inertia is constant + + # Distance vectors FROM the instantaneous motor CoM TO the component CoMs + prop_com_to_inst_com = self.center_of_propellant_mass - self.center_of_mass + dry_com_to_inst_com = ( + Function(lambda t: Vector([0, 0, self.center_of_dry_mass_position])) + - self.center_of_mass + ) - prop_I_11 = self.propellant_I_11 - dry_I_11 = self.dry_I_11 - - prop_to_cm = self.center_of_propellant_mass - self.center_of_mass - dry_to_cm = self.center_of_dry_mass_position - self.center_of_mass - - prop_I_11 = parallel_axis_theorem_from_com( - prop_I_11, self.propellant_mass, prop_to_cm + # Apply PAT relative to the instantaneous motor CoM + # Note: We need the negative of the distance vector for PAT formula (origin to point) + prop_I_11_cm = parallel_axis_theorem_I11( + prop_I_11_origin, self.propellant_mass, -prop_com_to_inst_com + ) + dry_I_11_cm = parallel_axis_theorem_I11( + dry_I_11_origin, self.dry_mass, -dry_com_to_inst_com ) - dry_I_11 = parallel_axis_theorem_from_com(dry_I_11, self.dry_mass, dry_to_cm) - return prop_I_11 + dry_I_11 + return prop_I_11_cm + dry_I_11_cm @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def I_22(self): @@ -859,7 +836,7 @@ def propellant_I_13(self): Returns ------- Function - Propellant inertia tensor 13 component at time t. + Propellant inertia tensor 13 components at time t. Notes ----- @@ -1231,9 +1208,16 @@ def get_attr_value(obj, attr_name, multiplier=1): # Write last line file.write(f"{self.thrust.source[-1, 0]:.4f} {0:.3f}\n") - def to_dict(self, **kwargs): + def to_dict(self, include_outputs=False): + thrust_source = self.thrust_source + + if isinstance(thrust_source, str): + thrust_source = self.thrust.source + elif callable(thrust_source) and not isinstance(thrust_source, Function): + thrust_source = Function(thrust_source) + data = { - "thrust_source": self.thrust, + "thrust_source": thrust_source, "dry_I_11": self.dry_I_11, "dry_I_22": self.dry_I_22, "dry_I_33": self.dry_I_33, @@ -1251,94 +1235,31 @@ def to_dict(self, **kwargs): "reference_pressure": self.reference_pressure, } - if kwargs.get("include_outputs", False): - total_mass = self.total_mass - propellant_mass = self.propellant_mass - mass_flow_rate = self.total_mass_flow_rate - center_of_mass = self.center_of_mass - center_of_propellant_mass = self.center_of_propellant_mass - exhaust_velocity = self.exhaust_velocity - I_11 = self.I_11 - I_22 = self.I_22 - I_33 = self.I_33 - I_12 = self.I_12 - I_13 = self.I_13 - I_23 = self.I_23 - propellant_I_11 = self.propellant_I_11 - propellant_I_22 = self.propellant_I_22 - propellant_I_33 = self.propellant_I_33 - propellant_I_12 = self.propellant_I_12 - propellant_I_13 = self.propellant_I_13 - propellant_I_23 = self.propellant_I_23 - if kwargs.get("discretize", False): - total_mass = total_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_mass = propellant_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - mass_flow_rate = mass_flow_rate.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - center_of_mass = center_of_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - center_of_propellant_mass = ( - center_of_propellant_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - ) - exhaust_velocity = exhaust_velocity.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - I_11 = I_11.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_22 = I_22.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_33 = I_33.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_12 = I_12.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_13 = I_13.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_23 = I_23.set_discrete_based_on_model(self.thrust, mutate_self=False) - propellant_I_11 = propellant_I_11.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_22 = propellant_I_22.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_33 = propellant_I_33.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_12 = propellant_I_12.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_13 = propellant_I_13.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_23 = propellant_I_23.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) + if include_outputs: data.update( { "vacuum_thrust": self.vacuum_thrust, - "total_mass": total_mass, - "propellant_mass": propellant_mass, - "mass_flow_rate": mass_flow_rate, - "center_of_mass": center_of_mass, - "center_of_propellant_mass": center_of_propellant_mass, + "total_mass": self.total_mass, + "propellant_mass": self.propellant_mass, + "mass_flow_rate": self.mass_flow_rate, + "center_of_mass": self.center_of_mass, + "center_of_propellant_mass": self.center_of_propellant_mass, "total_impulse": self.total_impulse, - "exhaust_velocity": exhaust_velocity, + "exhaust_velocity": self.exhaust_velocity, "propellant_initial_mass": self.propellant_initial_mass, "structural_mass_ratio": self.structural_mass_ratio, - "I_11": I_11, - "I_22": I_22, - "I_33": I_33, - "I_12": I_12, - "I_13": I_13, - "I_23": I_23, - "propellant_I_11": propellant_I_11, - "propellant_I_22": propellant_I_22, - "propellant_I_33": propellant_I_33, - "propellant_I_12": propellant_I_12, - "propellant_I_13": propellant_I_13, - "propellant_I_23": propellant_I_23, + "I_11": self.I_11, + "I_22": self.I_22, + "I_33": self.I_33, + "I_12": self.I_12, + "I_13": self.I_13, + "I_23": self.I_23, + "propellant_I_11": self.propellant_I_11, + "propellant_I_22": self.propellant_I_22, + "propellant_I_33": self.propellant_I_33, + "propellant_I_12": self.propellant_I_12, + "propellant_I_13": self.propellant_I_13, + "propellant_I_23": self.propellant_I_23, } ) @@ -1566,7 +1487,7 @@ def center_of_propellant_mass(self): Function Function representing the center of mass of the motor. """ - return Function(self.chamber_position).set_discrete_based_on_model(self.thrust) + return self.chamber_position @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def propellant_I_11(self): @@ -1588,11 +1509,11 @@ def propellant_I_11(self): ---------- https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - return Function( + return ( self.propellant_mass * (3 * self.chamber_radius**2 + self.chamber_height**2) / 12 - ).set_discrete_based_on_model(self.thrust) + ) @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): @@ -1636,21 +1557,19 @@ def propellant_I_33(self): ---------- https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - return Function( - self.propellant_mass * self.chamber_radius**2 / 2 - ).set_discrete_based_on_model(self.thrust) + return self.propellant_mass * self.chamber_radius**2 / 2 @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): - return Function(0).set_discrete_based_on_model(self.thrust) + return Function(0) @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): - return Function(0).set_discrete_based_on_model(self.thrust) + return Function(0) @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): - return Function(0).set_discrete_based_on_model(self.thrust) + return Function(0) @staticmethod def load_from_eng_file( @@ -1920,8 +1839,8 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "chamber_radius": self.chamber_radius, diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py index 8a00eeec9..04e0d3b25 100644 --- a/rocketpy/motors/solid_motor.py +++ b/rocketpy/motors/solid_motor.py @@ -7,8 +7,15 @@ from ..plots.solid_motor_plots import _SolidMotorPlots from ..prints.solid_motor_prints import _SolidMotorPrints from .motor import Motor - - +from ..tools import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I22, + parallel_axis_theorem_I33, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I23, + ) +from ..mathutils.vector_matrix import Vector class SolidMotor(Motor): """Class to specify characteristics and useful operations for solid motors. @@ -195,7 +202,6 @@ class SolidMotor(Motor): It will allow to obtain the net thrust in the Flight class. """ - # pylint: disable=too-many-arguments def __init__( self, thrust_source, @@ -319,6 +325,24 @@ class Function. Thrust units are Newtons. ------- None """ + self.rse_motor_data = None + self.description_eng_file = None + + if isinstance(thrust_source, str): + if thrust_source.endswith(".eng"): + + comments, description, thrust_source_data = Motor.import_eng( + thrust_source + ) + self.description_eng_file = description + self.comments_eng_file = comments + thrust_source = thrust_source_data + elif thrust_source.endswith(".rse"): + + rse_data, thrust_source_data = Motor.import_rse(thrust_source) + self.rse_motor_data = rse_data + thrust_source = thrust_source_data + super().__init__( thrust_source=thrust_source, dry_inertia=dry_inertia, @@ -332,11 +356,10 @@ class Function. Thrust units are Newtons. coordinate_system_orientation=coordinate_system_orientation, reference_pressure=reference_pressure, ) - # Nozzle parameters + self.throat_radius = throat_radius self.throat_area = np.pi * throat_radius**2 - # Grain parameters self.grains_center_of_mass_position = grains_center_of_mass_position self.grain_number = grain_number self.grain_separation = grain_separation @@ -345,7 +368,6 @@ class Function. Thrust units are Newtons. self.grain_initial_inner_radius = grain_initial_inner_radius self.grain_initial_height = grain_initial_height - # Grains initial geometrical parameters self.grain_initial_volume = ( self.grain_initial_height * np.pi @@ -355,9 +377,60 @@ class Function. Thrust units are Newtons. self.evaluate_geometry() - # Initialize plots and prints object self.prints = _SolidMotorPrints(self) self.plots = _SolidMotorPlots(self) + self.propellant_I_11_from_propellant_CM = self.propellant_I_11 + self.propellant_I_22_from_propellant_CM = self.propellant_I_22 + self.propellant_I_33_from_propellant_CM = self.propellant_I_33 + self.propellant_I_12_from_propellant_CM = self.propellant_I_12 + self.propellant_I_13_from_propellant_CM = self.propellant_I_13 + self.propellant_I_23_from_propellant_CM = self.propellant_I_23 + propellant_com_func = self.center_of_propellant_mass + + propellant_com_vector_func = Function( + lambda t: Vector([0, 0, propellant_com_func(t)]), + inputs="t", + outputs="Vector (m)", + ) + + # Utiliser les nouvelles fonctions PAT + self.propellant_I_11 = parallel_axis_theorem_I11( + self.propellant_I_11_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_22 = parallel_axis_theorem_I22( + self.propellant_I_22_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_33 = parallel_axis_theorem_I33( + self.propellant_I_33_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_12 = parallel_axis_theorem_I12( + self.propellant_I_12_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_13 = parallel_axis_theorem_I13( + self.propellant_I_13_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_23 = parallel_axis_theorem_I23( + self.propellant_I_23_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + + self.I_11 = Function(lambda t: self.dry_I_11) + self.propellant_I_11 + self.I_22 = Function(lambda t: self.dry_I_22) + self.propellant_I_22 + self.I_33 = Function(lambda t: self.dry_I_33) + self.propellant_I_33 + self.I_12 = Function(lambda t: self.dry_I_12) + self.propellant_I_12 + self.I_13 = Function(lambda t: self.dry_I_13) + self.propellant_I_13 + self.I_23 = Function(lambda t: self.dry_I_23) + self.propellant_I_23 @funcify_method("Time (s)", "Mass (kg)") def propellant_mass(self): @@ -765,8 +838,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "nozzle_radius": self.nozzle_radius, @@ -781,18 +854,13 @@ def to_dict(self, **kwargs): } ) - if kwargs.get("include_outputs", False): - burn_rate = self.burn_rate - if kwargs.get("discretize", False): - burn_rate = burn_rate.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) + if include_outputs: data.update( { "grain_inner_radius": self.grain_inner_radius, "grain_height": self.grain_height, "burn_area": self.burn_area, - "burn_rate": burn_rate, + "burn_rate": self.burn_rate, "Kn": self.Kn, } ) diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py index 8eaaded16..612d63f80 100644 --- a/rocketpy/plots/rocket_plots.py +++ b/rocketpy/plots/rocket_plots.py @@ -2,6 +2,7 @@ import numpy as np from rocketpy.motors import EmptyMotor, HybridMotor, LiquidMotor, SolidMotor +from rocketpy.motors.ClusterMotor import ClusterMotor from rocketpy.rocket.aero_surface import Fins, NoseCone, Tail from rocketpy.rocket.aero_surface.generic_surface import GenericSurface @@ -201,7 +202,7 @@ def draw(self, vis_args=None, plane="xz", *, filename=None): drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args, plane) last_radius, last_x = self._draw_tubes(ax, drawn_surfaces, vis_args) - self._draw_motor(last_radius, last_x, ax, vis_args) + self._draw_motor(last_radius, last_x, ax, vis_args, plane) self._draw_rail_buttons(ax, vis_args) self._draw_center_of_mass_and_pressure(ax) self._draw_sensors(ax, self.rocket.sensors, plane) @@ -224,15 +225,7 @@ def __validate_aerodynamic_surfaces(self): def _draw_aerodynamic_surfaces(self, ax, vis_args, plane): """Draws the aerodynamic surfaces and saves the position of the points of interest for the tubes.""" - # List of drawn surfaces with the position of points of interest - # and the radius of the rocket at that point drawn_surfaces = [] - # Idea is to get the shape of each aerodynamic surface in their own - # coordinate system and then plot them in the rocket coordinate system - # using the position of each surface - # For the tubes, the surfaces need to be checked in order to check for - # diameter changes. The final point of the last surface is the final - # point of the last tube for surface, position in self.rocket.aerodynamic_surfaces: if isinstance(surface, NoseCone): @@ -264,14 +257,14 @@ def _draw_nose_cone(self, ax, surface, position, drawn_surfaces, vis_args): color=vis_args["nose"], linewidth=vis_args["line_width"], ) - # close the nosecone + ax.plot( [x_nosecone[-1], x_nosecone[-1]], [y_nosecone[-1], -y_nosecone[-1]], color=vis_args["nose"], linewidth=vis_args["line_width"], ) - # Add the nosecone to the list of drawn surfaces + drawn_surfaces.append( (surface, x_nosecone[-1], surface.rocket_radius, x_nosecone[-1]) ) @@ -287,7 +280,7 @@ def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args): ax.plot( x_tail, -y_tail, color=vis_args["tail"], linewidth=vis_args["line_width"] ) - # close above and below the tail + ax.plot( [x_tail[-1], x_tail[-1]], [y_tail[-1], -y_tail[-1]], @@ -300,7 +293,7 @@ def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args): color=vis_args["tail"], linewidth=vis_args["line_width"], ) - # Add the tail to the list of drawn surfaces + drawn_surfaces.append((surface, position, surface.bottom_radius, x_tail[-1])) def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): @@ -312,16 +305,15 @@ def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): rotation_angles = [2 * np.pi * i / num_fins for i in range(num_fins)] for angle in rotation_angles: - # Create a rotation matrix for the current angle around the x-axis + rotation_matrix = np.array([[1, 0], [0, np.cos(angle)]]) - # Apply the rotation to the original fin points + rotated_points_2d = np.dot(rotation_matrix, np.vstack((x_fin, y_fin))) - # Extract x and y coordinates of the rotated points x_rotated, y_rotated = rotated_points_2d - # Project points above the XY plane back into the XY plane (set z-coordinate to 0) + x_rotated = np.where( rotated_points_2d[1] > 0, rotated_points_2d[0], x_rotated ) @@ -343,20 +335,20 @@ def _draw_generic_surface( surface, position, drawn_surfaces, - vis_args, # pylint: disable=unused-argument + vis_args, plane, ): """Draws the generic surface and saves the position of the points of interest for the tubes.""" if plane == "xz": - # z position of the sensor is the x position in the plot + x_pos = position[2] - # x position of the surface is the y position in the plot + y_pos = position[0] elif plane == "yz": - # z position of the surface is the x position in the plot + x_pos = position[2] - # y position of the surface is the y position in the plot + y_pos = position[1] else: # pragma: no cover raise ValueError("Plane must be 'xz' or 'yz'.") @@ -375,22 +367,19 @@ def _draw_tubes(self, ax, drawn_surfaces, vis_args): radius = 0 last_x = 0 for i, d_surface in enumerate(drawn_surfaces): - # Draw the tubes, from the end of the first surface to the beginning - # of the next surface, with the radius of the rocket at that point + surface, position, radius, last_x = d_surface if i == len(drawn_surfaces) - 1: - # If the last surface is a tail, do nothing + if isinstance(surface, Tail): continue - # Else goes to the end of the surface + x_tube = [position, last_x] y_tube = [radius, radius] y_tube_negated = [-radius, -radius] else: - # If it is not the last surface, the tube goes to the beginning - # of the next surface - # [next_surface, next_position, next_radius, next_last_x] + next_position = drawn_surfaces[i + 1][1] x_tube = [last_x, next_position] y_tube = [radius, radius] @@ -410,45 +399,128 @@ def _draw_tubes(self, ax, drawn_surfaces, vis_args): ) return radius, last_x - def _draw_motor(self, last_radius, last_x, ax, vis_args): + def _draw_motor(self, last_radius, last_x, ax, vis_args, plane="xz"): """Draws the motor from motor patches""" - total_csys = self.rocket._csys * self.rocket.motor._csys - nozzle_position = ( - self.rocket.motor_position + self.rocket.motor.nozzle_position * total_csys - ) - # Get motor patches translated to the correct position - motor_patches = self._generate_motor_patches(total_csys, ax) + motor_patches = self._generate_motor_patches(ax, plane) - # Draw patches if not isinstance(self.rocket.motor, EmptyMotor): - # Add nozzle last so it is in front of the other patches - nozzle = self.rocket.motor.plots._generate_nozzle( - translate=(nozzle_position, 0), csys=self.rocket._csys - ) - motor_patches += [nozzle] - outline = self.rocket.motor.plots._generate_motor_region( - list_of_patches=motor_patches - ) - # add outline first so it is behind the other patches - ax.add_patch(outline) - for patch in motor_patches: - ax.add_patch(patch) + if isinstance(self.rocket.motor, ClusterMotor): + + for patch in motor_patches: + ax.add_patch(patch) + + if self.rocket.motor.positions: - self._draw_nozzle_tube(last_radius, last_x, nozzle_position, ax, vis_args) + cluster_z_positions = [ + pos[2] for pos in self.rocket.motor.positions + ] + z_connector = ( + min(cluster_z_positions) + if self.rocket._csys == 1 + else max(cluster_z_positions) + ) + self._draw_nozzle_tube( + last_radius, last_x, z_connector, ax, vis_args + ) + + else: + + total_csys = self.rocket._csys * self.rocket.motor._csys + nozzle_position = ( + self.rocket.motor_position + + self.rocket.motor.nozzle_position * total_csys + ) + + nozzle = self.rocket.motor.plots._generate_nozzle( + translate=(nozzle_position, 0), csys=self.rocket._csys + ) + motor_patches += [nozzle] - def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-argument + outline = self.rocket.motor.plots._generate_motor_region( + list_of_patches=motor_patches + ) + + ax.add_patch(outline) + for patch in motor_patches: + ax.add_patch(patch) + + self._draw_nozzle_tube( + last_radius, last_x, nozzle_position, ax, vis_args + ) + + def _generate_motor_patches(self, ax, plane="xz"): """Generates motor patches for drawing""" motor_patches = [] + total_csys = self.rocket._csys * self.rocket.motor._csys - if isinstance(self.rocket.motor, SolidMotor): + if isinstance(self.rocket.motor, ClusterMotor): + cluster = self.rocket.motor + all_sub_patches = [] + + for sub_motor, sub_pos in zip(cluster.motors, cluster.positions): + + motor_longitudinal_pos = sub_pos[2] + + offset = 0 + if plane == "xz": + offset = sub_pos[0] + elif plane == "yz": + offset = sub_pos[1] + + sub_total_csys = self.rocket._csys * sub_motor._csys + + if isinstance(sub_motor, SolidMotor): + current_motor_patches = [] + + grains_cm_pos = ( + motor_longitudinal_pos + + sub_motor.grains_center_of_mass_position * sub_total_csys + ) + ax.scatter( + grains_cm_pos.real, + offset, + color="brown", + label=f"Grains CM ({sub_pos[0]:.2f}, {sub_pos[1]:.2f})", + s=8, + zorder=10, + ) + + chamber = sub_motor.plots._generate_combustion_chamber( + translate=(grains_cm_pos.real, offset), label=None + ) + grains = sub_motor.plots._generate_grains( + translate=(grains_cm_pos.real, offset) + ) + + nozzle_position = ( + motor_longitudinal_pos + + sub_motor.nozzle_position * sub_total_csys + ) + + nozzle = sub_motor.plots._generate_nozzle( + translate=(nozzle_position.real, offset) + ) + + current_motor_patches += [chamber, *grains, nozzle] + all_sub_patches.extend(current_motor_patches) + + if all_sub_patches: + + outline = self.rocket.motor.motors[0].plots._generate_motor_region( + list_of_patches=all_sub_patches + ) + motor_patches.append(outline) # Mettre l'outline en premier + motor_patches.extend(all_sub_patches) + + elif isinstance(self.rocket.motor, SolidMotor): grains_cm_position = ( self.rocket.motor_position + self.rocket.motor.grains_center_of_mass_position * total_csys ) ax.scatter( - grains_cm_position, + grains_cm_position.real, 0, color="brown", label="Grains Center of Mass", @@ -457,10 +529,10 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg ) chamber = self.rocket.motor.plots._generate_combustion_chamber( - translate=(grains_cm_position, 0), label=None + translate=(grains_cm_position.real, 0), label=None ) grains = self.rocket.motor.plots._generate_grains( - translate=(grains_cm_position, 0) + translate=(grains_cm_position.real, 0) ) motor_patches += [chamber, *grains] @@ -471,7 +543,7 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg + self.rocket.motor.grains_center_of_mass_position * total_csys ) ax.scatter( - grains_cm_position, + grains_cm_position.real, 0, color="brown", label="Grains Center of Mass", @@ -483,16 +555,16 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg translate=(self.rocket.motor_position, 0), csys=total_csys ) chamber = self.rocket.motor.plots._generate_combustion_chamber( - translate=(grains_cm_position, 0), label=None + translate=(grains_cm_position.real, 0), label=None ) grains = self.rocket.motor.plots._generate_grains( - translate=(grains_cm_position, 0) + translate=(grains_cm_position.real, 0) ) motor_patches += [chamber, *grains] for tank, center in tanks_and_centers: ax.scatter( - center[0], - center[1], + center[0].real, + center[1].real, color="black", alpha=0.2, s=5, @@ -506,8 +578,8 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg ) for tank, center in tanks_and_centers: ax.scatter( - center[0], - center[1], + center[0].real, + center[1].real, color="black", alpha=0.2, s=4, @@ -519,8 +591,6 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args): """Draws the tube from the last surface to the nozzle position.""" - # Check if nozzle is beyond the last surface, if so draw a tube - # to it, with the radius of the last surface if self.rocket._csys == 1: if nozzle_position < last_x: x_tube = [last_x, nozzle_position] @@ -539,7 +609,7 @@ def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args): color=vis_args["body"], linewidth=vis_args["line_width"], ) - else: # if self.rocket._csys == -1: + else: if nozzle_position > last_x: x_tube = [last_x, nozzle_position] y_tube = [last_radius, last_radius] @@ -575,13 +645,16 @@ def _draw_rail_buttons(self, ax, vis_args): def _draw_center_of_mass_and_pressure(self, ax): """Draws the center of mass and center of pressure of the rocket.""" - # Draw center of mass and center of pressure - cm = self.rocket.center_of_mass(0) - ax.scatter(cm, 0, color="#1565c0", label="Center of Mass", s=10) - cp = self.rocket.cp_position(0) + cm_vector = self.rocket.center_of_mass(0) + + cm_z = float(cm_vector.z.real) + + ax.scatter(cm_z, 0, color="#1565c0", label="Center of Mass", s=10) + + cp_z = float(self.rocket.cp_position(0).real) ax.scatter( - cp, 0, label="Static Center of Pressure", color="red", s=10, zorder=10 + cp_z, 0, label="Static Center of Pressure", color="red", s=10, zorder=10 ) def _draw_sensors(self, ax, sensors, plane): @@ -593,23 +666,22 @@ def _draw_sensors(self, ax, sensors, plane): sensor = sensor_pos[0] pos = sensor_pos[1] if plane == "xz": - # z position of the sensor is the x position in the plot + x_pos = pos[2] normal_x = sensor.normal_vector.z - # x position of the sensor is the y position in the plot + y_pos = pos[0] normal_y = sensor.normal_vector.x elif plane == "yz": - # z position of the sensor is the x position in the plot + x_pos = pos[2] normal_x = sensor.normal_vector.z - # y position of the sensor is the y position in the plot + y_pos = pos[1] normal_y = sensor.normal_vector.y - else: # pragma: no cover + else: raise ValueError("Plane must be 'xz' or 'yz'.") - # line length is 2/5 of the rocket radius line_length = self.rocket.radius / 2.5 ax.plot( @@ -644,34 +716,34 @@ def all(self): None """ - # Rocket draw + if len(self.rocket.aerodynamic_surfaces) > 0: print("\nRocket Draw") print("-" * 40) self.draw() - # Mass Plots + print("\nMass Plots") print("-" * 40) self.total_mass() self.reduced_mass() - # Aerodynamics Plots + print("\nAerodynamics Plots") print("-" * 40) # Drag Plots print("Drag Plots") - print("-" * 20) # Separator for Drag Plots + print("-" * 20) self.drag_curves() # Stability Plots print("\nStability Plots") - print("-" * 20) # Separator for Stability Plots + print("-" * 20) self.static_margin() self.stability_margin() - # Thrust-to-Weight Plot + print("\nThrust-to-Weight Plot") print("-" * 40) self.thrust_to_weight() diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 1112a98f3..1fe77570c 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1,6 +1,5 @@ import math import warnings -from typing import Iterable import numpy as np @@ -8,6 +7,7 @@ from rocketpy.mathutils.function import Function from rocketpy.mathutils.vector_matrix import Matrix, Vector from rocketpy.motors.empty_motor import EmptyMotor +from rocketpy.motors.ClusterMotor import ClusterMotor from rocketpy.plots.rocket_plots import _RocketPlots from rocketpy.prints.rocket_prints import _RocketPrints from rocketpy.rocket.aero_surface import ( @@ -24,9 +24,12 @@ from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute from rocketpy.tools import ( - deprecated, - find_obj_from_hash, - parallel_axis_theorem_from_com, + parallel_axis_theorem_I11, + parallel_axis_theorem_I22, + parallel_axis_theorem_I33, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I23, ) @@ -365,24 +368,48 @@ def __init__( # pylint: disable=too-many-statements self.reduced_mass = None self.total_mass = None self.dry_mass = None + self.propellant_initial_mass = 0 # Ajout pour la cohérence + + # Initialize plots and prints object + self.prints = _RocketPrints(self) + self.plots = _RocketPlots(self) + self._full_rocket_update() + + def _full_rocket_update(self): + """ + Centralized method to update all rocket properties after adding or modifying a motor. + """ + if hasattr(self.motor, "propellant_initial_mass"): + self.propellant_initial_mass = self.motor.propellant_initial_mass + else: + self.propellant_initial_mass = ( + self.motor.propellant_mass.get_value_opt(0) + if self.motor.propellant_mass + else 0 + ) - # calculate dynamic inertial quantities self.evaluate_dry_mass() - self.evaluate_structural_mass_ratio() self.evaluate_total_mass() self.evaluate_center_of_dry_mass() + + if hasattr(self.motor, "nozzle_position"): + self.nozzle_position = self.motor.nozzle_position + else: + self.nozzle_position = self.motor_position + + self.evaluate_nozzle_to_cdm() self.evaluate_center_of_mass() + self.evaluate_dry_inertias() + self.evaluate_inertias() self.evaluate_reduced_mass() self.evaluate_thrust_to_weight() - - # Evaluate stability (even though no aerodynamic surfaces are present yet) self.evaluate_center_of_pressure() + self.evaluate_surfaces_cp_to_cdm() self.evaluate_stability_margin() self.evaluate_static_margin() - - # Initialize plots and prints object - self.prints = _RocketPrints(self) - self.plots = _RocketPlots(self) + self.evaluate_com_to_cdm_function() + self.evaluate_nozzle_gyration_tensor() + self.evaluate_structural_mass_ratio() @property def nosecones(self): @@ -454,13 +481,23 @@ def evaluate_structural_mass_ratio(self): divided by total mass (Rocket + Motor + Propellant) (kg). """ try: - self.structural_mass_ratio = self.dry_mass / ( - self.dry_mass + self.motor.propellant_initial_mass + dry_mass = self.dry_mass if self.dry_mass is not None else self.mass + propellant_mass_initial = self.propellant_initial_mass + total_mass_initial = dry_mass + propellant_mass_initial + + if total_mass_initial <= 1e-9: + print( + "Warning: Total initial mass is near zero, setting structural mass ratio to 1.0" + ) + self.structural_mass_ratio = 1.0 + else: + self.structural_mass_ratio = dry_mass / total_mass_initial + except (AttributeError, TypeError) as e: + print( + f"Warning: Could not calculate structural mass ratio, setting to 1.0. Details: {e}" ) - except ZeroDivisionError as e: - raise ValueError( - "Total rocket mass (dry + propellant) cannot be zero" - ) from e + self.structural_mass_ratio = 1.0 + return self.structural_mass_ratio def evaluate_center_of_mass(self): @@ -475,10 +512,18 @@ def evaluate_center_of_mass(self): See :doc:`Positions and Coordinate Systems ` for more information. """ + + com_without_motor_vec = Function( + lambda t: Vector([0, 0, self.center_of_mass_without_motor]), + inputs="t", + outputs="Vector (m)", + ) + self.center_of_mass = ( - self.center_of_mass_without_motor * self.mass + com_without_motor_vec * self.mass + self.motor_center_of_mass_position * self.motor.total_mass ) / self.total_mass + self.center_of_mass.set_inputs("Time (s)") self.center_of_mass.set_outputs("Center of Mass Position (m)") self.center_of_mass.set_title( @@ -493,13 +538,19 @@ def evaluate_center_of_dry_mass(self): Returns ------- - self.center_of_dry_mass_position : int, float + self.center_of_dry_mass_position : Vector Rocket's center of dry mass position (with unloaded motor) """ + + com_without_motor_vec = Vector([0, 0, self.center_of_mass_without_motor]) + + # motor_center_of_dry_mass_position + motor_com_dry_vec = self.motor_center_of_dry_mass_position + self.center_of_dry_mass_position = ( - self.center_of_mass_without_motor * self.mass - + self.motor_center_of_dry_mass_position * self.motor.dry_mass + com_without_motor_vec * self.mass + motor_com_dry_vec * self.motor.dry_mass ) / self.dry_mass + return self.center_of_dry_mass_position def evaluate_reduced_mass(self): @@ -604,7 +655,7 @@ def __evaluate_single_surface_cp_to_cdm(self, surface, position): [ (position.x - self.cm_eccentricity_x) * self._csys - surface.cpx, (position.y - self.cm_eccentricity_y) - surface.cpy, - (position.z - self.center_of_dry_mass_position) * self._csys + (position.z - self.center_of_dry_mass_position.z) * self._csys - surface.cpz, ] ) @@ -613,21 +664,10 @@ def __evaluate_single_surface_cp_to_cdm(self, surface, position): def evaluate_stability_margin(self): """Calculates the stability margin of the rocket as a function of mach number and time. - - Returns - ------- - stability_margin : Function - Stability margin of the rocket, in calibers, as a function of mach - number and time. Stability margin is defined as the distance between - the center of pressure and the center of mass, divided by the - rocket's diameter. """ self.stability_margin.set_source( lambda mach, time: ( - ( - self.center_of_mass.get_value_opt(time) - - self.cp_position.get_value_opt(mach) - ) + (self.center_of_mass(time).z - self.cp_position.get_value_opt(mach)) / (2 * self.radius) ) * self._csys @@ -635,31 +675,21 @@ def evaluate_stability_margin(self): return self.stability_margin def evaluate_static_margin(self): - """Calculates the static margin of the rocket as a function of time. - - Returns - ------- - static_margin : Function - Static margin of the rocket, in calibers, as a function of time. - Static margin is defined as the distance between the center of - pressure and the center of mass, divided by the rocket's diameter. - """ - # Calculate static margin + """Calculates the static margin of the rocket as a function of time.""" self.static_margin.set_source( lambda time: ( - self.center_of_mass.get_value_opt(time) - - self.cp_position.get_value_opt(0) + self.center_of_mass(time).z - self.cp_position.get_value_opt(0) ) / (2 * self.radius) ) - # Change sign if coordinate system is upside down self.static_margin *= self._csys self.static_margin.set_inputs("Time (s)") self.static_margin.set_outputs("Static Margin (c)") self.static_margin.set_title("Static Margin") - self.static_margin.set_discrete( - lower=0, upper=self.motor.burn_out_time, samples=200 - ) + if hasattr(self.motor, "burn_out_time"): + self.static_margin.set_discrete( + lower=0, upper=self.motor.burn_out_time, samples=200 + ) return self.static_margin def evaluate_dry_inertias(self): @@ -667,76 +697,49 @@ def evaluate_dry_inertias(self): the rocket's center of dry mass. The inertias are saved and returned in units of kg*m². This does not consider propellant mass but does take into account the motor dry mass. - - Returns - ------- - self.dry_I_11 : float - Float value corresponding to rocket inertia tensor 11 - component, which corresponds to the inertia relative to the - e_1 axis, centered at the center of dry mass. - self.dry_I_22 : float - Float value corresponding to rocket inertia tensor 22 - component, which corresponds to the inertia relative to the - e_2 axis, centered at the center of dry mass. - self.dry_I_33 : float - Float value corresponding to rocket inertia tensor 33 - component, which corresponds to the inertia relative to the - e_3 axis, centered at the center of dry mass. - self.dry_I_12 : float - Float value corresponding to rocket inertia tensor 12 - component, which corresponds to the inertia relative to the - e_1 and e_2 axes, centered at the center of dry mass. - self.dry_I_13 : float - Float value corresponding to rocket inertia tensor 13 - component, which corresponds to the inertia relative to the - e_1 and e_3 axes, centered at the center of dry mass. - self.dry_I_23 : float - Float value corresponding to rocket inertia tensor 23 - component, which corresponds to the inertia relative to the - e_2 and e_3 axes, centered at the center of dry mass. - - Notes - ----- - #. The ``e_1`` and ``e_2`` directions are assumed to be the directions \ - perpendicular to the rocket axial direction. - #. The ``e_3`` direction is assumed to be the direction parallel to the \ - axis of symmetry of the rocket. - #. RocketPy follows the definition of the inertia tensor that includes \ - the minus sign for all products of inertia. - - See Also - -------- - `Inertia Tensor `_ """ # Get masses motor_dry_mass = self.motor.dry_mass mass = self.mass - # Compute axes distances (CDM: Center of Dry Mass) + center_of_mass_without_motor_vec = Vector( + [0, 0, self.center_of_mass_without_motor] + ) + + # Compute axes distances (CDM: Center of Dry Mass) using vector subtraction center_of_mass_without_motor_to_CDM = ( - self.center_of_mass_without_motor - self.center_of_dry_mass_position + center_of_mass_without_motor_vec - self.center_of_dry_mass_position ) motor_center_of_dry_mass_to_CDM = ( self.motor_center_of_dry_mass_position - self.center_of_dry_mass_position ) - # Compute dry inertias - self.dry_I_11 = parallel_axis_theorem_from_com( - self.I_11_without_motor, mass, center_of_mass_without_motor_to_CDM - ) + parallel_axis_theorem_from_com( - self.motor.dry_I_11, motor_dry_mass, motor_center_of_dry_mass_to_CDM - ) + d_rocket = center_of_mass_without_motor_to_CDM + d_motor = motor_center_of_dry_mass_to_CDM - self.dry_I_22 = parallel_axis_theorem_from_com( - self.I_22_without_motor, mass, center_of_mass_without_motor_to_CDM - ) + parallel_axis_theorem_from_com( - self.motor.dry_I_22, motor_dry_mass, motor_center_of_dry_mass_to_CDM - ) + self.dry_I_11 = parallel_axis_theorem_I11( + self.I_11_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I11(self.motor.dry_I_11, motor_dry_mass, d_motor) + + self.dry_I_22 = parallel_axis_theorem_I22( + self.I_22_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I22(self.motor.dry_I_22, motor_dry_mass, d_motor) + + self.dry_I_33 = parallel_axis_theorem_I33( + self.I_33_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I33(self.motor.dry_I_33, motor_dry_mass, d_motor) + + self.dry_I_12 = parallel_axis_theorem_I12( + self.I_12_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I12(self.motor.dry_I_12, motor_dry_mass, d_motor) + + self.dry_I_13 = parallel_axis_theorem_I13( + self.I_13_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I13(self.motor.dry_I_13, motor_dry_mass, d_motor) - self.dry_I_33 = self.I_33_without_motor + self.motor.dry_I_33 - self.dry_I_12 = self.I_12_without_motor + self.motor.dry_I_12 - self.dry_I_13 = self.I_13_without_motor + self.motor.dry_I_13 - self.dry_I_23 = self.I_23_without_motor + self.motor.dry_I_23 + self.dry_I_23 = parallel_axis_theorem_I23( + self.I_23_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I23(self.motor.dry_I_23, motor_dry_mass, d_motor) return ( self.dry_I_11, @@ -751,58 +754,43 @@ def evaluate_inertias(self): """Calculates and returns the rocket's inertias relative to the rocket's center of dry mass. The inertias are saved and returned in units of kg*m². - - Returns - ------- - self.I_11 : float - Float value corresponding to rocket inertia tensor 11 - component, which corresponds to the inertia relative to the - e_1 axis, centered at the center of dry mass. - self.I_22 : float - Float value corresponding to rocket inertia tensor 22 - component, which corresponds to the inertia relative to the - e_2 axis, centered at the center of dry mass. - self.I_33 : float - Float value corresponding to rocket inertia tensor 33 - component, which corresponds to the inertia relative to the - e_3 axis, centered at the center of dry mass. - - Notes - ----- - #. The ``e_1`` and ``e_2`` directions are assumed to be the directions \ - perpendicular to the rocket axial direction. - #. The ``e_3`` direction is assumed to be the direction parallel to the \ - axis of symmetry of the rocket. - #. RocketPy follows the definition of the inertia tensor that includes \ - the minus sign for all products of inertia. - - See Also - -------- - `Inertia Tensor `_ """ - # Get masses - prop_mass = self.motor.propellant_mass # Propellant mass as a function of time - # Compute axes distances - CDM_to_CPM = ( - self.center_of_dry_mass_position - self.center_of_propellant_position + prop_mass = self.motor.propellant_mass + + cdm_position_func = Function( + lambda t: self.center_of_dry_mass_position, inputs="t", outputs="Vector (m)" ) - # Compute inertias - self.I_11 = self.dry_I_11 + parallel_axis_theorem_from_com( + CDM_to_CPM = cdm_position_func - self.center_of_propellant_position + + propellant_inertia_term_11 = parallel_axis_theorem_I11( self.motor.propellant_I_11, prop_mass, CDM_to_CPM ) - - self.I_22 = self.dry_I_22 + parallel_axis_theorem_from_com( + propellant_inertia_term_22 = parallel_axis_theorem_I22( self.motor.propellant_I_22, prop_mass, CDM_to_CPM ) + propellant_inertia_term_33 = parallel_axis_theorem_I33( + self.motor.propellant_I_33, prop_mass, CDM_to_CPM + ) + propellant_inertia_term_12 = parallel_axis_theorem_I12( + self.motor.propellant_I_12, prop_mass, CDM_to_CPM + ) + propellant_inertia_term_13 = parallel_axis_theorem_I13( + self.motor.propellant_I_13, prop_mass, CDM_to_CPM + ) + propellant_inertia_term_23 = parallel_axis_theorem_I23( + self.motor.propellant_I_23, prop_mass, CDM_to_CPM + ) + + self.I_11 = Function(lambda t: self.dry_I_11) + propellant_inertia_term_11 + self.I_22 = Function(lambda t: self.dry_I_22) + propellant_inertia_term_22 + self.I_33 = Function(lambda t: self.dry_I_33) + propellant_inertia_term_33 - self.I_33 = self.dry_I_33 + self.motor.propellant_I_33 - self.I_12 = self.dry_I_12 + self.motor.propellant_I_12 - self.I_13 = self.dry_I_13 + self.motor.propellant_I_13 - self.I_23 = self.dry_I_23 + self.motor.propellant_I_23 + self.I_12 = Function(lambda t: self.dry_I_12) + propellant_inertia_term_12 + self.I_13 = Function(lambda t: self.dry_I_13) + propellant_inertia_term_13 + self.I_23 = Function(lambda t: self.dry_I_23) + propellant_inertia_term_23 - # Return inertias return ( self.I_11, self.I_22, @@ -815,15 +803,9 @@ def evaluate_inertias(self): def evaluate_nozzle_to_cdm(self): """Evaluates the distance between the nozzle exit and the rocket's center of dry mass. - - Returns - ------- - self.nozzle_to_cdm : float - Distance between the nozzle exit and the rocket's center of dry - mass position, in meters. """ self.nozzle_to_cdm = ( - -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys + -(self.nozzle_position - self.center_of_dry_mass_position.z) * self._csys ) return self.nozzle_to_cdm @@ -852,29 +834,21 @@ def evaluate_nozzle_gyration_tensor(self): def evaluate_com_to_cdm_function(self): """Evaluates the z-coordinate of the center of mass (COM) relative to the center of dry mass (CDM). + """ - Notes - ----- - 1. The `com_to_cdm_function` plus `center_of_mass` should be equal - to `center_of_dry_mass_position` at every time step. - 2. The `com_to_cdm_function` is a function of time and will usually - already be discretized. + cpm_z_func = Function( + lambda t: self.center_of_propellant_position(t).z, inputs="t" + ) + + cdm_z_scalar = self.center_of_dry_mass_position.z - Returns - ------- - self.com_to_cdm_function : Function - Function of time expressing the z-coordinate of the center of mass - relative to the center of dry mass. - """ self.com_to_cdm_function = ( -1 - * ( - (self.center_of_propellant_position - self.center_of_dry_mass_position) - * self._csys - ) + * ((cpm_z_func - cdm_z_scalar) * self._csys) * self.motor.propellant_mass / self.total_mass ) + self.com_to_cdm_function.set_inputs("Time (s)") self.com_to_cdm_function.set_outputs("Z Coordinate COM to CDM (m)") self.com_to_cdm_function.set_title("Z Coordinate COM to CDM") @@ -940,62 +914,79 @@ def get_inertia_tensor_derivative_at_time(self, t): ] ) - def add_motor(self, motor, position): # pylint: disable=too-many-statements - """Adds a motor to the rocket. + def add_motor(self, motor, position=0): + """ + Adds a motor or a motor cluster to the rocket. This method is universal + and handles standard motors, EmptyMotor, and ClusterMotor objects. + """ + if hasattr(self, "motor") and not isinstance(self.motor, EmptyMotor): + print("A motor is already attached. Overwriting previous motor.") - Parameters - ---------- - motor : Motor, SolidMotor, HybridMotor, LiquidMotor, GenericMotor - Motor to be added to the rocket. - position : int, float - Position, in meters, of the motor's coordinate system origin - relative to the user defined rocket coordinate system. + self.motor = motor + self.motor_position = position - See Also - -------- - :ref:`addsurface` + if isinstance(motor, EmptyMotor): + self.center_of_propellant_position = Function( + lambda t: Vector([0, 0, position]) + ) + self.motor_center_of_mass_position = Function( + lambda t: Vector([0, 0, position]) + ) + self.motor_center_of_dry_mass_position = Vector([0, 0, position]) + self.nozzle_position = position + self.total_mass_flow_rate = Function(0) + return - Returns - ------- - None + _ = self._csys * getattr(motor, "_csys", 1) + + is_cluster = hasattr(motor, "get_total_thrust_vector") + + if is_cluster: + + self.center_of_propellant_position = motor.center_of_propellant_mass + self.motor_center_of_mass_position = motor.center_of_mass + self.motor_center_of_dry_mass_position = motor.center_of_dry_mass_position + else: + + self.motor_center_of_mass_position = Function( + lambda t: Vector([0, 0, motor.center_of_mass(t) * _ + position]), + inputs="t", + outputs="Vector (m)", + ) + self.center_of_propellant_position = Function( + lambda t: Vector( + [0, 0, motor.center_of_propellant_mass(t) * _ + position] + ), + inputs="t", + outputs="Vector (m)", + ) + self.motor_center_of_dry_mass_position = Vector( + [0, 0, motor.center_of_dry_mass_position * _ + position] + ) + + self.nozzle_position = motor.nozzle_position * _ + position + self.total_mass_flow_rate = motor.total_mass_flow_rate + + self._full_rocket_update() + + def add_cluster_motor(self, motors, positions, orientations=None): """ - if hasattr(self, "motor"): - # pylint: disable=access-member-before-definition - if not isinstance(self.motor, EmptyMotor): - print( - "Only one motor per rocket is currently supported. " - + "Overwriting previous motor." - ) - self.motor = motor - self.motor_position = position - _ = self._csys * self.motor._csys - self.center_of_propellant_position = ( - self.motor.center_of_propellant_mass * _ + self.motor_position - ) - self.motor_center_of_mass_position = ( - self.motor.center_of_mass * _ + self.motor_position - ) - self.motor_center_of_dry_mass_position = ( - self.motor.center_of_dry_mass_position * _ + self.motor_position - ) - self.nozzle_position = self.motor.nozzle_position * _ + self.motor_position + Creates a ClusterMotor and adds it to the rocket. + """ + + cluster = ClusterMotor(motors, positions, orientations) + + self.motor = cluster + self.motor_position = 0 # Positions are handled within the cluster + + self.center_of_propellant_position = self.motor.center_of_propellant_mass + self.motor_center_of_mass_position = self.motor.center_of_mass + self.motor_center_of_dry_mass_position = self.motor.center_of_dry_mass_position self.total_mass_flow_rate = self.motor.total_mass_flow_rate - self.evaluate_dry_mass() - self.evaluate_structural_mass_ratio() - self.evaluate_total_mass() - self.evaluate_center_of_dry_mass() - self.evaluate_nozzle_to_cdm() - self.evaluate_center_of_mass() - self.evaluate_dry_inertias() - self.evaluate_inertias() - self.evaluate_reduced_mass() - self.evaluate_thrust_to_weight() - self.evaluate_center_of_pressure() - self.evaluate_surfaces_cp_to_cdm() - self.evaluate_stability_margin() - self.evaluate_static_margin() - self.evaluate_com_to_cdm_function() - self.evaluate_nozzle_gyration_tensor() + + # Call the centralized update method + self._full_rocket_update() + return cluster def __add_single_surface(self, surface, position): """Adds a single aerodynamic surface to the rocket. Makes checks for @@ -1178,16 +1169,16 @@ def add_nose( self.add_surfaces(nose, position) return nose - @deprecated( - reason="This method is set to be deprecated in version 1.0.0 and fully " - "removed by version 2.0.0", - alternative="Rocket.add_trapezoidal_fins", - ) def add_fins(self, *args, **kwargs): # pragma: no cover """See Rocket.add_trapezoidal_fins for documentation. This method is set to be deprecated in version 1.0.0 and fully removed by version 2.0.0. Use Rocket.add_trapezoidal_fins instead. It keeps the same arguments and signature.""" + warnings.warn( + "This method is set to be deprecated in version 1.0.0 and fully " + "removed by version 2.0.0. Use Rocket.add_trapezoidal_fins instead", + DeprecationWarning, + ) return self.add_trapezoidal_fins(*args, **kwargs) def add_trapezoidal_fins( @@ -1439,16 +1430,7 @@ def add_free_form_fins( return fin_set def add_parachute( - self, - name, - cd_s, - trigger, - sampling_rate=100, - lag=0, - noise=(0, 0, 0), - radius=1.5, - height=None, - porosity=0.0432, + self, name, cd_s, trigger, sampling_rate=100, lag=0, noise=(0, 0, 0) ): """Creates a new parachute, storing its parameters such as opening delay, drag coefficients and trigger function. @@ -1507,39 +1489,16 @@ def add_parachute( The values are used to add noise to the pressure signal which is passed to the trigger function. Default value is (0, 0, 0). Units are in pascal. - radius : float, optional - Length of the non-unique semi-axis (radius) of the inflated hemispheroid - parachute. Default value is 1.5. - Units are in meters. - height : float, optional - Length of the unique semi-axis (height) of the inflated hemispheroid - parachute. Default value is the radius of the parachute. - Units are in meters. - porosity : float, optional - Geometric porosity of the canopy (ratio of open area to total canopy area), - in [0, 1]. Affects only the added-mass scaling during descent; it does - not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass - of 1.0 (“neutral” behavior). Returns ------- parachute : Parachute - Parachute containing trigger, sampling_rate, lag, cd_s, noise, radius, - height, porosity and name. Furthermore, it stores clean_pressure_signal, + Parachute containing trigger, sampling_rate, lag, cd_s, noise + and name. Furthermore, it stores clean_pressure_signal, noise_signal and noisyPressureSignal which are filled in during Flight simulation. """ - parachute = Parachute( - name, - cd_s, - trigger, - sampling_rate, - lag, - noise, - radius, - height, - porosity, - ) + parachute = Parachute(name, cd_s, trigger, sampling_rate, lag, noise) self.parachutes.append(parachute) return self.parachutes[-1] @@ -1929,16 +1888,7 @@ def all_info(self): self.info() self.plots.all() - # pylint: disable=too-many-statements - def to_dict(self, **kwargs): - discretize = kwargs.get("discretize", False) - - power_off_drag = self.power_off_drag - power_on_drag = self.power_on_drag - if discretize: - power_off_drag = power_off_drag.set_discrete(0, 4, 50, mutate_self=False) - power_on_drag = power_on_drag.set_discrete(0, 4, 50, mutate_self=False) - + def to_dict(self, include_outputs=False): rocket_dict = { "radius": self.radius, "mass": self.mass, @@ -1948,8 +1898,8 @@ def to_dict(self, **kwargs): "I_12_without_motor": self.I_12_without_motor, "I_13_without_motor": self.I_13_without_motor, "I_23_without_motor": self.I_23_without_motor, - "power_off_drag": power_off_drag, - "power_on_drag": power_on_drag, + "power_off_drag": self.power_off_drag, + "power_on_drag": self.power_on_drag, "center_of_mass_without_motor": self.center_of_mass_without_motor, "coordinate_system_orientation": self.coordinate_system_orientation, "motor": self.motor, @@ -1962,51 +1912,7 @@ def to_dict(self, **kwargs): "sensors": self.sensors, } - if kwargs.get("include_outputs", False): - thrust_to_weight = self.thrust_to_weight - cp_position = self.cp_position - stability_margin = self.stability_margin - center_of_mass = self.center_of_mass - motor_center_of_mass_position = self.motor_center_of_mass_position - reduced_mass = self.reduced_mass - total_mass = self.total_mass - total_mass_flow_rate = self.total_mass_flow_rate - center_of_propellant_position = self.center_of_propellant_position - - if discretize: - thrust_to_weight = thrust_to_weight.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - cp_position = cp_position.set_discrete(0, 4, 25, mutate_self=False) - stability_margin = stability_margin.set_discrete( - (0, self.motor.burn_time[0]), - (2, self.motor.burn_time[1]), - (10, 10), - mutate_self=False, - ) - center_of_mass = center_of_mass.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - motor_center_of_mass_position = ( - motor_center_of_mass_position.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - ) - reduced_mass = reduced_mass.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - total_mass = total_mass.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - total_mass_flow_rate = total_mass_flow_rate.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - center_of_propellant_position = ( - center_of_propellant_position.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - ) - + if include_outputs: rocket_dict["area"] = self.area rocket_dict["center_of_dry_mass_position"] = ( self.center_of_dry_mass_position @@ -2014,26 +1920,30 @@ def to_dict(self, **kwargs): rocket_dict["center_of_mass_without_motor"] = ( self.center_of_mass_without_motor ) - rocket_dict["motor_center_of_mass_position"] = motor_center_of_mass_position + rocket_dict["motor_center_of_mass_position"] = ( + self.motor_center_of_mass_position + ) rocket_dict["motor_center_of_dry_mass_position"] = ( self.motor_center_of_dry_mass_position ) - rocket_dict["center_of_mass"] = center_of_mass - rocket_dict["reduced_mass"] = reduced_mass - rocket_dict["total_mass"] = total_mass - rocket_dict["total_mass_flow_rate"] = total_mass_flow_rate - rocket_dict["thrust_to_weight"] = thrust_to_weight + rocket_dict["center_of_mass"] = self.center_of_mass + rocket_dict["reduced_mass"] = self.reduced_mass + rocket_dict["total_mass"] = self.total_mass + rocket_dict["total_mass_flow_rate"] = self.total_mass_flow_rate + rocket_dict["thrust_to_weight"] = self.thrust_to_weight rocket_dict["cp_eccentricity_x"] = self.cp_eccentricity_x rocket_dict["cp_eccentricity_y"] = self.cp_eccentricity_y rocket_dict["thrust_eccentricity_x"] = self.thrust_eccentricity_x rocket_dict["thrust_eccentricity_y"] = self.thrust_eccentricity_y - rocket_dict["cp_position"] = cp_position - rocket_dict["stability_margin"] = stability_margin + rocket_dict["cp_position"] = self.cp_position + rocket_dict["stability_margin"] = self.stability_margin rocket_dict["static_margin"] = self.static_margin rocket_dict["nozzle_position"] = self.nozzle_position rocket_dict["nozzle_to_cdm"] = self.nozzle_to_cdm rocket_dict["nozzle_gyration_tensor"] = self.nozzle_gyration_tensor - rocket_dict["center_of_propellant_position"] = center_of_propellant_position + rocket_dict["center_of_propellant_position"] = ( + self.center_of_propellant_position + ) return rocket_dict @@ -2076,29 +1986,17 @@ def from_dict(cls, data): for parachute in data["parachutes"]: rocket.parachutes.append(parachute) - for sensor, position in data["sensors"]: - rocket.add_sensor(sensor, position) - - for air_brake in data["air_brakes"]: - rocket.air_brakes.append(air_brake) - - for controller in data["_controllers"]: - interactive_objects_hash = getattr(controller, "_interactive_objects_hash") - if interactive_objects_hash is not None: - is_iterable = isinstance(interactive_objects_hash, Iterable) - if not is_iterable: - interactive_objects_hash = [interactive_objects_hash] - for hash_ in interactive_objects_hash: - if (hashed_obj := find_obj_from_hash(data, hash_)) is not None: - if not is_iterable: - controller.interactive_objects = hashed_obj - else: - controller.interactive_objects.append(hashed_obj) - else: - warnings.warn( - "Could not find controller interactive objects." - "Deserialization will proceed, results may not be accurate." - ) - rocket._add_controllers(controller) + for air_brakes in data["air_brakes"]: + rocket.add_air_brakes( + drag_coefficient_curve=air_brakes["drag_coefficient_curve"], + controller_function=air_brakes["controller_function"], + sampling_rate=air_brakes["sampling_rate"], + clamp=air_brakes["clamp"], + reference_area=air_brakes["reference_area"], + initial_observed_variables=air_brakes["initial_observed_variables"], + override_rocket_drag=air_brakes["override_rocket_drag"], + name=air_brakes["name"], + controller_name=air_brakes["controller_name"], + ) return rocket diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a38be7d93..e7a02f57a 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1,21 +1,18 @@ # pylint: disable=too-many-lines +import json import math import warnings from copy import deepcopy from functools import cached_property - import numpy as np +import simplekml from scipy.integrate import BDF, DOP853, LSODA, RK23, RK45, OdeSolver, Radau - -from rocketpy.simulation.flight_data_exporter import FlightDataExporter - from ..mathutils.function import Function, funcify_method from ..mathutils.vector_matrix import Matrix, Vector from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints from ..tools import ( calculate_cubic_hermite_coefficients, - deprecated, euler313_to_quaternions, find_closest, find_root_linear_interpolation, @@ -144,7 +141,7 @@ class Flight: Flight.heading : int, float Launch heading angle relative to north given in degrees. Flight.initial_solution : list - List defines initial condition - [t_initial, x_init, + List defines initial condition - [tInit, x_init, y_init, z_init, vx_init, vy_init, vz_init, e0_init, e1_init, e2_init, e3_init, w1_init, w2_init, w3_init] Flight.t_initial : int, float @@ -156,6 +153,12 @@ class Flight: Current integration time. Flight.y : list Current integration state vector u. + Flight.post_processed : bool + Defines if solution data has been post processed. + Flight.initial_solution : list + List defines initial condition - [tInit, x_init, + y_init, z_init, vx_init, vy_init, vz_init, e0_init, e1_init, + e2_init, e3_init, w1_init, w2_init, w3_init] Flight.out_of_rail_time : int, float Time, in seconds, in which the rocket completely leaves the rail. @@ -541,7 +544,7 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements rtol : float, array, optional Maximum relative error tolerance to be tolerated in the integration scheme. Can be given as array for each - state space variable. Default is 1e-6. + state space variable. Default is 1e-3. atol : float, optional Maximum absolute error tolerance to be tolerated in the integration scheme. Can be given as array for each @@ -762,22 +765,7 @@ def __simulate(self, verbose): callbacks = [ lambda self, parachute_cd_s=parachute.cd_s: setattr( self, "parachute_cd_s", parachute_cd_s - ), - lambda self, parachute_radius=parachute.radius: setattr( - self, "parachute_radius", parachute_radius - ), - lambda self, parachute_height=parachute.height: setattr( - self, "parachute_height", parachute_height - ), - lambda self, parachute_porosity=parachute.porosity: setattr( - self, "parachute_porosity", parachute_porosity - ), - lambda self, - added_mass_coefficient=parachute.added_mass_coefficient: setattr( - self, - "parachute_added_mass_coefficient", - added_mass_coefficient, - ), + ) ] self.flight_phases.add_phase( node.t + parachute.lag, @@ -1020,34 +1008,9 @@ def __simulate(self, verbose): i += 1 # Create flight phase for time after inflation callbacks = [ - lambda self, - parachute_cd_s=parachute.cd_s: setattr( + lambda self, parachute_cd_s=parachute.cd_s: setattr( self, "parachute_cd_s", parachute_cd_s - ), - lambda self, - parachute_radius=parachute.radius: setattr( - self, - "parachute_radius", - parachute_radius, - ), - lambda self, - parachute_height=parachute.height: setattr( - self, - "parachute_height", - parachute_height, - ), - lambda self, - parachute_porosity=parachute.porosity: setattr( - self, - "parachute_porosity", - parachute_porosity, - ), - lambda self, - added_mass_coefficient=parachute.added_mass_coefficient: setattr( - self, - "parachute_added_mass_coefficient", - added_mass_coefficient, - ), + ) ] self.flight_phases.add_phase( overshootable_node.t + parachute.lag, @@ -1130,13 +1093,13 @@ def __init_solution_monitors(self): self.out_of_rail_time_index = 0 self.out_of_rail_state = np.array([0]) self.apogee_state = np.array([0]) - self.apogee = 0 self.apogee_time = 0 self.x_impact = 0 self.y_impact = 0 self.impact_velocity = 0 self.impact_state = np.array([0]) self.parachute_events = [] + self.post_processed = False self.__post_processed_variables = [] def __init_flight_state(self): @@ -1205,7 +1168,6 @@ def __init_flight_state(self): self.out_of_rail_state = self.initial_solution[1:] self.out_of_rail_time = self.initial_solution[0] self.out_of_rail_time_index = 0 - self.t_initial = self.initial_solution[0] self.initial_derivative = self.u_dot_generalized if self._controllers or self.sensors: # Handle post process during simulation, get initial accel/forces @@ -1439,7 +1401,9 @@ def udot_rail2(self, t, u, post_processing=False): # pragma: no cover # Hey! We will finish this function later, now we just can use u_dot return self.u_dot_generalized(t, u, post_processing=post_processing) - def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals,too-many-statements + def u_dot( + self, t, u, post_processing=False + ): # pylint: disable=too-many-locals,too-many-statements """Calculates derivative of u state vector with respect to time when rocket is flying in 6 DOF motion during ascent out of rail and descent without parachute. @@ -1714,12 +1678,6 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals ax, ay, az = K @ Vector(L) az -= self.env.gravity.get_value_opt(z) # Include gravity - # Coriolis acceleration - _, w_earth_y, w_earth_z = self.env.earth_rotation_vector - ax -= 2 * (vz * w_earth_y - vy * w_earth_z) - ay -= 2 * (vx * w_earth_z) - az -= 2 * (-vx * w_earth_y) - # Create u_dot u_dot = [ vx, @@ -1759,90 +1717,64 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals return u_dot - def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too-many-locals,too-many-statements - """Calculates derivative of u state vector with respect to time when the - rocket is flying in 6 DOF motion in space and significant mass variation - effects exist. Typical flight phases include powered ascent after launch - rail. - - Parameters - ---------- - t : float - Time in seconds - u : list - State vector defined by u = [x, y, z, vx, vy, vz, q0, q1, - q2, q3, omega1, omega2, omega3]. - post_processing : bool, optional - If True, adds flight data information directly to self variables - such as self.angle_of_attack, by default False. - - Returns - ------- - u_dot : list - State vector defined by u_dot = [vx, vy, vz, ax, ay, az, - e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. + def u_dot_generalized(self, t, u, post_processing=False): """ - # Retrieve integration data - _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u + Calculates the derivative of the state vector u (6-DOF) using a + generalized formalism suitable for motor clusters and including all + physical forces. + """ + # State variable extraction + _, _, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3 = u - # Create necessary vectors - # r = Vector([x, y, z]) # CDM position vector - v = Vector([vx, vy, vz]) # CDM velocity vector - e = [e0, e1, e2, e3] # Euler parameters/quaternions - w = Vector([omega1, omega2, omega3]) # Angular velocity vector + # Create vectors for calculations + velocity_vector = Vector([vx, vy, vz]) + euler_params = [e0, e1, e2, e3] + angular_velocity_vector = Vector([w1, w2, w3]) + w = angular_velocity_vector - # Retrieve necessary quantities - ## Rocket mass + # Rocket properties at time t total_mass = self.rocket.total_mass.get_value_opt(t) - total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) - total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t) - ## CM position vector and time derivatives relative to CDM in body frame - r_CM_z = self.rocket.com_to_cdm_function - r_CM_t = r_CM_z.get_value_opt(t) - r_CM = Vector([0, 0, r_CM_t]) - r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) - r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)]) - ## Nozzle position vector - r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm]) - ## Nozzle gyration tensor - S_nozzle = self.rocket.nozzle_gyration_tensor - ## Inertia tensor - inertia_tensor = self.rocket.get_inertia_tensor_at_time(t) - ## Inertia tensor time derivative in the body frame - I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) - - # Calculate the Inertia tensor relative to CM - H = (r_CM.cross_matrix @ -r_CM.cross_matrix) * total_mass - I_CM = inertia_tensor - H - # Prepare transformation matrices - K = Matrix.transformation(e) + # Transformation matrix from body frame to inertial frame + K = Matrix.transformation(euler_params) Kt = K.transpose - # Compute aerodynamic forces and moments + # --- Forces and Moments Calculation --- + + # Gravity in body frame + gravity_force_inertial = Vector( + [0, 0, -total_mass * self.env.gravity.get_value_opt(z)] + ) + gravity_force_body = Kt @ gravity_force_inertial + + # --- START OF INTEGRATED AERODYNAMIC CALCULATION BLOCK --- R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0 - ## Drag force - rho = self.env.density.get_value_opt(z) + # Get freestream speed wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) - wind_velocity = Vector([wind_velocity_x, wind_velocity_y, 0]) - free_stream_speed = abs((wind_velocity - Vector(v))) speed_of_sound = self.env.speed_of_sound.get_value_opt(z) + rho = self.env.density.get_value_opt(z) + + # Relative wind speed in inertial frame + stream_velocity_inertial = Vector( + [wind_velocity_x - vx, wind_velocity_y - vy, -vz] + ) + free_stream_speed = abs(stream_velocity_inertial) free_stream_mach = free_stream_speed / speed_of_sound - if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time: - pressure = self.env.pressure.get_value_opt(z) - net_thrust = max( - self.rocket.motor.thrust.get_value_opt(t) - + self.rocket.motor.pressure_thrust(pressure), - 0, - ) + # Determine the base drag of the rocket + if t < self.rocket.motor.burn_out_time: drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach) else: - net_thrust = 0 drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach) - R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff + + R3_base_drag = ( + -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff + ) + R3 += R3_base_drag + + # Add air brakes drag if active for air_brakes in self.rocket.air_brakes: if air_brakes.deployment_level > 0: air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( @@ -1856,36 +1788,40 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too * air_brakes_cd ) if air_brakes.override_rocket_drag: - R3 = air_brakes_force # Substitutes rocket drag coefficient + R3 = air_brakes_force else: R3 += air_brakes_force - # Get rocket velocity in body frame - velocity_in_body_frame = Kt @ v - # Calculate lift and moment for each component of the rocket + + # Moment due to CoP eccentricity + M1 += self.rocket.cp_eccentricity_y * R3 + M2 -= self.rocket.cp_eccentricity_x * R3 + + # Rocket velocity in body frame + velocity_body_frame = Kt @ velocity_vector + + # Calculate lift and moment for each aerodynamic surface for aero_surface, _ in self.rocket.aerodynamic_surfaces: - # Component cp relative to CDM in body frame comp_cp = self.rocket.surfaces_cp_to_cdm[aero_surface] # Component absolute velocity in body frame - comp_vb = velocity_in_body_frame + (w ^ comp_cp) + comp_vb = velocity_body_frame + (w ^ comp_cp) # Wind velocity at component altitude comp_z = z + (K @ comp_cp).z comp_wind_vx = self.env.wind_velocity_x.get_value_opt(comp_z) comp_wind_vy = self.env.wind_velocity_y.get_value_opt(comp_z) - # Component freestream velocity in body frame + # Wind velocity in body frame comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0]) comp_stream_velocity = comp_wind_vb - comp_vb comp_stream_speed = abs(comp_stream_velocity) comp_stream_mach = comp_stream_speed / speed_of_sound - # Reynolds at component altitude - # TODO: Reynolds is only used in generic surfaces. This calculation - # should be moved to the surface class for efficiency + comp_reynolds = ( self.env.density.get_value_opt(comp_z) * comp_stream_speed * aero_surface.reference_length / self.env.dynamic_viscosity.get_value_opt(comp_z) ) - # Forces and moments + + # Calculate forces and moments for the surface X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( comp_stream_velocity, comp_stream_speed, @@ -1902,70 +1838,86 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too M2 += N M3 += L - # Off center moment - M1 += ( - self.rocket.cp_eccentricity_y * R3 - + self.rocket.thrust_eccentricity_y * net_thrust - ) - M2 -= ( - self.rocket.cp_eccentricity_x * R3 - + self.rocket.thrust_eccentricity_x * net_thrust - ) + # Moment due to CoP eccentricity M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1 - weight_in_body_frame = Kt @ Vector( - [0, 0, -total_mass * self.env.gravity.get_value_opt(z)] - ) + # Create final vectors for aero forces and moments + R_aero_body = Vector([R1, R2, R3]) + M_aero_body = Vector([M1, M2, M3]) + # --- END OF AERODYNAMIC CALCULATION BLOCK --- + + # Thrust and moment from cluster/motor + if hasattr(self.rocket.motor, "get_total_thrust_vector"): # It's a ClusterMotor + thrust_vector_body = self.rocket.motor.get_total_thrust_vector(t) + rocket_cm = self.rocket.center_of_mass.get_value_opt(t) + moment_vector_body = self.rocket.motor.get_total_moment(t, rocket_cm) + net_thrust_scalar = abs(thrust_vector_body) + else: # It's a standard motor + pressure = self.env.pressure.get_value_opt(z) + net_thrust_scalar = max( + self.rocket.motor.thrust.get_value_opt(t) + + self.rocket.motor.pressure_thrust(pressure), + 0, + ) + thrust_vector_body = Vector([0, 0, net_thrust_scalar]) # Axial thrust + moment_vector_body = Vector( + [ # Moment due to eccentricity + self.rocket.thrust_eccentricity_y * net_thrust_scalar, + -self.rocket.thrust_eccentricity_x * net_thrust_scalar, + 0, + ] + ) - T00 = total_mass * r_CM - T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot - T04 = ( - Vector([0, 0, net_thrust]) - - total_mass * r_CM_ddot - - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (r_NOZ - r_CM) - ) - T05 = total_mass_dot * S_nozzle - I_dot - - T20 = ( - ((w ^ T00) ^ w) - + (w ^ T03) - + T04 - + weight_in_body_frame - + Vector([R1, R2, R3]) - ) + # Total force and moment in body frame + total_force_body = gravity_force_body + thrust_vector_body + R_aero_body + total_moment_body = moment_vector_body + M_aero_body - T21 = ( - ((inertia_tensor @ w) ^ w) - + T05 @ w - - (weight_in_body_frame ^ r_CM) - + Vector([M1, M2, M3]) - ) + # --- Equations of Motion --- - # Angular velocity derivative - w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) + # Linear acceleration + linear_acceleration_body = total_force_body / total_mass + linear_acceleration_inertial = K @ linear_acceleration_body - # Euler parameters derivative - e_dot = [ - 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3), - 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3), - 0.5 * (omega2 * e0 - omega3 * e1 + omega1 * e3), - 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2), - ] + # Angular acceleration (Euler's equation for rigid body with variable mass) + inertia_tensor = self.rocket.get_inertia_tensor_at_time(t) + I_w = inertia_tensor @ angular_velocity_vector + w_cross_Iw = angular_velocity_vector.cross(I_w) + I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) - # Velocity vector derivative + Coriolis acceleration - w_earth = Vector(self.env.earth_rotation_vector) - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v) + angular_acceleration_body = inertia_tensor.inverse @ ( + total_moment_body - w_cross_Iw - (I_dot @ angular_velocity_vector) + ) - # Position vector derivative - r_dot = [vx, vy, vz] + # Euler parameters derivatives + e0dot = 0.5 * (-w1 * e1 - w2 * e2 - w3 * e3) + e1dot = 0.5 * (w1 * e0 + w3 * e2 - w2 * e3) + e2dot = 0.5 * (w2 * e0 - w3 * e1 + w1 * e3) + e3dot = 0.5 * (w3 * e0 + w2 * e1 - w1 * e2) - # Create u_dot - u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] + # Assemble the derivative vector + u_dot = [ + vx, + vy, + vz, + *linear_acceleration_inertial, + e0dot, + e1dot, + e2dot, + e3dot, + *angular_acceleration_body, + ] + # Post-processing (if necessary) if post_processing: self.__post_processed_variables.append( - [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3, net_thrust] + [ + t, + *linear_acceleration_inertial, + *angular_acceleration_body, + *R_aero_body, + *total_moment_body, + net_thrust_scalar, + ] ) return u_dot @@ -2001,9 +1953,15 @@ def u_dot_parachute(self, t, u, post_processing=False): wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) + # Get Parachute data + cd_s = self.parachute_cd_s + # Get the mass of the rocket mp = self.rocket.dry_mass + # Define constants + # ka = 1 # Added mass coefficient (depends on parachute's porosity) + # R = 1.5 # Parachute radius # to = 1.2 # eta = 1 # Rdot = (6 * R * (1 - eta) / (1.2**6)) * ( @@ -2011,17 +1969,8 @@ def u_dot_parachute(self, t, u, post_processing=False): # ) # Rdot = 0 - # tf = 8 * nominal diameter / velocity at line stretch - # Calculate added mass - ma = ( - self.parachute_added_mass_coefficient - * rho - * (2 / 3) - * np.pi - * self.parachute_radius**2 - * self.parachute_height - ) + ma = 0 # Calculate freestream speed freestream_x = vx - wind_velocity_x @@ -2030,20 +1979,14 @@ def u_dot_parachute(self, t, u, post_processing=False): free_stream_speed = (freestream_x**2 + freestream_y**2 + freestream_z**2) ** 0.5 # Determine drag force - pseudo_drag = -0.5 * rho * self.parachute_cd_s * free_stream_speed + pseudo_drag = -0.5 * rho * cd_s * free_stream_speed # pseudo_drag = pseudo_drag - ka * rho * 4 * np.pi * (R**2) * Rdot - Dx = pseudo_drag * freestream_x # add eta efficiency for wake + Dx = pseudo_drag * freestream_x Dy = pseudo_drag * freestream_y Dz = pseudo_drag * freestream_z ax = Dx / (mp + ma) ay = Dy / (mp + ma) - az = (Dz - mp * self.env.gravity.get_value_opt(z)) / (mp + ma) - - # Add coriolis acceleration - _, w_earth_y, w_earth_z = self.env.earth_rotation_vector - ax -= 2 * (vz * w_earth_y - vy * w_earth_z) - ay -= 2 * (vx * w_earth_z) - az -= 2 * (-vx * w_earth_y) + az = (Dz - self.env.gravity.get_value_opt(z) * mp) / (mp + ma) if post_processing: self.__post_processed_variables.append( @@ -3106,16 +3049,14 @@ def __calculate_rail_button_forces(self): # TODO: complex method. null_force = Function(0) if self.out_of_rail_time_index == 0: # No rail phase, no rail button forces warnings.warn( - "Trying to calculate rail button forces without a rail phase defined. " - + "The rail button forces will be set to zero.", - UserWarning, + "Trying to calculate rail button forces without a rail phase defined." + + "The rail button forces will be set to zero." ) return null_force, null_force, null_force, null_force if len(self.rocket.rail_buttons) == 0: warnings.warn( - "Trying to calculate rail button forces without rail buttons defined. " - + "The rail button forces will be set to zero.", - UserWarning, + "Trying to calculate rail button forces without rail buttons defined." + + "The rail button forces will be set to zero." ) return null_force, null_force, null_force, null_force @@ -3227,6 +3168,28 @@ def __evaluate_post_process(self): 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 + computed just in time. + + Post-process all Flight information produced during + simulation. Includes the calculation of maximum values, + calculation of secondary values such as energy and conversion + of lists to Function objects to facilitate plotting. + + Returns + ------- + None + """ + # pylint: disable=unused-argument + warnings.warn( + "The method post_process is deprecated and will be removed in v1.10. " + "All attributes that need to be post processed are computed just in time.", + DeprecationWarning, + ) + self.post_processed = True + 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. @@ -3266,53 +3229,191 @@ def calculate_stall_wind_velocity(self, stall_angle): # TODO: move to utilities + f" of attack exceeds {stall_angle:.3f}°: {w_v:.3f} m/s" ) - @deprecated( - reason="Moved to FlightDataExporter.export_pressures()", - version="v1.12.0", - alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_pressures", - ) - def export_pressures(self, file_name, time_step): - """ - .. deprecated:: 1.11 - Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` - and call ``.export_pressures(...)``. + 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 + parachutes, although it is not possible to get a noisy pressure signal + if no parachute is added. + + If a parachute is added, the file will contain 3 columns: time in + seconds, clean pressure in Pascals and noisy pressure in Pascals. + For flights without parachutes, the third column will be discarded + + This function was created especially for the 'Projeto Jupiter' + Electronics Subsystems team and aims to help in configuring + micro-controllers. + + Parameters + ---------- + file_name : string + The final file name, + time_step : float + Time step desired for the final file + + Return + ------ + None """ - return FlightDataExporter(self).export_pressures(file_name, time_step) + time_points = np.arange(0, self.t_final, time_step) + # pylint: disable=W1514, E1121 + with open(file_name, "w") as file: + 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.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.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 - @deprecated( - reason="Moved to FlightDataExporter.export_data()", - version="v1.12.0", - alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_data", - ) def export_data(self, file_name, *variables, time_step=None): + """Exports flight data to a comma separated value file (.csv). + + Data is exported in columns, with the first column representing time + steps. The first line of the file is a header line, specifying the + meaning of each column and its units. + + Parameters + ---------- + file_name : string + The file name or path of the exported file. Example: flight_data.csv + Do not use forbidden characters, such as / in Linux/Unix and + `<, >, :, ", /, \\, | ?, *` in Windows. + variables : strings, optional + Names of the data variables which shall be exported. Must be Flight + class attributes which are instances of the Function class. Usage + example: test_flight.export_data('test.csv', 'z', 'angle_of_attack', + 'mach_number'). + time_step : float, optional + Time step desired for the data. If None, all integration time steps + will be exported. Otherwise, linear interpolation is carried out to + calculate values at the desired time steps. Example: 0.001. """ - .. deprecated:: 1.11 - Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` - and call ``.export_data(...)``. - """ - return FlightDataExporter(self).export_data( - file_name, *variables, time_step=time_step + # 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: + np.savetxt( + file_name, + self.solution, + fmt="%.6f", + delimiter=",", + header="" + "Time (s)," + "X (m)," + "Y (m)," + "Z (m)," + "E0," + "E1," + "E2," + "E3," + "W1 (rad/s)," + "W2 (rad/s)," + "W3 (rad/s)", + ) + return + + # Not so fast evaluation for general case + if variables is None: + variables = [ + "x", + "y", + "z", + "vx", + "vy", + "vz", + "e0", + "e1", + "e2", + "e3", + "w1", + "w2", + "w3", + ] + + if time_step is None: + time_points = self.time + else: + time_points = np.arange(self.t_initial, self.t_final, time_step) + + exported_matrix = [time_points] + exported_header = "Time (s)" + + # Loop through variables, get points and names (for the header) + for variable in variables: + if variable in self.__dict__: + variable_function = self.__dict__[variable] + # Deal with decorated Flight methods + else: + try: + obj = getattr(self.__class__, variable) + variable_function = obj.__get__(self, self.__class__) + except AttributeError as exc: + raise AttributeError( + f"Variable '{variable}' not found in Flight class" + ) from exc + variable_points = variable_function(time_points) + exported_matrix += [variable_points] + exported_header += f", {variable_function.__outputs__[0]}" + + exported_matrix = np.array(exported_matrix).T # Fix matrix orientation + + np.savetxt( + file_name, + exported_matrix, + fmt="%.6f", + delimiter=",", + header=exported_header, + encoding="utf-8", ) - @deprecated( - reason="Moved to FlightDataExporter.export_sensor_data()", - version="v1.12.0", - alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_sensor_data", - ) def export_sensor_data(self, file_name, sensor=None): + """Exports sensors data to a file. The file format can be either .csv or + .json. + + Parameters + ---------- + file_name : str + The file name or path of the exported file. Example: flight_data.csv + Do not use forbidden characters, such as / in Linux/Unix and + `<, >, :, ", /, \\, | ?, *` in Windows. + sensor : Sensor, string, optional + The sensor to export data from. Can be given as a Sensor object or + as a string with the sensor name. If None, all sensors data will be + exported. Default is None. """ - .. deprecated:: 1.11 - Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` - and call ``.export_sensor_data(...)``. - """ - return FlightDataExporter(self).export_sensor_data(file_name, sensor=sensor) + if sensor is None: + data_dict = {} + for used_sensor, measured_data in self.sensor_data.items(): + data_dict[used_sensor.name] = measured_data + else: + # export data of only that sensor + data_dict = {} + + if not isinstance(sensor, str): + data_dict[sensor.name] = self.sensor_data[sensor] + else: # sensor is a string + matching_sensors = [s for s in self.sensor_data if s.name == sensor] + + if len(matching_sensors) > 1: + data_dict[sensor] = [] + for s in matching_sensors: + data_dict[s.name].append(self.sensor_data[s]) + elif len(matching_sensors) == 1: + data_dict[sensor] = self.sensor_data[matching_sensors[0]] + else: + raise ValueError("Sensor not found in the Flight.sensor_data.") - @deprecated( - reason="Moved to FlightDataExporter.export_kml()", - version="v1.12.0", - alternative="rocketpy.simulation.flight_data_exporter.FlightDataExporter.export_kml", - ) - def export_kml( + with open(file_name, "w") as file: + json.dump(data_dict, file) + print("Sensor data exported to: ", file_name) + + def export_kml( # TODO: should be moved out of this class. self, file_name="trajectory.kml", time_step=None, @@ -3320,18 +3421,78 @@ def export_kml( color="641400F0", altitude_mode="absolute", ): + """Exports flight data to a .kml file, which can be opened with Google + Earth to display the rocket's trajectory. + + Parameters + ---------- + file_name : string + The file name or path of the exported file. Example: flight_data.csv + time_step : float, optional + Time step desired for the data. If None, all integration time steps + will be exported. Otherwise, linear interpolation is carried out to + calculate values at the desired time steps. Example: 0.001. + extrude: bool, optional + To be used if you want to project the path over ground by using an + extruded polygon. In case False only the linestring containing the + flight path will be created. Default is True. + color : str, optional + Color of your trajectory path, need to be used in specific kml + format. Refer to http://www.zonums.com/gmaps/kml_color/ for more + info. + altitude_mode: str + Select elevation values format to be used on the kml file. Use + 'relativetoground' if you want use Above Ground Level elevation, or + 'absolute' if you want to parse elevation using Above Sea Level. + Default is 'relativetoground'. Only works properly if the ground + level is flat. Change to 'absolute' if the terrain is to irregular + or contains mountains. """ - .. deprecated:: 1.11 - Use :class:`rocketpy.simulation.flight_data_exporter.FlightDataExporter` - and call ``.export_kml(...)``. - """ - return FlightDataExporter(self).export_kml( - file_name=file_name, - time_step=time_step, - extrude=extrude, - color=color, - altitude_mode=altitude_mode, - ) + # Define time points vector + if time_step is None: + time_points = self.time + else: + time_points = np.arange(self.t_initial, self.t_final + time_step, time_step) + # Open kml file with simplekml library + kml = simplekml.Kml(open=1) + trajectory = kml.newlinestring(name="Rocket Trajectory - Powered by RocketPy") + + 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 + 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 + 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 + trajectory.style.linestyle.color = color + trajectory.style.polystyle.color = color + if extrude: + trajectory.extrude = 1 + # Save the KML + kml.save(file_name) + print("File ", file_name, " saved with success!") def info(self): """Prints out a summary of the data available about the Flight.""" @@ -3348,7 +3509,7 @@ def time_iterator(self, node_list): yield i, node_list[i] i += 1 - def to_dict(self, **kwargs): + def to_dict(self, include_outputs=False, **_kwargs): data = { "rocket": self.rocket, "env": self.env, @@ -3377,6 +3538,7 @@ def to_dict(self, **kwargs): "x_impact": self.x_impact, "y_impact": self.y_impact, "t_final": self.t_final, + "flight_phases": self.flight_phases, "function_evaluations": self.function_evaluations, "ax": self.ax, "ay": self.ay, @@ -3390,10 +3552,9 @@ def to_dict(self, **kwargs): "M1": self.M1, "M2": self.M2, "M3": self.M3, - "net_thrust": self.net_thrust, } - if kwargs.get("include_outputs", False): + if include_outputs: data.update( { "time": self.time, @@ -3571,9 +3732,7 @@ def add(self, flight_phase, index=None): # TODO: quite complex method new_index = ( index - 1 if flight_phase.t < previous_phase.t - else index + 1 - if flight_phase.t > next_phase.t - else index + else index + 1 if flight_phase.t > next_phase.t else index ) flight_phase.t += adjust self.add(flight_phase, new_index) diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 0e80152a7..a1f7af4eb 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -14,7 +14,6 @@ import math import re import time -import warnings from bisect import bisect_left import dill @@ -24,71 +23,13 @@ from cftime import num2pydate from matplotlib.patches import Ellipse from packaging import version as packaging_version +from rocketpy.mathutils.function import Function +from rocketpy.mathutils.vector_matrix import Vector # Mapping of module name and the name of the package that should be installed INSTALL_MAPPING = {"IPython": "ipython"} -def deprecated(reason=None, version=None, alternative=None): - """ - Decorator to mark functions or methods as deprecated. - - This decorator issues a DeprecationWarning when the decorated function - is called, indicating that it will be removed in future versions. - - Parameters - ---------- - reason : str, optional - Custom deprecation message. If not provided, a default message will be used. - version : str, optional - Version when the function will be removed. If provided, it will be - included in the warning message. - alternative : str, optional - Name of the alternative function/method that should be used instead. - If provided, it will be included in the warning message. - - Returns - ------- - callable - The decorated function with deprecation warning functionality. - - Examples - -------- - >>> @deprecated(reason="This function is obsolete", version="v2.0.0", - ... alternative="new_function") - ... def old_function(): - ... return "old result" - - >>> @deprecated() - ... def another_old_function(): - ... return "result" - """ - - def decorator(func): - @functools.wraps(func) - def wrapper(*args, **kwargs): - # Build the deprecation message - if reason: - message = reason - else: - message = f"The function `{func.__name__}` is deprecated" - - if version: - message += f" and will be removed in {version}" - - if alternative: - message += f". Use `{alternative}` instead" - - message += "." - - warnings.warn(message, DeprecationWarning, stacklevel=2) - return func(*args, **kwargs) - - return wrapper - - return decorator - - def tuple_handler(value): """Transforms the input value into a tuple that represents a range. If the input is an int or float, the output is a tuple from zero to the input @@ -182,13 +123,12 @@ def find_roots_cubic_function(a, b, c, d): Examples -------- >>> from rocketpy.tools import find_roots_cubic_function - >>> import cmath 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 = find_roots_cubic_function(a, b, c, d) - >>> cmath.isclose(x1, (-1+0j)) - True + >>> x1 + (-1+0j) To get the real part of the roots, use the real attribute of the complex number. @@ -1081,31 +1021,113 @@ def wrapper(*args, **kwargs): return decorator -def parallel_axis_theorem_from_com(com_inertia_moment, mass, distance): - """Calculates the moment of inertia of a object relative to a new axis using - the parallel axis theorem. The new axis is parallel to and at a distance - 'distance' from the original axis, which *must* passes through the object's - center of mass. - Parameters - ---------- - com_inertia_moment : float - Moment of inertia relative to the center of mass of the object. - mass : float - Mass of the object. - distance : float - Perpendicular distance between the original and new axis. - Returns - ------- - float - Moment of inertia relative to the new axis. +def _pat_dynamic_helper(com_inertia_moment, mass, distance_vec_3d, axes_term_lambda): + - References - ---------- - https://en.wikipedia.org/wiki/Parallel_axis_theorem - """ - return com_inertia_moment + mass * distance**2 + + is_dynamic = ( + isinstance(com_inertia_moment, Function) + or isinstance(mass, Function) + or isinstance(distance_vec_3d, Function) + ) + + def get_val(arg, t): + + return arg(t) if isinstance(arg, Function) else arg + + if not is_dynamic: + + d_vec = Vector(distance_vec_3d) + mass_term = mass * axes_term_lambda(d_vec) + return com_inertia_moment + mass_term + else: + + def new_source(t): + d_vec_t = get_val(distance_vec_3d, t) + mass_t = get_val(mass, t) + inertia_t = get_val(com_inertia_moment, t) + + mass_term = mass_t * axes_term_lambda(d_vec_t) + return inertia_t + mass_term + + return Function(new_source, inputs="t", outputs="Inertia (kg*m^2)") + + +def _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, product_term_lambda +): + + + is_dynamic = ( + isinstance(com_inertia_product, Function) + or isinstance(mass, Function) + or isinstance(distance_vec_3d, Function) + ) + + def get_val(arg, t): + + return arg(t) if isinstance(arg, Function) else arg + + if not is_dynamic: + + d_vec = Vector(distance_vec_3d) + mass_term = mass * product_term_lambda(d_vec) + return com_inertia_product + mass_term + else: + + def new_source(t): + d_vec_t = get_val(distance_vec_3d, t) + mass_t = get_val(mass, t) + inertia_t = get_val(com_inertia_product, t) + + mass_term = mass_t * product_term_lambda(d_vec_t) + return inertia_t + mass_term + + return Function(new_source, inputs="t", outputs="Inertia (kg*m^2)") + + +def parallel_axis_theorem_I11(com_inertia_moment, mass, distance_vec_3d): + + return _pat_dynamic_helper( + com_inertia_moment, mass, distance_vec_3d, lambda d_vec: d_vec.y**2 + d_vec.z**2 + ) + + +def parallel_axis_theorem_I22(com_inertia_moment, mass, distance_vec_3d): + + return _pat_dynamic_helper( + com_inertia_moment, mass, distance_vec_3d, lambda d_vec: d_vec.x**2 + d_vec.z**2 + ) + + +def parallel_axis_theorem_I33(com_inertia_moment, mass, distance_vec_3d): + + return _pat_dynamic_helper( + com_inertia_moment, mass, distance_vec_3d, lambda d_vec: d_vec.x**2 + d_vec.y**2 + ) + + +def parallel_axis_theorem_I12(com_inertia_product, mass, distance_vec_3d): + + return _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, lambda d_vec: d_vec.x * d_vec.y + ) + + +def parallel_axis_theorem_I13(com_inertia_product, mass, distance_vec_3d): + + return _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, lambda d_vec: d_vec.x * d_vec.z + ) + + +def parallel_axis_theorem_I23(com_inertia_product, mass, distance_vec_3d): + + return _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, lambda d_vec: d_vec.y * d_vec.z + ) # Flight @@ -1296,43 +1318,6 @@ def from_hex_decode(obj_bytes, decoder=base64.b85decode): return dill.loads(decoder(bytes.fromhex(obj_bytes))) -def find_obj_from_hash(obj, hash_, depth_limit=None): - """Searches the object (and its children) for - an object whose '__rpy_hash' field has a particular hash value. - - Parameters - ---------- - obj : object - Object to search. - hash_ : int - Hash value to search for in the '__rpy_hash' field. - depth_limit : int, optional - Maximum depth to search recursively. If None, no limit. - - Returns - ------- - object - The object whose '__rpy_hash' matches hash_, or None if not found. - """ - - stack = [(obj, 0)] - while stack: - current_obj, current_depth = stack.pop() - if depth_limit is not None and current_depth > depth_limit: - continue - - if getattr(current_obj, "__rpy_hash", None) == hash_: - return current_obj - - if isinstance(current_obj, dict): - stack.extend((v, current_depth + 1) for v in current_obj.values()) - - elif isinstance(current_obj, (list, tuple, set)): - stack.extend((item, current_depth + 1) for item in current_obj) - - return None - - if __name__ == "__main__": # pragma: no cover import doctest