Anyone knows why this keeps running forever

function TWOSTGRQ2
clc
clear variables
%%Variables
g0 = 9.81; % Gravity
m01 = 120000;% Lift-off mass
n = 15; % Mass Ratio
Isp1 = 320; Isp2 =350; % Specific Impulse (s)
c1 = Isp1*g0; c2 = Isp2*g0;%
T2W1 = 1.4; T2W2 = 1.5; % Thrust to weight ratio
p0 = 1.225; % Sea level density of atmosphere (kg/m^3)
CD = 0.5; % Drag Coefficient (assumed constant)
h0 = 7.5*10^3; % Density scale height (m)
hTower = 110; % Tower Height (m)
Y_ini = 89.85*pi/180; % Flight Path before gravity turn maneuovre
RE = 6378000;
p0 = 1.225;
D = 5; %Rocket Diameter
A = pi*(D/2)^2; %Rocket Cross Sectional Area
%%Calculated Parameters
mf1 = m01/n; m02 = mf1; % Mass of rocket after engine burn stage 1- m02
mf2 = m02/n; % Mass of rocket after engine burn stage 1 - m02
T1 = T2W1*m01*g0; T2 = T2W2*m02*g0; % Thrust
rme1 = T1/c1; rme2 = T2/c2;% Rate at which exhaust mass flows accross the nozzle exit plane
mp1 =m01-mf1; mp2 = m02-mf2; % Propellant mass
tburn1 = mp1/rme1; tburn2 = mp2/rme2; % Burnout times
t0 = 0; tf = tburn1+tburn2;
%%Differential Equations - ode45 - actual- 2b
init_params = [0 Y_ini 0 0 0 0];
options = odeset('RelTol', 1e-4);
[t,x] = ode45(@(t,x) rocket_analysis_actual(t,x,T1,T2,m01,m02,rme1,rme2,mp1,mp2,g0,RE,Y_ini,h0,CD,p0,A), [t0 tf],init_params,options);
%%Vectorization for outputs - Actual
v = x(:,1);
Y = x(:,2);
h = x(:,3);
x = x(:,4);
VelLostDrag = x(:,5);
VelLostGravity = x(:,6);
Y_a_deg = radtodeg(Y);
%%Dynamic Pressure & Karman Line
density = p0*exp(-(h)./h0);
dynamicpressure = 0.5*density.*v.^2;
j=100000; escapevel = 11200; timek = 167;  % J is karman line height
[dynamicpressuremax, Index_maxdynamicpressure] = max(double(dynamicpressure));
HeightMaxPos = h(Index_maxdynamicpressure);
fprintf("Max Dynamic Pressure is: %g\nHeight at Max Dynamic Pressure is: %g",dynamicpressuremax, HeightMaxPos);
%
%%Plots 2b and 2c
figure(1); plot(t,v,'r-'); hold on
escape = line([t0 timek],[escapevel escapevel]); tl = line([timek timek],[0 escapevel]); 
tl.Color = 'green'; escape.Color = 'magenta';
xlabel('Time (s)'); ylabel('Velocity (m/s)');
title('Graph showing Velocity against Time (Actual)')
legend('Velocity','Escape Velocity')
figure(2); plot(t,Y_a_deg,'m-')
xlabel ('Time (s)'); ylabel('Flight Path Angle (degrees)')
title('Graph showing Flight Path Angle against Time (Actual)')
figure(3); plot(t,h,'b-'); hold on
karman = line([t0 tburn],[j,j]); karman.Color = 'red';
legend('Altitude','Karman Line');
xlabel('Time (s)'); ylabel('Altitude (m)')
title('Graph showing Altitude against Time (Actual)')
figure(4); plot(t,x,'g-')
xlabel('Time (s)'); ylabel('Horizontal Travel (m)')
title('Graph showing Horizontal Travel against Time (Approximate)')
figure(5); plot(t,VelLostDrag,'b-')
xlabel('Time (s)'); ylabel('Velocity Lost Due to Drag (m/s)')
title('Graph showing Velocity Lost to due Drag against Time (Cumulative)')
figure(6); plot(t,VelLostGravity,'g-')
xlabel('Time (s)'); ylabel('Velocity Lost due to Gravity (m/s)')
title('Graph showing Velocity Lost due to Gravity against Time (Cumalitive)')
figure(7); plot(t,dynamicpressure,'r-'); hold on
dynpress = line([0 tburn],[dynamicpressuremax dynamicpressuremax]); dynpress.Color = 'green';
xlabel('Time (s)'); ylabel('Dynamic Pressure')
figure(8); plot(h,dynamicpressure,'r-'); hold on
dynpressalt = line([0 HeightMaxPos],[dynamicpressuremax dynamicpressuremax]); dynpressalt.Color = 'green';
dynpressaltvert = line([HeightMaxPos HeightMaxPos],[0 dynamicpressuremax]); dynpressaltvert.Color = 'green';
xlabel('Height (m)'); ylabel('Dynamic Pressure')
end
%%Rocket Analysis Function Before Tower Cleared
function f_dot = rocket_analysis_actual(t,x,T1,T2,m01,m02,rme1,rme2,mp1,mp2,g0,RE,Y_ini,h0,CD,p0,A)
tburn1 = mp1/rme1; tburn2 = mp2/rme2; % Burnout times
tf = tburn1+tburn2;
while t<=tburn1 
    m0 = m01;
    T = T1;
    rme = rme1;     
end
while t>tburn1 && t<=tf
    m0 = m02;
    T = T2;
    rme = rme2;
end
hTower = 110; % Tower Height (m)
gc = g0/(1+(x(3)/RE))^2; % Gravity change with altitude
m = m0-rme.*t;
p = p0*exp(-x(3)./h0); % Density change with height
q = 0.5*p.*x(1).^2;% Dynamic Pressure as it changes with density
drag = q*A*CD; % Drag
f_dot = [0 0 0 0 0 0]';
if x(3) < hTower
    f_dot(1) = (T-drag-(m*gc*sin(Y_ini)))/m; % Acceleration in vertical
    f_dot(2) = 0; % Flight Path angle with time - zero before tower cleared
    f_dot(3) = x(1)*sin(Y_ini); % Altitude change with time
    f_dot(4) = (RE/(RE+x(3)))*(x(1)*cos(Y_ini)); % Change in horizontal
    f_dot(5) = drag/m; % Velocity Lost due to drag
    f_dot(6) = gc*sin(Y_ini); % Velocity lost due to gravity
else
    f_dot(1) = (T-drag-m*gc*sin(x(2)))/m; % Acceleration in vertical
    f_dot(2) = (-1/x(1))*(gc-(((x(1))^2)/(RE+x(3))))*cos(x(2)); % Flight Path angle with time - zero before tower cleared
    f_dot(3) = x(1)*sin(x(2)); % Altitude change with time
    f_dot(4) = (RE/(RE+x(3)))*(x(1)*cos(x(2))); % Change in horizontal
    f_dot(5) = drag/m;
    f_dot(6) = gc*sin(x(2)); 
end
end

[per isakson]

Check these while-loops

    while t<=tburn1
        m0 = m01;
        T = T1;
        rme = rme1;
    end
    while t>tburn1 && t<=tf
        m0 = m02;
        T = T2;
        rme = rme2;
    end

[David Oshidero]

Thing is I know it started happening when I added that but i dont know how to get that bit to work or what can be used as a substitute?

[David Fletcher]

The variables in the conditions aren't being updated in the loop, so if the condition is true when it enters the loop, it will always be true

[per isakson]

"Thing is I know it started happening when I added" Why didn't you tell us that?

• See Debug a MATLAB Program
• See Examine Values While Debugging

[David Oshidero]

yeah my bad for that, I thought it was something else upper in the code conflicting with it, If i wanted to make it so that the variable changes after tburn1 ends how would i go about doing it?

[David Fletcher]

What are the loops actually doing? It just seems like assignments that will be the same on every iteration.

[David Oshidero]

the loops are just telling the functions below it, that if the time is less than tburn assign these variables to the main ones and if it goes past tburn1 assign the next set. Im trying to do an ode for a twostagerocket and draw the graphs but this was the only way i could think of doing it

[David Fletcher]

I think there is something fundamental that I'm missing - I just don't see why you'd use loops there rather than a decision block(s). Nothing in the loop varies from one iteration to the next.

[David Oshidero]

Ive never used that. What is a decision block(s)? Im a student and lecturers arent really doing too well teaching us matlab

[David Oshidero]

or do you mean if statements?

[David Fletcher]

Yes, if statements

[David Oshidero]

i tried if statements and i now get this error, while loops were replaced with if statements

Warning: Failure at t=2.333333e+02.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.547474e-13) at time t. 
> In ode45 (line 308)
  In TWOSTGRQ2 (line 31) 
Index exceeds matrix dimensions.
Error in TWOSTGRQ2 (line 37)
VelLostDrag = x(:,5);

Answer 1

[Abraham Boayue]

It would be much better to use a single for loop to implement this function. The two while loops and the if statement do not seem to be working together as they should.