In Projects

MATLAB Rocket Trajectory Simulation

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

Results are shown for a multi-stage rocket with variable mass, thrust, and drag.

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


\begin{equation} c_d = {F_d \over {{1 \over 2} \rho v^2 A}} \end{equation}
\begin{equation} F_d = {{1 \over 2} \rho v^2 c_d A} \end{equation}
  • #
    Figure 1
    Rocket Trajectory Results

    Plots of rocket trajectory simulation results.

MATLAB Code Download

Download the MATLAB program.

MATLAB 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'});