From c999891f923f68b265089b3feba61b1e9c7a5e5a Mon Sep 17 00:00:00 2001 From: GuilhermeAsura Date: Sat, 6 Dec 2025 00:29:05 -0300 Subject: [PATCH 01/12] build: add 'vedo' as optional dependency for animation --- pyproject.toml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 4f1ecced4..16cc479e8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,7 +68,11 @@ monte-carlo = [ "contextily>=1.0.0; python_version < '3.14'", ] -all = ["rocketpy[env-analysis]", "rocketpy[monte-carlo]"] +animation = [ + "vedo>=2024.5.1" +] + +all = ["rocketpy[env-analysis]", "rocketpy[monte-carlo]", "rocketpy[animation]"] [tool.coverage.report] From 61e526a141a4e351df07117e74d8820010b466c8 Mon Sep 17 00:00:00 2001 From: GuilhermeAsura Date: Sat, 6 Dec 2025 00:33:38 -0300 Subject: [PATCH 02/12] feat: add animate_trajectory and animate_rotate to Flight class - Ported methods from legacy 'animate_flight' branch - Adapted to new Flight class structure (removed postProcess) - Added vedo as optional dependency Co-authored-by: Patrick Sampaio --- rocketpy/simulation/flight.py | 196 ++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 2293d9706..afa95d596 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1,5 +1,6 @@ # pylint: disable=too-many-lines import math +import time import warnings from copy import deepcopy from functools import cached_property @@ -4122,7 +4123,202 @@ def export_kml( color=color, altitude_mode=altitude_mode, ) + def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwargs): + """ + 6-DOF Animation of the flight trajectory using Vedo. + + Parameters + ---------- + file_name : str + 3D object file representing your rocket, usually in .stl format. + Example: "rocket.stl" + start : int, float, optional + Time for starting animation, in seconds. Default is 0. + stop : int, float, optional + Time for ending animation, in seconds. If None, uses self.t_final. + Default is None. + time_step : float, optional + Time step for data interpolation in the animation. Default is 0.1. + **kwargs : dict, optional + Additional keyword arguments to be passed to vedo.Plotter.show(). + Common arguments: + - azimuth (float): Rotation in degrees around the vertical axis. + - elevation (float): Rotation in degrees above the horizon. + - roll (float): Rotation in degrees around the view axis. + - zoom (float): Zoom level (default 1). + + Returns + ------- + None + + Raises + ------ + ImportError + If the 'vedo' package is not installed. + + Notes + ----- + This feature requires the 'vedo' package. Install it with: + pip install rocketpy[animation] + """ + try: + from vedo import Box, Line, Mesh, Plotter, settings + except ImportError as e: + raise ImportError( + "The animation feature requires the 'vedo' package. " + "Install it with:\n" + " pip install rocketpy[animation]\n" + "Or directly:\n" + " pip install vedo>=2024.5.1" + ) from e + + # Enable interaction if needed + try: + settings.allow_interaction = True + except AttributeError: + pass # Not available in newer versions of vedo + + # Handle stop time + if stop is None: + stop = self.t_final + + # Define the world bounds based on trajectory + max_x = max(self.x[:, 1]) + max_y = max(self.y[:, 1]) + # Use simple logic for bounds + world = Box( + pos=[self.x(start), self.y(start), self.apogee], + length=max_x * 2 if max_x != 0 else 1000, + width=max_y * 2 if max_y != 0 else 1000, + height=self.apogee, + ).wireframe() + + # Load rocket mesh + rocket = Mesh(file_name).c("green") + rocket.pos(self.x(start), self.y(start), 0).add_trail(n=len(self.x[:, 1])) + # Create trail + trail_points = [[self.x(t), self.y(t), self.z(t) - self.env.elevation] + for t in np.arange(start, stop, time_step)] + trail = Line(trail_points, c="k", alpha=0.5) + # Setup Plotter + plt = Plotter(axes=1, interactive=False) + plt.show(world, rocket, __doc__, viewup="z", **kwargs) + + # Animation Loop + for t in np.arange(start, stop, time_step): + # Calculate rotation angle and vector from quaternions + # Note: This simple rotation logic mimics the old branch. + # Ideally, vedo handles orientation via matrix, but we stick + # to the provided logic for now. + + # e0 is the scalar part of the quaternion + angle = np.arccos(2 * self.e0(t)**2 - 1) + k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 + + # Update position and rotation + # Adjusting for ground elevation + rocket.pos(self.x(t), self.y(t), self.z(t) - self.env.elevation) + rocket.rotate_x(self.e1(t) / k) + rocket.rotate_y(self.e2(t) / k) + rocket.rotate_z(self.e3(t) / k) + + # update the scene + plt.show(world, rocket, trail) + + # slow down to make animation visible + start_pause = time.time() + while time.time() - start_pause < time_step: + plt.render() + + if getattr(plt, 'escaped', False): + break + + plt.interactive().close() + return None + + def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs): + """ + Animation of rocket attitude (rotation) during the flight. + + Parameters + ---------- + file_name : str + 3D object file representing your rocket, usually in .stl format. + start : int, float, optional + Time for starting animation, in seconds. Default is 0. + stop : int, float, optional + Time for ending animation, in seconds. If None, uses self.t_final. + Default is None. + time_step : float, optional + Time step for data interpolation. Default is 0.1. + **kwargs : dict, optional + Additional keyword arguments to be passed to vedo.Plotter.show(). + Common arguments: + - azimuth (float): Rotation in degrees around the vertical axis. + - elevation (float): Rotation in degrees above the horizon. + - roll (float): Rotation in degrees around the view axis. + - zoom (float): Zoom level (default 1). + + Returns + ------- + None + Raises + ------ + ImportError + If the 'vedo' package is not installed. + """ + try: + from vedo import Box, Mesh, Plotter, settings + except ImportError as e: + raise ImportError( + "The animation feature requires the 'vedo' package. " + "Install it with:\n" + " pip install rocketpy[animation]\n" + ) from e + + # Enable interaction if needed + try: + settings.allow_interaction = True + except AttributeError: + pass # Not available in newer versions of vedo + + if stop is None: + stop = self.t_final + + # Smaller box for rotation view + world = Box( + pos=[self.x(start), self.y(start), self.apogee], + length=max(self.x[:, 1]) * 0.2, + width=max(self.y[:, 1]) * 0.2, + height=self.apogee * 0.1, + ).wireframe() + + rocket = Mesh(file_name).c("green") + # Initialize at origin relative to view + rocket.pos(self.x(start), self.y(start), 0).add_trail(n=len(self.x[:, 1])) + + plt = Plotter(axes=1, interactive=False) + plt.show(world, rocket, __doc__, viewup="z", **kwargs) + + for t in np.arange(start, stop, time_step): + angle = np.arccos(2 * self.e0(t)**2 - 1) + k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 + + # Keep position static (relative start) to observe only rotation + rocket.pos(self.x(start), self.y(start), 0) + rocket.rotate_x(self.e1(t) / k) + rocket.rotate_y(self.e2(t) / k) + rocket.rotate_z(self.e3(t) / k) + + plt.show(world, rocket) + + if getattr(plt, 'escaped', False): + break + + plt.interactive().close() + return None + def info(self): """Prints out a summary of the data available about the Flight.""" self.prints.all() From 51ca422f6ac96289c967d043af5a9dd2adc9868a Mon Sep 17 00:00:00 2001 From: GuilhermeAsura Date: Sat, 6 Dec 2025 00:34:09 -0300 Subject: [PATCH 03/12] test: add modular verification suite for 3D animation methods --- tests/animation_verification/__init__.py | 0 tests/animation_verification/main.py | 71 ++++++ tests/animation_verification/rocket_setup.py | 82 +++++++ tests/animation_verification/rocket_stl.py | 220 +++++++++++++++++++ 4 files changed, 373 insertions(+) create mode 100644 tests/animation_verification/__init__.py create mode 100644 tests/animation_verification/main.py create mode 100644 tests/animation_verification/rocket_setup.py create mode 100644 tests/animation_verification/rocket_stl.py diff --git a/tests/animation_verification/__init__.py b/tests/animation_verification/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/animation_verification/main.py b/tests/animation_verification/main.py new file mode 100644 index 000000000..7779b3f76 --- /dev/null +++ b/tests/animation_verification/main.py @@ -0,0 +1,71 @@ +import os +import traceback +from rocketpy import Environment, Flight +from rocket_stl import create_rocket_stl +from rocket_setup import get_calisto_rocket + +def run_simulation_and_test_animation(): + print("🚀 Setting up simulation (Calisto Example)...") + + # 1. Setup Environment + env = Environment(latitude=32.990254, longitude=-106.974998, elevation=1400) + env.set_date((2025, 12, 5, 12)) + env.set_atmospheric_model(type="standard_atmosphere") + + # 2. Get Rocket + try: + calisto = get_calisto_rocket() + except Exception as e: + print(f"❌ Failed to configure rocket: {e}") + return + + # 3. Simulate Flight + test_flight = Flight( + rocket=calisto, environment=env, rail_length=5.2, inclination=85, heading=0 + ) + + print(f"✅ Flight simulated successfully! Apogee: {test_flight.apogee:.2f} m") + + # 4. Test Animation Methods + stl_file = "rocket_model.stl" + # Note: Depending on where you run this, you might need to adjust imports + # or ensure create_rocket_stl is available in scope. + create_rocket_stl(stl_file, length=300, radius=50) + + print("\n🎥 Testing animate_trajectory()...") + + try: + test_flight.animate_trajectory( + file_name=stl_file, + stop=15.0, + time_step=0.05, + azimuth=-45, # Rotates view 45 degrees left + elevation=30, # Tilts view 30 degrees up + zoom=1.2 + ) + print("✅ animate_trajectory() executed successfully.") + except Exception as e: + print(f"❌ animate_trajectory() Failed: {e}") + traceback.print_exc() + + print("\n🔄 Testing animate_rotate()...") + + try: + test_flight.animate_rotate( + file_name=stl_file, + time_step=1.0, + azimuth=-45, # Rotates view 45 degrees left + elevation=30, # Tilts view 30 degrees up + zoom=1.2 + ) + print("✅ animate_rotate() executed successfully.") + except Exception as e: + print(f"❌ animate_rotate() Failed: {e}") + traceback.print_exc() + + # Cleanup + if os.path.exists(stl_file): + os.remove(stl_file) + +if __name__ == "__main__": + run_simulation_and_test_animation() \ No newline at end of file diff --git a/tests/animation_verification/rocket_setup.py b/tests/animation_verification/rocket_setup.py new file mode 100644 index 000000000..0dc337728 --- /dev/null +++ b/tests/animation_verification/rocket_setup.py @@ -0,0 +1,82 @@ +import os +from rocketpy import SolidMotor, Rocket + +CURRENT_DIR = os.path.dirname(os.path.abspath(__file__)) +ROOT_DIR = os.path.dirname(os.path.dirname(CURRENT_DIR)) +DATA_DIR = os.path.join(ROOT_DIR, "data") + +OFF_DRAG_PATH = os.path.join(DATA_DIR, "rockets/calisto/powerOffDragCurve.csv") +ON_DRAG_PATH = os.path.join(DATA_DIR, "rockets/calisto/powerOnDragCurve.csv") +AIRFOIL_PATH = os.path.join(DATA_DIR, "airfoils/NACA0012-radians.txt") +MOTOR_PATH = os.path.join(DATA_DIR, "motors/cesaroni/Cesaroni_M1670.eng") + +def get_motor(): + """Locates the motor file and returns a configured SolidMotor object.""" + # We can now point directly to the file without searching + if not os.path.exists(MOTOR_PATH): + raise FileNotFoundError(f"Could not find Cesaroni_M1670.eng at: {MOTOR_PATH}") + + return SolidMotor( + thrust_source=MOTOR_PATH, + dry_mass=1.815, + dry_inertia=(0.125, 0.125, 0.002), + nozzle_radius=33 / 1000, + grain_number=5, + grain_density=1815, + grain_outer_radius=33 / 1000, + grain_initial_inner_radius=15 / 1000, + grain_initial_height=120 / 1000, + grain_separation=5 / 1000, + grains_center_of_mass_position=0.397, + center_of_dry_mass_position=0.317, + nozzle_position=0, + burn_time=3.9, + throat_radius=11 / 1000, + coordinate_system_orientation="nozzle_to_combustion_chamber", + ) + +def get_calisto_rocket(): + """Configures and returns the Calisto Rocket object.""" + motor = get_motor() + + calisto = Rocket( + radius=127 / 2000, + mass=14.426, + inertia=(6.321, 6.321, 0.034), + power_off_drag=OFF_DRAG_PATH, + power_on_drag=ON_DRAG_PATH, + center_of_mass_without_motor=0, + coordinate_system_orientation="tail_to_nose", + ) + + calisto.add_motor(motor, position=-1.255) + calisto.set_rail_buttons( + upper_button_position=0.0818, + lower_button_position=-0.618, + angular_position=45, + ) + + # Aerodynamic surfaces + calisto.add_nose(length=0.55829, kind="vonKarman", position=1.27) + calisto.add_trapezoidal_fins( + n=4, + root_chord=0.120, + tip_chord=0.060, + span=0.110, + position=-1.04956, + cant_angle=0, + airfoil=(AIRFOIL_PATH, "radians"), + ) + calisto.add_tail( + top_radius=0.0635, bottom_radius=0.0435, length=0.060, position=-1.194656 + ) + + # Parachutes + calisto.add_parachute( + name="Main", cd_s=10.0, trigger=800, sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5) + ) + calisto.add_parachute( + name="Drogue", cd_s=1.0, trigger="apogee", sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5) + ) + + return calisto \ No newline at end of file diff --git a/tests/animation_verification/rocket_stl.py b/tests/animation_verification/rocket_stl.py new file mode 100644 index 000000000..7f91678c0 --- /dev/null +++ b/tests/animation_verification/rocket_stl.py @@ -0,0 +1,220 @@ +import math +import numpy as np +from stl import mesh # Requires numpy-stl package: pip install numpy-stl +import struct + +def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): + """ + Creates a detailed rocket-shaped STL file with proper scale. + + Parameters + ---------- + filename : str + Output filename + length : float + Rocket length in meters + radius : float + Rocket base radius in meters + """ + + # More realistic proportions + nose_length = length * 0.25 + body_length = length * 0.65 + engine_length = length * 0.1 + + # Fin parameters + fin_height = radius * 3 + fin_width = radius * 2 + fin_thickness = radius * 0.2 + num_fins = 4 + + # Calculate number of vertices for smooth curves + num_segments = 36 + vertices = [] + faces = [] + + # Helper function to add triangle + def add_triangle(v1, v2, v3): + nonlocal faces + idx = len(vertices) + vertices.extend([v1, v2, v3]) + faces.append([idx, idx+1, idx+2]) + + # Helper function to add quadrilateral as two triangles + def add_quad(v1, v2, v3, v4): + add_triangle(v1, v2, v3) + add_triangle(v1, v3, v4) + + # 1. Create nose cone (pointed ogive) + nose_vertices = [] + for i in range(num_segments + 1): + angle = 2 * math.pi * i / num_segments + x = math.cos(angle) * radius + y = math.sin(angle) * radius + z = length # Tip at full length + nose_vertices.append([x, y, z]) + + # Create nose cone triangles + nose_tip = [0, 0, length] + for i in range(num_segments): + v1 = nose_tip + v2 = nose_vertices[i] + v3 = nose_vertices[(i + 1) % num_segments] + add_triangle(v1, v2, v3) + + # 2. Create main body cylinder + body_z_start = length - nose_length + body_z_end = engine_length + + # Create vertices for top and bottom circles of body + top_circle = [] + bottom_circle = [] + + for i in range(num_segments): + angle = 2 * math.pi * i / num_segments + x = math.cos(angle) * radius + y = math.sin(angle) * radius + + top_circle.append([x, y, body_z_start]) + bottom_circle.append([x, y, body_z_end]) + + # Create body cylinder triangles + for i in range(num_segments): + next_i = (i + 1) % num_segments + + # Side quad + v1 = top_circle[i] + v2 = top_circle[next_i] + v3 = bottom_circle[next_i] + v4 = bottom_circle[i] + add_quad(v1, v2, v3, v4) + + # Top circle (connecting to nose) + v1 = top_circle[i] + v2 = top_circle[next_i] + v3 = nose_vertices[i] + add_triangle(v1, v2, v3) + + # 3. Create engine nozzle (truncated cone) + engine_radius = radius * 0.7 + engine_z_end = 0 + + # Create vertices for engine circles + engine_top_circle = [] + engine_bottom_circle = [] + + for i in range(num_segments): + angle = 2 * math.pi * i / num_segments + x_top = math.cos(angle) * radius + y_top = math.sin(angle) * radius + x_bottom = math.cos(angle) * engine_radius + y_bottom = math.sin(angle) * engine_radius + + engine_top_circle.append([x_top, y_top, body_z_end]) + engine_bottom_circle.append([x_bottom, y_bottom, engine_z_end]) + + # Create engine nozzle triangles + for i in range(num_segments): + next_i = (i + 1) % num_segments + + # Side quad + v1 = engine_top_circle[i] + v2 = engine_top_circle[next_i] + v3 = engine_bottom_circle[next_i] + v4 = engine_bottom_circle[i] + add_quad(v1, v2, v3, v4) + + # Connect to body + v1 = engine_top_circle[i] + v2 = engine_top_circle[next_i] + v3 = bottom_circle[i] + add_triangle(v1, v2, v3) + + # 4. Create fins + for fin_num in range(num_fins): + fin_angle = 2 * math.pi * fin_num / num_fins + + # Fin vertices + fin_base_center = [0, 0, body_z_end * 0.5] + + # Inner fin edge (attached to rocket) + inner_x = math.cos(fin_angle) * radius + inner_y = math.sin(fin_angle) * radius + inner_top = [inner_x, inner_y, body_z_end + fin_height * 0.3] + inner_bottom = [inner_x, inner_y, body_z_end - fin_height * 0.7] + + # Outer fin edge + outer_x = math.cos(fin_angle) * (radius + fin_width) + outer_y = math.sin(fin_angle) * (radius + fin_width) + outer_top = [outer_x, outer_y, body_z_end + fin_height * 0.3] + outer_bottom = [outer_x, outer_y, body_z_end - fin_height * 0.7] + + # Fin side that attaches to rocket + add_quad(inner_top, inner_bottom, bottom_circle[fin_num * num_segments // num_fins], + bottom_circle[(fin_num * num_segments // num_fins + 1) % num_segments]) + + # Fin surfaces + # Front face + add_quad(inner_top, outer_top, outer_bottom, inner_bottom) + + # Top face + add_triangle(inner_top, outer_top, + [inner_top[0] + math.cos(fin_angle + math.pi/2) * fin_thickness * 0.5, + inner_top[1] + math.sin(fin_angle + math.pi/2) * fin_thickness * 0.5, + inner_top[2]]) + + # Bottom face + add_triangle(inner_bottom, outer_bottom, + [inner_bottom[0] + math.cos(fin_angle + math.pi/2) * fin_thickness * 0.5, + inner_bottom[1] + math.sin(fin_angle + math.pi/2) * fin_thickness * 0.5, + inner_bottom[2]]) + + # 5. Create bottom cap + center_bottom = [0, 0, engine_z_end] + for i in range(num_segments): + next_i = (i + 1) % num_segments + add_triangle(center_bottom, engine_bottom_circle[i], engine_bottom_circle[next_i]) + + # Convert to numpy arrays + vertices_array = np.array(vertices) + faces_array = np.array(faces) + + # Create the mesh + rocket_mesh = mesh.Mesh(np.zeros(faces_array.shape[0], dtype=mesh.Mesh.dtype)) + for i, face in enumerate(faces_array): + for j in range(3): + rocket_mesh.vectors[i][j] = vertices_array[face[j]] + + # Write the mesh to file + rocket_mesh.save(filename) + print(f"Rocket STL saved to {filename}") + print(f"Total triangles: {len(faces)}") + +# Alternative version using pure Python STL generation +def create_rocket_stl_simple(filename="rocket_simple.stl", length=100, radius=10): + """Simpler version using pure Python without numpy-stl dependency""" + + def write_stl(faces, filename): + with open(filename, 'wb') as f: + # Write 80 byte header + f.write(b'\x00' * 80) + # Write number of faces + f.write(struct.pack(' Date: Sat, 6 Dec 2025 15:33:50 -0300 Subject: [PATCH 04/12] maint: fix unused variables and pin requests>=2.25.0 for CI stability --- requirements.txt | 2 +- tests/animation_verification/rocket_stl.py | 8 -------- 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/requirements.txt b/requirements.txt index 61a594320..a4fe9cf26 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ numpy>=1.13 scipy>=1.0 matplotlib>=3.9.0 # Released May 15th 2024 netCDF4>=1.6.4 -requests +requests>=2.25.0 pytz simplekml dill diff --git a/tests/animation_verification/rocket_stl.py b/tests/animation_verification/rocket_stl.py index 7f91678c0..5b10ccf1b 100644 --- a/tests/animation_verification/rocket_stl.py +++ b/tests/animation_verification/rocket_stl.py @@ -19,7 +19,6 @@ def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): # More realistic proportions nose_length = length * 0.25 - body_length = length * 0.65 engine_length = length * 0.1 # Fin parameters @@ -134,9 +133,6 @@ def add_quad(v1, v2, v3, v4): for fin_num in range(num_fins): fin_angle = 2 * math.pi * fin_num / num_fins - # Fin vertices - fin_base_center = [0, 0, body_z_end * 0.5] - # Inner fin edge (attached to rocket) inner_x = math.cos(fin_angle) * radius inner_y = math.sin(fin_angle) * radius @@ -211,10 +207,6 @@ def write_stl(faces, filename): f.write(struct.pack(' Date: Sat, 6 Dec 2025 23:36:04 -0300 Subject: [PATCH 05/12] style: format verification scripts and pin pytz>=2020.1 for CI stability --- requirements.txt | 2 +- tests/animation_verification/main.py | 28 ++-- tests/animation_verification/rocket_setup.py | 22 +++- tests/animation_verification/rocket_stl.py | 130 +++++++++++-------- 4 files changed, 108 insertions(+), 74 deletions(-) diff --git a/requirements.txt b/requirements.txt index a4fe9cf26..e08c9a81d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,6 +3,6 @@ scipy>=1.0 matplotlib>=3.9.0 # Released May 15th 2024 netCDF4>=1.6.4 requests>=2.25.0 -pytz +pytz>=2020.1 simplekml dill diff --git a/tests/animation_verification/main.py b/tests/animation_verification/main.py index 7779b3f76..5538bd7fe 100644 --- a/tests/animation_verification/main.py +++ b/tests/animation_verification/main.py @@ -4,9 +4,10 @@ from rocket_stl import create_rocket_stl from rocket_setup import get_calisto_rocket + def run_simulation_and_test_animation(): print("🚀 Setting up simulation (Calisto Example)...") - + # 1. Setup Environment env = Environment(latitude=32.990254, longitude=-106.974998, elevation=1400) env.set_date((2025, 12, 5, 12)) @@ -23,40 +24,40 @@ def run_simulation_and_test_animation(): test_flight = Flight( rocket=calisto, environment=env, rail_length=5.2, inclination=85, heading=0 ) - + print(f"✅ Flight simulated successfully! Apogee: {test_flight.apogee:.2f} m") # 4. Test Animation Methods stl_file = "rocket_model.stl" - # Note: Depending on where you run this, you might need to adjust imports + # Note: Depending on where you run this, you might need to adjust imports # or ensure create_rocket_stl is available in scope. create_rocket_stl(stl_file, length=300, radius=50) print("\n🎥 Testing animate_trajectory()...") - + try: test_flight.animate_trajectory( - file_name=stl_file, + file_name=stl_file, stop=15.0, time_step=0.05, - azimuth=-45, # Rotates view 45 degrees left + azimuth=-45, # Rotates view 45 degrees left elevation=30, # Tilts view 30 degrees up - zoom=1.2 - ) + zoom=1.2, + ) print("✅ animate_trajectory() executed successfully.") except Exception as e: print(f"❌ animate_trajectory() Failed: {e}") traceback.print_exc() print("\n🔄 Testing animate_rotate()...") - + try: test_flight.animate_rotate( - file_name=stl_file, + file_name=stl_file, time_step=1.0, - azimuth=-45, # Rotates view 45 degrees left + azimuth=-45, # Rotates view 45 degrees left elevation=30, # Tilts view 30 degrees up - zoom=1.2 + zoom=1.2, ) print("✅ animate_rotate() executed successfully.") except Exception as e: @@ -67,5 +68,6 @@ def run_simulation_and_test_animation(): if os.path.exists(stl_file): os.remove(stl_file) + if __name__ == "__main__": - run_simulation_and_test_animation() \ No newline at end of file + run_simulation_and_test_animation() diff --git a/tests/animation_verification/rocket_setup.py b/tests/animation_verification/rocket_setup.py index 0dc337728..4f3c11ba1 100644 --- a/tests/animation_verification/rocket_setup.py +++ b/tests/animation_verification/rocket_setup.py @@ -10,6 +10,7 @@ AIRFOIL_PATH = os.path.join(DATA_DIR, "airfoils/NACA0012-radians.txt") MOTOR_PATH = os.path.join(DATA_DIR, "motors/cesaroni/Cesaroni_M1670.eng") + def get_motor(): """Locates the motor file and returns a configured SolidMotor object.""" # We can now point directly to the file without searching @@ -35,10 +36,11 @@ def get_motor(): coordinate_system_orientation="nozzle_to_combustion_chamber", ) + def get_calisto_rocket(): """Configures and returns the Calisto Rocket object.""" motor = get_motor() - + calisto = Rocket( radius=127 / 2000, mass=14.426, @@ -73,10 +75,20 @@ def get_calisto_rocket(): # Parachutes calisto.add_parachute( - name="Main", cd_s=10.0, trigger=800, sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5) + name="Main", + cd_s=10.0, + trigger=800, + sampling_rate=105, + lag=1.5, + noise=(0, 8.3, 0.5), ) calisto.add_parachute( - name="Drogue", cd_s=1.0, trigger="apogee", sampling_rate=105, lag=1.5, noise=(0, 8.3, 0.5) + name="Drogue", + cd_s=1.0, + trigger="apogee", + sampling_rate=105, + lag=1.5, + noise=(0, 8.3, 0.5), ) - - return calisto \ No newline at end of file + + return calisto diff --git a/tests/animation_verification/rocket_stl.py b/tests/animation_verification/rocket_stl.py index 5b10ccf1b..700285d07 100644 --- a/tests/animation_verification/rocket_stl.py +++ b/tests/animation_verification/rocket_stl.py @@ -3,10 +3,11 @@ from stl import mesh # Requires numpy-stl package: pip install numpy-stl import struct + def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): """ Creates a detailed rocket-shaped STL file with proper scale. - + Parameters ---------- filename : str @@ -16,34 +17,34 @@ def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): radius : float Rocket base radius in meters """ - + # More realistic proportions nose_length = length * 0.25 engine_length = length * 0.1 - + # Fin parameters fin_height = radius * 3 fin_width = radius * 2 fin_thickness = radius * 0.2 num_fins = 4 - + # Calculate number of vertices for smooth curves num_segments = 36 vertices = [] faces = [] - + # Helper function to add triangle def add_triangle(v1, v2, v3): nonlocal faces idx = len(vertices) vertices.extend([v1, v2, v3]) - faces.append([idx, idx+1, idx+2]) - + faces.append([idx, idx + 1, idx + 2]) + # Helper function to add quadrilateral as two triangles def add_quad(v1, v2, v3, v4): add_triangle(v1, v2, v3) add_triangle(v1, v3, v4) - + # 1. Create nose cone (pointed ogive) nose_vertices = [] for i in range(num_segments + 1): @@ -52,7 +53,7 @@ def add_quad(v1, v2, v3, v4): y = math.sin(angle) * radius z = length # Tip at full length nose_vertices.append([x, y, z]) - + # Create nose cone triangles nose_tip = [0, 0, length] for i in range(num_segments): @@ -60,153 +61,172 @@ def add_quad(v1, v2, v3, v4): v2 = nose_vertices[i] v3 = nose_vertices[(i + 1) % num_segments] add_triangle(v1, v2, v3) - + # 2. Create main body cylinder body_z_start = length - nose_length body_z_end = engine_length - + # Create vertices for top and bottom circles of body top_circle = [] bottom_circle = [] - + for i in range(num_segments): angle = 2 * math.pi * i / num_segments x = math.cos(angle) * radius y = math.sin(angle) * radius - + top_circle.append([x, y, body_z_start]) bottom_circle.append([x, y, body_z_end]) - + # Create body cylinder triangles for i in range(num_segments): next_i = (i + 1) % num_segments - + # Side quad v1 = top_circle[i] v2 = top_circle[next_i] v3 = bottom_circle[next_i] v4 = bottom_circle[i] add_quad(v1, v2, v3, v4) - + # Top circle (connecting to nose) v1 = top_circle[i] v2 = top_circle[next_i] v3 = nose_vertices[i] add_triangle(v1, v2, v3) - + # 3. Create engine nozzle (truncated cone) engine_radius = radius * 0.7 engine_z_end = 0 - + # Create vertices for engine circles engine_top_circle = [] engine_bottom_circle = [] - + for i in range(num_segments): angle = 2 * math.pi * i / num_segments x_top = math.cos(angle) * radius y_top = math.sin(angle) * radius x_bottom = math.cos(angle) * engine_radius y_bottom = math.sin(angle) * engine_radius - + engine_top_circle.append([x_top, y_top, body_z_end]) engine_bottom_circle.append([x_bottom, y_bottom, engine_z_end]) - + # Create engine nozzle triangles for i in range(num_segments): next_i = (i + 1) % num_segments - + # Side quad v1 = engine_top_circle[i] v2 = engine_top_circle[next_i] v3 = engine_bottom_circle[next_i] v4 = engine_bottom_circle[i] add_quad(v1, v2, v3, v4) - + # Connect to body v1 = engine_top_circle[i] v2 = engine_top_circle[next_i] v3 = bottom_circle[i] add_triangle(v1, v2, v3) - + # 4. Create fins for fin_num in range(num_fins): fin_angle = 2 * math.pi * fin_num / num_fins - + # Inner fin edge (attached to rocket) inner_x = math.cos(fin_angle) * radius inner_y = math.sin(fin_angle) * radius inner_top = [inner_x, inner_y, body_z_end + fin_height * 0.3] inner_bottom = [inner_x, inner_y, body_z_end - fin_height * 0.7] - + # Outer fin edge outer_x = math.cos(fin_angle) * (radius + fin_width) outer_y = math.sin(fin_angle) * (radius + fin_width) outer_top = [outer_x, outer_y, body_z_end + fin_height * 0.3] outer_bottom = [outer_x, outer_y, body_z_end - fin_height * 0.7] - + # Fin side that attaches to rocket - add_quad(inner_top, inner_bottom, bottom_circle[fin_num * num_segments // num_fins], - bottom_circle[(fin_num * num_segments // num_fins + 1) % num_segments]) - + add_quad( + inner_top, + inner_bottom, + bottom_circle[fin_num * num_segments // num_fins], + bottom_circle[(fin_num * num_segments // num_fins + 1) % num_segments], + ) + # Fin surfaces # Front face add_quad(inner_top, outer_top, outer_bottom, inner_bottom) - + # Top face - add_triangle(inner_top, outer_top, - [inner_top[0] + math.cos(fin_angle + math.pi/2) * fin_thickness * 0.5, - inner_top[1] + math.sin(fin_angle + math.pi/2) * fin_thickness * 0.5, - inner_top[2]]) - + add_triangle( + inner_top, + outer_top, + [ + inner_top[0] + math.cos(fin_angle + math.pi / 2) * fin_thickness * 0.5, + inner_top[1] + math.sin(fin_angle + math.pi / 2) * fin_thickness * 0.5, + inner_top[2], + ], + ) + # Bottom face - add_triangle(inner_bottom, outer_bottom, - [inner_bottom[0] + math.cos(fin_angle + math.pi/2) * fin_thickness * 0.5, - inner_bottom[1] + math.sin(fin_angle + math.pi/2) * fin_thickness * 0.5, - inner_bottom[2]]) - + add_triangle( + inner_bottom, + outer_bottom, + [ + inner_bottom[0] + + math.cos(fin_angle + math.pi / 2) * fin_thickness * 0.5, + inner_bottom[1] + + math.sin(fin_angle + math.pi / 2) * fin_thickness * 0.5, + inner_bottom[2], + ], + ) + # 5. Create bottom cap center_bottom = [0, 0, engine_z_end] for i in range(num_segments): next_i = (i + 1) % num_segments - add_triangle(center_bottom, engine_bottom_circle[i], engine_bottom_circle[next_i]) - + add_triangle( + center_bottom, engine_bottom_circle[i], engine_bottom_circle[next_i] + ) + # Convert to numpy arrays vertices_array = np.array(vertices) faces_array = np.array(faces) - + # Create the mesh rocket_mesh = mesh.Mesh(np.zeros(faces_array.shape[0], dtype=mesh.Mesh.dtype)) for i, face in enumerate(faces_array): for j in range(3): rocket_mesh.vectors[i][j] = vertices_array[face[j]] - + # Write the mesh to file rocket_mesh.save(filename) print(f"Rocket STL saved to {filename}") print(f"Total triangles: {len(faces)}") + # Alternative version using pure Python STL generation def create_rocket_stl_simple(filename="rocket_simple.stl", length=100, radius=10): """Simpler version using pure Python without numpy-stl dependency""" - + def write_stl(faces, filename): - with open(filename, 'wb') as f: + with open(filename, "wb") as f: # Write 80 byte header - f.write(b'\x00' * 80) + f.write(b"\x00" * 80) # Write number of faces - f.write(struct.pack(' Date: Sat, 10 Jan 2026 22:21:40 -0300 Subject: [PATCH 06/12] make format --- rocketpy/simulation/flight.py | 41 +++++++++++--------- tests/animation_verification/main.py | 6 ++- tests/animation_verification/rocket_setup.py | 3 +- tests/animation_verification/rocket_stl.py | 3 +- 4 files changed, 31 insertions(+), 22 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index afa95d596..b63d3be13 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -4123,7 +4123,10 @@ def export_kml( color=color, altitude_mode=altitude_mode, ) - def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwargs): + + def animate_trajectory( + self, file_name, start=0, stop=None, time_step=0.1, **kwargs + ): """ 6-DOF Animation of the flight trajectory using Vedo. @@ -4155,7 +4158,7 @@ def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwa ------ ImportError If the 'vedo' package is not installed. - + Notes ----- This feature requires the 'vedo' package. Install it with: @@ -4197,8 +4200,10 @@ def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwa rocket = Mesh(file_name).c("green") rocket.pos(self.x(start), self.y(start), 0).add_trail(n=len(self.x[:, 1])) # Create trail - trail_points = [[self.x(t), self.y(t), self.z(t) - self.env.elevation] - for t in np.arange(start, stop, time_step)] + trail_points = [ + [self.x(t), self.y(t), self.z(t) - self.env.elevation] + for t in np.arange(start, stop, time_step) + ] trail = Line(trail_points, c="k", alpha=0.5) # Setup Plotter plt = Plotter(axes=1, interactive=False) @@ -4207,21 +4212,21 @@ def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwa # Animation Loop for t in np.arange(start, stop, time_step): # Calculate rotation angle and vector from quaternions - # Note: This simple rotation logic mimics the old branch. - # Ideally, vedo handles orientation via matrix, but we stick + # Note: This simple rotation logic mimics the old branch. + # Ideally, vedo handles orientation via matrix, but we stick # to the provided logic for now. - + # e0 is the scalar part of the quaternion - angle = np.arccos(2 * self.e0(t)**2 - 1) + angle = np.arccos(2 * self.e0(t) ** 2 - 1) k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 - + # Update position and rotation # Adjusting for ground elevation rocket.pos(self.x(t), self.y(t), self.z(t) - self.env.elevation) rocket.rotate_x(self.e1(t) / k) rocket.rotate_y(self.e2(t) / k) rocket.rotate_z(self.e3(t) / k) - + # update the scene plt.show(world, rocket, trail) @@ -4230,7 +4235,7 @@ def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwa while time.time() - start_pause < time_step: plt.render() - if getattr(plt, 'escaped', False): + if getattr(plt, "escaped", False): break plt.interactive().close() @@ -4258,7 +4263,7 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) - elevation (float): Rotation in degrees above the horizon. - roll (float): Rotation in degrees around the view axis. - zoom (float): Zoom level (default 1). - + Returns ------- None @@ -4276,13 +4281,13 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) "Install it with:\n" " pip install rocketpy[animation]\n" ) from e - + # Enable interaction if needed try: settings.allow_interaction = True except AttributeError: pass # Not available in newer versions of vedo - + if stop is None: stop = self.t_final @@ -4302,9 +4307,9 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) plt.show(world, rocket, __doc__, viewup="z", **kwargs) for t in np.arange(start, stop, time_step): - angle = np.arccos(2 * self.e0(t)**2 - 1) + angle = np.arccos(2 * self.e0(t) ** 2 - 1) k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 - + # Keep position static (relative start) to observe only rotation rocket.pos(self.x(start), self.y(start), 0) rocket.rotate_x(self.e1(t) / k) @@ -4313,12 +4318,12 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) plt.show(world, rocket) - if getattr(plt, 'escaped', False): + if getattr(plt, "escaped", False): break plt.interactive().close() return None - + def info(self): """Prints out a summary of the data available about the Flight.""" self.prints.all() diff --git a/tests/animation_verification/main.py b/tests/animation_verification/main.py index 5538bd7fe..15c3d42d2 100644 --- a/tests/animation_verification/main.py +++ b/tests/animation_verification/main.py @@ -1,8 +1,10 @@ import os import traceback -from rocketpy import Environment, Flight -from rocket_stl import create_rocket_stl + from rocket_setup import get_calisto_rocket +from rocket_stl import create_rocket_stl + +from rocketpy import Environment, Flight def run_simulation_and_test_animation(): diff --git a/tests/animation_verification/rocket_setup.py b/tests/animation_verification/rocket_setup.py index 4f3c11ba1..26a5640a7 100644 --- a/tests/animation_verification/rocket_setup.py +++ b/tests/animation_verification/rocket_setup.py @@ -1,5 +1,6 @@ import os -from rocketpy import SolidMotor, Rocket + +from rocketpy import Rocket, SolidMotor CURRENT_DIR = os.path.dirname(os.path.abspath(__file__)) ROOT_DIR = os.path.dirname(os.path.dirname(CURRENT_DIR)) diff --git a/tests/animation_verification/rocket_stl.py b/tests/animation_verification/rocket_stl.py index 700285d07..11654cc42 100644 --- a/tests/animation_verification/rocket_stl.py +++ b/tests/animation_verification/rocket_stl.py @@ -1,7 +1,8 @@ import math +import struct + import numpy as np from stl import mesh # Requires numpy-stl package: pip install numpy-stl -import struct def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): From 6b4cab83bb0c3bfb35bebfdfbca4179daa9a9fea Mon Sep 17 00:00:00 2001 From: GuilhermeAsura Date: Sat, 6 Dec 2025 00:33:38 -0300 Subject: [PATCH 07/12] feat: add animate_trajectory and animate_rotate to Flight class - Ported methods from legacy 'animate_flight' branch - Adapted to new Flight class structure (removed postProcess) - Added vedo as optional dependency Co-authored-by: Patrick Sampaio --- rocketpy/simulation/flight.py | 680 ++++++++++++++++++---------------- 1 file changed, 366 insertions(+), 314 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index b63d3be13..82f8c99f4 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -505,6 +505,7 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements equations_of_motion="standard", ode_solver="LSODA", simulation_mode="6 DOF", + weathercock_coeff=0.0, ): """Run a trajectory simulation. @@ -588,6 +589,16 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements A custom ``scipy.integrate.OdeSolver`` can be passed as well. For more information on the integration methods, see the scipy documentation [1]_. + weathercock_coeff : float, optional + Proportionality coefficient (rate coefficient) for the alignment rate of the rocket's body axis + with the relative wind direction in 3-DOF simulations, in rad/s. The actual angular velocity + applied to align the rocket is calculated as ``weathercock_coeff * sin(angle)``, where ``angle`` + is the angle between the rocket's axis and the wind direction. A higher value means faster alignment + (quasi-static weathercocking). This parameter is only used when simulation_mode is '3 DOF'. + Default is 0.0 to mimic a pure 3-DOF simulation without any weathercocking (fixed attitude). + Set to a positive value to enable quasi-static weathercocking behaviour. + + Returns ------- None @@ -617,6 +628,7 @@ def __init__( # pylint: disable=too-many-arguments,too-many-statements self.equations_of_motion = equations_of_motion self.simulation_mode = simulation_mode self.ode_solver = ode_solver + self.weathercock_coeff = weathercock_coeff # Controller initialization self.__init_controllers() @@ -699,6 +711,14 @@ def __simulate(self, verbose): self.__process_sensors_and_controllers_at_current_node(node, phase) + for controller in node._controllers: + controller( + self.t, + self.y_sol, + self.solution, + self.sensors, + ) + for parachute in node.parachutes: # Calculate and save pressure signal ( @@ -740,12 +760,11 @@ def __simulate(self, verbose): 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, - ) + lambda self, + added_mass_coefficient=parachute.added_mass_coefficient: setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, ), ] self.flight_phases.add_phase( @@ -845,7 +864,7 @@ def __process_sensors_and_controllers_at_current_node(self, node, phase): phase : FlightPhase The current flight phase. """ - if node._component_sensors: + if self.sensors: u_dot = phase.derivative(self.t, self.y_sol) self.__measure_sensors(node._component_sensors, u_dot) @@ -855,7 +874,6 @@ def __process_sensors_and_controllers_at_current_node(self, node, phase): self.y_sol, self.solution, self.sensors, - self.env, ) def __measure_sensors(self, component_sensors, u_dot, t=None, y_sol=None): @@ -964,12 +982,11 @@ def __check_and_handle_parachute_triggers( 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, - ) + lambda self, + added_mass_coefficient=parachute.added_mass_coefficient: setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, ), ] self.flight_phases.add_phase( @@ -1277,7 +1294,6 @@ def __process_overshootable_nodes(self, phase, phase_index, node_index): overshootable_node.y_sol, self.solution, self.sensors, - self.env, ) # Process sensors at overshootable node @@ -1372,12 +1388,11 @@ def __check_overshootable_parachute_triggers( 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, - ) + lambda self, + added_mass_coefficient=parachute.added_mass_coefficient: setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, ), ] self.flight_phases.add_phase( @@ -1683,28 +1698,89 @@ def lateral_surface_wind(self): return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad) - def __compute_drag_7d_inputs( - self, - stream_velocity_body, - stream_speed, - stream_mach, - density, - dynamic_viscosity, - ): - """Build drag-model inputs in the 7D order used by Rocket drag functions.""" - aerodynamic_stream_velocity = -stream_velocity_body - alpha = np.arctan2( - aerodynamic_stream_velocity[1], aerodynamic_stream_velocity[2] - ) - beta = np.arctan2( - aerodynamic_stream_velocity[0], aerodynamic_stream_velocity[2] - ) - reynolds = ( - density * stream_speed * (2 * self.rocket.radius) / dynamic_viscosity - if dynamic_viscosity > 0 - else 0 - ) - return alpha, beta, stream_mach, reynolds + def __get_drag_coefficient(self, drag_function, mach, z, freestream_velocity_body): + """Calculate drag coefficient, handling both 1D and multi-dimensional functions. + + Parameters + ---------- + drag_function : Function + The drag coefficient function (power_on_drag or power_off_drag) + mach : float + Mach number + z : float + Altitude in meters + freestream_velocity_body : Vector or array-like + Freestream velocity in body frame [stream_vx_b, stream_vy_b, stream_vz_b] + + Returns + ------- + float + Drag coefficient value + """ + # Early return for 1D drag functions (only mach number) + if not isinstance(drag_function, Function) or not getattr( + drag_function, "is_multidimensional", False + ): + return drag_function.get_value_opt(mach) + + # Multi-dimensional drag function - calculate additional parameters + + # Calculate Reynolds number: Re = rho * V * L / mu + # where L is characteristic length (rocket diameter) + rho = self.env.density.get_value_opt(z) + mu = self.env.dynamic_viscosity.get_value_opt(z) + freestream_speed = np.linalg.norm(freestream_velocity_body) + characteristic_length = 2 * self.rocket.radius # Diameter + # Defensive: avoid division by zero or non-finite viscosity values. + # Use a small epsilon fallback if `mu` is zero, negative, NaN or infinite. + try: + mu_val = float(mu) + except (TypeError, ValueError, OverflowError): + # Only catch errors related to invalid numeric conversion. + # Avoid catching broad Exception to satisfy linters and + # allow other unexpected errors to surface. + mu_val = 0.0 + if not np.isfinite(mu_val) or mu_val <= 0.0: + mu_safe = 1e-10 + else: + mu_safe = mu_val + + reynolds = rho * freestream_speed * characteristic_length / mu_safe + + # Calculate angle of attack + # Angle between freestream velocity and rocket axis (z-axis in body frame) + # The z component of freestream velocity in body frame + if hasattr(freestream_velocity_body, "z"): + stream_vz_b = -freestream_velocity_body.z + else: + stream_vz_b = -freestream_velocity_body[2] + + # Normalize and calculate angle + if freestream_speed > 1e-6: + cos_alpha = stream_vz_b / freestream_speed + # Clamp to [-1, 1] to avoid numerical issues + cos_alpha = np.clip(cos_alpha, -1.0, 1.0) + alpha_rad = np.arccos(cos_alpha) + alpha_deg = np.rad2deg(alpha_rad) + else: + alpha_deg = 0.0 + + # Determine which parameters to pass based on input names + input_names = [name.lower() for name in drag_function.__inputs__] + args = [] + + for name in input_names: + if "mach" in name or name == "m": + args.append(mach) + elif "reynolds" in name or name == "re": + args.append(reynolds) + elif "alpha" in name or name == "a" or "attack" in name: + args.append(alpha_deg) + else: + # Unknown parameter, default to mach + args.append(mach) + + return drag_function.get_value_opt(*args) def udot_rail1(self, t, u, post_processing=False): """Calculates derivative of u state vector with respect to time @@ -1736,28 +1812,38 @@ def udot_rail1(self, t, u, post_processing=False): total_mass_at_t = self.rocket.total_mass.get_value_opt(t) # Get freestream speed - free_stream_velocity = Vector( - [ - self.env.wind_velocity_x.get_value_opt(z) - vx, - self.env.wind_velocity_y.get_value_opt(z) - vy, - -vz, - ] - ) - free_stream_speed = abs(free_stream_velocity) + free_stream_speed = ( + (self.env.wind_velocity_x.get_value_opt(z) - vx) ** 2 + + (self.env.wind_velocity_y.get_value_opt(z) - vy) ** 2 + + (vz) ** 2 + ) ** 0.5 free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z) - rho = self.env.density.get_value_opt(z) - stream_velocity_body = ( - Matrix.transformation([e0, e1, e2, e3]).transpose @ free_stream_velocity - ) - dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) - alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( - stream_velocity_body, - free_stream_speed, + + # For rail motion, rocket is constrained - velocity mostly along z-axis in body frame + # Calculate velocity in body frame (simplified for rail) + a11 = 1 - 2 * (e2**2 + e3**2) + a12 = 2 * (e1 * e2 - e0 * e3) + a13 = 2 * (e1 * e3 + e0 * e2) + a21 = 2 * (e1 * e2 + e0 * e3) + a22 = 1 - 2 * (e1**2 + e3**2) + a23 = 2 * (e2 * e3 - e0 * e1) + a31 = 2 * (e1 * e3 - e0 * e2) + a32 = 2 * (e2 * e3 + e0 * e1) + a33 = 1 - 2 * (e1**2 + e2**2) + + # Freestream velocity in body frame + wind_vx = self.env.wind_velocity_x.get_value_opt(z) + wind_vy = self.env.wind_velocity_y.get_value_opt(z) + stream_vx_b = a11 * (wind_vx - vx) + a21 * (wind_vy - vy) + a31 * (-vz) + stream_vy_b = a12 * (wind_vx - vx) + a22 * (wind_vy - vy) + a32 * (-vz) + stream_vz_b = a13 * (wind_vx - vx) + a23 * (wind_vy - vy) + a33 * (-vz) + + drag_coeff = self.__get_drag_coefficient( + self.rocket.power_on_drag, free_stream_mach, - rho, - dynamic_viscosity, + z, + [stream_vx_b, stream_vy_b, stream_vz_b], ) - drag_coeff = self.rocket.power_on_drag_7d(alpha, beta, mach, reynolds, 0, 0, 0) # Calculate Forces pressure = self.env.pressure.get_value_opt(z) @@ -1766,6 +1852,7 @@ def udot_rail1(self, t, u, post_processing=False): + self.rocket.motor.pressure_thrust(pressure), 0, ) + rho = self.env.density.get_value_opt(z) R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff) # Calculate Linear acceleration @@ -1921,42 +2008,44 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) speed_of_sound = self.env.speed_of_sound.get_value_opt(z) - free_stream_velocity = Vector([wind_velocity_x - vx, wind_velocity_y - vy, -vz]) - free_stream_speed = abs(free_stream_velocity) + free_stream_speed = ( + (wind_velocity_x - vx) ** 2 + (wind_velocity_y - vy) ** 2 + (vz) ** 2 + ) ** 0.5 free_stream_mach = free_stream_speed / speed_of_sound - stream_velocity_body = Kt @ free_stream_velocity + + # Get rocket velocity in body frame (needed for drag calculation) + vx_b = a11 * vx + a21 * vy + a31 * vz + vy_b = a12 * vx + a22 * vy + a32 * vz + vz_b = a13 * vx + a23 * vy + a33 * vz + + # Calculate freestream velocity in body frame + stream_vx_b = ( + a11 * (wind_velocity_x - vx) + a21 * (wind_velocity_y - vy) + a31 * (-vz) + ) + stream_vy_b = ( + a12 * (wind_velocity_x - vx) + a22 * (wind_velocity_y - vy) + a32 * (-vz) + ) + stream_vz_b = ( + a13 * (wind_velocity_x - vx) + a23 * (wind_velocity_y - vy) + a33 * (-vz) + ) # Determine aerodynamics forces # Determine Drag Force - rho = self.env.density.get_value_opt(z) - dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) - alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( - stream_velocity_body, - free_stream_speed, - free_stream_mach, - rho, - dynamic_viscosity, - ) if t < self.rocket.motor.burn_out_time: - drag_coeff = self.rocket.power_on_drag_7d( - alpha, - beta, - mach, - reynolds, - omega1, - omega2, - omega3, + drag_coeff = self.__get_drag_coefficient( + self.rocket.power_on_drag, + free_stream_mach, + z, + [stream_vx_b, stream_vy_b, stream_vz_b], ) else: - drag_coeff = self.rocket.power_off_drag_7d( - alpha, - beta, - mach, - reynolds, - omega1, - omega2, - omega3, + drag_coeff = self.__get_drag_coefficient( + self.rocket.power_off_drag, + free_stream_mach, + z, + [stream_vx_b, stream_vy_b, stream_vz_b], ) + rho = self.env.density.get_value_opt(z) R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff for air_brakes in self.rocket.air_brakes: if air_brakes.deployment_level > 0: @@ -1977,10 +2066,6 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Off center moment M1 += self.rocket.cp_eccentricity_y * R3 M2 -= self.rocket.cp_eccentricity_x * R3 - # Get rocket velocity in body frame - vx_b = a11 * vx + a21 * vy + a31 * vz - vy_b = a12 * vx + a22 * vy + a32 * vz - vz_b = a13 * vx + a23 * vy + a33 * vz # Calculate lift and moment for each component of the rocket velocity_in_body_frame = Vector([vx_b, vy_b, vz_b]) w = Vector([omega1, omega2, omega3]) @@ -2001,13 +2086,11 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Reynolds at component altitude # TODO: Reynolds is only used in generic surfaces. This calculation # should be moved to the surface class for efficiency - comp_density = self.env.density.get_value_opt(comp_z) - comp_dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(comp_z) comp_reynolds = ( - comp_density + self.env.density.get_value_opt(comp_z) * comp_stream_speed * aero_surface.reference_length - / comp_dynamic_viscosity + / self.env.dynamic_viscosity.get_value_opt(comp_z) ) # Forces and moments X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( @@ -2201,25 +2284,12 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): 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 - stream_velocity_body = Kt @ free_stream_velocity - dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) - alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( - stream_velocity_body, - free_stream_speed, - mach, - rho, - dynamic_viscosity, - ) # Drag computation if t < self.rocket.motor.burn_out_time: - cd = self.rocket.power_on_drag_7d( - alpha, beta, mach, reynolds, omega1, omega2, omega3 - ) + cd = self.rocket.power_on_drag.get_value_opt(mach) else: - cd = self.rocket.power_off_drag_7d( - alpha, beta, mach, reynolds, omega1, omega2, omega3 - ) + 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 @@ -2253,13 +2323,11 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): rel_speed = abs(rel_velocity) rel_mach = rel_speed / speed_of_sound - comp_density = self.env.density.get_value_opt(comp_z) - comp_dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(comp_z) reynolds = ( - comp_density + self.env.density.get_value_opt(comp_z) * rel_speed * surface.reference_length - / comp_dynamic_viscosity + / self.env.dynamic_viscosity.get_value_opt(comp_z) ) fx, fy, fz, *_ = surface.compute_forces_and_moments( @@ -2290,8 +2358,7 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): r_dot = [vx, vy, vz] # Weathercocking: evolve body axis direction toward relative wind # The body z-axis (attitude vector) should align with -freestream_velocity - weathercock_coeff = getattr(self.rocket, "weathercock_coeff", 0.0) - if weathercock_coeff > 0 and free_stream_speed > 1e-6: + if self.weathercock_coeff > 0 and free_stream_speed > 1e-6: # Current body z-axis in inertial frame (attitude vector) # From rotation matrix: column 3 gives the body z-axis in inertial frame body_z_inertial = Vector( @@ -2323,7 +2390,7 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): sin_angle = min(1.0, max(-1.0, rotation_axis_mag)) # Angular velocity magnitude proportional to misalignment angle - omega_mag = weathercock_coeff * sin_angle + omega_mag = self.weathercock_coeff * sin_angle # Angular velocity in inertial frame, then transform to body frame omega_body = Kt @ (rotation_axis * omega_mag) @@ -2344,7 +2411,7 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): ) rotation_axis = perp_axis.unit_vector # 180 degree rotation: sin(angle) = 1 - omega_mag = weathercock_coeff * 1.0 + omega_mag = self.weathercock_coeff * 1.0 omega_body = Kt @ (rotation_axis * omega_mag) # else: aligned (dot > 0.999) - no rotation needed, omega_body stays None @@ -2445,19 +2512,15 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too 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_velocity = wind_velocity - v - free_stream_speed = abs(free_stream_velocity) + free_stream_speed = abs((wind_velocity - Vector(v))) speed_of_sound = self.env.speed_of_sound.get_value_opt(z) free_stream_mach = free_stream_speed / speed_of_sound - stream_velocity_body = Kt @ free_stream_velocity - dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) - alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( - stream_velocity_body, - free_stream_speed, - free_stream_mach, - rho, - dynamic_viscosity, - ) + + # Get rocket velocity in body frame (needed for drag calculation) + velocity_in_body_frame = Kt @ v + # Calculate freestream velocity in body frame + freestream_velocity = wind_velocity - v + freestream_velocity_body = Kt @ freestream_velocity if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time: pressure = self.env.pressure.get_value_opt(z) @@ -2466,25 +2529,19 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too + self.rocket.motor.pressure_thrust(pressure), 0, ) - drag_coeff = self.rocket.power_on_drag_7d( - alpha, - beta, - mach, - reynolds, - omega1, - omega2, - omega3, + drag_coeff = self.__get_drag_coefficient( + self.rocket.power_on_drag, + free_stream_mach, + z, + freestream_velocity_body, ) else: net_thrust = 0 - drag_coeff = self.rocket.power_off_drag_7d( - alpha, - beta, - mach, - reynolds, - omega1, - omega2, - omega3, + drag_coeff = self.__get_drag_coefficient( + self.rocket.power_off_drag, + free_stream_mach, + z, + freestream_velocity_body, ) R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff for air_brakes in self.rocket.air_brakes: @@ -2503,8 +2560,6 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too R3 = air_brakes_force # Substitutes rocket drag coefficient else: R3 += air_brakes_force - # Get rocket velocity in body frame - velocity_in_body_frame = Kt @ v # Calculate lift and moment for each component of the rocket for aero_surface, _ in self.rocket.aerodynamic_surfaces: # Component cp relative to CDM in body frame @@ -2523,13 +2578,11 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too # Reynolds at component altitude # TODO: Reynolds is only used in generic surfaces. This calculation # should be moved to the surface class for efficiency - comp_density = self.env.density.get_value_opt(comp_z) - comp_dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(comp_z) comp_reynolds = ( - comp_density + self.env.density.get_value_opt(comp_z) * comp_stream_speed * aero_surface.reference_length - / comp_dynamic_viscosity + / self.env.dynamic_viscosity.get_value_opt(comp_z) ) # Forces and moments X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( @@ -3107,7 +3160,7 @@ def max_acceleration_power_off_time(self): max_acceleration_time_index = np.argmax( self.acceleration[burn_out_time_index:, 1] ) - return self.acceleration[burn_out_time_index + max_acceleration_time_index, 0] + return self.acceleration[max_acceleration_time_index, 0] @cached_property def max_acceleration_power_off(self): @@ -3668,142 +3721,6 @@ def max_rail_button2_shear_force(self): """Maximum lower rail button shear force, in Newtons.""" return np.abs(self.rail_button2_shear_force.y_array).max() - @cached_property - def calculate_rail_button_bending_moments(self): - """Calculate internal bending moments at rail button attachment points. - - Uses beam theory to determine the internal structural moments for - stress analysis of the rail button attachments (fasteners and airframe). - - The bending moment at each button attachment consists of: - - 1. Normal force moment: $M = N \\times d$, where $N$ is the normal - reaction force and $d$ is the distance from button to center of - dry mass. - 2. Shear force cantilever moment: $M = S \\times h$, where $S$ is the - shear (tangential) force and $h$ is the button standoff height. - - Returns - ------- - tuple - rail_button1_bending_moment : Function - Bending moment at upper rail button as a function of time (N·m). - max_rail_button1_bending_moment : float - Maximum upper rail button bending moment (N·m). - rail_button2_bending_moment : Function - Bending moment at lower rail button as a function of time (N·m). - max_rail_button2_bending_moment : float - Maximum lower rail button bending moment (N·m). - - Notes - ----- - - Calculated only during the rail phase of flight - - Maximum values use absolute values for worst-case stress analysis - - The bending moments represent internal stresses in the rocket - airframe at the rail button attachment points - - **Assumptions:** - - - Rail buttons act as simple supports: provide reaction forces (normal - and shear) but no moment reaction at the rail contact point - - The rocket acts as a beam supported at two points (rail buttons) - - Bending moments arise from the lever arm effect of reaction forces - and the cantilever moment from button standoff height - """ - # Check if rail buttons exist - null_moment = Function(0) - if len(self.rocket.rail_buttons) == 0: - warnings.warn( - "Trying to calculate rail button bending moments without " - "rail buttons defined. Setting moments to zero.", - UserWarning, - ) - return (null_moment, 0.0, null_moment, 0.0) - - # Get rail button geometry - rail_buttons_tuple = self.rocket.rail_buttons[0] - # Rail button standoff height - h_button = rail_buttons_tuple.component.button_height - if h_button is None: - warnings.warn( - "Rail button height not defined. Bending moments cannot be " - "calculated. Setting moments to zero.", - UserWarning, - ) - return (null_moment, 0.0, null_moment, 0.0) - upper_button_position = ( - rail_buttons_tuple.component.buttons_distance - + rail_buttons_tuple.position.z - ) - lower_button_position = rail_buttons_tuple.position.z - - # Get center of dry mass (handle both callable and property) - if callable(self.rocket.center_of_dry_mass_position): - cdm = self.rocket.center_of_dry_mass_position(self.rocket._csys) - else: - cdm = self.rocket.center_of_dry_mass_position - - # Distances from buttons to center of dry mass - d1 = abs(upper_button_position - cdm) - d2 = abs(lower_button_position - cdm) - - # forces - N1 = self.rail_button1_normal_force - N2 = self.rail_button2_normal_force - S1 = self.rail_button1_shear_force - S2 = self.rail_button2_shear_force - t = N1.source[:, 0] - - # Calculate bending moments at attachment points - # Primary contribution from shear force acting at button height - # Secondary contribution from normal force creating moment about attachment - m1_values = N2.source[:, 1] * d2 + S1.source[:, 1] * h_button - m2_values = N1.source[:, 1] * d1 + S2.source[:, 1] * h_button - - rail_button1_bending_moment = Function( - np.column_stack([t, m1_values]), - inputs="Time (s)", - outputs="Bending Moment (N·m)", - interpolation="linear", - ) - rail_button2_bending_moment = Function( - np.column_stack([t, m2_values]), - inputs="Time (s)", - outputs="Bending Moment (N·m)", - interpolation="linear", - ) - - # Maximum bending moments (absolute value for stress calculations) - max_rail_button1_bending_moment = float(np.max(np.abs(m1_values))) - max_rail_button2_bending_moment = float(np.max(np.abs(m2_values))) - - return ( - rail_button1_bending_moment, - max_rail_button1_bending_moment, - rail_button2_bending_moment, - max_rail_button2_bending_moment, - ) - - @property - def rail_button1_bending_moment(self): - """Upper rail button bending moment as a Function of time.""" - return self.calculate_rail_button_bending_moments[0] - - @property - def max_rail_button1_bending_moment(self): - """Maximum upper rail button bending moment, in N·m.""" - return self.calculate_rail_button_bending_moments[1] - - @property - def rail_button2_bending_moment(self): - """Lower rail button bending moment as a Function of time.""" - return self.calculate_rail_button_bending_moments[2] - - @property - def max_rail_button2_bending_moment(self): - """Maximum lower rail button bending moment, in N·m.""" - return self.calculate_rail_button_bending_moments[3] - @funcify_method( "Time (s)", "Horizontal Distance to Launch Point (m)", "spline", "constant" ) @@ -4123,10 +4040,7 @@ def export_kml( color=color, altitude_mode=altitude_mode, ) - - def animate_trajectory( - self, file_name, start=0, stop=None, time_step=0.1, **kwargs - ): + def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwargs): """ 6-DOF Animation of the flight trajectory using Vedo. @@ -4158,7 +4072,7 @@ def animate_trajectory( ------ ImportError If the 'vedo' package is not installed. - + Notes ----- This feature requires the 'vedo' package. Install it with: @@ -4200,10 +4114,8 @@ def animate_trajectory( rocket = Mesh(file_name).c("green") rocket.pos(self.x(start), self.y(start), 0).add_trail(n=len(self.x[:, 1])) # Create trail - trail_points = [ - [self.x(t), self.y(t), self.z(t) - self.env.elevation] - for t in np.arange(start, stop, time_step) - ] + trail_points = [[self.x(t), self.y(t), self.z(t) - self.env.elevation] + for t in np.arange(start, stop, time_step)] trail = Line(trail_points, c="k", alpha=0.5) # Setup Plotter plt = Plotter(axes=1, interactive=False) @@ -4212,21 +4124,21 @@ def animate_trajectory( # Animation Loop for t in np.arange(start, stop, time_step): # Calculate rotation angle and vector from quaternions - # Note: This simple rotation logic mimics the old branch. - # Ideally, vedo handles orientation via matrix, but we stick + # Note: This simple rotation logic mimics the old branch. + # Ideally, vedo handles orientation via matrix, but we stick # to the provided logic for now. - + # e0 is the scalar part of the quaternion - angle = np.arccos(2 * self.e0(t) ** 2 - 1) + angle = np.arccos(2 * self.e0(t)**2 - 1) k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 - + # Update position and rotation # Adjusting for ground elevation rocket.pos(self.x(t), self.y(t), self.z(t) - self.env.elevation) rocket.rotate_x(self.e1(t) / k) rocket.rotate_y(self.e2(t) / k) rocket.rotate_z(self.e3(t) / k) - + # update the scene plt.show(world, rocket, trail) @@ -4235,7 +4147,7 @@ def animate_trajectory( while time.time() - start_pause < time_step: plt.render() - if getattr(plt, "escaped", False): + if getattr(plt, 'escaped', False): break plt.interactive().close() @@ -4263,7 +4175,7 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) - elevation (float): Rotation in degrees above the horizon. - roll (float): Rotation in degrees around the view axis. - zoom (float): Zoom level (default 1). - + Returns ------- None @@ -4281,13 +4193,13 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) "Install it with:\n" " pip install rocketpy[animation]\n" ) from e - + # Enable interaction if needed try: settings.allow_interaction = True except AttributeError: pass # Not available in newer versions of vedo - + if stop is None: stop = self.t_final @@ -4307,9 +4219,9 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) plt.show(world, rocket, __doc__, viewup="z", **kwargs) for t in np.arange(start, stop, time_step): - angle = np.arccos(2 * self.e0(t) ** 2 - 1) + angle = np.arccos(2 * self.e0(t)**2 - 1) k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 - + # Keep position static (relative start) to observe only rotation rocket.pos(self.x(start), self.y(start), 0) rocket.rotate_x(self.e1(t) / k) @@ -4318,12 +4230,12 @@ def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs) plt.show(world, rocket) - if getattr(plt, "escaped", False): + if getattr(plt, 'escaped', False): break plt.interactive().close() return None - + def info(self): """Prints out a summary of the data available about the Flight.""" self.prints.all() @@ -4734,6 +4646,7 @@ def merge(self): tmp_dict[time]._controllers += node._controllers tmp_dict[time].callbacks += node.callbacks tmp_dict[time]._component_sensors += node._component_sensors + tmp_dict[time]._controllers += node._controllers except KeyError: # If the node does not exist, add it to the dictionary tmp_dict[time] = node @@ -4799,3 +4712,142 @@ def __lt__(self, other): otherwise. """ return self.t < other.t + + @cached_property + def calculate_rail_button_bending_moments(self): + """ + Calculate internal bending moments at rail button attachment points. + + Uses beam theory to determine internal structural moments for stress + analysis of the rail button attachments (fasteners and airframe). + + The bending moment at each button attachment consists of: + 1. Bending from shear force at button contact point: M = S × h + where S is the shear (tangential) force and h is button height + 2. Direct moment contribution from the button's reaction forces + + Assumptions + ----------- + - Rail buttons act as simple supports: provide reaction forces (normal + and shear) but no moment reaction at the rail contact point. + - The rocket acts as a beam supported at two points (rail buttons). + - Bending moments arise from the lever arm effect of reaction forces + and the cantilever moment from button standoff height. + + The bending moment at each button attachment consists of: + 1. Normal force moment: M = N x d, where N is normal reaction force + and d is distance from button to center of dry mass + 2. Shear force cantilever moment: M = S x h, where S is shear force + and h is button standoff height + + Notes + ----- + - Calculated only during the rail phase of flight + - Maximum values use absolute values for worst-case stress analysis + - The bending moments represent internal stresses in the rocket + airframe at the rail button attachment points + + Returns + ------- + tuple + (rail_button1_bending_moment : Function, + max_rail_button1_bending_moment : float, + rail_button2_bending_moment : Function, + max_rail_button2_bending_moment : float) + + Where rail_button1/2_bending_moment are Function objects of time + in N·m, and max values are floats in N·m. + """ + # Check if rail buttons exist + null_moment = Function(0) + if len(self.rocket.rail_buttons) == 0: + warnings.warn( + "Trying to calculate rail button bending moments without " + "rail buttons defined. Setting moments to zero.", + UserWarning, + ) + return (null_moment, 0.0, null_moment, 0.0) + + # Get rail button geometry + rail_buttons_tuple = self.rocket.rail_buttons[0] + # Rail button standoff height + h_button = rail_buttons_tuple.component.button_height + if h_button is None: + warnings.warn( + "Rail button height not defined. Bending moments cannot be " + "calculated. Setting moments to zero.", + UserWarning, + ) + return (null_moment, 0.0, null_moment, 0.0) + upper_button_position = ( + rail_buttons_tuple.component.buttons_distance + + rail_buttons_tuple.position.z + ) + lower_button_position = rail_buttons_tuple.position.z + + # Get center of dry mass (handle both callable and property) + if callable(self.rocket.center_of_dry_mass_position): + cdm = self.rocket.center_of_dry_mass_position(self.rocket._csys) + else: + cdm = self.rocket.center_of_dry_mass_position + + # Distances from buttons to center of dry mass + d1 = abs(upper_button_position - cdm) + d2 = abs(lower_button_position - cdm) + + # forces + N1 = self.rail_button1_normal_force + N2 = self.rail_button2_normal_force + S1 = self.rail_button1_shear_force + S2 = self.rail_button2_shear_force + t = N1.source[:, 0] + + # Calculate bending moments at attachment points + # Primary contribution from shear force acting at button height + # Secondary contribution from normal force creating moment about attachment + m1_values = N2.source[:, 1] * d2 + S1.source[:, 1] * h_button + m2_values = N1.source[:, 1] * d1 + S2.source[:, 1] * h_button + + rail_button1_bending_moment = Function( + np.column_stack([t, m1_values]), + inputs="Time (s)", + outputs="Bending Moment (N·m)", + interpolation="linear", + ) + rail_button2_bending_moment = Function( + np.column_stack([t, m2_values]), + inputs="Time (s)", + outputs="Bending Moment (N·m)", + interpolation="linear", + ) + + # Maximum bending moments (absolute value for stress calculations) + max_rail_button1_bending_moment = float(np.max(np.abs(m1_values))) + max_rail_button2_bending_moment = float(np.max(np.abs(m2_values))) + + return ( + rail_button1_bending_moment, + max_rail_button1_bending_moment, + rail_button2_bending_moment, + max_rail_button2_bending_moment, + ) + + @property + def rail_button1_bending_moment(self): + """Upper rail button bending moment as a Function of time.""" + return self.calculate_rail_button_bending_moments[0] + + @property + def max_rail_button1_bending_moment(self): + """Maximum upper rail button bending moment, in N·m.""" + return self.calculate_rail_button_bending_moments[1] + + @property + def rail_button2_bending_moment(self): + """Lower rail button bending moment as a Function of time.""" + return self.calculate_rail_button_bending_moments[2] + + @property + def max_rail_button2_bending_moment(self): + """Maximum lower rail button bending moment, in N·m.""" + return self.calculate_rail_button_bending_moments[3] From a3b37d746bbff97b594d9923bab01b7a0443b713 Mon Sep 17 00:00:00 2001 From: GuilhermeAsura Date: Sat, 6 Dec 2025 23:36:04 -0300 Subject: [PATCH 08/12] style: format verification scripts and pin pytz>=2020.1 for CI stability --- tests/animation_verification/main.py | 1 + tests/animation_verification/rocket_stl.py | 1 + 2 files changed, 2 insertions(+) diff --git a/tests/animation_verification/main.py b/tests/animation_verification/main.py index 15c3d42d2..c224f9b25 100644 --- a/tests/animation_verification/main.py +++ b/tests/animation_verification/main.py @@ -7,6 +7,7 @@ from rocketpy import Environment, Flight + def run_simulation_and_test_animation(): print("🚀 Setting up simulation (Calisto Example)...") diff --git a/tests/animation_verification/rocket_stl.py b/tests/animation_verification/rocket_stl.py index 11654cc42..6aecc6971 100644 --- a/tests/animation_verification/rocket_stl.py +++ b/tests/animation_verification/rocket_stl.py @@ -5,6 +5,7 @@ from stl import mesh # Requires numpy-stl package: pip install numpy-stl + def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): """ Creates a detailed rocket-shaped STL file with proper scale. From 822891b735b31bb6a0a09f36f2de87c4fa6ef6b1 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sat, 14 Mar 2026 14:38:25 -0300 Subject: [PATCH 09/12] Refactors the PR solving all reviewer's comments - Implemented 3D flight trajectory and attitude animations in the Flight plots layer. - Added methods `animate_trajectory` and `animate_rotate` to visualize rocket flight. - Included validation for animation inputs and error handling for missing STL files. - Updated documentation to reflect new animation features and installation requirements. - Added optional dependency for `vedo` in `requirements-optional.txt`. - Created a default STL model for the rocket. - Removed outdated animation verification tests and replaced them with unit tests for new animation methods. --- CHANGELOG.md | 37 +- docs/user/flight.rst | 49 +- pyproject.toml | 5 +- requirements-optional.txt | 3 +- rocketpy/plots/assets/default_rocket.stl | 16 + rocketpy/plots/flight_plots.py | 254 ++++++ rocketpy/simulation/flight.py | 807 +++++++------------ tests/animation_verification/__init__.py | 0 tests/animation_verification/main.py | 76 -- tests/animation_verification/rocket_setup.py | 95 --- tests/animation_verification/rocket_stl.py | 234 ------ tests/unit/test_plots.py | 200 +++++ 12 files changed, 811 insertions(+), 965 deletions(-) create mode 100644 rocketpy/plots/assets/default_rocket.stl delete mode 100644 tests/animation_verification/__init__.py delete mode 100644 tests/animation_verification/main.py delete mode 100644 tests/animation_verification/rocket_setup.py delete mode 100644 tests/animation_verification/rocket_stl.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 2183005f8..ae43d86d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,37 +32,7 @@ Attention: The newest changes should be on top --> ### Added -- ENH: Monte Carlo Formatting Options [#947](https://github.com/RocketPy-Team/RocketPy/pull/947) -- ENH: ENH: Auto-Detection of Pressure Conversion Factor [#966](https://github.com/RocketPy-Team/RocketPy/pull/966) -- ENH: Auto-Detection of Pressure Conversion Factor [#966](https://github.com/RocketPy-Team/RocketPy/pull/966) -- ENH: MNT: introduce pressure unit conversion when using forecast/reanalysis/ensemble data [#955](https://github.com/RocketPy-Team/RocketPy/pull/955) -- ENH: Auto Populate Changelog [#919](https://github.com/RocketPy-Team/RocketPy/pull/919) -- ENH: Adaptive Monte Carlo via Convergence Criteria [#922](https://github.com/RocketPy-Team/RocketPy/pull/922) -- TST: Add acceptance tests for 3DOF flight simulation based on Bella Lui rocket [#914](https://github.com/RocketPy-Team/RocketPy/pull/914) - -### Changed - -- - -### Deprecated - -- MNT: Rename `radius` to `radius_function` in `CylindricalTank` and `SphericalTank`; old `radius=` keyword argument now raises `DeprecationWarning` [#957](https://github.com/RocketPy-Team/RocketPy/pull/957) - -### Fixed - -- BUG: fix NaN in ND linear interpolation outside convex hull [#926](https://github.com/RocketPy-Team/RocketPy/issues/926) -- BUG: Add wraparound logic for wind direction in environment plots [#939](https://github.com/RocketPy-Team/RocketPy/pull/939) - -## [v1.12.1] - 2026-04-03 - -### Fixed - -- HOTFIX: Migrate Forecasts to UCAR THREDDS [#943](https://github.com/RocketPy-Team/RocketPy/pull/943) - -## [v1.12.0] - 2026-03-08 - -### Added - +- ENH: Add 3D flight trajectory and attitude animations in Flight plots layer [#909](https://github.com/RocketPy-Team/RocketPy/pull/909) - ENH: Air brakes controller functions now support 8-parameter signature [#854](https://github.com/RocketPy-Team/RocketPy/pull/854) - TST: Add acceptance tests for 3DOF flight simulation based on Bella Lui rocket [#914] (https://github.com/RocketPy-Team/RocketPy/pull/914_ - ENH: Add background map auto download functionality to Monte Carlo plots [#896](https://github.com/RocketPy-Team/RocketPy/pull/896) @@ -89,8 +59,6 @@ Attention: The newest changes should be on top --> ### Fixed -- BUG: Restore `Rocket.power_off_drag` and `Rocket.power_on_drag` as `Function` objects while preserving raw inputs in `power_off_drag_input` and `power_on_drag_input` [#941](https://github.com/RocketPy-Team/RocketPy/pull/941) -- BUG: Add explicit timeouts to ThrustCurve API requests [#935](https://github.com/RocketPy-Team/RocketPy/pull/935) - BUG: Fix hard-coded radius value for parachute added mass calculation [#889](https://github.com/RocketPy-Team/RocketPy/pull/889) - DOC: Fix documentation build [#908](https://github.com/RocketPy-Team/RocketPy/pull/908) - BUG: energy_data plot not working for 3 dof sims [[#906](https://github.com/RocketPy-Team/RocketPy/issues/906)] @@ -98,8 +66,6 @@ Attention: The newest changes should be on top --> - BUG: Fix parallel Monte Carlo simulation showing incorrect iteration count [#806](https://github.com/RocketPy-Team/RocketPy/pull/806) - BUG: Fix missing titles in roll parameter plots for fin sets [#934](https://github.com/RocketPy-Team/RocketPy/pull/934) - BUG: Duplicate _controllers in Flight.TimeNodes.merge() [#931](https://github.com/RocketPy-Team/RocketPy/pull/931) -- BUG: Fix incorrect Jacobian in `only_radial_burn` branch of `SolidMotor.evaluate_geometry` [#935](https://github.com/RocketPy-Team/RocketPy/pull/935) -- BUG: Add explicit timeouts to ThrustCurve API requests [#935](https://github.com/RocketPy-Team/RocketPy/pull/935) ## [v1.11.0] - 2025-11-01 @@ -349,7 +315,6 @@ You can install this version by running `pip install rocketpy==1.4.0` - ENH: CP and Thrust Eccentricity Effects Generate Roll Moment [#617](https://github.com/RocketPy-Team/RocketPy/pull/617) - ENH: Add Prandtl-Gauss transformation to NoseCone and Tail [#609](https://github.com/RocketPy-Team/RocketPy/pull/609) - ENH: Implement power series nose cones [#603](https://github.com/RocketPy-Team/RocketPy/pull/603) -- DOC: Valkyrie-Bisky Team Flight Example [#xxx] (xxx) ### Changed diff --git a/docs/user/flight.rst b/docs/user/flight.rst index 31e7ab588..b005477db 100644 --- a/docs/user/flight.rst +++ b/docs/user/flight.rst @@ -274,7 +274,7 @@ During the rail launch phase, RocketPy calculates reaction forces and internal b **Rail Button Forces (N):** - ``rail_button1_normal_force`` : Normal reaction force at upper rail button -- ``rail_button1_shear_force`` : Shear (tangential) reaction force at upper rail button +- ``rail_button1_shear_force`` : Shear (tangential) reaction force at upper rail button - ``rail_button2_normal_force`` : Normal reaction force at lower rail button - ``rail_button2_shear_force`` : Shear (tangential) reaction force at lower rail button @@ -282,7 +282,7 @@ During the rail launch phase, RocketPy calculates reaction forces and internal b - ``rail_button1_bending_moment`` : Time-dependent bending moment at upper rail button attachment - ``max_rail_button1_bending_moment`` : Maximum absolute bending moment at upper rail button -- ``rail_button2_bending_moment`` : Time-dependent bending moment at lower rail button attachment +- ``rail_button2_bending_moment`` : Time-dependent bending moment at lower rail button attachment - ``max_rail_button2_bending_moment`` : Maximum absolute bending moment at lower rail button **Calculation Method:** @@ -454,6 +454,51 @@ Flight Data Plots # Flight path and orientation flight.plots.flight_path_angle_data() +3D Flight Animation +~~~~~~~~~~~~~~~~~~~ + +RocketPy can animate the simulated flight trajectory and attitude through the +Flight plots layer. + +.. note:: + + Install optional animation dependencies first: + + .. code-block:: bash + + pip install rocketpy[animation] + +.. code-block:: python + + # Fast start using RocketPy's built-in default STL model + flight.plots.animate_trajectory( + start=0.0, + stop=min(flight.t_final, 20.0), + time_step=0.05, + ) + + # Or provide your own STL model file + flight.plots.animate_trajectory( + file_name="rocket.stl", + start=0.0, + stop=flight.t_final, + time_step=0.05, + azimuth=45, + elevation=20, + ) + + # Keep rocket centered and animate only attitude changes + flight.plots.animate_rotate( + file_name="rocket.stl", + start=0.0, + stop=min(flight.t_final, 20.0), + time_step=0.05, + ) + +Both methods validate the selected time interval and STL path before rendering. +If ``vedo`` is not installed, RocketPy raises an informative ``ImportError`` +with installation instructions. + Forces and Moments ~~~~~~~~~~~~~~~~~~ diff --git a/pyproject.toml b/pyproject.toml index 16cc479e8..09ce86273 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,9 @@ build-backend = "setuptools.build_meta" [tool.setuptools] packages = { find = { where = ["."], include = ["rocketpy*"] } } +[tool.setuptools.package-data] +"rocketpy.plots" = ["assets/*.stl"] + [tool.setuptools.dynamic] dependencies = { file = ["requirements.txt"] } @@ -61,7 +64,7 @@ env-analysis = [ ] monte-carlo = [ - "imageio", + "imageio", "multiprocess>=0.70", "statsmodels", "prettytable", diff --git a/requirements-optional.txt b/requirements-optional.txt index 58ed1030b..2961946ca 100644 --- a/requirements-optional.txt +++ b/requirements-optional.txt @@ -6,4 +6,5 @@ timezonefinder imageio multiprocess>=0.70 statsmodels -prettytable \ No newline at end of file +prettytable +vedo>=2024.5.1 \ No newline at end of file diff --git a/rocketpy/plots/assets/default_rocket.stl b/rocketpy/plots/assets/default_rocket.stl new file mode 100644 index 000000000..e3889fe36 --- /dev/null +++ b/rocketpy/plots/assets/default_rocket.stl @@ -0,0 +1,16 @@ +solid default_rocket + facet normal 0 0 1 + outer loop + vertex 0 0 0 + vertex 1 0 0 + vertex 0 1 0 + endloop + endfacet + facet normal 0 0 -1 + outer loop + vertex 0 0 0 + vertex 0 1 0 + vertex 1 0 0 + endloop + endfacet +endsolid default_rocket diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py index 7eb41a8b2..86118037f 100644 --- a/rocketpy/plots/flight_plots.py +++ b/rocketpy/plots/flight_plots.py @@ -1,8 +1,12 @@ +import os +import time from functools import cached_property +from importlib import resources import matplotlib.pyplot as plt import numpy as np +from ..tools import import_optional_dependency from .plot_helpers import show_or_save_plot @@ -133,6 +137,256 @@ def trajectory_3d(self, *, filename=None): # pylint: disable=too-many-statement ax1.set_box_aspect(None, zoom=0.95) # 95% for label adjustment show_or_save_plot(filename) + def _resolve_animation_model_path(self, file_name): + """Resolve model path, defaulting to the built-in STL when omitted.""" + if file_name is not None: + return file_name + + return str( + resources.files("rocketpy.plots").joinpath("assets/default_rocket.stl") + ) + + def _validate_animation_inputs(self, file_name, start, stop, time_step): + """Validate shared input parameters for 3D animation methods.""" + if time_step <= 0: + raise ValueError( + f"Invalid time_step: {time_step}. It must be greater than 0." + ) + + if stop is None: + stop = self.flight.t_final + + if ( + start < 0 + or stop < 0 + or start > self.flight.t_final + or stop > self.flight.t_final + or start >= stop + ): + raise ValueError( + f"Invalid animation time range: start={start}, stop={stop}. " + f"Both must be within [0, {self.flight.t_final}] and start < stop." + ) + + if not os.path.isfile(file_name): + raise FileNotFoundError( + f"Could not find the 3D model file: '{file_name}'. " + "Provide a valid .stl file path." + ) + + return stop + + @staticmethod + def _rotation_from_quaternion(q0, q1, q2, q3): + """Convert unit quaternion to axis-angle representation in degrees.""" + norm = np.sqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3) + if norm == 0: + return 0.0, (1.0, 0.0, 0.0) + + q0 = q0 / norm + q1 = q1 / norm + q2 = q2 / norm + q3 = q3 / norm + + # q and -q represent the same orientation. Keep q0 non-negative to + # reduce discontinuities in axis-angle interpolation across frames. + if q0 < 0: + q0 = -q0 + q1 = -q1 + q2 = -q2 + q3 = -q3 + + q0 = np.clip(q0, -1.0, 1.0) + angle = 2 * np.arccos(q0) + sin_half = np.sqrt(max(1 - q0 * q0, 0.0)) + + if sin_half < 1e-12: + return 0.0, (1.0, 0.0, 0.0) + + axis = (q1 / sin_half, q2 / sin_half, q3 / sin_half) + return np.degrees(angle), axis + + def _create_animation_box(self, start, scale=1.0): + """Create a world box with minimum visible dimensions.""" + min_box_dim = 10.0 + x_values = self.flight.x[:, 1] + y_values = self.flight.y[:, 1] + z_values = self.flight.z[:, 1] - self.flight.env.elevation + + center_x = 0.5 * (np.max(x_values) + np.min(x_values)) + center_y = 0.5 * (np.max(y_values) + np.min(y_values)) + center_z = max(self.flight.z(start) - self.flight.env.elevation, 0.0) + + length = max(np.ptp(x_values) * scale, min_box_dim) + width = max(np.ptp(y_values) * scale, min_box_dim) + height = max(np.ptp(z_values) * scale, min_box_dim) + + # Keep z center inside visible space while preserving minimum box size. + center_z = max(center_z, 0.5 * min_box_dim) + + vedo = import_optional_dependency("vedo") + Box = vedo.Box + + return Box( + pos=[center_x, center_y, center_z], + length=length, + width=width, + height=height, + ).wireframe() + + def animate_trajectory( + self, file_name=None, start=0, stop=None, time_step=0.1, **kwargs + ): + """Animate 6-DOF trajectory and attitude using vedo. + + Parameters + ---------- + file_name : str | None, optional + Path to a 3D model file representing the rocket, usually ``.stl``. + If None, RocketPy uses a built-in default STL model. + Default is None. + start : int, float, optional + Animation start time in seconds. Default is 0. + stop : int, float | None, optional + Animation end time in seconds. If None, uses ``flight.t_final``. + Default is None. + time_step : float, optional + Animation frame step in seconds. Must be greater than 0. + Default is 0.1. + **kwargs : dict, optional + Additional keyword arguments passed to ``vedo.Plotter.show``. + """ + + vedo = import_optional_dependency("vedo") + + Line = vedo.Line + Mesh = vedo.Mesh + Plotter = vedo.Plotter + settings = vedo.settings + + file_name = self._resolve_animation_model_path(file_name) + stop = self._validate_animation_inputs(file_name, start, stop, time_step) + + try: + settings.allow_interaction = True + except AttributeError: + pass + + world = self._create_animation_box(start, scale=1.2) + base_rocket = Mesh(file_name).c("green") + time_steps = np.arange(start, stop, time_step) + trajectory_points = [] + + plt = Plotter(axes=1, interactive=False) + plt.show(world, "Rocket Trajectory Animation", viewup="z", **kwargs) + + for t in time_steps: + rocket = base_rocket.clone() + x_position = self.flight.x(t) + y_position = self.flight.y(t) + z_position = self.flight.z(t) - self.flight.env.elevation + + angle_deg, axis = self._rotation_from_quaternion( + self.flight.e0(t), + self.flight.e1(t), + self.flight.e2(t), + self.flight.e3(t), + ) + + rocket.pos(x_position, y_position, z_position) + if angle_deg != 0.0: + rocket.rotate(angle_deg, axis=axis) + + trajectory_points.append([x_position, y_position, z_position]) + actors = [world, rocket] + if len(trajectory_points) > 1: + actors.append(Line(trajectory_points, c="k", alpha=0.5)) + + plt.show(*actors, resetcam=False) + + start_pause = time.time() + while time.time() - start_pause < time_step: + plt.render() + + if getattr(plt, "escaped", False): + break + + plt.interactive().close() + return None + + def animate_rotate( + self, file_name=None, start=0, stop=None, time_step=0.1, **kwargs + ): + """Animate rocket attitude (rotation-only view) using vedo. + + Parameters + ---------- + file_name : str | None, optional + Path to a 3D model file representing the rocket, usually ``.stl``. + If None, RocketPy uses a built-in default STL model. + Default is None. + start : int, float, optional + Animation start time in seconds. Default is 0. + stop : int, float | None, optional + Animation end time in seconds. If None, uses ``flight.t_final``. + Default is None. + time_step : float, optional + Animation frame step in seconds. Must be greater than 0. + Default is 0.1. + **kwargs : dict, optional + Additional keyword arguments passed to ``vedo.Plotter.show``. + """ + + vedo = import_optional_dependency("vedo") + + Mesh = vedo.Mesh + Plotter = vedo.Plotter + settings = vedo.settings + + file_name = self._resolve_animation_model_path(file_name) + stop = self._validate_animation_inputs(file_name, start, stop, time_step) + + try: + settings.allow_interaction = True + except AttributeError: + pass + + world = self._create_animation_box(start, scale=0.3) + base_rocket = Mesh(file_name).c("green") + time_steps = np.arange(start, stop, time_step) + + x_start = self.flight.x(start) + y_start = self.flight.y(start) + z_start = self.flight.z(start) - self.flight.env.elevation + + plt = Plotter(axes=1, interactive=False) + plt.show(world, "Rocket Rotation Animation", viewup="z", **kwargs) + + for t in time_steps: + rocket = base_rocket.clone() + angle_deg, axis = self._rotation_from_quaternion( + self.flight.e0(t), + self.flight.e1(t), + self.flight.e2(t), + self.flight.e3(t), + ) + + rocket.pos(x_start, y_start, z_start) + if angle_deg != 0.0: + rocket.rotate(angle_deg, axis=axis) + + plt.show(world, rocket, resetcam=False) + + start_pause = time.time() + while time.time() - start_pause < time_step: + plt.render() + + if getattr(plt, "escaped", False): + break + + plt.interactive().close() + return None + def linear_kinematics_data(self, *, filename=None): # pylint: disable=too-many-statements """Prints out all Kinematics graphs available about the Flight diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 82f8c99f4..eb82b1998 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1,6 +1,5 @@ # pylint: disable=too-many-lines import math -import time import warnings from copy import deepcopy from functools import cached_property @@ -717,6 +716,7 @@ def __simulate(self, verbose): self.y_sol, self.solution, self.sensors, + self.env, ) for parachute in node.parachutes: @@ -760,11 +760,12 @@ def __simulate(self, verbose): 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, + lambda self, added_mass_coefficient=parachute.added_mass_coefficient: ( + setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, + ) ), ] self.flight_phases.add_phase( @@ -874,6 +875,7 @@ def __process_sensors_and_controllers_at_current_node(self, node, phase): self.y_sol, self.solution, self.sensors, + self.env, ) def __measure_sensors(self, component_sensors, u_dot, t=None, y_sol=None): @@ -982,11 +984,12 @@ def __check_and_handle_parachute_triggers( 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, + lambda self, added_mass_coefficient=parachute.added_mass_coefficient: ( + setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, + ) ), ] self.flight_phases.add_phase( @@ -1294,6 +1297,7 @@ def __process_overshootable_nodes(self, phase, phase_index, node_index): overshootable_node.y_sol, self.solution, self.sensors, + self.env, ) # Process sensors at overshootable node @@ -1388,11 +1392,12 @@ def __check_overshootable_parachute_triggers( 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, + lambda self, added_mass_coefficient=parachute.added_mass_coefficient: ( + setattr( + self, + "parachute_added_mass_coefficient", + added_mass_coefficient, + ) ), ] self.flight_phases.add_phase( @@ -1698,89 +1703,28 @@ def lateral_surface_wind(self): return -wind_u * np.cos(heading_rad) + wind_v * np.sin(heading_rad) - def __get_drag_coefficient(self, drag_function, mach, z, freestream_velocity_body): - """Calculate drag coefficient, handling both 1D and multi-dimensional functions. - - Parameters - ---------- - drag_function : Function - The drag coefficient function (power_on_drag or power_off_drag) - mach : float - Mach number - z : float - Altitude in meters - freestream_velocity_body : Vector or array-like - Freestream velocity in body frame [stream_vx_b, stream_vy_b, stream_vz_b] - - Returns - ------- - float - Drag coefficient value - """ - # Early return for 1D drag functions (only mach number) - if not isinstance(drag_function, Function) or not getattr( - drag_function, "is_multidimensional", False - ): - return drag_function.get_value_opt(mach) - - # Multi-dimensional drag function - calculate additional parameters - - # Calculate Reynolds number: Re = rho * V * L / mu - # where L is characteristic length (rocket diameter) - rho = self.env.density.get_value_opt(z) - mu = self.env.dynamic_viscosity.get_value_opt(z) - freestream_speed = np.linalg.norm(freestream_velocity_body) - characteristic_length = 2 * self.rocket.radius # Diameter - # Defensive: avoid division by zero or non-finite viscosity values. - # Use a small epsilon fallback if `mu` is zero, negative, NaN or infinite. - try: - mu_val = float(mu) - except (TypeError, ValueError, OverflowError): - # Only catch errors related to invalid numeric conversion. - # Avoid catching broad Exception to satisfy linters and - # allow other unexpected errors to surface. - mu_val = 0.0 - if not np.isfinite(mu_val) or mu_val <= 0.0: - mu_safe = 1e-10 - else: - mu_safe = mu_val - - reynolds = rho * freestream_speed * characteristic_length / mu_safe - - # Calculate angle of attack - # Angle between freestream velocity and rocket axis (z-axis in body frame) - # The z component of freestream velocity in body frame - if hasattr(freestream_velocity_body, "z"): - stream_vz_b = -freestream_velocity_body.z - else: - stream_vz_b = -freestream_velocity_body[2] - - # Normalize and calculate angle - if freestream_speed > 1e-6: - cos_alpha = stream_vz_b / freestream_speed - # Clamp to [-1, 1] to avoid numerical issues - cos_alpha = np.clip(cos_alpha, -1.0, 1.0) - alpha_rad = np.arccos(cos_alpha) - alpha_deg = np.rad2deg(alpha_rad) - else: - alpha_deg = 0.0 - - # Determine which parameters to pass based on input names - input_names = [name.lower() for name in drag_function.__inputs__] - args = [] - - for name in input_names: - if "mach" in name or name == "m": - args.append(mach) - elif "reynolds" in name or name == "re": - args.append(reynolds) - elif "alpha" in name or name == "a" or "attack" in name: - args.append(alpha_deg) - else: - # Unknown parameter, default to mach - args.append(mach) - - return drag_function.get_value_opt(*args) + def __compute_drag_7d_inputs( + self, + stream_velocity_body, + stream_speed, + stream_mach, + density, + dynamic_viscosity, + ): + """Build drag-model inputs in the 7D order used by Rocket drag functions.""" + aerodynamic_stream_velocity = -stream_velocity_body + alpha = np.arctan2( + aerodynamic_stream_velocity[1], aerodynamic_stream_velocity[2] + ) + beta = np.arctan2( + aerodynamic_stream_velocity[0], aerodynamic_stream_velocity[2] + ) + reynolds = ( + density * stream_speed * (2 * self.rocket.radius) / dynamic_viscosity + if dynamic_viscosity > 0 + else 0 + ) + return alpha, beta, stream_mach, reynolds def udot_rail1(self, t, u, post_processing=False): """Calculates derivative of u state vector with respect to time @@ -1812,38 +1756,28 @@ def udot_rail1(self, t, u, post_processing=False): total_mass_at_t = self.rocket.total_mass.get_value_opt(t) # Get freestream speed - free_stream_speed = ( - (self.env.wind_velocity_x.get_value_opt(z) - vx) ** 2 - + (self.env.wind_velocity_y.get_value_opt(z) - vy) ** 2 - + (vz) ** 2 - ) ** 0.5 + free_stream_velocity = Vector( + [ + self.env.wind_velocity_x.get_value_opt(z) - vx, + self.env.wind_velocity_y.get_value_opt(z) - vy, + -vz, + ] + ) + free_stream_speed = abs(free_stream_velocity) free_stream_mach = free_stream_speed / self.env.speed_of_sound.get_value_opt(z) - - # For rail motion, rocket is constrained - velocity mostly along z-axis in body frame - # Calculate velocity in body frame (simplified for rail) - a11 = 1 - 2 * (e2**2 + e3**2) - a12 = 2 * (e1 * e2 - e0 * e3) - a13 = 2 * (e1 * e3 + e0 * e2) - a21 = 2 * (e1 * e2 + e0 * e3) - a22 = 1 - 2 * (e1**2 + e3**2) - a23 = 2 * (e2 * e3 - e0 * e1) - a31 = 2 * (e1 * e3 - e0 * e2) - a32 = 2 * (e2 * e3 + e0 * e1) - a33 = 1 - 2 * (e1**2 + e2**2) - - # Freestream velocity in body frame - wind_vx = self.env.wind_velocity_x.get_value_opt(z) - wind_vy = self.env.wind_velocity_y.get_value_opt(z) - stream_vx_b = a11 * (wind_vx - vx) + a21 * (wind_vy - vy) + a31 * (-vz) - stream_vy_b = a12 * (wind_vx - vx) + a22 * (wind_vy - vy) + a32 * (-vz) - stream_vz_b = a13 * (wind_vx - vx) + a23 * (wind_vy - vy) + a33 * (-vz) - - drag_coeff = self.__get_drag_coefficient( - self.rocket.power_on_drag, + rho = self.env.density.get_value_opt(z) + stream_velocity_body = ( + Matrix.transformation([e0, e1, e2, e3]).transpose @ free_stream_velocity + ) + dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) + alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( + stream_velocity_body, + free_stream_speed, free_stream_mach, - z, - [stream_vx_b, stream_vy_b, stream_vz_b], + rho, + dynamic_viscosity, ) + drag_coeff = self.rocket.power_on_drag_7d(alpha, beta, mach, reynolds, 0, 0, 0) # Calculate Forces pressure = self.env.pressure.get_value_opt(z) @@ -1852,7 +1786,6 @@ def udot_rail1(self, t, u, post_processing=False): + self.rocket.motor.pressure_thrust(pressure), 0, ) - rho = self.env.density.get_value_opt(z) R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff) # Calculate Linear acceleration @@ -2008,44 +1941,42 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) speed_of_sound = self.env.speed_of_sound.get_value_opt(z) - free_stream_speed = ( - (wind_velocity_x - vx) ** 2 + (wind_velocity_y - vy) ** 2 + (vz) ** 2 - ) ** 0.5 + free_stream_velocity = Vector([wind_velocity_x - vx, wind_velocity_y - vy, -vz]) + free_stream_speed = abs(free_stream_velocity) free_stream_mach = free_stream_speed / speed_of_sound - - # Get rocket velocity in body frame (needed for drag calculation) - vx_b = a11 * vx + a21 * vy + a31 * vz - vy_b = a12 * vx + a22 * vy + a32 * vz - vz_b = a13 * vx + a23 * vy + a33 * vz - - # Calculate freestream velocity in body frame - stream_vx_b = ( - a11 * (wind_velocity_x - vx) + a21 * (wind_velocity_y - vy) + a31 * (-vz) - ) - stream_vy_b = ( - a12 * (wind_velocity_x - vx) + a22 * (wind_velocity_y - vy) + a32 * (-vz) - ) - stream_vz_b = ( - a13 * (wind_velocity_x - vx) + a23 * (wind_velocity_y - vy) + a33 * (-vz) - ) + stream_velocity_body = Kt @ free_stream_velocity # Determine aerodynamics forces # Determine Drag Force + rho = self.env.density.get_value_opt(z) + dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) + alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( + stream_velocity_body, + free_stream_speed, + free_stream_mach, + rho, + dynamic_viscosity, + ) if t < self.rocket.motor.burn_out_time: - drag_coeff = self.__get_drag_coefficient( - self.rocket.power_on_drag, - free_stream_mach, - z, - [stream_vx_b, stream_vy_b, stream_vz_b], + drag_coeff = self.rocket.power_on_drag_7d( + alpha, + beta, + mach, + reynolds, + omega1, + omega2, + omega3, ) else: - drag_coeff = self.__get_drag_coefficient( - self.rocket.power_off_drag, - free_stream_mach, - z, - [stream_vx_b, stream_vy_b, stream_vz_b], + drag_coeff = self.rocket.power_off_drag_7d( + alpha, + beta, + mach, + reynolds, + omega1, + omega2, + omega3, ) - rho = self.env.density.get_value_opt(z) R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff for air_brakes in self.rocket.air_brakes: if air_brakes.deployment_level > 0: @@ -2066,6 +1997,10 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Off center moment M1 += self.rocket.cp_eccentricity_y * R3 M2 -= self.rocket.cp_eccentricity_x * R3 + # Get rocket velocity in body frame + vx_b = a11 * vx + a21 * vy + a31 * vz + vy_b = a12 * vx + a22 * vy + a32 * vz + vz_b = a13 * vx + a23 * vy + a33 * vz # Calculate lift and moment for each component of the rocket velocity_in_body_frame = Vector([vx_b, vy_b, vz_b]) w = Vector([omega1, omega2, omega3]) @@ -2086,11 +2021,13 @@ def u_dot(self, t, u, post_processing=False): # pylint: disable=too-many-locals # Reynolds at component altitude # TODO: Reynolds is only used in generic surfaces. This calculation # should be moved to the surface class for efficiency + comp_density = self.env.density.get_value_opt(comp_z) + comp_dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(comp_z) comp_reynolds = ( - self.env.density.get_value_opt(comp_z) + comp_density * comp_stream_speed * aero_surface.reference_length - / self.env.dynamic_viscosity.get_value_opt(comp_z) + / comp_dynamic_viscosity ) # Forces and moments X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( @@ -2284,12 +2221,25 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): 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 + stream_velocity_body = Kt @ free_stream_velocity + dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) + alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( + stream_velocity_body, + free_stream_speed, + mach, + rho, + dynamic_viscosity, + ) # Drag computation if t < self.rocket.motor.burn_out_time: - cd = self.rocket.power_on_drag.get_value_opt(mach) + cd = self.rocket.power_on_drag_7d( + alpha, beta, mach, reynolds, omega1, omega2, omega3 + ) else: - cd = self.rocket.power_off_drag.get_value_opt(mach) + cd = self.rocket.power_off_drag_7d( + alpha, beta, mach, reynolds, omega1, omega2, omega3 + ) R1, R2 = 0, 0 R3 = -0.5 * rho * free_stream_speed**2 * self.rocket.area * cd @@ -2323,11 +2273,13 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): rel_speed = abs(rel_velocity) rel_mach = rel_speed / speed_of_sound + comp_density = self.env.density.get_value_opt(comp_z) + comp_dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(comp_z) reynolds = ( - self.env.density.get_value_opt(comp_z) + comp_density * rel_speed * surface.reference_length - / self.env.dynamic_viscosity.get_value_opt(comp_z) + / comp_dynamic_viscosity ) fx, fy, fz, *_ = surface.compute_forces_and_moments( @@ -2512,15 +2464,19 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too 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))) + 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) free_stream_mach = free_stream_speed / speed_of_sound - - # Get rocket velocity in body frame (needed for drag calculation) - velocity_in_body_frame = Kt @ v - # Calculate freestream velocity in body frame - freestream_velocity = wind_velocity - v - freestream_velocity_body = Kt @ freestream_velocity + stream_velocity_body = Kt @ free_stream_velocity + dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(z) + alpha, beta, mach, reynolds = self.__compute_drag_7d_inputs( + stream_velocity_body, + free_stream_speed, + free_stream_mach, + rho, + dynamic_viscosity, + ) if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time: pressure = self.env.pressure.get_value_opt(z) @@ -2529,19 +2485,25 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too + self.rocket.motor.pressure_thrust(pressure), 0, ) - drag_coeff = self.__get_drag_coefficient( - self.rocket.power_on_drag, - free_stream_mach, - z, - freestream_velocity_body, + drag_coeff = self.rocket.power_on_drag_7d( + alpha, + beta, + mach, + reynolds, + omega1, + omega2, + omega3, ) else: net_thrust = 0 - drag_coeff = self.__get_drag_coefficient( - self.rocket.power_off_drag, - free_stream_mach, - z, - freestream_velocity_body, + drag_coeff = self.rocket.power_off_drag_7d( + alpha, + beta, + mach, + reynolds, + omega1, + omega2, + omega3, ) R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff for air_brakes in self.rocket.air_brakes: @@ -2560,6 +2522,8 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too R3 = air_brakes_force # Substitutes rocket drag coefficient else: R3 += air_brakes_force + # Get rocket velocity in body frame + velocity_in_body_frame = Kt @ v # Calculate lift and moment for each component of the rocket for aero_surface, _ in self.rocket.aerodynamic_surfaces: # Component cp relative to CDM in body frame @@ -2578,11 +2542,13 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too # Reynolds at component altitude # TODO: Reynolds is only used in generic surfaces. This calculation # should be moved to the surface class for efficiency + comp_density = self.env.density.get_value_opt(comp_z) + comp_dynamic_viscosity = self.env.dynamic_viscosity.get_value_opt(comp_z) comp_reynolds = ( - self.env.density.get_value_opt(comp_z) + comp_density * comp_stream_speed * aero_surface.reference_length - / self.env.dynamic_viscosity.get_value_opt(comp_z) + / comp_dynamic_viscosity ) # Forces and moments X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( @@ -3160,7 +3126,7 @@ def max_acceleration_power_off_time(self): max_acceleration_time_index = np.argmax( self.acceleration[burn_out_time_index:, 1] ) - return self.acceleration[max_acceleration_time_index, 0] + return self.acceleration[burn_out_time_index + max_acceleration_time_index, 0] @cached_property def max_acceleration_power_off(self): @@ -3721,6 +3687,142 @@ def max_rail_button2_shear_force(self): """Maximum lower rail button shear force, in Newtons.""" return np.abs(self.rail_button2_shear_force.y_array).max() + @cached_property + def calculate_rail_button_bending_moments(self): + """Calculate internal bending moments at rail button attachment points. + + Uses beam theory to determine the internal structural moments for + stress analysis of the rail button attachments (fasteners and airframe). + + The bending moment at each button attachment consists of: + + 1. Normal force moment: $M = N \\times d$, where $N$ is the normal + reaction force and $d$ is the distance from button to center of + dry mass. + 2. Shear force cantilever moment: $M = S \\times h$, where $S$ is the + shear (tangential) force and $h$ is the button standoff height. + + Returns + ------- + tuple + rail_button1_bending_moment : Function + Bending moment at upper rail button as a function of time (N·m). + max_rail_button1_bending_moment : float + Maximum upper rail button bending moment (N·m). + rail_button2_bending_moment : Function + Bending moment at lower rail button as a function of time (N·m). + max_rail_button2_bending_moment : float + Maximum lower rail button bending moment (N·m). + + Notes + ----- + - Calculated only during the rail phase of flight + - Maximum values use absolute values for worst-case stress analysis + - The bending moments represent internal stresses in the rocket + airframe at the rail button attachment points + + **Assumptions:** + + - Rail buttons act as simple supports: provide reaction forces (normal + and shear) but no moment reaction at the rail contact point + - The rocket acts as a beam supported at two points (rail buttons) + - Bending moments arise from the lever arm effect of reaction forces + and the cantilever moment from button standoff height + """ + # Check if rail buttons exist + null_moment = Function(0) + if len(self.rocket.rail_buttons) == 0: + warnings.warn( + "Trying to calculate rail button bending moments without " + "rail buttons defined. Setting moments to zero.", + UserWarning, + ) + return (null_moment, 0.0, null_moment, 0.0) + + # Get rail button geometry + rail_buttons_tuple = self.rocket.rail_buttons[0] + # Rail button standoff height + h_button = rail_buttons_tuple.component.button_height + if h_button is None: + warnings.warn( + "Rail button height not defined. Bending moments cannot be " + "calculated. Setting moments to zero.", + UserWarning, + ) + return (null_moment, 0.0, null_moment, 0.0) + upper_button_position = ( + rail_buttons_tuple.component.buttons_distance + + rail_buttons_tuple.position.z + ) + lower_button_position = rail_buttons_tuple.position.z + + # Get center of dry mass (handle both callable and property) + if callable(self.rocket.center_of_dry_mass_position): + cdm = self.rocket.center_of_dry_mass_position(self.rocket._csys) + else: + cdm = self.rocket.center_of_dry_mass_position + + # Distances from buttons to center of dry mass + d1 = abs(upper_button_position - cdm) + d2 = abs(lower_button_position - cdm) + + # forces + N1 = self.rail_button1_normal_force + N2 = self.rail_button2_normal_force + S1 = self.rail_button1_shear_force + S2 = self.rail_button2_shear_force + t = N1.source[:, 0] + + # Calculate bending moments at attachment points + # Primary contribution from shear force acting at button height + # Secondary contribution from normal force creating moment about attachment + m1_values = N2.source[:, 1] * d2 + S1.source[:, 1] * h_button + m2_values = N1.source[:, 1] * d1 + S2.source[:, 1] * h_button + + rail_button1_bending_moment = Function( + np.column_stack([t, m1_values]), + inputs="Time (s)", + outputs="Bending Moment (N·m)", + interpolation="linear", + ) + rail_button2_bending_moment = Function( + np.column_stack([t, m2_values]), + inputs="Time (s)", + outputs="Bending Moment (N·m)", + interpolation="linear", + ) + + # Maximum bending moments (absolute value for stress calculations) + max_rail_button1_bending_moment = float(np.max(np.abs(m1_values))) + max_rail_button2_bending_moment = float(np.max(np.abs(m2_values))) + + return ( + rail_button1_bending_moment, + max_rail_button1_bending_moment, + rail_button2_bending_moment, + max_rail_button2_bending_moment, + ) + + @property + def rail_button1_bending_moment(self): + """Upper rail button bending moment as a Function of time.""" + return self.calculate_rail_button_bending_moments[0] + + @property + def max_rail_button1_bending_moment(self): + """Maximum upper rail button bending moment, in N·m.""" + return self.calculate_rail_button_bending_moments[1] + + @property + def rail_button2_bending_moment(self): + """Lower rail button bending moment as a Function of time.""" + return self.calculate_rail_button_bending_moments[2] + + @property + def max_rail_button2_bending_moment(self): + """Maximum lower rail button bending moment, in N·m.""" + return self.calculate_rail_button_bending_moments[3] + @funcify_method( "Time (s)", "Horizontal Distance to Launch Point (m)", "spline", "constant" ) @@ -4040,202 +4142,7 @@ def export_kml( color=color, altitude_mode=altitude_mode, ) - def animate_trajectory(self, file_name, start=0, stop=None, time_step=0.1, **kwargs): - """ - 6-DOF Animation of the flight trajectory using Vedo. - - Parameters - ---------- - file_name : str - 3D object file representing your rocket, usually in .stl format. - Example: "rocket.stl" - start : int, float, optional - Time for starting animation, in seconds. Default is 0. - stop : int, float, optional - Time for ending animation, in seconds. If None, uses self.t_final. - Default is None. - time_step : float, optional - Time step for data interpolation in the animation. Default is 0.1. - **kwargs : dict, optional - Additional keyword arguments to be passed to vedo.Plotter.show(). - Common arguments: - - azimuth (float): Rotation in degrees around the vertical axis. - - elevation (float): Rotation in degrees above the horizon. - - roll (float): Rotation in degrees around the view axis. - - zoom (float): Zoom level (default 1). - - Returns - ------- - None - Raises - ------ - ImportError - If the 'vedo' package is not installed. - - Notes - ----- - This feature requires the 'vedo' package. Install it with: - pip install rocketpy[animation] - """ - try: - from vedo import Box, Line, Mesh, Plotter, settings - except ImportError as e: - raise ImportError( - "The animation feature requires the 'vedo' package. " - "Install it with:\n" - " pip install rocketpy[animation]\n" - "Or directly:\n" - " pip install vedo>=2024.5.1" - ) from e - - # Enable interaction if needed - try: - settings.allow_interaction = True - except AttributeError: - pass # Not available in newer versions of vedo - - # Handle stop time - if stop is None: - stop = self.t_final - - # Define the world bounds based on trajectory - max_x = max(self.x[:, 1]) - max_y = max(self.y[:, 1]) - # Use simple logic for bounds - world = Box( - pos=[self.x(start), self.y(start), self.apogee], - length=max_x * 2 if max_x != 0 else 1000, - width=max_y * 2 if max_y != 0 else 1000, - height=self.apogee, - ).wireframe() - - # Load rocket mesh - rocket = Mesh(file_name).c("green") - rocket.pos(self.x(start), self.y(start), 0).add_trail(n=len(self.x[:, 1])) - # Create trail - trail_points = [[self.x(t), self.y(t), self.z(t) - self.env.elevation] - for t in np.arange(start, stop, time_step)] - trail = Line(trail_points, c="k", alpha=0.5) - # Setup Plotter - plt = Plotter(axes=1, interactive=False) - plt.show(world, rocket, __doc__, viewup="z", **kwargs) - - # Animation Loop - for t in np.arange(start, stop, time_step): - # Calculate rotation angle and vector from quaternions - # Note: This simple rotation logic mimics the old branch. - # Ideally, vedo handles orientation via matrix, but we stick - # to the provided logic for now. - - # e0 is the scalar part of the quaternion - angle = np.arccos(2 * self.e0(t)**2 - 1) - k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 - - # Update position and rotation - # Adjusting for ground elevation - rocket.pos(self.x(t), self.y(t), self.z(t) - self.env.elevation) - rocket.rotate_x(self.e1(t) / k) - rocket.rotate_y(self.e2(t) / k) - rocket.rotate_z(self.e3(t) / k) - - # update the scene - plt.show(world, rocket, trail) - - # slow down to make animation visible - start_pause = time.time() - while time.time() - start_pause < time_step: - plt.render() - - if getattr(plt, 'escaped', False): - break - - plt.interactive().close() - return None - - def animate_rotate(self, file_name, start=0, stop=None, time_step=0.1, **kwargs): - """ - Animation of rocket attitude (rotation) during the flight. - - Parameters - ---------- - file_name : str - 3D object file representing your rocket, usually in .stl format. - start : int, float, optional - Time for starting animation, in seconds. Default is 0. - stop : int, float, optional - Time for ending animation, in seconds. If None, uses self.t_final. - Default is None. - time_step : float, optional - Time step for data interpolation. Default is 0.1. - **kwargs : dict, optional - Additional keyword arguments to be passed to vedo.Plotter.show(). - Common arguments: - - azimuth (float): Rotation in degrees around the vertical axis. - - elevation (float): Rotation in degrees above the horizon. - - roll (float): Rotation in degrees around the view axis. - - zoom (float): Zoom level (default 1). - - Returns - ------- - None - - Raises - ------ - ImportError - If the 'vedo' package is not installed. - """ - try: - from vedo import Box, Mesh, Plotter, settings - except ImportError as e: - raise ImportError( - "The animation feature requires the 'vedo' package. " - "Install it with:\n" - " pip install rocketpy[animation]\n" - ) from e - - # Enable interaction if needed - try: - settings.allow_interaction = True - except AttributeError: - pass # Not available in newer versions of vedo - - if stop is None: - stop = self.t_final - - # Smaller box for rotation view - world = Box( - pos=[self.x(start), self.y(start), self.apogee], - length=max(self.x[:, 1]) * 0.2, - width=max(self.y[:, 1]) * 0.2, - height=self.apogee * 0.1, - ).wireframe() - - rocket = Mesh(file_name).c("green") - # Initialize at origin relative to view - rocket.pos(self.x(start), self.y(start), 0).add_trail(n=len(self.x[:, 1])) - - plt = Plotter(axes=1, interactive=False) - plt.show(world, rocket, __doc__, viewup="z", **kwargs) - - for t in np.arange(start, stop, time_step): - angle = np.arccos(2 * self.e0(t)**2 - 1) - k = np.sin(angle / 2) if np.sin(angle / 2) != 0 else 1 - - # Keep position static (relative start) to observe only rotation - rocket.pos(self.x(start), self.y(start), 0) - rocket.rotate_x(self.e1(t) / k) - rocket.rotate_y(self.e2(t) / k) - rocket.rotate_z(self.e3(t) / k) - - plt.show(world, rocket) - - if getattr(plt, 'escaped', False): - break - - plt.interactive().close() - return None - def info(self): """Prints out a summary of the data available about the Flight.""" self.prints.all() @@ -4646,7 +4553,6 @@ def merge(self): tmp_dict[time]._controllers += node._controllers tmp_dict[time].callbacks += node.callbacks tmp_dict[time]._component_sensors += node._component_sensors - tmp_dict[time]._controllers += node._controllers except KeyError: # If the node does not exist, add it to the dictionary tmp_dict[time] = node @@ -4712,142 +4618,3 @@ def __lt__(self, other): otherwise. """ return self.t < other.t - - @cached_property - def calculate_rail_button_bending_moments(self): - """ - Calculate internal bending moments at rail button attachment points. - - Uses beam theory to determine internal structural moments for stress - analysis of the rail button attachments (fasteners and airframe). - - The bending moment at each button attachment consists of: - 1. Bending from shear force at button contact point: M = S × h - where S is the shear (tangential) force and h is button height - 2. Direct moment contribution from the button's reaction forces - - Assumptions - ----------- - - Rail buttons act as simple supports: provide reaction forces (normal - and shear) but no moment reaction at the rail contact point. - - The rocket acts as a beam supported at two points (rail buttons). - - Bending moments arise from the lever arm effect of reaction forces - and the cantilever moment from button standoff height. - - The bending moment at each button attachment consists of: - 1. Normal force moment: M = N x d, where N is normal reaction force - and d is distance from button to center of dry mass - 2. Shear force cantilever moment: M = S x h, where S is shear force - and h is button standoff height - - Notes - ----- - - Calculated only during the rail phase of flight - - Maximum values use absolute values for worst-case stress analysis - - The bending moments represent internal stresses in the rocket - airframe at the rail button attachment points - - Returns - ------- - tuple - (rail_button1_bending_moment : Function, - max_rail_button1_bending_moment : float, - rail_button2_bending_moment : Function, - max_rail_button2_bending_moment : float) - - Where rail_button1/2_bending_moment are Function objects of time - in N·m, and max values are floats in N·m. - """ - # Check if rail buttons exist - null_moment = Function(0) - if len(self.rocket.rail_buttons) == 0: - warnings.warn( - "Trying to calculate rail button bending moments without " - "rail buttons defined. Setting moments to zero.", - UserWarning, - ) - return (null_moment, 0.0, null_moment, 0.0) - - # Get rail button geometry - rail_buttons_tuple = self.rocket.rail_buttons[0] - # Rail button standoff height - h_button = rail_buttons_tuple.component.button_height - if h_button is None: - warnings.warn( - "Rail button height not defined. Bending moments cannot be " - "calculated. Setting moments to zero.", - UserWarning, - ) - return (null_moment, 0.0, null_moment, 0.0) - upper_button_position = ( - rail_buttons_tuple.component.buttons_distance - + rail_buttons_tuple.position.z - ) - lower_button_position = rail_buttons_tuple.position.z - - # Get center of dry mass (handle both callable and property) - if callable(self.rocket.center_of_dry_mass_position): - cdm = self.rocket.center_of_dry_mass_position(self.rocket._csys) - else: - cdm = self.rocket.center_of_dry_mass_position - - # Distances from buttons to center of dry mass - d1 = abs(upper_button_position - cdm) - d2 = abs(lower_button_position - cdm) - - # forces - N1 = self.rail_button1_normal_force - N2 = self.rail_button2_normal_force - S1 = self.rail_button1_shear_force - S2 = self.rail_button2_shear_force - t = N1.source[:, 0] - - # Calculate bending moments at attachment points - # Primary contribution from shear force acting at button height - # Secondary contribution from normal force creating moment about attachment - m1_values = N2.source[:, 1] * d2 + S1.source[:, 1] * h_button - m2_values = N1.source[:, 1] * d1 + S2.source[:, 1] * h_button - - rail_button1_bending_moment = Function( - np.column_stack([t, m1_values]), - inputs="Time (s)", - outputs="Bending Moment (N·m)", - interpolation="linear", - ) - rail_button2_bending_moment = Function( - np.column_stack([t, m2_values]), - inputs="Time (s)", - outputs="Bending Moment (N·m)", - interpolation="linear", - ) - - # Maximum bending moments (absolute value for stress calculations) - max_rail_button1_bending_moment = float(np.max(np.abs(m1_values))) - max_rail_button2_bending_moment = float(np.max(np.abs(m2_values))) - - return ( - rail_button1_bending_moment, - max_rail_button1_bending_moment, - rail_button2_bending_moment, - max_rail_button2_bending_moment, - ) - - @property - def rail_button1_bending_moment(self): - """Upper rail button bending moment as a Function of time.""" - return self.calculate_rail_button_bending_moments[0] - - @property - def max_rail_button1_bending_moment(self): - """Maximum upper rail button bending moment, in N·m.""" - return self.calculate_rail_button_bending_moments[1] - - @property - def rail_button2_bending_moment(self): - """Lower rail button bending moment as a Function of time.""" - return self.calculate_rail_button_bending_moments[2] - - @property - def max_rail_button2_bending_moment(self): - """Maximum lower rail button bending moment, in N·m.""" - return self.calculate_rail_button_bending_moments[3] diff --git a/tests/animation_verification/__init__.py b/tests/animation_verification/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/tests/animation_verification/main.py b/tests/animation_verification/main.py deleted file mode 100644 index c224f9b25..000000000 --- a/tests/animation_verification/main.py +++ /dev/null @@ -1,76 +0,0 @@ -import os -import traceback - -from rocket_setup import get_calisto_rocket -from rocket_stl import create_rocket_stl - -from rocketpy import Environment, Flight - - - -def run_simulation_and_test_animation(): - print("🚀 Setting up simulation (Calisto Example)...") - - # 1. Setup Environment - env = Environment(latitude=32.990254, longitude=-106.974998, elevation=1400) - env.set_date((2025, 12, 5, 12)) - env.set_atmospheric_model(type="standard_atmosphere") - - # 2. Get Rocket - try: - calisto = get_calisto_rocket() - except Exception as e: - print(f"❌ Failed to configure rocket: {e}") - return - - # 3. Simulate Flight - test_flight = Flight( - rocket=calisto, environment=env, rail_length=5.2, inclination=85, heading=0 - ) - - print(f"✅ Flight simulated successfully! Apogee: {test_flight.apogee:.2f} m") - - # 4. Test Animation Methods - stl_file = "rocket_model.stl" - # Note: Depending on where you run this, you might need to adjust imports - # or ensure create_rocket_stl is available in scope. - create_rocket_stl(stl_file, length=300, radius=50) - - print("\n🎥 Testing animate_trajectory()...") - - try: - test_flight.animate_trajectory( - file_name=stl_file, - stop=15.0, - time_step=0.05, - azimuth=-45, # Rotates view 45 degrees left - elevation=30, # Tilts view 30 degrees up - zoom=1.2, - ) - print("✅ animate_trajectory() executed successfully.") - except Exception as e: - print(f"❌ animate_trajectory() Failed: {e}") - traceback.print_exc() - - print("\n🔄 Testing animate_rotate()...") - - try: - test_flight.animate_rotate( - file_name=stl_file, - time_step=1.0, - azimuth=-45, # Rotates view 45 degrees left - elevation=30, # Tilts view 30 degrees up - zoom=1.2, - ) - print("✅ animate_rotate() executed successfully.") - except Exception as e: - print(f"❌ animate_rotate() Failed: {e}") - traceback.print_exc() - - # Cleanup - if os.path.exists(stl_file): - os.remove(stl_file) - - -if __name__ == "__main__": - run_simulation_and_test_animation() diff --git a/tests/animation_verification/rocket_setup.py b/tests/animation_verification/rocket_setup.py deleted file mode 100644 index 26a5640a7..000000000 --- a/tests/animation_verification/rocket_setup.py +++ /dev/null @@ -1,95 +0,0 @@ -import os - -from rocketpy import Rocket, SolidMotor - -CURRENT_DIR = os.path.dirname(os.path.abspath(__file__)) -ROOT_DIR = os.path.dirname(os.path.dirname(CURRENT_DIR)) -DATA_DIR = os.path.join(ROOT_DIR, "data") - -OFF_DRAG_PATH = os.path.join(DATA_DIR, "rockets/calisto/powerOffDragCurve.csv") -ON_DRAG_PATH = os.path.join(DATA_DIR, "rockets/calisto/powerOnDragCurve.csv") -AIRFOIL_PATH = os.path.join(DATA_DIR, "airfoils/NACA0012-radians.txt") -MOTOR_PATH = os.path.join(DATA_DIR, "motors/cesaroni/Cesaroni_M1670.eng") - - -def get_motor(): - """Locates the motor file and returns a configured SolidMotor object.""" - # We can now point directly to the file without searching - if not os.path.exists(MOTOR_PATH): - raise FileNotFoundError(f"Could not find Cesaroni_M1670.eng at: {MOTOR_PATH}") - - return SolidMotor( - thrust_source=MOTOR_PATH, - dry_mass=1.815, - dry_inertia=(0.125, 0.125, 0.002), - nozzle_radius=33 / 1000, - grain_number=5, - grain_density=1815, - grain_outer_radius=33 / 1000, - grain_initial_inner_radius=15 / 1000, - grain_initial_height=120 / 1000, - grain_separation=5 / 1000, - grains_center_of_mass_position=0.397, - center_of_dry_mass_position=0.317, - nozzle_position=0, - burn_time=3.9, - throat_radius=11 / 1000, - coordinate_system_orientation="nozzle_to_combustion_chamber", - ) - - -def get_calisto_rocket(): - """Configures and returns the Calisto Rocket object.""" - motor = get_motor() - - calisto = Rocket( - radius=127 / 2000, - mass=14.426, - inertia=(6.321, 6.321, 0.034), - power_off_drag=OFF_DRAG_PATH, - power_on_drag=ON_DRAG_PATH, - center_of_mass_without_motor=0, - coordinate_system_orientation="tail_to_nose", - ) - - calisto.add_motor(motor, position=-1.255) - calisto.set_rail_buttons( - upper_button_position=0.0818, - lower_button_position=-0.618, - angular_position=45, - ) - - # Aerodynamic surfaces - calisto.add_nose(length=0.55829, kind="vonKarman", position=1.27) - calisto.add_trapezoidal_fins( - n=4, - root_chord=0.120, - tip_chord=0.060, - span=0.110, - position=-1.04956, - cant_angle=0, - airfoil=(AIRFOIL_PATH, "radians"), - ) - calisto.add_tail( - top_radius=0.0635, bottom_radius=0.0435, length=0.060, position=-1.194656 - ) - - # Parachutes - calisto.add_parachute( - name="Main", - cd_s=10.0, - trigger=800, - sampling_rate=105, - lag=1.5, - noise=(0, 8.3, 0.5), - ) - calisto.add_parachute( - name="Drogue", - cd_s=1.0, - trigger="apogee", - sampling_rate=105, - lag=1.5, - noise=(0, 8.3, 0.5), - ) - - return calisto diff --git a/tests/animation_verification/rocket_stl.py b/tests/animation_verification/rocket_stl.py deleted file mode 100644 index 6aecc6971..000000000 --- a/tests/animation_verification/rocket_stl.py +++ /dev/null @@ -1,234 +0,0 @@ -import math -import struct - -import numpy as np -from stl import mesh # Requires numpy-stl package: pip install numpy-stl - - - -def create_rocket_stl(filename="rocket_model.stl", length=100, radius=10): - """ - Creates a detailed rocket-shaped STL file with proper scale. - - Parameters - ---------- - filename : str - Output filename - length : float - Rocket length in meters - radius : float - Rocket base radius in meters - """ - - # More realistic proportions - nose_length = length * 0.25 - engine_length = length * 0.1 - - # Fin parameters - fin_height = radius * 3 - fin_width = radius * 2 - fin_thickness = radius * 0.2 - num_fins = 4 - - # Calculate number of vertices for smooth curves - num_segments = 36 - vertices = [] - faces = [] - - # Helper function to add triangle - def add_triangle(v1, v2, v3): - nonlocal faces - idx = len(vertices) - vertices.extend([v1, v2, v3]) - faces.append([idx, idx + 1, idx + 2]) - - # Helper function to add quadrilateral as two triangles - def add_quad(v1, v2, v3, v4): - add_triangle(v1, v2, v3) - add_triangle(v1, v3, v4) - - # 1. Create nose cone (pointed ogive) - nose_vertices = [] - for i in range(num_segments + 1): - angle = 2 * math.pi * i / num_segments - x = math.cos(angle) * radius - y = math.sin(angle) * radius - z = length # Tip at full length - nose_vertices.append([x, y, z]) - - # Create nose cone triangles - nose_tip = [0, 0, length] - for i in range(num_segments): - v1 = nose_tip - v2 = nose_vertices[i] - v3 = nose_vertices[(i + 1) % num_segments] - add_triangle(v1, v2, v3) - - # 2. Create main body cylinder - body_z_start = length - nose_length - body_z_end = engine_length - - # Create vertices for top and bottom circles of body - top_circle = [] - bottom_circle = [] - - for i in range(num_segments): - angle = 2 * math.pi * i / num_segments - x = math.cos(angle) * radius - y = math.sin(angle) * radius - - top_circle.append([x, y, body_z_start]) - bottom_circle.append([x, y, body_z_end]) - - # Create body cylinder triangles - for i in range(num_segments): - next_i = (i + 1) % num_segments - - # Side quad - v1 = top_circle[i] - v2 = top_circle[next_i] - v3 = bottom_circle[next_i] - v4 = bottom_circle[i] - add_quad(v1, v2, v3, v4) - - # Top circle (connecting to nose) - v1 = top_circle[i] - v2 = top_circle[next_i] - v3 = nose_vertices[i] - add_triangle(v1, v2, v3) - - # 3. Create engine nozzle (truncated cone) - engine_radius = radius * 0.7 - engine_z_end = 0 - - # Create vertices for engine circles - engine_top_circle = [] - engine_bottom_circle = [] - - for i in range(num_segments): - angle = 2 * math.pi * i / num_segments - x_top = math.cos(angle) * radius - y_top = math.sin(angle) * radius - x_bottom = math.cos(angle) * engine_radius - y_bottom = math.sin(angle) * engine_radius - - engine_top_circle.append([x_top, y_top, body_z_end]) - engine_bottom_circle.append([x_bottom, y_bottom, engine_z_end]) - - # Create engine nozzle triangles - for i in range(num_segments): - next_i = (i + 1) % num_segments - - # Side quad - v1 = engine_top_circle[i] - v2 = engine_top_circle[next_i] - v3 = engine_bottom_circle[next_i] - v4 = engine_bottom_circle[i] - add_quad(v1, v2, v3, v4) - - # Connect to body - v1 = engine_top_circle[i] - v2 = engine_top_circle[next_i] - v3 = bottom_circle[i] - add_triangle(v1, v2, v3) - - # 4. Create fins - for fin_num in range(num_fins): - fin_angle = 2 * math.pi * fin_num / num_fins - - # Inner fin edge (attached to rocket) - inner_x = math.cos(fin_angle) * radius - inner_y = math.sin(fin_angle) * radius - inner_top = [inner_x, inner_y, body_z_end + fin_height * 0.3] - inner_bottom = [inner_x, inner_y, body_z_end - fin_height * 0.7] - - # Outer fin edge - outer_x = math.cos(fin_angle) * (radius + fin_width) - outer_y = math.sin(fin_angle) * (radius + fin_width) - outer_top = [outer_x, outer_y, body_z_end + fin_height * 0.3] - outer_bottom = [outer_x, outer_y, body_z_end - fin_height * 0.7] - - # Fin side that attaches to rocket - add_quad( - inner_top, - inner_bottom, - bottom_circle[fin_num * num_segments // num_fins], - bottom_circle[(fin_num * num_segments // num_fins + 1) % num_segments], - ) - - # Fin surfaces - # Front face - add_quad(inner_top, outer_top, outer_bottom, inner_bottom) - - # Top face - add_triangle( - inner_top, - outer_top, - [ - inner_top[0] + math.cos(fin_angle + math.pi / 2) * fin_thickness * 0.5, - inner_top[1] + math.sin(fin_angle + math.pi / 2) * fin_thickness * 0.5, - inner_top[2], - ], - ) - - # Bottom face - add_triangle( - inner_bottom, - outer_bottom, - [ - inner_bottom[0] - + math.cos(fin_angle + math.pi / 2) * fin_thickness * 0.5, - inner_bottom[1] - + math.sin(fin_angle + math.pi / 2) * fin_thickness * 0.5, - inner_bottom[2], - ], - ) - - # 5. Create bottom cap - center_bottom = [0, 0, engine_z_end] - for i in range(num_segments): - next_i = (i + 1) % num_segments - add_triangle( - center_bottom, engine_bottom_circle[i], engine_bottom_circle[next_i] - ) - - # Convert to numpy arrays - vertices_array = np.array(vertices) - faces_array = np.array(faces) - - # Create the mesh - rocket_mesh = mesh.Mesh(np.zeros(faces_array.shape[0], dtype=mesh.Mesh.dtype)) - for i, face in enumerate(faces_array): - for j in range(3): - rocket_mesh.vectors[i][j] = vertices_array[face[j]] - - # Write the mesh to file - rocket_mesh.save(filename) - print(f"Rocket STL saved to {filename}") - print(f"Total triangles: {len(faces)}") - - -# Alternative version using pure Python STL generation -def create_rocket_stl_simple(filename="rocket_simple.stl", length=100, radius=10): - """Simpler version using pure Python without numpy-stl dependency""" - - def write_stl(faces, filename): - with open(filename, "wb") as f: - # Write 80 byte header - f.write(b"\x00" * 80) - # Write number of faces - f.write(struct.pack(" Date: Tue, 26 May 2026 21:57:19 -0300 Subject: [PATCH 10/12] MNT: fix pylint failures in flight animation methods Reference vedo classes directly instead of capitalized local aliases (invalid-name), drop useless returns, and silence too-many-statements on the animate methods to match the file convention. Lint now scores 10/10 and ruff is clean. Co-Authored-By: Claude Opus 4.7 (1M context) --- rocketpy/plots/flight_plots.py | 32 ++++++++++---------------------- 1 file changed, 10 insertions(+), 22 deletions(-) diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py index 86118037f..a7c750dad 100644 --- a/rocketpy/plots/flight_plots.py +++ b/rocketpy/plots/flight_plots.py @@ -225,16 +225,15 @@ def _create_animation_box(self, start, scale=1.0): center_z = max(center_z, 0.5 * min_box_dim) vedo = import_optional_dependency("vedo") - Box = vedo.Box - return Box( + return vedo.Box( pos=[center_x, center_y, center_z], length=length, width=width, height=height, ).wireframe() - def animate_trajectory( + def animate_trajectory( # pylint: disable=too-many-statements self, file_name=None, start=0, stop=None, time_step=0.1, **kwargs ): """Animate 6-DOF trajectory and attitude using vedo. @@ -259,25 +258,20 @@ def animate_trajectory( vedo = import_optional_dependency("vedo") - Line = vedo.Line - Mesh = vedo.Mesh - Plotter = vedo.Plotter - settings = vedo.settings - file_name = self._resolve_animation_model_path(file_name) stop = self._validate_animation_inputs(file_name, start, stop, time_step) try: - settings.allow_interaction = True + vedo.settings.allow_interaction = True except AttributeError: pass world = self._create_animation_box(start, scale=1.2) - base_rocket = Mesh(file_name).c("green") + base_rocket = vedo.Mesh(file_name).c("green") time_steps = np.arange(start, stop, time_step) trajectory_points = [] - plt = Plotter(axes=1, interactive=False) + plt = vedo.Plotter(axes=1, interactive=False) plt.show(world, "Rocket Trajectory Animation", viewup="z", **kwargs) for t in time_steps: @@ -300,7 +294,7 @@ def animate_trajectory( trajectory_points.append([x_position, y_position, z_position]) actors = [world, rocket] if len(trajectory_points) > 1: - actors.append(Line(trajectory_points, c="k", alpha=0.5)) + actors.append(vedo.Line(trajectory_points, c="k", alpha=0.5)) plt.show(*actors, resetcam=False) @@ -312,9 +306,8 @@ def animate_trajectory( break plt.interactive().close() - return None - def animate_rotate( + def animate_rotate( # pylint: disable=too-many-statements self, file_name=None, start=0, stop=None, time_step=0.1, **kwargs ): """Animate rocket attitude (rotation-only view) using vedo. @@ -339,27 +332,23 @@ def animate_rotate( vedo = import_optional_dependency("vedo") - Mesh = vedo.Mesh - Plotter = vedo.Plotter - settings = vedo.settings - file_name = self._resolve_animation_model_path(file_name) stop = self._validate_animation_inputs(file_name, start, stop, time_step) try: - settings.allow_interaction = True + vedo.settings.allow_interaction = True except AttributeError: pass world = self._create_animation_box(start, scale=0.3) - base_rocket = Mesh(file_name).c("green") + base_rocket = vedo.Mesh(file_name).c("green") time_steps = np.arange(start, stop, time_step) x_start = self.flight.x(start) y_start = self.flight.y(start) z_start = self.flight.z(start) - self.flight.env.elevation - plt = Plotter(axes=1, interactive=False) + plt = vedo.Plotter(axes=1, interactive=False) plt.show(world, "Rocket Rotation Animation", viewup="z", **kwargs) for t in time_steps: @@ -385,7 +374,6 @@ def animate_rotate( break plt.interactive().close() - return None def linear_kinematics_data(self, *, filename=None): # pylint: disable=too-many-statements """Prints out all Kinematics graphs available about the Flight From 8136706be3c825c7b28880ef52e1f42c23c93fb9 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 May 2026 21:58:07 -0300 Subject: [PATCH 11/12] MNT: revert unrelated requirements.txt version pins Keep the PR scoped to the animation feature; requests/pytz pinning was unrelated. Co-Authored-By: Claude Opus 4.7 (1M context) --- requirements.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index e08c9a81d..61a594320 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ numpy>=1.13 scipy>=1.0 matplotlib>=3.9.0 # Released May 15th 2024 netCDF4>=1.6.4 -requests>=2.25.0 -pytz>=2020.1 +requests +pytz simplekml dill From 292eadf25c58ecd1ff3539fd525dff1775d2b759 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 26 May 2026 22:37:52 -0300 Subject: [PATCH 12/12] add feature to the documentation --- docs/notebooks/getting_started_colab.ipynb | 12 +++ docs/user/flight.rst | 109 +++++++++++++++++---- docs/user/installation.rst | 34 +++++-- 3 files changed, 128 insertions(+), 27 deletions(-) diff --git a/docs/notebooks/getting_started_colab.ipynb b/docs/notebooks/getting_started_colab.ipynb index ef0f7e711..154a8cd4c 100644 --- a/docs/notebooks/getting_started_colab.ipynb +++ b/docs/notebooks/getting_started_colab.ipynb @@ -596,6 +596,18 @@ ")" ] }, + { + "cell_type": "markdown", + "source": "## 3D Flight Animation\n\nRocketPy can render an interactive 3D animation of the rocket's trajectory and attitude using [vedo](https://vedo.embl.es/). This feature requires the optional `animation` extra:\n\n```bash\npip install rocketpy[animation]\n```\n\n> **Note:** The interactive animation window opens in a desktop environment. It will not display inside Google Colab or other headless notebook servers — run it locally for the best experience.\n\nTwo animation modes are available:\n\n| Method | What it shows |\n|---|---|\n| `flight.plots.animate_trajectory()` | Rocket moves through 3D space following the simulated trajectory |\n| `flight.plots.animate_rotate()` | Rocket stays centred; only attitude (rotation) is animated |\n\nBoth methods accept an optional `file_name` argument pointing to a custom `.stl` model. When omitted, RocketPy uses a built-in default rocket shape.", + "metadata": {} + }, + { + "cell_type": "code", + "source": "# Install the optional animation dependency (skip if already installed)\n!pip install \"rocketpy[animation]\"\n\n# Animate the full trajectory — rocket moves through 3D space\n# Press Escape or close the window to exit the animation\ntest_flight.plots.animate_trajectory(\n start=0,\n stop=test_flight.t_final,\n time_step=0.05,\n)\n\n# Alternatively, animate only the attitude changes (rocket stays centred)\n# test_flight.plots.animate_rotate(\n# start=0,\n# stop=test_flight.t_final,\n# time_step=0.05,\n# )\n\n# To use your own 3D model, pass its path via file_name:\n# test_flight.plots.animate_trajectory(file_name=\"my_rocket.stl\")", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, { "attachments": {}, "cell_type": "markdown", diff --git a/docs/user/flight.rst b/docs/user/flight.rst index b005477db..a220554ab 100644 --- a/docs/user/flight.rst +++ b/docs/user/flight.rst @@ -457,47 +457,118 @@ Flight Data Plots 3D Flight Animation ~~~~~~~~~~~~~~~~~~~ -RocketPy can animate the simulated flight trajectory and attitude through the -Flight plots layer. +RocketPy can produce real-time interactive 3D animations of the simulated +flight using `vedo `_, a scientific visualization +library built on top of VTK. Two complementary animation modes are provided: + +.. list-table:: + :header-rows: 1 + :widths: 30 70 + + * - Method + - What it shows + * - ``flight.plots.animate_trajectory()`` + - The rocket 3D model moves through space following the simulated + trajectory; a black trail line is drawn behind it. + * - ``flight.plots.animate_rotate()`` + - The rocket 3D model stays centred in the scene; only its attitude + (orientation derived from the quaternion solution) is animated. .. note:: - Install optional animation dependencies first: + The animation window opens on the desktop via VTK. It will **not** render + inside headless environments such as Google Colab. For notebook use, run + the cell on a local Jupyter server or JupyterLab installation. - .. code-block:: bash +**Installation** - pip install rocketpy[animation] +The ``vedo`` dependency is not installed by default. Add the optional extra +before calling either animation method: + +.. code-block:: bash + + pip install rocketpy[animation] + +If ``vedo`` is not available when an animation method is called, RocketPy +raises an :class:`ImportError` with the above install command embedded in the +message. + +**animate_trajectory — full 6-DOF trajectory animation** .. code-block:: python - # Fast start using RocketPy's built-in default STL model + # Quickstart: uses RocketPy's built-in default STL rocket model + flight.plots.animate_trajectory() + + # Customise the time window and frame rate flight.plots.animate_trajectory( - start=0.0, - stop=min(flight.t_final, 20.0), - time_step=0.05, + start=0.0, # start time in seconds (default: 0) + stop=flight.t_final, # end time in seconds (default: t_final) + time_step=0.05, # seconds per frame (default: 0.1) ) - # Or provide your own STL model file + # Provide your own 3D model (any STL file) flight.plots.animate_trajectory( - file_name="rocket.stl", + file_name="my_rocket.stl", start=0.0, stop=flight.t_final, time_step=0.05, - azimuth=45, - elevation=20, ) - # Keep rocket centered and animate only attitude changes +**animate_rotate — attitude-only animation** + +Useful for inspecting roll, pitch, and yaw behaviour without the distraction of +the trajectory translation. The rocket mesh is fixed at its position at +``start`` and only rotated according to the quaternion solution. + +.. code-block:: python + flight.plots.animate_rotate( - file_name="rocket.stl", start=0.0, - stop=min(flight.t_final, 20.0), + stop=flight.t_final, time_step=0.05, ) -Both methods validate the selected time interval and STL path before rendering. -If ``vedo`` is not installed, RocketPy raises an informative ``ImportError`` -with installation instructions. +**Parameters** + +Both methods share the same signature: + +.. list-table:: + :header-rows: 1 + :widths: 20 15 65 + + * - Parameter + - Default + - Description + * - ``file_name`` + - ``None`` + - Path to a ``.stl`` model file. When ``None``, the built-in default + rocket shape packaged with RocketPy is used. + * - ``start`` + - ``0`` + - Animation start time in seconds. Must be within + ``[0, flight.t_final]``. + * - ``stop`` + - ``None`` + - Animation end time in seconds. ``None`` defaults to + ``flight.t_final``. + * - ``time_step`` + - ``0.1`` + - Duration of each frame in seconds. Smaller values produce smoother + animations at the cost of longer render times. Must be > 0. + * - ``**kwargs`` + - — + - Additional keyword arguments forwarded to ``vedo.Plotter.show`` + (e.g. ``viewup``, ``azimuth``, ``elevation``). + +**Tips** + +- A ``time_step`` of ``0.05`` (20 fps) is a good balance between smoothness + and performance for flights lasting tens of seconds. +- Press **Escape** or close the vedo window to exit the animation loop early. +- Both methods validate ``start``, ``stop``, ``time_step``, and the STL path + before any rendering begins, raising a :class:`ValueError` with a descriptive + message on invalid input. Forces and Moments ~~~~~~~~~~~~~~~~~~ diff --git a/docs/user/installation.rst b/docs/user/installation.rst index 72ba1f42d..a08c65aeb 100644 --- a/docs/user/installation.rst +++ b/docs/user/installation.rst @@ -138,22 +138,40 @@ To update Scipy and install netCDF4 using Conda, the following code is used: Optional Packages ^^^^^^^^^^^^^^^^^ -The EnvironmentAnalysis class requires a few extra packages to be installed. -In case you want to use this class, you will need to install the following packages: +RocketPy has several optional feature sets that can be installed individually. -- `timezonefinder` : to allow for automatic timezone detection, -- `windrose` : to allow for windrose plots, -- `ipywidgets` : to allow for GIFs generation, -- `jsonpickle` : to allow for saving and loading of class instances. +**Environment Analysis** — extra plots and tools for the +:class:`rocketpy.EnvironmentAnalysis` class: -You can install all these packages by simply running the following lines in your preferred terminal: +- `timezonefinder` : automatic timezone detection +- `windrose` : windrose plots +- `ipywidgets` : GIF generation +- `jsonpickle` : saving and loading class instances .. code-block:: shell pip install rocketpy[env_analysis] +**3D Flight Animation** — interactive 3D animations of rocket trajectory and +attitude using `vedo `_ (requires a desktop environment): -Alternatively, you can instal all extra packages by running the following line in your preferred terminal: +.. code-block:: shell + + pip install rocketpy[animation] + +Once installed, you can render animations from a :class:`rocketpy.Flight` object: + +.. code-block:: python + + # Animate rocket moving through 3D space + flight.plots.animate_trajectory(start=0, stop=flight.t_final, time_step=0.05) + + # Animate attitude changes only (rocket stays centred) + flight.plots.animate_rotate(start=0, stop=flight.t_final, time_step=0.05) + +See :ref:`flightusage` for full details and parameter descriptions. + +**All extras** — install every optional dependency at once: .. code-block:: shell