MATLAB® Tutorial#

RocketPy in MATLAB®#

Did you know you can run RocketPy inside MATLAB®?

MATLAB® provides a nice method to run Python libraries, allowing everyone familiar with MATLAB® to run RocketPy, getting the best of both worlds.

Below, you will find a quick tutorial created using Live Script which reproduces RocketPy’s classic Getting Started tutorial.

Download Files#

You can download the Live Script or the MATLAB® script file to follow along.

Tutorial#

Getting Started with RocketPy in MATLAB®

In this Live Script, you will learn how to run RocketPy using the MATLAB® environment.

First things first: clone/download RocketPy's repository and set it as your MATLAB® working directory so that you can run this live script without any issues.

After that, we start by configuring our Python environment. You can do so by following the guidelines presented in the MATLAB® documentation: Configure Your System to Use Python - MATLAB & Simulink (mathworks.com).

Once the Python environment is configured, RocketPy needs to installed using pip as outlined in RocketPy's documentation: Installation — RocketPy documentation.

Finally, all the prerequisites are complete and you can comeback to MATLAB®! We just need to set the execution mode as out of process and start working. MATLAB® can run Python scripts and functions in a separate process. Running Python in a separate process enables you to:

Use some third-party libraries in the Python code that are not compatible with MATLAB®.
Isolate the MATLAB® process from crashes in the Python code.

pyenv('ExecutionMode','OutOfProcess');

Note: if MATLAB® is not able to find Python automatically, you may have to run the command line above including a path to the Python exectuable installed on your computer:

% pyenv('ExecutionMode','OutOfProcess', 'Version', '/path/to/python/executable');

Now, we will go through a simplified rocket trajectory simulation to get you started. Let's start by importing the rocketpy module and its classes.

rocketpy = py.importlib.import_module('rocketpy');
Environment = rocketpy.environment.Environment;
SolidMotor = rocketpy.motors.solid_motor.SolidMotor;
Rocket = rocketpy.rocket.Rocket;
Flight = rocketpy.simulation.flight.Flight;

Setting Up a Simulation

Creating an Environment for Spaceport America

% rocketpy.Environment(latitude, longitude, elevation);
env = Environment(pyargs(...
    '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.

Tomorrow = datetime('tomorrow');
env.set_date({int32(Tomorrow.Year), int32(Tomorrow.Month), int32(Tomorrow.Day), int32(12)}) % Hour given in UTC time (noon UTC)

Now, we tell our Environment object to retrieve a weather forecast for our specified location and date using GFS:

env.set_atmospheric_model(pyargs( ...
    'type', "Forecast", ...
    'file', "GFS" ...
))

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

env.info()

Plots will open in a separate window, so be sure to run this last cell to see them!

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(pyargs( ...
    'thrust_source', "../../data/motors/Cesaroni_M1670.eng", ...
    'dry_mass', 1.815, ...
    'dry_inertia', py.tuple({0.125, 0.125, 0.002}), ...
    'nozzle_radius', 33 / 1000, ...
    'grain_number', int32(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.85704, ...
    'center_of_dry_mass_position', 0.317, ...
    'nozzle_position', 0, ...
    'burn_time', 3.9, ...
    'throat_radius', 11 / 1000, ...
    'interpolation_method', "linear", ...
    'coordinate_system_orientation', "nozzle_to_combustion_chamber" ...
));

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

Plots will open in a separate window, so be sure to run this last cell to see them!

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, supplying it with the Pro75M1670 engine, entering its inertia properties, some dimensions and also its drag curves.

calisto = Rocket(pyargs( ...
    'radius', 127 / 2000, ...
    'mass', 14.426, ...
    'inertia', py.tuple({6.321, 6.321, 0.0334}), ...
    '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" ...
));
 
Rocket.set_rail_buttons(pyargs( ...
    'self', calisto, ...
    'upper_button_position', 0.0818, ...
    'lower_button_position', -0.618, ...
    'angular_position', 45 ...
));
Rocket.add_motor(pyargs( ...
    'self', calisto, ...
    'motor', Pro75M1670, ...
    'position', -1.255 ...
));
% Adding Aerodynamic Surfaces
% Now we define the aerodynamic surfaces. They are really straight forward.
 
nosecone = Rocket.add_nose(pyargs( ...
    'self', calisto, ...
    'length', 0.55829, ...
    'kind', "von karman", ...
    'position', 1.278 ...
));
 
fin_set = Rocket.add_trapezoidal_fins(pyargs( ...
    'self', calisto, ...
    'n', int32(4), ...
    'span', 0.110, ...
    'root_chord', 0.120, ...
    'tip_chord', 0.060, ...
    'position', -1.04956 ...
));
 
tail = Rocket.add_tail(pyargs( ...
    'self', calisto, ...
    'top_radius', 0.0635, ...
    'bottom_radius', 0.0435, ...
    'length', 0.060, ...
    'position', -1.194656 ...
));
% 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 
 
% py.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. Unfortunately, defining these with MATLAB® code is not yet possible.

% Drogue parachute is triggered when vertical velocity is negative, i.e. rocket is falling past apogee
drogue_trigger = py.eval("lambda p, y: y[5] < 0", py.dict);
 
% Main parachute is triggered when vertical velocity is negative and altitude is below 800 AGL 
main_trigger = py.eval("lambda p, y: (y[5] < 0) and (y[2] < 800 + 1400)", py.dict); 

Now we add both the drogue and the main parachute to our rocket.

Main = Rocket.add_parachute(pyargs( ...
    'self', calisto, ...
    'name', "Main", ...
    'cd_s', 10.0, ...
    'trigger', main_trigger, ...
    'sampling_rate', 105, ...
    'lag', 1.5, ...
    'noise', py.tuple({0, 8.3, 0.5}) ...
));
 
Drogue = Rocket.add_parachute(pyargs( ...
    'self', calisto, ...
    'name', "Drogue", ...
    'cd_s', 1.0, ...
    'trigger', drogue_trigger, ...
    'sampling_rate', 105, ...
    'lag', 1.5, ...
    'noise', py.tuple({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(pyargs(...
    '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()

Plots will open in a separate window, so be sure to run this last cell to see them!

Working with Data Generated by RocketPy in MATLAB®

You can access the entire trajectory solution matrix with the following line of code. The returned matrix contain the following columns: time t (s), x (m), y (m), z (m),

(m/s),

(rad/s),

(rad/s).

solution_matrix = double(py.numpy.array(test_flight.solution))

Support for accessing secondary values calculated during post processing, such as energy, mach number, and angle of attack, is also available for all versions of RocketPy greater than or equal to version 0.11.0.

To showcase this, let's get the angle of attack of the rocket and plot it using MATLAB®:

angle_of_attack = double(test_flight.angle_of_attack.source)
plot(angle_of_attack(:,1), angle_of_attack(:,2)) % First column is time (s), second is the angle of attack
ylabel('Angle of Attack (Deg)')
xlabel('Time(s)')
xlim([test_flight.out_of_rail_time, 10])

You can also convert data from other objects besides the Flight class, such as from the Environment. For example, let's say you want to get the wind velocity, both the x and y component:

wind_velocity_x = double(env.wind_velocity_x.source) % First column is altitude (ASL m), while the second one is the speed (m/s)
wind_velocity_y = double(env.wind_velocity_y.source) % First column is altitude (ASL m), while the second one is the speed (m/s)

Time to Fly!

This is all you need to get started using RocketPy in MATLAB®! Now it is time to play around and create amazing rockets. Have a great launch!