# ENGINEERING

Home Engineering
• Engineering
• Senior Project
• Projects
• 3 Doors Game Simulation
• Woodworking: Bedstand
• Woodworking: Desk
• Numbers Game
• Path Planning
• Snake Video Game AI
• Minimum Spanning Tree
• Heat Map
• MATLAB: Rocket Simulation
• MATLAB: Bike Suspension
• MATLAB: Thermal FEA SS
• MATLAB: Thermal FEA Tr
• Notes
• Statistics
• Steoscopic Estimates
• Cubed Roots
• Other
• Free Fall Distance
• Speed of Sound
• Travel
• Solar System Scale
• Videos
• Yellowstone Geysers
• Geysers and Waterfalls
• Yellowstone and Glacier
• Grinnell Glacier
• Yellowstone Geysers
• Photography

## MATLAB Rocket Trajectory Simulation

This program simulates the trajectory of a rocket given initial conditions and phyiscal parameters. Numerical analysis is used to simulate the flight of the rocket, including effects due to thrust, drag, change in mass, and gravity.

The drag coefficient, Equation 1, is used as input to the dynamics model. The drag force of the rocket is calculated using Equation 2.

$$c_d = {F_d \over {{1 \over 2} \rho v^2 A}}$$
$$F_d = {{1 \over 2} \rho v^2 c_d A}$$

## MATLAB Source Code:

% Program:
% Rocket_Trajectory_Simulation.m
% Multi-stage rocket dynamics and trajectory simulation.
%
% Description:
% Predicts multi-stage rocket dynamics and trajectories based on the given
% rocket mass, engine thrust, launch parameters, and drag coefficient.
%
% Variable List:
% Delta = Time step (s)
% t = Time (s)
% Thrust = Thrust (N)
% Mass = Mass (kg)
% Mass_Rocket_With_Motor = Mass with motor (kg)
% Mass_Rocket_Without_Motor = Mass without motor (kg)
% Theta = Angle (deg)
% C = Drag coefficient
% Rho = Air density (kg/m^3)
% A = Rocket projected area (m^2)
% Gravity = Gravity (m/s^2)
% Launch_Rod_Length = Length of launch rod (m)
% n = Counter
% Fn = Normal force (N)
% Drag = Drag force (N)
% Fx = Sum of forces in the horizontal direction (N)
% Fy = Sum of forces in the vertical direction (N)
% Vx = Velocity in the horizontal direction (m/s)
% Vy = Velocity in the vertical direction (m/s)
% Ax = Acceleration in the horizontal direction (m/s^2)
% Ay = Acceleration in the vertical direction (m/s^2)
% x = Horizontal position (m)
% y = Vertical position (m)
% Distance_x = Horizontal distance travelled (m)
% Distance_y = Vertical travelled (m)
% Distance = Total distance travelled (m)
% Memory_Allocation = Maximum number of time steps expected

clear, clc      % Clear command window and workspace

% Parameters
Delta = 0.001;                  % Time step
Memory_Allocation = 30000;      % Maximum number of time steps expected

% Preallocate memory for arrays
t = zeros(1, Memory_Allocation);
Thrust = zeros(1, Memory_Allocation);
Mass = zeros(1, Memory_Allocation);
Theta = zeros(1, Memory_Allocation);
Fn = zeros(1, Memory_Allocation);
Drag = zeros(1, Memory_Allocation);
Fx = zeros(1, Memory_Allocation);
Fy = zeros(1, Memory_Allocation);
Ax = zeros(1, Memory_Allocation);
Ay = zeros(1, Memory_Allocation);
Vx = zeros(1, Memory_Allocation);
Vy = zeros(1, Memory_Allocation);
x = zeros(1, Memory_Allocation);
y = zeros(1, Memory_Allocation);
Distance_x = zeros(1, Memory_Allocation);
Distance_y = zeros(1, Memory_Allocation);
Distance = zeros(1, Memory_Allocation);

C = 0.4;                                % Drag coefficient
Rho = 1.2;                              % Air density (kg/m^3)
A = 4.9*10^-4;                          % Rocket projected area (m^2)
Gravity = 9.81;                         % Gravity (m/s^2)
Launch_Rod_Length = 1;                  % Length of launch rod (m)
Mass_Rocket_With_Motor = 0.01546;       % Mass with motor (kg)
Mass_Rocket_Without_Motor = 0.0117;     % Mass without motor (kg)

Theta(1) = 60;                  % Initial angle (deg)
Vx(1) = 0;                      % Initial horizontal speed (m/s)
Vy(1) = 0;                      % Initial vertical speed (m/s)
x(1) = 0;                       % Initial horizontal position (m)
y(1) = 0.1;                     % Initial vertical position (m)
Distance_x(1) = 0;              % Initial horizontal distance travelled (m)
Distance_y(1) = 0;              % Initial vertical distance travelled (m)
Distance(1) = 0;                % Initial  distance travelled (m)
Mass(1) = Mass_Rocket_With_Motor;       % Initial rocket mass (kg)

n = 1;                          % Initial time step

while y(n) > 0                  % Run until rocket hits the ground
n = n+1;                    % Increment time step

t(n)= (n-1)*Delta;          % Elapsed time

% Determine rocket thrust and mass based on launch phase
if t(n) <= 0.1                              % Launch phase 1
Thrust(n) = 56*t(n);
Mass(n) = Mass_Rocket_With_Motor;
elseif t(n) > 0.1 && t(n) < 0.5            % Launch phase 2
Thrust(n) = 5.6;
Mass(n) = Mass_Rocket_With_Motor;
elseif t(n) >= 0.5 && t(n) < 3.5            % Launch phase 3
Thrust(n) = 0;
Mass(n) = Mass_Rocket_With_Motor;
elseif t(n) >= 3.5                          % Launch phase 4
Thrust(n) = 0;
Mass(n) = Mass_Rocket_Without_Motor;    % Rocket motor ejects
end

% Normal force calculations
if Distance(n-1) <= Launch_Rod_Length       % Launch rod normal force
Fn(n) = Mass(n)*Gravity*cosd(Theta(1));
else
Fn(n) = 0;                              % No longer on launch rod
end

% Drag force calculation
Drag(n)= 0.5*C*Rho*A*(Vx(n-1)^2+Vy(n-1)^2); % Calculate drag force

% Sum of forces calculations
Fx(n)= Thrust(n)*cosd(Theta(n-1))-Drag(n)*cosd(Theta(n-1))...
-Fn(n)*sind(Theta(n-1));                            % Sum x forces
Fy(n)= Thrust(n)*sind(Theta(n-1))-(Mass(n)*Gravity)-...
Drag(n)*sind(Theta(n-1))+Fn(n)*cosd(Theta(n-1));    % Sum y forces

% Acceleration calculations
Ax(n)= Fx(n)/Mass(n);                       % Net accel in x direction
Ay(n)= Fy(n)/Mass(n);                       % Net accel in y direction

% Velocity calculations
Vx(n)= Vx(n-1)+Ax(n)*Delta;                 % Velocity in x direction
Vy(n)= Vy(n-1)+Ay(n)*Delta;                 % Velocity in y direction

% Position calculations
x(n)= x(n-1)+Vx(n)*Delta;                   % Position in x direction
y(n)= y(n-1)+Vy(n)*Delta;                   % Position in y direction

% Distance calculations
Distance_x(n) = Distance_x(n-1)+abs(Vx(n)*Delta);      % Distance in x
Distance_y(n) = Distance_y(n-1)+abs(Vy(n)*Delta);      % Distance in y
Distance(n) = (Distance_x(n)^2+Distance_y(n)^2)^(1/2); % Total distance

% Rocket angle calculation
Theta(n)= atand(Vy(n)/Vx(n));      % Angle defined by velocity vector
end

figure('units','normalized','outerposition',[0 0 1 1]) % Maximize plot window

% Figure 1
subplot(3,3,1)
plot(x(1:n),y(1:n));
xlim([0 inf]);
ylim([0 inf]);
xlabel({'Range (m)'});
ylabel({'Altitude (m)'});
title({'Trajectory'});

% Figure 2
subplot(3,3,2)
plot(t(1:n),Vx(1:n));
xlabel({'Time (s)'});
ylabel({'Vx (m/s)'});
title({'Vertical Velocity'});

% Figure 3
subplot(3,3,3)
plot(t(1:n),Vy(1:n));
xlabel({'Time (s)'});
ylabel({'Vy (m/s)'});
title({'Horizontal Velocity'});

% Figure 4
subplot(3,3,4)
plot(t(1:n),Theta(1:n));
xlabel({'Time (s)'});
ylabel({'Theta (Deg)'});
title({'Theta'});

% Figure 5
subplot(3,3,5)
plot(Distance(1:n),Theta(1:n));
xlim([0 2]);
ylim([59 61]);
xlabel({'Distance (m)'});
ylabel({'Theta (Deg)'});
title({'Theta at Launch'});

% Figure 6
subplot(3,3,6)
plot(t(1:n),Mass(1:n));
ylim([.0017 .02546]);
xlabel({'Time (s)'});
ylabel({'Mass (kg)'});
title({'Rocket Mass'});

% Figure 7
subplot(3,3,7)
plot(t(1:n),Thrust(1:n));
xlim([0 0.8]);
xlabel({'Time (s)'});
ylabel({'Thrust (N)'});
title({'Thrust'});

% Figure 8
subplot(3,3,8)
plot(t(1:n),Drag(1:n));
xlabel({'Time (s)'});
ylabel({'Drag (N)'});
title({'Drag Force'});

% Figure 9
subplot(3,3,9)
plot(Distance(1:n),Fn(1:n));
xlim([0 2]);
xlabel({'Distance (m)'});
ylabel({'Normal Force (N)'});
title({'Normal Force'});