How do you appropriately integrate time steps to a while loop program for a single stage model rocket?

When you run the code and open the ModelRocket array you immediately notice the values don't change after the first pass. It has only iterated for one time step, when I need it for the full 153s duration. I'm assuming its due to an error with my while loop and usage of t with dt. If this is the case please share your thoughts on how to adjust my mistakes.

Equationsdata = fopen('ModelRocket.dat','w');
dt = 1;
while dt>=1
    
    % Earth and Launch Site Data
    g0 = 9.81;
    Re = 6378000;
    hs = 6000;
    
    %LV Configuration 
    CD = 0.2;
    Sref = 451.229;
    T = 3.3e+07; 
    m_dot = 250;
    ff = T/(g0*m_dot);
    L = 110.98;
    ve = 2.58e3;
    
    % initial parameters 
    m0 = 2951000;
    pitch0 = 90;
    pitch_rate = 0;
    v0 = 0;
    x0 = 0;
    x_dot = 0;
    h_dot = 0;
    h0 = 0;
    D = 0;
    t0 = 0;
    t = t0;
    rho0 = 1.225;
    a = T/m0-g0;
    fprintf(Equationsdata,'%f %f %f %f %f %f\n',t,h0,pitch0,v0,a,m0);
    t = t+dt;
    
    while t<153 % Stage 1
    m = m0 - ff*dt;
    v = v0 + a*dt;
    h_dot = v*sind(pitch0);
    h = h0 + h_dot*dt;
    rho = rho0*exp(-h/hs);
    D = CD*Sref*rho.*0.5*v^2;
    g_local = g0/(1+h/Re).^2;
    Pitch_rate = -(g_local/v)+(v/Re+h0)*cosd(pitch0);
    pitch = pitch0 + Pitch_rate*dt;
    x_dot = v*cosd(pitch0)*(Re/(Re+h));
    x = x0 + x_dot*dt;
    a = ((T-D)/m)-g_local*sind(pitch);
    fprintf(Equationsdata,'%f %f %f %f %f %f\n',t,h,pitch,v,a,m);
    t = t+dt;
    end
    dt = dt*.1; 
end
fclose(Equationsdata);
load ModelRocket.dat

Answer 1

[Sam Chak]

Your code runs, but I suspect that something is wrong in your Rocket Model.

Does the line return the pitch angle in degree or radian?

pitch = pitch0 + Pitch_rate*dt

Your pitch0 is described in degree.

It probably the best to use ode45() command to integrate all the governing differential equations and put all angles in radian unit. Can you list out the ODEs? I see mdot, xdot, and hdot. There should be a pitchrate_dot too.

dt = 1;
while dt >= 1
    
    % Earth and Launch Site Data
    g0 = 9.81;
    Re = 6378000;
    hs = 6000;
    
    %LV Configuration 
    CD    = 0.2;
    Sref  = 451.229;
    T     = 3.3e+07; 
    m_dot = 250;
    ff    = T/(g0*m_dot);
    L     = 110.98;
    ve    = 2.58e3;
    
    % initial parameters 
    m0         = 2951000;
    pitch0     = pi/2;  % 90 deg
    pitch_rate = 0;
    v0         = 0;
    x0         = 0;
    x_dot      = 0;
    h_dot      = 0;
    h0         = 0;
    D          = 0;
    t0         = 0;
    t          = t0;
    rho0       = 1.225;
    a          = T/m0 - g0;
    % fprintf(Equationsdata,'%f %f %f %f %f %f\n',t,h0,pitch0,v0,a,m0);
    t          = t + dt;
    
    while t < 153 % Stage 1
    m          = m0 - ff*dt;
    v          = v0 + a*dt;
    h_dot      = v*sin(pitch0);
    h          = h0 + h_dot*dt;
    rho        = rho0*exp(- h/hs);
    D          = CD*Sref*rho.*0.5*v^2;
    g_local    = g0/(1 + h/Re).^2;
    Pitch_rate = -(g_local/v) + (v/Re + h0)*cos(pitch0);
    pitch      = pitch0 + Pitch_rate*dt;
    x_dot      = v*cos(pitch0)*(Re/(Re + h));
    x          = x0 + x_dot*dt;
    a          = ((T - D)/m) - g_local*sin(pitch);
    % fprintf(Equationsdata,'%f %f %f %f %f %f\n',t,h,pitch,v,a,m);
    
    % just added this line to see what happen to h
    plot(t, h, 'b.'), hold on
    
    t          = t + dt;
    end
    
    dt = dt*.1;
    
end