For Loop Graph Not Plotting

The code maps the change in flight angle as it rotates about the centre of mass between 0 and pi/2 as shown by the diagram:

The first part of the code finds how the centre of mass of the body changes as the fuel rod is used up, the second is a for loop that maps the change in theta. However the graphs do not generate any values as shown below:

Why does theta_double_dot (angular acceleration) stay zero?

clear

clc

%time

timestep = 1;

t_end = 1200;

t_start = 0;

t = t_start:timestep:t_end;

%mass

m_payload = 1; %chnage this value

m_fuel_initial = 100*300;

dmdt = 3;

m_fuel = m_fuel_initial - dmdt*t;

m_engine = 100;

m_to = m_payload + m_fuel_initial + m_engine;

m(1) = m_to;

g_0 = 9.81;

%m = m_payload + m_fuel + m_engine;

v_e = 3E3; %exhaust velocity

T_max = m_to*10;

dxdt = 1.5E-3; % change of fuel rod length m/s

l_payload = 2.5-1.5;

l_fuel_initial = 8;

l_fuel = l_fuel_initial - dxdt*t;

l_engine = 2.4;

l_rocket_initial = l_payload + l_fuel_initial + l_engine;

l_rocket = l_rocket_initial - dxdt*t;

x_payload = l_payload/2; %distance from top of rocket

m_payload = 1;

x_0 = l_fuel/2;

x_fuel = (l_rocket -x_0) - dxdt*t;

m_fuel = m_fuel_initial - dmdt*t;

x_engine_intial = l_rocket - 2.5/2;

x_engine = x_engine_intial - dxdt*t; %distance from top of rocket

m_x_payload = m_payload*x_payload;

m_x_fuel = m_fuel.*x_fuel;

m_x_engine = m_engine.*x_engine;

m_x = m_x_payload + m_x_fuel + m_x_engine;

moment(1) = m_x(1);

x_COM = l_rocket - (m_x./(m_payload + m_fuel + m_engine));

%COM = centre of mass taken from the bottom of the rocket

%absolutre distance of each componet from the COM

x_p = abs(x_payload - x_COM);

x_b = abs(x_fuel - x_COM);

x_e = abs(x_engine - x_COM);

diameter = 15;

radius = diameter/2;

%moment of inertia

I_p = m_payload*radius^2*83/320;

I_e = m_engine*l_engine^3;

I_fuel = (radius^2/2)*m_fuel;

I = I_p + I_fuel + I_e + m_payload*x_p.^2+ m_fuel.*x_b.^2 + m_engine*x_e.^2;

%sum of moments

m_x_I = g_0*m_x.*I.^-1;

m_x_I_initial = m_x_I(1);

thetastep = 0.5*pi/t_end ; %change in theta

theta(1) = thetastep;

theta_double_dot(1) = m_x_I(1)*sin(thetastep);

theta_dot(1) = theta_double_dot(1)*timestep;

%something is going wrong here

for i = 2:length(t)

if (theta <= 0) & (theta >= pi/2)

m(i) = m(i-1) - timestep*dmdt;

rho(i)= 0;

theta_double_dot(i) = moment(i)/I(i)*g*sin(theta(i-1));

theta_dot(i) = theta_dot(i-1) + timestep*thetastep(i-1);

theta(i) = theta_dot(i-1)+ timestep*thetadot(i-1);

end

end

subplot(3,1,1)

plot(t,theta_double_dot,'o')

xlabel('time s')

ylabel('theta rad/s^2')

subplot(3,1,2)

plot(t,theta_dot,'x')

xlabel('time s')

ylabel('theta rad/s')

subplot(3,1,3)

plot(t,theta,'o')

xlabel('time s')

ylabel('theta rad')

Answer 1

[Star Strider]

I could be missing something, however, how could this:

if (theta <= 0) & (theta >= pi/2)

ever be true?

It would likely work better if the test were on ‘theta(i-1)’ however that would still not resolve the logical conflict!

.

[Asit Rahman]

That works, sorry, I didn't realise my inequalities were impossible

[Star Strider]

No worries!

The inequalities likely need to be reversed, so the first is >= and the second is <=. Using | (or) instead of & would also work, if the intent was to exclude the region between 0 and.pi/2.