Getting Started#

Here we go through a complete rocket trajectory simulation to get you started.

Also, a more in detail complete description of the simulation can also be found in the First Simulation Section of the RocketPy Documentation.

If you want, you can open this notebook in Google Colab by clicking the badge below.

Open In Colab

Let’s start by importing the rocketpy module.

[1]:
%load_ext autoreload
%autoreload 2
[2]:
from rocketpy import Environment, SolidMotor, Rocket, Flight

If you are using Jupyter Notebooks, it is recommended to run the following line to make matplotlib plots which will be shown later interactive and higher quality.

[ ]:
%matplotlib widget

Setting Up a Simulation#

Creating an Environment for Spaceport America#

The Environment class is used to define the atmosphere, the winds, and the gravity models.

You can find more information about the Environment class in the Environment Class Usage Docs.

[3]:
env = Environment(latitude=32.990254, longitude=-106.974998, elevation=1400)

To get weather data from the GFS forecast, available online, we run the following lines.

First, we set tomorrow’s date.

[4]:
import datetime

tomorrow = datetime.date.today() + datetime.timedelta(days=1)

env.set_date(
    (tomorrow.year, tomorrow.month, tomorrow.day, 12)
)  # Hour given in UTC time

Then, we tell env to use a GFS forecast to get the atmospheric conditions for flight.

Don’t mind the warning, it just means that not all variables, such as wind speed or atmospheric temperature, are available at all altitudes given by the forecast.

[5]:
env.set_atmospheric_model(type="Forecast", file="GFS")

We can see what the weather will look like by calling the info method!

[6]:
env.max_expected_height = 5000  # adjust the plots to this height
env.info()

Gravity Details

Acceleration of gravity at surface level:    9.7911 m/s²
Acceleration of gravity at   5.000 km (ASL): 9.7802 m/s²


Launch Site Details

Launch Date: 2023-10-10 12:00:00 UTC
Launch Site Latitude: 32.99025°
Launch Site Longitude: -106.97500°
Reference Datum: SIRGAS2000
Launch Site UTM coordinates: 315468.64 W    3651938.65 N
Launch Site UTM zone: 13S
Launch Site Surface Elevation: 1471.5 m


Atmospheric Model Details

Atmospheric Model Type: Forecast
Forecast Maximum Height: 5.000 km
Forecast Time Period: From  2023-10-09 12:00:00  to  2023-10-25 12:00:00  UTC
Forecast Hour Interval: 3  hrs
Forecast Latitude Range: From  -90.0 ° To  90.0 °
Forecast Longitude Range: From  0.0 ° To  359.75 °


Surface Atmospheric Conditions

Surface Wind Speed: 1.74 m/s
Surface Wind Direction: 168.62°
Surface Wind Heading: 335.65°
Surface Pressure: 848.20 hPa
Surface Temperature: 291.05 K
Surface Air Density: 1.015 kg/m³
Surface Speed of Sound: 342.00 m/s


Earth Model Details

Earth Radius at Launch site: 6371.83 km
Semi-major Axis: 6378.14 km
Semi-minor Axis: 6356.75 km
Flattening: 0.0034


Atmospheric Model Plots

../_images/notebooks_getting_started_14_1.png

Creating a Motor#

A solid rocket motor is used in this case. To create a motor, the SolidMotor class is used and the required arguments are given.

The SolidMotor class requires the user to have a thrust curve ready. This can come either from a .eng file for a commercial motor, such as below, or a .csv file from a static test measurement.

Besides the thrust curve, other parameters such as grain properties and nozzle dimensions must also be given.

See Solid Motor Class Usage Docs for more information.

[7]:
Pro75M1670 = SolidMotor(
    thrust_source="../../data/motors/Cesaroni_M1670.eng",
    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",
)

Pay special attention to position related parameters: More details on Positions and Coordinate Systems

To see what our thrust curve looks like, along with other import properties, we invoke the info method yet again. You may try the all_info method if you want more information all at once!

[8]:
Pro75M1670.info()
Nozzle Details
Nozzle Radius: 0.033 m
Nozzle Throat Radius: 0.011 m

Grain Details
Number of Grains: 5
Grain Spacing: 0.005 m
Grain Density: 1815 kg/m3
Grain Outer Radius: 0.033 m
Grain Inner Radius: 0.015 m
Grain Height: 0.12 m
Grain Volume: 0.000 m3
Grain Mass: 0.591 kg

Motor Details
Total Burning Time: 3.9 s
Total Propellant Mass: 2.956 kg
Average Propellant Exhaust Velocity: 2038.745 m/s
Average Thrust: 1545.218 N
Maximum Thrust: 2200.0 N at 0.15 s after ignition.
Total Impulse: 6026.350 Ns

../_images/notebooks_getting_started_19_1.png

Creating a Rocket#

A rocket is composed of several components. Namely, we must have a motor (good thing we have the Pro75M1670 ready), a couple of aerodynamic surfaces (nose cone, fins and tail) and parachutes (if we are not launching a missile).

You can find more information about the Rocket class in the Rocket Class Usage Docs.

Let’s start by initializing our rocket, named Calisto, entering inertia properties, some dimensions and drag curves.

Pay special attention to position related parameters: More details on Positions and Coordinate Systems

[9]:
calisto = Rocket(
    radius=127 / 2000,
    mass=14.426,
    inertia=(6.321, 6.321, 0.034),
    power_off_drag="../../data/calisto/powerOffDragCurve.csv",
    power_on_drag="../../data/calisto/powerOnDragCurve.csv",
    center_of_mass_without_motor=0,
    coordinate_system_orientation="tail_to_nose",
)

rail_buttons = calisto.set_rail_buttons(
    upper_button_position=0.0818,
    lower_button_position=-0.618,
    angular_position=45,
)

To add the motor to our rocket we need only inform what motor we are adding (Pro75M1670) and inform the position, in meters, of the motor’s nozzle exit area relative to the previously defined coordinate system.

[10]:
calisto.add_motor(Pro75M1670, position=-1.255)

Adding Aerodynamic Surfaces#

Now we define the aerodynamic surfaces. They are really straight forward with special attention needed only for the position values. Here is a quick guide:

  • The positions given must be relative to the same coordinate system as the rockets;

  • Position of the Nosecone refers to the tip of the nose;

  • Position of fins refers to the point belonging to the root chord which is highest in the rocket coordinate system;

  • Position of the tail the point belonging to the tail which is highest in the rocket coordinate system.

See more details in Positions and Coordinate Systems

[11]:
nose_cone = calisto.add_nose(length=0.55829, kind="vonKarman", position=1.278)

fin_set = calisto.add_trapezoidal_fins(
    n=4,
    root_chord=0.120,
    tip_chord=0.060,
    span=0.110,
    position=-1.04956,
    cant_angle=0.5,
    airfoil=("../../data/calisto/NACA0012-radians.csv", "radians"),
)

tail = calisto.add_tail(
    top_radius=0.0635, bottom_radius=0.0435, length=0.060, position=-1.194656
)

To see all information regarding the rocket we just defined we run:

[12]:
calisto.all_info()

Inertia Details

Rocket Mass: 14.426 kg
Rocket Dry Mass: 16.241 kg (With Motor)
Rocket Mass: 19.197 kg (With Propellant)
Rocket Inertia (with motor, but without propellant) 11: 7.864 kg*m2
Rocket Inertia (with motor, but without propellant) 22: 7.864 kg*m2
Rocket Inertia (with motor, but without propellant) 33: 0.036 kg*m2
Rocket Inertia (with motor, but without propellant) 12: 0.000 kg*m2
Rocket Inertia (with motor, but without propellant) 13: 0.000 kg*m2
Rocket Inertia (with motor, but without propellant) 23: 0.000 kg*m2


Geometrical Parameters

Rocket Maximum Radius: 0.0635 m
Rocket Frontal Area: 0.012668 m2

Rocket Distances
Rocket Center of Dry Mass - Center of Mass without Motor: 0.105 m
Rocket Center of Dry Mass - Nozzle Exit: 1.150 m
Rocket Center of Dry Mass - Center of Propellant Mass: 0.753 m
Rocket Center of Mass - Rocket Loaded Center of Mass: 0.116 m


Aerodynamics Lift Coefficient Derivatives

Nose Cone Lift Coefficient Derivative: 2.000/rad
Fins Lift Coefficient Derivative: 6.280/rad
Tail Lift Coefficient Derivative: -1.061/rad

Center of Pressure

Nose Cone Center of Pressure position: 0.999 m
Fins Center of Pressure position: -1.100 m
Tail Center of Pressure position: -1.223 m

Stability

Center of Mass position (time=0): -0.221 m
Initial Static Margin (mach=0, time=0): 2.199 c
Final Static Margin (mach=0, time=burn_out): 3.112 c
Rocket Center of Mass (time=0) - Center of Pressure (mach=0): 0.279 m


Mass Plots
----------------------------------------
../_images/notebooks_getting_started_29_1.png
../_images/notebooks_getting_started_29_2.png

Aerodynamics Plots
----------------------------------------
Drag Plots
--------------------
../_images/notebooks_getting_started_29_4.png
../_images/notebooks_getting_started_29_5.png

Stability Plots
--------------------
../_images/notebooks_getting_started_29_7.png
../_images/notebooks_getting_started_29_8.png

Thrust-to-Weight Plot
----------------------------------------
../_images/notebooks_getting_started_29_10.png

Adding Parachutes#

Finally, we have parachutes! Calisto will have two parachutes, Drogue and Main. The Drogue parachute will open at apogee while the Main parachute will open at 800m above ground level

For more details see Adding Parachutes

[13]:
Main = calisto.add_parachute(
    "Main",
    cd_s=10.0,
    trigger=800,
    sampling_rate=105,
    lag=1.5,
    noise=(0, 8.3, 0.5),
)

Drogue = calisto.add_parachute(
    "Drogue",
    cd_s=1.0,
    trigger="apogee",
    sampling_rate=105,
    lag=1.5,
    noise=(0, 8.3, 0.5),
)

Just be careful if you run this last cell multiple times! If you do so, your rocket will end up with lots of parachutes which activate together, which may cause problems during the flight simulation. We advise you to re-run all cells which define our rocket before running this, preventing unwanted old parachutes. Alternatively, you can run the following lines to remove parachutes.

Calisto.parachutes.remove(Drogue)
Calisto.parachutes.remove(Main)

Simulating a Flight#

Simulating a flight trajectory is as simple as initializing a Flight class object givin the rocket and environnement set up above as inputs. The launch rail inclination and heading are also given here.

[14]:
test_flight = Flight(
    rocket=calisto, environment=env, rail_length=5.2, inclination=85, heading=0
)

Analyzing the Results#

RocketPy gives you many plots, thats for sure! They are divided into sections to keep them organized. Alternatively, see the Flight class documentation to see how to get plots for specific variables only, instead of all of them at once.

[15]:
test_flight.all_info()

Initial Conditions

Position - x: 0.00 m | y: 0.00 m | z: 1471.47 m
Velocity - Vx: 0.00 m/s | Vy: 0.00 m/s | Vz: 0.00 m/s
Attitude - e0: 0.999 | e1: -0.044 | e2: -0.000 | e3: 0.000
Euler Angles - Spin φ : 0.00° | Nutation θ: -5.00° | Precession ψ: 0.00°
Angular Velocity - ω1: 0.00 rad/s | ω2: 0.00 rad/s| ω3: 0.00 rad/s


Surface Wind Conditions

Frontal Surface Wind Speed: 1.68 m/s
Lateral Surface Wind Speed: 0.32 m/s


Launch Rail

Launch Rail Length: 5.2  m
Launch Rail Inclination: 85.00°
Launch Rail Heading: 0.00°


Rail Departure State

Rail Departure Time: 0.368 s
Rail Departure Velocity: 26.207 m/s
Rail Departure Stability Margin: 2.276 c
Rail Departure Angle of Attack: 3.738°
Rail Departure Thrust-Weight Ratio: 10.152
Rail Departure Reynolds Number: 1.867e+05


Burn out State

Burn out time: 3.900 s
Altitude at burn out: 659.936 m (AGL)
Rocket velocity at burn out: 280.945 m/s
Freestream velocity at burn out: 280.804 m/s
Mach Number at burn out: 0.825
Kinetic energy at burn out: 6.410e+05 J


Apogee State

Apogee Altitude: 4851.421 m (ASL) | 3379.955 m (AGL)
Apogee Time: 26.211 s
Apogee Freestream Speed: 17.793 m/s


Parachute Events

Drogue Ejection Triggered at: 26.219 s
Drogue Parachute Inflated at: 27.719 s
Drogue Parachute Inflated with Freestream Speed of: 22.764 m/s
Drogue Parachute Inflated at Height of: 3369.148 m (AGL)
Main Ejection Triggered at: 167.590 s
Main Parachute Inflated at: 169.090 s
Main Parachute Inflated with Freestream Speed of: 17.391 m/s
Main Parachute Inflated at Height of: 775.506 m (AGL)


Impact Conditions

X Impact: 1107.036 m
Y Impact: 875.712 m
Latitude: 32.9981279°
Longitude: -106.9631289°
Time of Impact: 312.356 s
Velocity at Impact: -5.278 m/s


Stability Margin

Maximum Stability Margin: 3.682 c at 3.91 s
Minimum Stability Margin: 2.199 c at 0.00 s


Maximum Values

Maximum Speed: 286.800 m/s at 3.39 s
Maximum Mach Number: 0.841 Mach at 3.39 s
Maximum Reynolds Number: 1.957e+06 at 3.32 s
Maximum Dynamic Pressure: 3.958e+04 Pa at 3.35 s
Maximum Acceleration During Motor Burn: 105.255 m/s² at 0.15 s
Maximum Gs During Motor Burn: 10.733 g at 0.15 s
Maximum Acceleration After Motor Burn: 10.202 m/s² at 20.14 s
Maximum Gs After Motor Burn: 1.040 g at 20.14 s
Maximum Stability Margin: 3.682 c at 3.91 s
Maximum Upper Rail Button Normal Force: 0.813 N
Maximum Upper Rail Button Shear Force: 1.163 N
Maximum Lower Rail Button Normal Force: 0.358 N
Maximum Lower Rail Button Shear Force: 0.541 N


Numerical Integration Settings

Maximum Allowed Flight Time: 600.000000 s
Maximum Allowed Time Step: inf s
Minimum Allowed Time Step: 0.000000e+00 s
Relative Error Tolerance:  1e-06
Absolute Error Tolerance:  [0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 1e-06, 1e-06, 1e-06, 1e-06, 0.001, 0.001, 0.001]
Allow Event Overshoot:  True
Terminate Simulation on Apogee:  False
Number of Time Steps Used:  1343
Number of Derivative Functions Evaluation:  2814
Average Function Evaluations per Time Step: 2.095309



Trajectory 3d Plot

../_images/notebooks_getting_started_37_1.png


Trajectory Kinematic Plots

../_images/notebooks_getting_started_37_3.png


Angular Position Plots

../_images/notebooks_getting_started_37_5.png


Path, Attitude and Lateral Attitude Angle plots

../_images/notebooks_getting_started_37_7.png


Trajectory Angular Velocity and Acceleration Plots

../_images/notebooks_getting_started_37_9.png


Aerodynamic Forces Plots

../_images/notebooks_getting_started_37_11.png


Rail Buttons Forces Plots

../_images/notebooks_getting_started_37_13.png


Trajectory Energy Plots

../_images/notebooks_getting_started_37_15.png


Trajectory Fluid Mechanics Plots

../_images/notebooks_getting_started_37_17.png


Trajectory Stability and Control Plots

../_images/notebooks_getting_started_37_19.png


Rocket and Parachute Pressure Plots

../_images/notebooks_getting_started_37_21.png

Parachute:  Main
../_images/notebooks_getting_started_37_23.png
../_images/notebooks_getting_started_37_24.png
../_images/notebooks_getting_started_37_25.png

Parachute:  Drogue
../_images/notebooks_getting_started_37_27.png
../_images/notebooks_getting_started_37_28.png
../_images/notebooks_getting_started_37_29.png

Export Flight Trajectory to a .kml file so it can be opened on Google Earth

[16]:
test_flight.export_kml(
    file_name="trajectory.kml",
    extrude=True,
    altitude_mode="relative_to_ground",
)
File  trajectory.kml  saved with success!

Using Simulation for Design#

Here, we go through a couple of examples which make use of RocketPy in cool ways to help us design our rocket.

Apogee as a Function of Mass#

This one is a classic one! We always need to know how much our rocket’s apogee will change when our payload gets heavier.

[17]:
from rocketpy.utilities import apogee_by_mass

apogee_by_mass(flight=test_flight, min_mass=5, max_mass=20, points=10, plot=True)
../_images/notebooks_getting_started_42_0.png
[17]:
'Function from R1 to R1 : (Rocket Mass without motor (kg)) → (Apogee AGL (m))'

Out of Rail Speed as a Function of Mass#

Lets make a really important plot. Out of rail speed is the speed our rocket has when it is leaving the launch rail. This is crucial to make sure it can fly safely after leaving the rail. A common rule of thumb is that our rocket’s out of rail speed should be 4 times the wind speed so that it does not stall and become unstable.

[18]:
from rocketpy.utilities import liftoff_speed_by_mass

liftoff_speed_by_mass(flight=test_flight, min_mass=5, max_mass=20, points=10, plot=True)
../_images/notebooks_getting_started_44_0.png
[18]:
'Function from R1 to R1 : (Rocket Mass without motor (kg)) → (Out of Rail Speed (m/s))'

Dynamic Stability Analysis#

Ever wondered how static stability translates into dynamic stability? Different static margins result in different dynamic behavior, which also depends on the rocket’s rotational inertial.

Let’s make use of RocketPy’s helper class called Function to explore how the dynamic stability of Calisto varies if we change the fins span by a certain factor.

[19]:
# Helper class
from rocketpy import Function
import copy

# Prepare a copy of the rocket
calisto2 = copy.deepcopy(calisto)

# Prepare Environment Class
custom_env = Environment()
custom_env.set_atmospheric_model(type="custom_atmosphere", wind_v=-5)

# Simulate Different Static Margins by Varying Fin Position
simulation_results = []

for factor in [-0.5, -0.2, 0.1, 0.4, 0.7]:
    # Modify rocket fin set by removing previous one and adding new one
    calisto2.aerodynamic_surfaces.pop(-1)

    fin_set = calisto2.add_trapezoidal_fins(
        n=4,
        root_chord=0.120,
        tip_chord=0.040,
        span=0.100,
        position=-1.04956 * factor,
    )
    # Simulate
    print(
        "Simulating Rocket with Static Margin of {:1.3f}->{:1.3f} c".format(
            calisto2.static_margin(0),
            calisto2.static_margin(calisto2.motor.burn_out_time),
        )
    )
    test_flight = Flight(
        rocket=calisto2,
        environment=custom_env,
        rail_length=5.2,
        inclination=90,
        heading=0,
        max_time_step=0.01,
        max_time=5,
        terminate_on_apogee=True,
        verbose=True,
    )
    # Store Results
    static_margin_at_ignition = calisto2.static_margin(0)
    static_margin_at_out_of_rail = calisto2.static_margin(test_flight.out_of_rail_time)
    static_margin_at_steady_state = calisto2.static_margin(test_flight.t_final)
    simulation_results += [
        (
            test_flight.attitude_angle,
            "{:1.2f} c | {:1.2f} c | {:1.2f} c".format(
                static_margin_at_ignition,
                static_margin_at_out_of_rail,
                static_margin_at_steady_state,
            ),
        )
    ]

Function.compare_plots(
    simulation_results,
    lower=0,
    upper=1.5,
    xlabel="Time (s)",
    ylabel="Attitude Angle (deg)",
)
Simulating Rocket with Static Margin of 0.113->0.846 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 1.064->1.796 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 2.014->2.747 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 2.964->3.697 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 3.914->4.647 c
Simulation Completed at Time: 5.0000 s
../_images/notebooks_getting_started_46_1.png

Characteristic Frequency Calculation#

Here we analyse the characteristic frequency of oscillation of our rocket just as it leaves the launch rail. Note that when we ran test_flight.all_info(), one of the plots already showed us the frequency spectrum of our flight. Here, however, we have more control of what we are plotting.

[20]:
import numpy as np
import matplotlib.pyplot as plt

# Simulate first 5 seconds of Flight
flight = Flight(
    rocket=calisto,
    environment=env,
    rail_length=5.2,
    inclination=90,
    heading=0,
    max_time_step=0.01,
    max_time=5,
)

# Perform a Fourier Analysis
Fs = 100.0
# sampling rate
Ts = 1.0 / Fs
# sampling interval
t = np.arange(1, 400, Ts)  # time vector
ff = 5
# frequency of the signal
y = flight.attitude_angle(t) - np.mean(flight.attitude_angle(t))
n = len(y)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n // 2)]  # one side frequency range
Y = np.fft.fft(y) / n  # fft computing and normalization
Y = Y[range(n // 2)]

# Create the plot
fig, ax = plt.subplots(2, 1)
ax[0].plot(t, y)
ax[0].set_xlabel("Time")
ax[0].set_ylabel("Signal")
ax[0].set_xlim((0, 5))
ax[0].grid()
ax[1].plot(frq, abs(Y), "r")  # plotting the spectrum
ax[1].set_xlabel("Freq (Hz)")
ax[1].set_ylabel("|Y(freq)|")
ax[1].set_xlim((0, 5))
ax[1].grid()
plt.subplots_adjust(hspace=0.5)
plt.show()
../_images/notebooks_getting_started_48_0.png