Getting Started - RocketPy in Colab

We start by setting up our environment. To run this notebook, we will need:

• RocketPy

• netCDF4 (to get weather forecasts)

• Data files (we will clone RocketPy’s repository for these)

Therefore, let’s run the following lines of code:

%pip install rocketpy netCDF4
%git clone https://github.com/giovaniceotto/RocketPy.git

import os

os.chdir("RocketPy/docs/notebooks")

Now we can start!

Here we go through a simplified rocket trajectory simulation to get you started. Let’s start by importing the rocketpy module.

%load_ext autoreload
%autoreload 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.

%config InlineBackend.figure_formats = ['svg']
%matplotlib inline

Setting Up a Simulation

Creating an Environment for Spaceport America

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.

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.

env.set_atmospheric_model(type="Forecast", file="GFS")

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

env.info()

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.

Pro75M1670 = SolidMotor(
    thrust_source="../../data/motors/Cesaroni_M1670.eng",
    dry_mass=1.815,
    dry_inertia=(0.125, 0.125, 0.002),
    center_of_dry_mass=0.317,
    grains_center_of_mass_position=0.397,
    burn_time=3.9,
    grain_number=5,
    grain_separation=5 / 1000,
    grain_density=1815,
    grain_outer_radius=33 / 1000,
    grain_initial_inner_radius=15 / 1000,
    grain_initial_height=120 / 1000,
    nozzle_radius=33 / 1000,
    throat_radius=11 / 1000,
    interpolation_method="linear",
    nozzle_position=0,
    coordinate_system_orientation="nozzle_to_combustion_chamber",
)

The arguments ``nozzle_position`` and ``coordinate_system_orientation`` need to be handled with care since the coordinate system origin is chosen by the user. This means that ``nozzle_position`` is given with respect to an arbitrary reference. The definition of these arguments are really helpful:

• coordinate_system_orientation : string, optional

  Orientation of the motor's coordinate system. The coordinate system is
  defined by the motor's axis of symmetry. The origin of the coordinate
  system may be placed anywhere along such axis, such as at the nozzle area,
  and must be kept the same for all other positions specified.
  Options are "nozzle_to_combustion_chamber" and "combustion_chamber_to_nozzle".
  Default is "nozzle_to_combustion_chamber".
  

• nozzle_position : float, optional

  Motor's nozzle outlet position in meters. More specifically, the coordinate
  of the nozzle outlet specified in the motor's coordinate system.
  Default is 0, in which case the origin of the motor's coordinate system
  is placed at the motor's nozzle outlet.
  

• grains_center_of_mass_position : float

  Position of the center of mass of the grains in meters. More specifically,
  the coordinate of the center of mass specified in the motor's coordinate
  

Here is a useful schematic for explanation:

In the Calisto example, the coordinate system origin we chose is at the center of dry mass of the rocket. Meaning the nozzle is 1.255 meters away from it. The minus sign comes from the coordinate_system_orientation argument, which is responsible for defining the positive direction of the coordinate system. Here is a representation of the reference used in Calisto:

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!

Pro75M1670.info()

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).

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

Here are some definitions that might prove useful. Check out the documentation for more.

• ``mass`` : int, float

  Unloaded rocket total mass (without propellant) in kg.
  

• ``inertia_i`` : int, float

  Unloaded rocket lateral (perpendicular to axis of symmetry)
  moment of inertia (without propellant) in kg m^2.
  

• ``inertia_z`` : int, float

  Unloaded rocket axial moment of inertia (without propellant) in kg m^2.
  

• ``power_off_drag`` : int, float, callable, string, array

  Rocket's drag coefficient when the motor is off.  If int or float
  is given, it is assumed constant. If callable, string or array is
  given, it must be a function of Mach number only.
  

• ``power_on_drag`` : int, float, callable, string, array

  Rocket's drag coefficient when the motor is on. If int or float is
  given, it is assumed constant. If callable, string or array is
  given, it must be a function of Mach number only.
  

calisto = Rocket(
    radius=127 / 2000,
    mass=19.197 - 2.956 - 1.815,
    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.2 - 0.1182359460624346,
    lower_button_position=-0.5 - 0.1182359460624346,
    angular_position=45,
)

Similar to the motor, the last two arguments, ``center_of_dry_mass_position`` and ``coordinate_system_orientation``, need special care. Here, the coordinate system origin is chosen by the user. This means that ``center_of_dry_mass_position`` is given with respect to an arbitrary reference. ``coordinate_system_orientation``. The definitions are also helpful:

• coordinate_system_orientation : string, optional

  String defining the orientation of the rocket's coordinate system. The
  coordinate system is defined by the rocket's axis of symmetry. The system's
  origin may be placed anywhere along such axis and must be kept the same for
  all other positions specified.
  The two options available are: "tail_to_nose" and "nose_to_tail". The first
  defines the coordinate system with the rocket's axis of symmetry pointing
  from the rocket's tail to the rocket's nose cone. The second option defines
  the coordinate system with the rocket's axis of symmetry pointing from the
  rocket's nose cone to the rocket's tail. Default is "tail_to_nose".
  

• center_of_dry_mass_position : int, float, optional

  Position, in m, of the rocket's center of dry mass (i.e. center of mass
  without propellant) relative to the rocket's coordinate system.
  Default is 0, which means the center of dry mass is chosen as the origin.
  

And a schematic for explanation:

In the Calisto example, the coordinate system origin we chose is at the center of dry mass of the rocket. Meaning the nozzle is 1.255 meters away from it. The minus sign comes from the coordinate_system_orientation argument, which is responsible for defining the positive direction of the coordinate system. The positions of the motor, the fins, the nose cone, and the tail must also be given by

It is important to note that the position of the origin of the coordinate system can be anywhere, as long as all positions are given coherently. This means that you could input all positions based on the nosecone’s tip, or based on the coordinates of a CAD model.

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.

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.

NoseCone = calisto.add_nose(
    length=0.55829, kind="vonKarman", position=0.71971 + 0.55829
)

fin_set = calisto.add_trapezoidal_fins(
    n=4,
    root_chord=0.120,
    tip_chord=0.040,
    span=0.100,
    position=-1.04956,
    cant_angle=0,
    radius=None,
    airfoil=None,
)

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:

calisto.all_info()

Adding Parachutes

Finally, we have parachutes! Calisto will have two parachutes, Drogue and Main.

Both parachutes are activated by some special algorithm, which is usually really complex and a trade secret. Most algorithms are based on pressure sampling only, while some also use acceleration info.

RocketPy allows you to define a trigger function which will decide when to activate the ejection event for each parachute. This trigger function is supplied with pressure measurement at a predefined sampling rate. This pressure signal is usually noisy, so artificial noise parameters can be given. Call help(Rocket.add_parachute) for more details. Furthermore, the trigger function also receives the complete state vector of the rocket, allowing us to use velocity, acceleration or even attitude to decide when the parachute event should be triggered.

Here, we define our trigger functions rather simply using Python. However, you can call the exact code which will fly inside your rocket as well.

def drogue_trigger(p, h, y):
    # p = pressure considering parachute noise signal
    # h = height above ground level considering parachute noise signal
    # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]

    # activate drogue when vz < 0 m/s.
    return True if y[5] < 0 else False


def main_trigger(p, h, y):
    # p = pressure considering parachute noise signal
    # h = height above ground level considering parachute noise signal
    # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]

    # activate main when vz < 0 m/s and z < 800 m
    return True if y[5] < 0 and h < 800 else False


Main = calisto.add_parachute(
    "Main",
    cd_s=10.0,
    trigger=main_trigger,
    sampling_rate=105,
    lag=1.5,
    noise=(0, 8.3, 0.5),
)

Drogue = calisto.add_parachute(
    "Drogue",
    cd_s=1.0,
    trigger=drogue_trigger,
    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.

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.

test_flight.all_info()

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

test_flight.export_kml(
    file_name="trajectory.kml",
    extrude=True,
    altitude_mode="relative_to_ground",
)

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.

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.

# 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)",
)

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.

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[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))
plt.subplots_adjust(hspace=0.5)
plt.show()

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.

def apogee(mass):
    # Prepare a copy of the rocket to avoid modifying the original
    mock_rocket = copy.deepcopy(calisto)

    # Modify the mass
    mock_rocket.mass = mass
    mock_rocket.evaluate_dry_mass()
    mock_rocket.evaluate_total_mass()
    mock_rocket.evaluate_center_of_dry_mass()
    mock_rocket.evaluate_center_of_mass()
    mock_rocket.evaluate_reduced_mass()
    mock_rocket.evaluate_thrust_to_weight()
    mock_rocket.evaluate_static_margin()

    # Simulate Flight until Apogee
    mock_flight = Flight(
        rocket=mock_rocket,
        environment=env,
        rail_length=5.2,
        inclination=85,
        heading=0,
        terminate_on_apogee=True,
        # equations_of_motion="solid_propulsion",
    )
    return mock_flight.apogee - env.elevation


apogeebymass = Function(apogee, inputs="Mass (kg)", outputs="Estimated Apogee (m)")
apogeebymass.plot(8, 20, 5)

Out of Rail Speed as a Function of Mass

To finish off, 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.

def speed(mass):
    # Prepare a copy of the rocket to avoid modifying the original
    mock_rocket = copy.deepcopy(calisto)

    # Modify the mass
    mock_rocket.mass = mass
    mock_rocket.evaluate_dry_mass()
    mock_rocket.evaluate_total_mass()
    mock_rocket.evaluate_center_of_dry_mass()
    mock_rocket.evaluate_center_of_mass()
    mock_rocket.evaluate_reduced_mass()
    mock_rocket.evaluate_thrust_to_weight()
    mock_rocket.evaluate_static_margin()

    # Simulate Flight until Apogee
    test_flight = Flight(
        rocket=mock_rocket,
        environment=env,
        rail_length=5.2,
        inclination=85,
        heading=0,
        terminate_on_apogee=True,
        # equations_of_motion="solid_propulsion",
    )
    return test_flight.out_of_rail_velocity


speedbymass = Function(speed, inputs="Mass (kg)", outputs="Out of Rail Speed (m/s)")
speedbymass.plot(8, 20, 5)