diff --git a/.vscode/settings.json b/.vscode/settings.json index 6399d8142..5d8f887bf 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -241,6 +241,7 @@ "outerboundaryis", "overshootable", "planform", + "pointmassmotor", "polystyle", "powerseries", "Prandtl", diff --git a/docs/user/index.rst b/docs/user/index.rst index 8620e3ad5..113afc3aa 100644 --- a/docs/user/index.rst +++ b/docs/user/index.rst @@ -7,10 +7,11 @@ RocketPy's User Guide Installation and Requirements First Simulation + 3-DOF Simulations .. toctree:: :maxdepth: 1 - :caption: Main Classes Usage + :caption: Main Classes Usage Positions and Coordinate Systems Motors diff --git a/docs/user/three_dof_simulation.rst b/docs/user/three_dof_simulation.rst new file mode 100644 index 000000000..5375a9eb7 --- /dev/null +++ b/docs/user/three_dof_simulation.rst @@ -0,0 +1,424 @@ +.. _threedofsimulation: + +3-DOF Rocket Simulations +========================= + +RocketPy supports simplified 3-DOF (3 Degrees of Freedom) trajectory simulations, +where the rocket is modeled as a point mass. This mode is useful for quick +analyses, educational purposes, or when rotational dynamics are negligible. + +.. contents:: Table of Contents + :local: + :depth: 2 + +Overview +-------- + +In a 3-DOF simulation, the rocket's motion is described by three translational +degrees of freedom (x, y, z positions), ignoring all rotational dynamics. This +simplification: + +- **Reduces computational complexity** - Faster simulations for initial design studies +- **Focuses on trajectory** - Ideal for apogee predictions and flight path analysis +- **Simplifies model setup** - Requires fewer input parameters than full 6-DOF + +When to Use 3-DOF Simulations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +3-DOF simulations are appropriate when: + +- You need quick trajectory estimates during preliminary design +- Rotational stability is not a concern (e.g., highly stable rockets) +- You're performing educational demonstrations +- You want to validate basic flight performance before detailed analysis + +.. warning:: + + 3-DOF simulations **do not** account for: + + - Rocket rotation and attitude dynamics + - Stability margin and center of pressure effects + - Aerodynamic moments and angular motion + - Fin effectiveness and control surfaces + + For complete flight analysis including stability, use standard 6-DOF simulations. + +Setting Up a 3-DOF Simulation +------------------------------ + +A 3-DOF simulation requires three main components: + +1. :class:`rocketpy.PointMassMotor` - Motor without rotational inertia +2. :class:`rocketpy.PointMassRocket` - Rocket without rotational properties +3. :class:`rocketpy.Flight` with ``simulation_mode="3 DOF"`` + +Step 1: Define the Environment +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The environment setup is identical to standard simulations: + +.. jupyter-execute:: + + from rocketpy import Environment + + env = Environment( + latitude=32.990254, + longitude=-106.974998, + elevation=1400 + ) + + env.set_atmospheric_model(type="StandardAtmosphere") + +Step 2: Create a PointMassMotor +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The :class:`rocketpy.PointMassMotor` class represents a motor as a point mass, +without rotational inertia or grain geometry: + +.. jupyter-execute:: + + from rocketpy import PointMassMotor + + # Using a thrust curve file + motor = PointMassMotor( + thrust_source="../data/motors/cesaroni/Cesaroni_M1670.eng", + dry_mass=1.815, + propellant_initial_mass=2.5, + ) + +You can also define a constant thrust profile: + +.. jupyter-execute:: + + # Constant thrust of 250 N for 3 seconds + motor_constant = PointMassMotor( + thrust_source=250, + dry_mass=1.0, + propellant_initial_mass=0.5, + burn_time=3.0, + ) + +Or use a custom thrust function: + +.. jupyter-execute:: + + def custom_thrust(t): + """Custom thrust profile: ramps up, plateaus, then ramps down""" + if t < 0.5: + return 500 * t / 0.5 # Ramp up + elif t < 2.5: + return 500 # Plateau + elif t < 3.0: + return 500 * (3.0 - t) / 0.5 # Ramp down + else: + return 0 + + motor_custom = PointMassMotor( + thrust_source=custom_thrust, + dry_mass=1.2, + propellant_initial_mass=0.6, + burn_time=3.0, + ) + +.. seealso:: + + For detailed information about :class:`rocketpy.PointMassMotor` parameters, + see the :ref:`API reference `. + +Step 3: Create a PointMassRocket +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The :class:`rocketpy.PointMassRocket` class represents a rocket as a point mass: + +.. jupyter-execute:: + + from rocketpy import PointMassRocket + + rocket = PointMassRocket( + radius=0.0635, # meters + mass=5.0, # kg (dry mass without motor) + center_of_mass_without_motor=0.0, + power_off_drag=0.5, # Constant drag coefficient + power_on_drag=0.5, + ) + + # Add the motor + rocket.add_motor(motor, position=0) + +You can also specify drag as a function of Mach number: + +.. jupyter-execute:: + + # Drag coefficient vs Mach number + drag_curve = [ + [0.0, 0.50], + [0.5, 0.48], + [0.9, 0.52], + [1.1, 0.65], + [2.0, 0.55], + [3.0, 0.50], + ] + + rocket_with_drag_curve = PointMassRocket( + radius=0.0635, + mass=5.0, + center_of_mass_without_motor=0.0, + power_off_drag=drag_curve, + power_on_drag=drag_curve, + ) + +.. note:: + + Unlike the standard :class:`rocketpy.Rocket` class, :class:`rocketpy.PointMassRocket` + does **not** support: + + - Aerodynamic surfaces (fins, nose cones) + - Inertia tensors + - Center of pressure calculations + - Stability margin analysis + +Step 4: Run the Simulation +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Create a :class:`rocketpy.Flight` object with ``simulation_mode="3 DOF"``: + +.. jupyter-execute:: + + from rocketpy import Flight + + flight = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + inclination=85, # degrees from horizontal + heading=0, # degrees (0 = North, 90 = East) + simulation_mode="3 DOF", + max_time=100, + terminate_on_apogee=False, + ) + +.. important:: + + The ``simulation_mode="3 DOF"`` parameter is **required** to enable 3-DOF mode. + Without it, RocketPy will attempt a full 6-DOF simulation and may fail with + :class:`rocketpy.PointMassRocket`. + +Analyzing Results +----------------- + +Once the simulation is complete, you can access trajectory data and generate plots. + +Trajectory Information +^^^^^^^^^^^^^^^^^^^^^^ + +View key flight metrics: + +.. jupyter-execute:: + + flight.info() + +This will display: + +- Apogee altitude and time +- Maximum velocity +- Flight time +- Landing position + +Plotting Trajectory +^^^^^^^^^^^^^^^^^^^ + +Visualize the 3D flight path: + +.. jupyter-execute:: + + flight.plots.trajectory_3d() + +.. note:: + + In 3-DOF mode, the rocket maintains a fixed orientation (no pitch, yaw, or roll), + so attitude plots are not meaningful. + +Available Plots +^^^^^^^^^^^^^^^ + +The following plots are available for 3-DOF simulations: + +.. jupyter-execute:: + + # Altitude vs time + flight.z.plot() + + # Velocity components + flight.vx.plot() + flight.vy.plot() + flight.vz.plot() + + # Total velocity + flight.speed.plot() + + # Acceleration + flight.ax.plot() + +Export Data +^^^^^^^^^^^ + +Export trajectory data to CSV: + +.. jupyter-execute:: + + from rocketpy.simulation import FlightDataExporter + + exporter = FlightDataExporter(flight) + exporter.export_data( + "trajectory_3dof.csv", + "x", + "y", + "z", + "vx", + "vy", + "vz", + ) + +Complete Example +---------------- + +Here's a complete 3-DOF simulation from start to finish: + +.. jupyter-execute:: + + from rocketpy import Environment, PointMassMotor, PointMassRocket, Flight + + # 1. Environment + env = Environment( + latitude=39.3897, + longitude=-8.2889, + elevation=100 + ) + env.set_atmospheric_model(type="StandardAtmosphere") + + # 2. Motor + motor = PointMassMotor( + thrust_source=1500, # Constant 1500 N thrust + dry_mass=2.0, + propellant_initial_mass=3.0, + burn_time=4.0, + ) + + # 3. Rocket + rocket = PointMassRocket( + radius=0.0635, + mass=8.0, + center_of_mass_without_motor=0.0, + power_off_drag=0.45, + power_on_drag=0.45, + ) + rocket.add_motor(motor, position=0) + + # 4. Simulate + flight = Flight( + rocket=rocket, + environment=env, + rail_length=5.2, + inclination=85, + heading=0, + simulation_mode="3 DOF", + max_time=120, + ) + + # 5. Results + print(f"Apogee: {flight.apogee:.2f} m") + print(f"Max velocity: {flight.max_speed:.2f} m/s") + print(f"Flight time: {flight.t_final:.2f} s") + + flight.plots.trajectory_3d() + +Comparison: 3-DOF vs 6-DOF +--------------------------- + +Understanding the differences between simulation modes: + +.. list-table:: 3-DOF vs 6-DOF Comparison + :header-rows: 1 + :widths: 30 35 35 + + * - Feature + - 3-DOF + - 6-DOF + * - Computational Speed + - Fast + - Slower + * - Rocket Orientation + - Fixed (no rotation) + - Full attitude dynamics + * - Stability Analysis + - ❌ Not available + - ✅ Full stability margin + * - Aerodynamic Surfaces + - ❌ Not modeled + - ✅ Fins, nose, tail + * - Center of Pressure + - ❌ Not computed + - ✅ Computed + * - Moments of Inertia + - ❌ Not needed + - ✅ Required + * - Use Cases + - Quick estimates, education + - Detailed design, stability + * - Trajectory Accuracy + - Good for stable rockets + - Highly accurate + +Best Practices +-------------- + +1. **Validate with 6-DOF**: After getting initial results with 3-DOF, validate + critical designs with full 6-DOF simulations. + +2. **Check Drag Coefficient**: Ensure your drag coefficient is realistic for + your rocket's geometry. Use wind tunnel data or CFD if available. + +3. **Use Realistic Launch Conditions**: Even in 3-DOF mode, wind conditions + and rail length affect trajectory. + +4. **Document Assumptions**: Clearly document that your analysis uses 3-DOF + and its limitations. + +Limitations and Warnings +------------------------ + +.. danger:: + + **Critical Limitations:** + + - **No stability checking** - The simulation cannot detect unstable rockets + - **No attitude control** - Air brakes and thrust vectoring are not supported + - **Assumes perfect alignment** - Rocket always points along velocity vector + - **No wind weathercocking** - Wind effects on orientation are ignored + +.. warning:: + + 3-DOF simulations should **not** be used for: + + - Final design verification + - Stability margin analysis + - Control system design + - Fin sizing and optimization + - Safety-critical trajectory predictions + +See Also +-------- + +- :ref:`First Simulation ` - Standard 6-DOF simulation tutorial +- :ref:`Rocket Class Usage ` - Full rocket modeling capabilities +- :ref:`Flight Class Usage ` - Complete flight simulation options +- :doc:`../examples/3_dof_trial_sim` - Jupyter notebook example + +Further Reading +--------------- + +For more information about point mass trajectory simulations: + +- `Trajectory Optimization `_ +- `Equations of Motion `_ +- `Point Mass Model `_ diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index f99a70f28..852a16aef 100644 --- a/rocketpy/__init__.py +++ b/rocketpy/__init__.py @@ -17,6 +17,7 @@ MassBasedTank, MassFlowRateBasedTank, Motor, + PointMassMotor, SolidMotor, SphericalTank, Tank, @@ -35,6 +36,7 @@ LinearGenericSurface, NoseCone, Parachute, + PointMassRocket, RailButtons, Rocket, Tail, diff --git a/rocketpy/motors/__init__.py b/rocketpy/motors/__init__.py index 73f49633c..b13ff9392 100644 --- a/rocketpy/motors/__init__.py +++ b/rocketpy/motors/__init__.py @@ -3,6 +3,7 @@ from .hybrid_motor import HybridMotor from .liquid_motor import LiquidMotor from .motor import GenericMotor, Motor +from .point_mass_motor import PointMassMotor from .solid_motor import SolidMotor from .tank import ( LevelBasedTank, diff --git a/rocketpy/motors/point_mass_motor.py b/rocketpy/motors/point_mass_motor.py new file mode 100644 index 000000000..7cd42636f --- /dev/null +++ b/rocketpy/motors/point_mass_motor.py @@ -0,0 +1,106 @@ +from functools import cached_property +from typing import Callable + +import numpy as np + +from rocketpy.mathutils.function import Function, funcify_method + +from .motor import Motor + + +class PointMassMotor(Motor): + """Motor modeled as a point mass for 3-DOF simulations.""" + + def __init__( + self, + thrust_source, + dry_mass, + propellant_initial_mass, + burn_time=None, + propellant_final_mass=None, + reshape_thrust_curve=False, + interpolation_method="linear", + ): + if isinstance(thrust_source, (int, float, Callable)): + if propellant_initial_mass is None: + raise ValueError( + "For constant or callable thrust, 'propellant_initial_mass' is required." + ) + if burn_time is None and propellant_final_mass is None: + raise ValueError( + "For constant or callable thrust, either 'burn_time' or " + "'propellant_final_mass' must be provided." + ) + elif isinstance(thrust_source, (Function, np.ndarray, str)): + if propellant_initial_mass is None: + raise ValueError( + "For thrust from a Function, NumPy array, or CSV, 'propellant_initial_mass' is required." + ) + else: + raise TypeError( + "Invalid 'thrust_source' type. Must be int, float, callable, str, numpy.ndarray, or Function." + ) + + self._propellant_initial_mass = propellant_initial_mass + self.propellant_final_mass = propellant_final_mass + + super().__init__( + thrust_source=thrust_source, + dry_inertia=(0, 0, 0), + nozzle_radius=0, + center_of_dry_mass_position=0, + dry_mass=dry_mass, + nozzle_position=0, + burn_time=burn_time, + reshape_thrust_curve=reshape_thrust_curve, + interpolation_method=interpolation_method, + coordinate_system_orientation="nozzle_to_combustion_chamber", + ) + + @property + def propellant_initial_mass(self): + return self._propellant_initial_mass + + @funcify_method("Time (s)", "Exhaust velocity (m/s)") + def exhaust_velocity(self): + """Assume constant exhaust velocity: total impulse / propellant mass""" + v_e = self.total_impulse / self.propellant_initial_mass + return Function(v_e).set_discrete_based_on_model(self.thrust) + + @cached_property + def total_mass_flow_rate(self) -> Function: + """Mass flow rate: -thrust / exhaust_velocity""" + return -self.thrust / self.exhaust_velocity + + @cached_property + def center_of_propellant_mass(self): + """Center of propellant mass is always zero""" + return Function(0.0) + + # Propellant inertias: always zero, but return as Function objects + def _zero_inertia_func(self): + return Function(0.0) + + @cached_property + def propellant_I_11(self): + return self._zero_inertia_func() + + @cached_property + def propellant_I_12(self): + return self._zero_inertia_func() + + @cached_property + def propellant_I_13(self): + return self._zero_inertia_func() + + @cached_property + def propellant_I_22(self): + return self._zero_inertia_func() + + @cached_property + def propellant_I_23(self): + return self._zero_inertia_func() + + @cached_property + def propellant_I_33(self): + return self._zero_inertia_func() diff --git a/rocketpy/rocket/__init__.py b/rocketpy/rocket/__init__.py index 0687b5ee5..463cbe3b3 100644 --- a/rocketpy/rocket/__init__.py +++ b/rocketpy/rocket/__init__.py @@ -14,4 +14,5 @@ ) from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute +from rocketpy.rocket.point_mass_rocket import PointMassRocket from rocketpy.rocket.rocket import Rocket diff --git a/rocketpy/rocket/point_mass_rocket.py b/rocketpy/rocket/point_mass_rocket.py new file mode 100644 index 000000000..9cdf86d47 --- /dev/null +++ b/rocketpy/rocket/point_mass_rocket.py @@ -0,0 +1,103 @@ +"""Point mass rocket class for simplified 3-DOF trajectory simulations.""" + +from rocketpy.rocket.rocket import Rocket + + +class PointMassRocket(Rocket): + """A simplified rocket class for trajectory simulations where the rocket + is modeled as a point mass. + + This class omits rotational dynamics and complex inertial properties, + focusing solely on translational (3-DOF) motion based on mass and + aerodynamics. Appropriate for educational use, quick analyses, or when + rotational effects are negligible. + + Parameters + ---------- + radius : float + Rocket's largest radius in meters. + mass : float + Rocket's mass without motor in kg. + center_of_mass_without_motor : float + Position, in meters, of the rocket's center of mass without motor + relative to the rocket's coordinate system. + power_off_drag : float, callable, array, string, Function + Drag coefficient as a function of Mach number when the motor is off. + power_on_drag : float, callable, array, string, Function + Drag coefficient as a function of Mach number when the motor is on. + + Attributes + ---------- + radius : float + Rocket's largest radius in meters. + mass : float + Rocket's mass without motor in kg. + center_of_mass_without_motor : float + Position, in meters, of the rocket's center of mass without motor + relative to the rocket's coordinate system. + power_off_drag : Function + Drag coefficient as a function of Mach number when the motor is off. + power_on_drag : Function + Drag coefficient as a function of Mach number when the motor is on. + """ + + def __init__( + self, + radius: float, + mass: float, + center_of_mass_without_motor: float, + power_off_drag, + power_on_drag, + ): + self._center_of_mass_without_motor_pointmass = center_of_mass_without_motor + self._center_of_dry_mass_position = center_of_mass_without_motor + self._center_of_mass = center_of_mass_without_motor + # Dry inertias are zero for point mass + self.dry_I_11 = 0.0 + self.dry_I_22 = 0.0 + self.dry_I_33 = 0.0 + self.dry_I_12 = 0.0 + self.dry_I_13 = 0.0 + self.dry_I_23 = 0.0 + + # Call base init with safe defaults + super().__init__( + radius=radius, + mass=mass, + inertia=(0, 0, 0), + power_off_drag=power_off_drag, + power_on_drag=power_on_drag, + center_of_mass_without_motor=center_of_mass_without_motor, + ) + + def evaluate_dry_inertias(self): + """Override to ensure inertias remain zero for point mass model. + + Returns + ------- + tuple + All inertia components as zeros. + """ + self.dry_I_11 = 0.0 + self.dry_I_22 = 0.0 + self.dry_I_33 = 0.0 + self.dry_I_12 = 0.0 + self.dry_I_13 = 0.0 + self.dry_I_23 = 0.0 + return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + + def evaluate_inertias(self): + """Override to ensure inertias remain zero for point mass model. + + Returns + ------- + tuple + All inertia components as zeros. + """ + self.I_11 = 0.0 + self.I_22 = 0.0 + self.I_33 = 0.0 + self.I_12 = 0.0 + self.I_13 = 0.0 + self.I_23 = 0.0 + return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a38be7d93..30ea66466 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -11,8 +11,10 @@ from ..mathutils.function import Function, funcify_method from ..mathutils.vector_matrix import Matrix, Vector +from ..motors.point_mass_motor import PointMassMotor from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints +from ..rocket import PointMassRocket from ..tools import ( calculate_cubic_hermite_coefficients, deprecated, @@ -466,6 +468,8 @@ class Flight: Defined as the minimum angle between the attitude vector and the freestream velocity vector. Can be called or accessed as array. + Flight.simulation_mode : str + Simulation mode for the flight. Can be "6 DOF" or "3 DOF". """ def __init__( # pylint: disable=too-many-arguments,too-many-statements @@ -487,6 +491,7 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements name="Flight", equations_of_motion="standard", ode_solver="LSODA", + simulation_mode="6 DOF", ): """Run a trajectory simulation. @@ -599,6 +604,7 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements self.terminate_on_apogee = terminate_on_apogee self.name = name self.equations_of_motion = equations_of_motion + self.simulation_mode = simulation_mode self.ode_solver = ode_solver # Controller initialization @@ -1229,9 +1235,34 @@ def __init_solver_monitors(self): def __init_equations_of_motion(self): """Initialize equations of motion.""" - if self.equations_of_motion == "solid_propulsion": - # NOTE: The u_dot is faster, but only works for solid propulsion - self.u_dot_generalized = self.u_dot + # Determine if a point-mass model is used. + is_point_mass = isinstance(self.rocket, PointMassRocket) or ( + hasattr(self.rocket, "motor") + and isinstance(self.rocket.motor, PointMassMotor) + ) + # Set simulation mode based on model type. + if is_point_mass: + if self.simulation_mode != "3 DOF": + warnings.warn( + "A point-mass model was detected. Simulation mode should be '3 DOF'.", + UserWarning, + ) + self.simulation_mode = "3 DOF" + + # Set the equations of motion based on the final simulation mode. + if self.simulation_mode == "3 DOF": + self.u_dot_generalized = self.u_dot_generalized_3dof + elif self.simulation_mode == "6 DOF": + self.u_dot_generalized = ( + self.u_dot + if self.equations_of_motion == "solid_propulsion" + else self.u_dot_generalized + ) + else: + raise ValueError( + f"Invalid simulation_mode: {self.simulation_mode}. " + "Must be '3 DOF' or '6 DOF'." + ) def __init_controllers(self): """Initialize controllers and sensors""" @@ -1465,6 +1496,7 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Retrieve integration data _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u # Determine lift force and moment + omega1, omega2, omega3 = 0, 0, 0 R1, R2, M1, M2, M3 = 0, 0, 0, 0, 0 # Thrust correction parameters pressure = self.env.pressure.get_value_opt(z) @@ -1485,6 +1517,7 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals mass_flow_rate_at_t = self.rocket.motor.mass_flow_rate.get_value_opt(t) propellant_mass_at_t = self.rocket.motor.propellant_mass.get_value_opt(t) # Thrust + net_thrust = max( self.rocket.motor.thrust.get_value_opt(t) + self.rocket.motor.pressure_thrust(pressure), @@ -1693,7 +1726,6 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals e1dot = 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3) e2dot = 0.5 * (omega2 * e0 - omega3 * e1 + omega1 * e3) e3dot = 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2) - # Linear acceleration L = [ ( @@ -1756,6 +1788,123 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals net_thrust, ] ) + return u_dot + + def u_dot_generalized_3dof(self, t, u, post_processing=False): + """Calculates derivative of u state vector with respect to time when the + rocket is flying in 3 DOF motion in space and significant mass variation + effects exist. + + Parameters + ---------- + t : float + Time in seconds. + u : list + State vector: [x, y, z, vx, vy, vz, q0, q1, q2, q3, omega1, omega2, omega3]. + post_processing : bool, optional + If True, adds flight data to self variables like self.angle_of_attack. + + Returns + ------- + list + Derivative state vector: [vx, vy, vz, ax, ay, az, + e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. + """ + # Unpack state + _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u + + # Define vectors + v = Vector([vx, vy, vz]) + e = [e0, e1, e2, e3] + w = Vector([omega1, omega2, omega3]) + + # Mass and transformation + total_mass = self.rocket.total_mass.get_value_opt(t) + K = Matrix.transformation(e) + Kt = K.transpose + + # Atmospheric and wind data + rho = self.env.density.get_value_opt(z) + wind_vx = self.env.wind_velocity_x.get_value_opt(z) + wind_vy = self.env.wind_velocity_y.get_value_opt(z) + wind_velocity = Vector([wind_vx, wind_vy, 0]) + + free_stream_velocity = wind_velocity - v + free_stream_speed = abs(free_stream_velocity) + speed_of_sound = self.env.speed_of_sound.get_value_opt(z) + mach = free_stream_speed / speed_of_sound + + # Drag computation + if t < self.rocket.motor.burn_out_time: + cd = self.rocket.power_on_drag.get_value_opt(mach) + else: + cd = self.rocket.power_off_drag.get_value_opt(mach) + + R1, R2 = 0, 0 + R3 = -0.5 * rho * free_stream_speed**2 * self.rocket.area * cd + + for air_brake in self.rocket.air_brakes: + if air_brake.deployment_level > 0: + ab_cd = air_brake.drag_coefficient.get_value_opt( + air_brake.deployment_level, mach + ) + ab_force = ( + -0.5 * rho * free_stream_speed**2 * air_brake.reference_area * ab_cd + ) + if air_brake.override_rocket_drag: + R3 = ab_force + else: + R3 += ab_force + + # Velocity in body frame + vb_body = Kt @ v + + for surface, _ in self.rocket.aerodynamic_surfaces: + cp = self.rocket.surfaces_cp_to_cdm[surface] + vb_component = vb_body + (w ^ cp) + + comp_z = z + (K @ cp).z + wind_cx = self.env.wind_velocity_x.get_value_opt(comp_z) + wind_cy = self.env.wind_velocity_y.get_value_opt(comp_z) + wind_body = Kt @ Vector([wind_cx, wind_cy, 0]) + + rel_velocity = wind_body - vb_component + rel_speed = abs(rel_velocity) + rel_mach = rel_speed / speed_of_sound + + reynolds = ( + self.env.density.get_value_opt(comp_z) + * rel_speed + * surface.reference_length + / self.env.dynamic_viscosity.get_value_opt(comp_z) + ) + + fx, fy, fz, *_ = surface.compute_forces_and_moments( + rel_velocity, rel_speed, rel_mach, rho, cp, w, reynolds + ) + R1 += fx + R2 += fy + R3 += fz + + # Thrust and weight + thrust = self.rocket.motor.thrust.get_value_opt(t) + gravity = self.env.gravity.get_value_opt(z) + weight_body = Kt @ Vector([0, 0, -total_mass * gravity]) + + total_force = Vector([0, 0, thrust]) + weight_body + Vector([R1, R2, R3]) + + # Dynamics + v_dot = K @ (total_force / total_mass) + e_dot = [0, 0, 0, 0] # Euler derivatives unused in 3DOF + w_dot = [0, 0, 0] # No angular dynamics in 3DOF + r_dot = [vx, vy, vz] + + u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] + + if post_processing: + self.__post_processed_variables.append( + [t, *v_dot, *w_dot, R1, R2, R3, 0, 0, 0] + ) return u_dot diff --git a/tests/conftest.py b/tests/conftest.py index 2f6124900..976e177db 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,5 +1,9 @@ +import matplotlib import pytest +# Configure matplotlib to use non-interactive backend for tests +matplotlib.use("Agg") + # Pytest configuration pytest_plugins = [ "tests.fixtures.environment.environment_fixtures", diff --git a/tests/unit/motors/test_point_mass_motor.py b/tests/unit/motors/test_point_mass_motor.py new file mode 100644 index 000000000..1eb248723 --- /dev/null +++ b/tests/unit/motors/test_point_mass_motor.py @@ -0,0 +1,79 @@ +import pytest + +from rocketpy.motors.point_mass_motor import PointMassMotor + + +def test_init_required_args(): + """Tests that PointMassMotor initializes correctly with required arguments. + + Verifies that the motor is properly instantiated and that basic properties + like dry_mass, propellant_initial_mass, and burn_time are correctly set. + """ + m = PointMassMotor( + thrust_source=10, dry_mass=1.0, propellant_initial_mass=0.5, burn_time=1.2 + ) + assert isinstance(m, PointMassMotor) + assert m.dry_mass == 1.0 + assert m.propellant_initial_mass == 0.5 + assert m.burn_time == (0, 1.2) # burn_time is always a tuple (start, end) + assert m.burn_duration == 1.2 + + +def test_missing_required_args_raises(): + """Tests that PointMassMotor raises errors for missing required arguments. + + Verifies that ValueError is raised when propellant_initial_mass is None + or when burn_time is not provided with constant thrust. Also verifies + TypeError is raised for invalid thrust_source types. + """ + # TODO: in the future we would like to capture specific RocketPy Exceptions + with pytest.raises(ValueError): + PointMassMotor(10, 1.0, None, 1) + with pytest.raises(ValueError): + PointMassMotor(10, 1.0, 1.2) + with pytest.raises(TypeError): + PointMassMotor([], 1.0, 1.2, 1) + + +def test_exhaustvelocity_and_totalmassflowrate(): + """Tests that exhaust_velocity and total_mass_flow_rate return Function objects. + + Verifies that both properties are callable functions with get_value method, + which is required for numerical evaluation during simulation. + """ + m = PointMassMotor( + thrust_source=10, dry_mass=1.0, propellant_initial_mass=1.0, burn_time=2.0 + ) + ev_func = m.exhaust_velocity + assert hasattr(ev_func, "get_value") + tmf_func = m.total_mass_flow_rate + assert hasattr(tmf_func, "get_value") + + +def test_zero_inertias(): + """Tests that all propellant inertia components are zero for point mass motor. + + Verifies that propellant_I_11, propellant_I_22, and propellant_I_33 all + return zero, as expected for a point mass model without rotational dynamics. + """ + m = PointMassMotor( + thrust_source=10, dry_mass=1.0, propellant_initial_mass=1.0, burn_time=2.0 + ) + assert m.propellant_I_11.get_value(0) == 0 + assert m.propellant_I_22.get_value(0) == 0 + assert m.propellant_I_33.get_value(0) == 0 + + +def test_callable_thrust(): + """Tests that PointMassMotor accepts a callable function as thrust_source. + + Verifies that when a lambda function is used for thrust_source, the motor + correctly evaluates thrust values at different times. + """ + m = PointMassMotor( + thrust_source=lambda t: 100 * t, + dry_mass=2, + propellant_initial_mass=2, + burn_time=4, + ) + assert m.thrust_source(0.5) == 50 diff --git a/tests/unit/rocket/test_point_mass_rocket.py b/tests/unit/rocket/test_point_mass_rocket.py new file mode 100644 index 000000000..c9ad709c9 --- /dev/null +++ b/tests/unit/rocket/test_point_mass_rocket.py @@ -0,0 +1,59 @@ +from rocketpy.motors.point_mass_motor import PointMassMotor +from rocketpy.rocket.point_mass_rocket import PointMassRocket + + +def test_init_sets_basic_properties_correctly(): + """Tests that PointMassRocket initializes with correct basic properties. + + Verifies that radius, mass, motor assignment, and zero inertias are + correctly set when a PointMassRocket is created and a motor is added. + """ + motor = PointMassMotor(10, 1.0, 0.5, 1.0) + rocket = PointMassRocket( + radius=0.05, + mass=2.0, + center_of_mass_without_motor=0.1, + power_off_drag=0.7, + power_on_drag=0.8, + ) + rocket.add_motor(motor, 0) + assert rocket.radius == 0.05 + assert rocket.mass == 2.0 + assert rocket.motor is motor + assert rocket.dry_I_11 == 0.0 # 3-DOF: inertias are forced zero + + +def test_structural_and_total_mass(): + """Tests structural and total mass properties of PointMassRocket. + + Verifies that the rocket's structural mass and total mass (including motor + dry mass and propellant) are calculated correctly at time t=0. + """ + motor = PointMassMotor(10, 1.0, 1.1, 2.0) + rocket = PointMassRocket( + radius=0.03, + mass=2.5, + center_of_mass_without_motor=0, + power_off_drag=0.3, + power_on_drag=0.4, + ) + rocket.add_motor(motor, 0) + + # Test that structural mass and total mass are calculated correctly + assert rocket.mass == 2.5 + expected_total_mass = rocket.mass + motor.dry_mass + motor.propellant_initial_mass + assert abs(rocket.total_mass(0) - expected_total_mass) < 1e-6 + + +def test_add_motor_overwrites(): + """Tests that adding a new motor to PointMassRocket overwrites the previous motor. + + Verifies that when add_motor is called multiple times, only the most + recently added motor is retained. + """ + motor1 = PointMassMotor(10, 1, 1.1, 2.0) + motor2 = PointMassMotor(20, 2, 1.5, 3.0) + rocket = PointMassRocket(0.02, 1.0, 0.0, 0.2, 0.5) + rocket.add_motor(motor1, position=0) + rocket.add_motor(motor2, position=0) + assert rocket.motor == motor2 diff --git a/tests/unit/simulation/test_flight_3dof.py b/tests/unit/simulation/test_flight_3dof.py new file mode 100644 index 000000000..a44448baa --- /dev/null +++ b/tests/unit/simulation/test_flight_3dof.py @@ -0,0 +1,139 @@ +import numpy as np +import pytest + +from rocketpy.motors.point_mass_motor import PointMassMotor +from rocketpy.rocket.point_mass_rocket import PointMassRocket +from rocketpy.simulation.flight import Flight + + +@pytest.fixture +def point_mass_motor(): + """Creates a simple PointMassMotor for 3-DOF flight tests. + + Returns + ------- + rocketpy.PointMassMotor + A point mass motor with constant thrust of 10 N, 1.0 kg dry mass, + 0.5 kg propellant mass, and 2.2 s burn time. + """ + return PointMassMotor( + thrust_source=10, + dry_mass=1.0, + propellant_initial_mass=0.5, + burn_time=2.2, + ) + + +@pytest.fixture +def point_mass_rocket(point_mass_motor): + """Creates a simple PointMassRocket for 3-DOF flight tests. + + Parameters + ---------- + point_mass_motor : rocketpy.PointMassMotor + The motor to be added to the rocket. + + Returns + ------- + rocketpy.PointMassRocket + A point mass rocket with 0.05 m radius, 2.0 kg mass, and the + provided motor attached at position 0. + """ + rocket = PointMassRocket( + radius=0.05, + mass=2.0, + center_of_mass_without_motor=0.1, + power_off_drag=0.5, + power_on_drag=0.6, + ) + rocket.add_motor(point_mass_motor, position=0) + return rocket + + +def test_simulation_mode_sets_3dof_with_point_mass_rocket( + example_plain_env, point_mass_rocket +): + """Tests that simulation mode is correctly set to 3 DOF for PointMassRocket. + + Parameters + ---------- + example_plain_env : rocketpy.Environment + A basic environment fixture for flight simulation. + point_mass_rocket : rocketpy.PointMassRocket + A point mass rocket fixture for 3-DOF simulation. + """ + flight = Flight( + rocket=point_mass_rocket, + environment=example_plain_env, + rail_length=1, + simulation_mode="3 DOF", + ) + assert flight.simulation_mode == "3 DOF" + + +def test_3dof_simulation_mode_warning(example_plain_env, point_mass_rocket): + """Tests that a warning is issued when 6 DOF mode is requested with PointMassRocket. + + When a PointMassRocket is used with simulation_mode="6 DOF", the Flight + class should emit a UserWarning and automatically switch to 3 DOF mode. + + Parameters + ---------- + example_plain_env : rocketpy.Environment + A basic environment fixture for flight simulation. + point_mass_rocket : rocketpy.PointMassRocket + A point mass rocket fixture for 3-DOF simulation. + """ + with pytest.warns(UserWarning): + flight = Flight( + rocket=point_mass_rocket, + environment=example_plain_env, + rail_length=1, + simulation_mode="6 DOF", + ) + assert flight.simulation_mode == "3 DOF" + + +def test_u_dot_generalized_3dof_returns_valid_result( + example_plain_env, point_mass_rocket +): + """Tests that 3-DOF equations of motion return valid derivative results. + + Verifies that the u_dot_generalized_3dof method returns a list or numpy + array representing the state derivative vector. + + Parameters + ---------- + example_plain_env : rocketpy.Environment + A basic environment fixture for flight simulation. + point_mass_rocket : rocketpy.PointMassRocket + A point mass rocket fixture for 3-DOF simulation. + """ + flight = Flight( + rocket=point_mass_rocket, + environment=example_plain_env, + rail_length=1, + simulation_mode="3 DOF", + ) + u = [0] * 13 # Generalized state vector size + result = flight.u_dot_generalized_3dof(0, u) + assert isinstance(result, (list, np.ndarray)) + + +def test_invalid_simulation_mode(example_plain_env, calisto): + """Tests that invalid simulation mode raises ValueError. + + Parameters + ---------- + example_plain_env : rocketpy.Environment + A basic environment fixture for flight simulation. + calisto : rocketpy.Rocket + The Calisto rocket fixture from the test suite. + """ + with pytest.raises(ValueError): + Flight( + rocket=calisto, + environment=example_plain_env, + rail_length=1, + simulation_mode="2 DOF", + )