Plotting range for a model rocket

Im trying to plot the range (x-axis) and altitude (y-axis) of a model rocket at a 5 degree launch angle but cannot get a correct plot. So far ive gto the correct trajectory for the burn ohase and the coast phase of the rocket, however it looks like the rocket drifts quite significantly (~7000 m) during the descent phase, and i cannot figure out why that is. It should drift only slightly at a 5 degree launch angle (roughly a couple meters or so), so im not sure why this is the case. Any help is appreciated. Thanks!

clear,clc;

launch_angle = 5;
dt = 0.001;
t_delay = 6; %delay time [s]

d_v = 0.0532; % vehicle diameter [m]
A_v = (pi/4)*(d_v).^2; % vehicle area [m^2]
d_p = 0.771144; % parachute diameter [m]
A_p = (pi/4)*(d_p)^2; % parachute area [m^2]

CD_V = 0.72; % vehicle drag coefficient (Will confirm after CFD)
CD_P = 1.5; % parachute drag coefficient
g = 9.81; %[m/s^2]
density = 1.225; %[kg/m^3]
I_tot = 62.2; % total impulse [Ns]
m_prop = 0.0431; %propellant mass [kg]
w_prop = m_prop*g; % propellant weight [N]
Isp_avg = I_tot/w_prop; % average specific impulse [s]

for j = 1

    time = 0;
    mass = 0.456; % initial mass [kg]
    if thrust(time) > 0
        t_burn = 2.3; % burn time [s] *provided by manufacturer*
    end
    V_0 = 0; %velocity
    Y_0 = 0; %altitude
    Vx_0 = 0;
    X_0 =0;

    i = 1;
   
    cont = true;
    while cont
        if time == 0
            Y(i) = Y_0;
            V(i) = V_0;
            Vx(i) = Vx_0;
            X(i) = X_0;
        end

        % burn phase
        if thrust(time)>0
            mdot = thrust(time)/(9.81*Isp_avg); %propellant mass flow rate [kg/s]
            mass = 0.456 - mdot*dt; %updating mass
            D = CD_V(j)*0.5*density*(V(i)^2)*A_v; % drag [N]
            dv = (thrust(time)/mass)*dt - (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
            dvx = 9.81*sind(launch_angle)*dt;
            V(i + 1) = V(i) + dv; 
            Vx(i + 1) = Vx(i) + dvx;
            Y(i + 1) = Y(i) + V(i)*dt;
            X(i + 1) = X(i) + Vx(i)*dt;
        end

        % coast phase
        if thrust(time)<=0 && (time < t_delay+t_burn)
            D = CD_V(j)*0.5*density*(V(i)^2)*A_v;
            if(V(i) > 0 )
                D = -D;
            end
            mass = 0.456-m_prop; % initial mass - prop mass [kg]
            dv = (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
            dvx = 9.81*sind(launch_angle)*dt;
            V(i + 1) = V(i) + dv; 
            Vx(i + 1) = Vx(i) + dvx;
            Y(i + 1) = Y(i) + V(i)*dt;
            X(i + 1) = X(i) + Vx(i)*dt;
        end

        % descent phase
        if (time > t_delay+t_burn)
            D = CD_P(j)*0.5*density*(V(i)^2)*A_p;
            if(V(i) > 0 )
                D = -D;
            end
            dv = (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
            dvx = 9.81*sind(launch_angle)*dt;
            V(i + 1) = V(i) + dv; 
            Vx(i + 1) = Vx(i) + dvx;
            Y(i + 1) = Y(i) + V(i + 1)*dt;
            X(i + 1) = X(i) + Vx(i + 1)*dt;
        end

        if(Y(i)<0)
            cont = false;
        end

        i = i + 1;
        time = time + dt;

    end
    t = 0:dt:(i - 1)*dt;
end

    figure(1)
    sgtitle ('Simulated Rocket Range')
    set(gcf,'position',[500 200 1000 600])
    plot(X,Y)
    %xlim([-50 50])
    set(gca,'XMinorTick','on','YMinorTick','on')
    grid on,xlabel('Range [m]')
    ylabel('Altitude [m]')

% Data on AeroTech F26-6FJ
function T = thrust(time)

   temp = [ 0.041	38.289
0.114	36.318
0.293	34.347
0.497	32.939
0.774	32.376
1	31.25
1.254	28.716
1.498	25.338
1.743	22.241
2.003	17.737
2.077	15.484
2.304	5.349
2.484	1.689
2.61	0];

T = interp1(temp(:,1),temp(:,2),time,'linear','extrap'); 

end

Answer 1

[Pratik]

Hi Garrett,

As per my understanding, you want to plot the range vs altitude plot of a model rocket. However, there is a huge drift in the range of the model during the descent phase of the rocket.

Upon reviewing the code, I noticed that there is an issue with how the force is being applied to the model rocket. To ensure the range of the rocket is accurate, you may revisit the following points.

1. Gravity should only influence the vertical component of the rocket's velocity.
2. The thrust should impact both the vertical and horizontal components of the rocket's velocity.
3. During the descent phase, the horizontal component of the velocity should remain constant as there are no forces acting in that direction.

Please consider the modifications to the 'dv' and 'dvx' variables for different phases, in the code snippet below.

% burn phase
dv = (thrust(time)/mass)*cosd(launch_angle)*dt - (D/mass)*dt - 9.81*dt;
dvx = (thrust(time)/mass)*sind(launch_angle)*dt;

% coast phase
dv = (D/mass)*dt - 9.81*dt;
dvx = 0;

% descent phase
dv = (D/mass)*dt - 9.81*dt;
dvx = 0;

Hope this helps!