I have this code for Octave and I cannot figure out how to get it to run and giv
ID: 3834641 • Letter: I
Question
I have this code for Octave and I cannot figure out how to get it to run and give me the plots I need. It keeps telling me division by zero inside the for loop for the F_G2 line. PLEASE HELP!
clc
clear
T_init = 39; % initial temperature
T_init = (T_init-32)*5/9; % initial temperature in celsius
p_init = 10332.29101; % kg per square meter
ts=1;
t=2:ts:21600;
A = 6.835; % square meters
G = 6.67408*10^-11; % gravity constant
h = zeros(1,305); % height
acc = zeros(1305); % acceleration
v = zeros(1305); % velocity
h(1) = 6371090; % meters of launch pad from center of earth
%h(t-1) = 6371090;
m_1 = 5.972*10^24; % mass of earth (kg)
% stage I mass in kg
m_boosters_drymass = 3.8*10^3; % mass of bossters in kg
m_boosters_propel = 39.6*10^3; % mass of booster propellant
% stage II mass in kg
m_core_drymass = 6.55*10^3;
m_core_propel = 92.95*10^3;
% stage III mass in kg
m3_drymass = 2.41*10^3;
m3_propel = 22.79*10^3;
%% soyuz(spacecraft) mass in kg
m_soyuz_drymass = 6.965*10^3;
m_soyuz_propel = .99*10^3;
%% force acting on soyuz at any time = gravitational force,air drag, rocket thrust
%% F_total = rocket thrust-gravitational force - atmospheric pressure force
%% F_airdrag = r*A*v^2
F_thrust_boosters = 838500 + 792480; % newtons
% 423*10*10^3; % in newtons
F_thrust_core = 792480;
% 79*10*10^3;
F_thrust_m = 297930;
% 3.92*10^3;
% h(1) = 6371090; % meters of launch pad from center of earth
m_1 = 5.972*10^24; % mass of earth (kg)
%% force acting on soyuz at any time = gravitational force,air drag, rocket thrust
%% F_total = rocket thrust-gravitational force - atmospheric pressure force
%% F_airdrag = r*A*v^2
F_thrust_boosters = 838500 + 792480; % newtons
% 423*10*10^3; % in newtons
F_thrust_core = 792480;
% 79*10*10^3;
F_thrust_m = 297930;
% 3.92*10^3;
% h(1) = 6371090; % meters of launch pad from center of earth
m_1 = 5.972*10^24; % mass of earth (kg)
% acceleration(k)= -k_s*pos(k-1);%+delta(k-1);
% velocity(k)= velocity(k-1)+TS*acceleration(k-1);
% pos(k)=pos(k-1)+TS*velocity(k-1);
for t = 2:1:305 % 5.08 minutes %% should maybe be t
if t<=118 % stage I(boosters and core)
m_1 = 4*(m_boosters_drymass + m_boosters_propel) + m_core_drymass + m_core_propel + m3_drymass + m3_propel + m_soyuz_drymass + m_soyuz_propel; % total mass
propel_dec = ((4*m_boosters_propel) + m_core_propel)/118; % decay of propellant mass
m_1new = m_1-t*propel_dec; % current total mass
F_G = G*m_1*m_1new/(h(t-1)^2); % current force due to gravity
if (h(t-1)-6371090)<= 11000 % troposphere
T_tropo = 15.04 - .00649 * (h(1)-6371090); % current temperature of troposphere
p_tropo = 101.29 * [(T_tropo + 273.1)/288.08]^5.256; % current pressure of troposphere
r_tropo = p_tropo / [.2869 * (T_tropo + 273.1)]; % current densityof troposphere
F_airdrag = r_tropo*A*v(t-1)^2; % current air reistnace of troposphere
F_current= F_thrust_boosters - F_G - F_airdrag; % current forces
acc(t) = (6*F_current/(2*m_1new));
v(t) = v(t-1) + (acc(t)-acc(t-1))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
if (h(t-1)-6371090)>11000 && (h(t-1)-6371090)<=25000 % lower stratosphere
T_lowstrato = -56.46;
p_lowstrato = 22.65 * exp(1.73 - .000157 * (h(1)-6371090));
r_lowstrato = p_lowstrato / [.2869 * (T_lowstrato + 273.1)];
F_airdrag2 = r_lowstrato*A*v(t-1)^2;
F_current2= F_thrust_boosters - F_G - F_airdrag;
acc(t) = (F_current2/m_2new);
v(t) = v(t-1) + (acc(t)-acc(t-1))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
end
if t>118 && t<=160 % stage II(core only)
m_2 = m_core_drymass + m_core_propel + m3_drymass + m3_propel + m_soyuz_drymass + m_soyuz_propel;
propel_dec2 = m_core_propel/(160-118);
m_2new = m_2-(t-117)*propel_dec2;
F_G2 = G*m_2*m_2new/(h(t-1)^2);
if (h(t-1)-6371090)>11000 && (h(t-1)-6371090)<=25000 % lower stratosphere
T_lowstrato = -56.46;
p_lowstrato = 22.65 * exp(1.73 - .000157 * (h(1)-6371090));
r_lowstrato = p_lowstrato / [.2869 * (T_lowstrato + 273.1)];
F_airdrag3 = (10e-6)*r_lowstrato*A*v(t-1)^2;
F_current3= F_thrust_core - F_G2 - F_airdrag3;
acc(t) = F_current3*10^3/m_2new;
v(t) = v(t-1) + (acc(t)-acc(t-1))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
if (h(t-1)-6371090)>25000 % upper stratosphere
T_upstrato = -131.21 + .00299 * (h(1)-6371090);
p_upstrato = 2.488 * [(T_upstrato + 273.1)/ 216.6]^-11.388;
r_upstrato = p_upstrato / [.2869 * (T_upstrato + 273.1)];
F_airdrag4 = (10e-6)*r_upstrato*A*v(t-1)^2;
F_current4= F_thrust_core - F_G2 - F_airdrag4;
acc(t) = F_current4*10^3/m_2new;
v(t) = v(t-1) + ((acc(t)-acc(t-1)))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
end
if t>160 && t<=295 % stage III - all upper stratosphere
m_3 = m3_drymass + m3_propel + m_soyuz_drymass + m_soyuz_propel;
propel_dec3 = m3_propel/(160-118);
m_3new = m_3-(t-117)*propel_dec3;
F_G3 = G*m_3*m_3new/(h(t-1)^2);
T_upstrato = -131.21 + .00299 * (h(1)-6371090);
p_upstrato = 2.488 * [(T_upstrato + 273.1)/ 216.6]^-11.388;
r_upstrato = p_upstrato / [.2869 * (T_upstrato + 273.1)];
F_airdrag5 = r_upstrato*A*v(t-1)^2;
F_current5= F_thrust_m - F_G3 - F_airdrag5;
acc(t) = abs(F_current5/m_3new);
v(t) = v(t-1) + abs((acc(t)-acc(t-1)))/2;
h(t) = h(t-1)+ abs((v(t)-v(t-1)))/2;
end
if t>295 % spacecraft - all upper stratosphere
m_4 = m_soyuz_drymass + m_soyuz_propel;
propel_dec4 = m_soyuz_propel/(160-118);
m_4new = m_4-(t-117)*propel_dec4;
F_G4 = G*m_4*m_4new/(h(t-1)^2);
T_upstrato = -131.21 + .00299 * (h(1)-6371090);
p_upstrato = 2.488 * [(T_upstrato + 273.1)/ 216.6]^-11.388;
r_upstrato = p_upstrato / [.2869 * (T_upstrato + 273.1)];
F_airdrag6 = r_upstrato*A*v(t-1)^2;
F_current6= F_thrust_m - F_G4 - F_airdrag6;
acc(t) = F_current6/m_4new;
v(t) = v(t-1) + abs((acc(t)-acc(t-1)))/2;
h(t) = h(t-1)+ abs((v(t)-v(t-1)))/2;
end
end
figure(1)
plot(t,h)
figure(2)
plot(t,v)
figure(3)
plot(t,acc)
Explanation / Answer
clc
clear
T_init = 39; % initial temperature
T_init = (T_init-32)*5/9; % initial temperature in celsius
p_init = 10332.29101; % kg per square meter
ts=1;
t=2:ts:21600;
A = 6.835; % square meters
G = 6.67408*10^-11; % gravity constant
h = zeros(1,305); % height
acc = zeros(1305); % acceleration
v = zeros(1305); % velocity
h(1) = 6371090; % meters of launch pad from center of earth
%h(t-1) = 6371090;
m_1 = 5.972*10^24; % mass of earth (kg)
% stage I mass in kg
m_boosters_drymass = 3.8*10^3; % mass of bossters in kg
m_boosters_propel = 39.6*10^3; % mass of booster propellant
% stage II mass in kg
m_core_drymass = 6.55*10^3;
m_core_propel = 92.95*10^3;
% stage III mass in kg
m3_drymass = 2.41*10^3;
m3_propel = 22.79*10^3;
%% soyuz(spacecraft) mass in kg
m_soyuz_drymass = 6.965*10^3;
m_soyuz_propel = .99*10^3;
%% force acting on soyuz at any time = gravitational force,air drag, rocket thrust
%% F_total = rocket thrust-gravitational force - atmospheric pressure force
%% F_airdrag = r*A*v^2
F_thrust_boosters = 838500 + 792480; % newtons
% 423*10*10^3; % in newtons
F_thrust_core = 792480;
% 79*10*10^3;
F_thrust_m = 297930;
% 3.92*10^3;
% h(1) = 6371090; % meters of launch pad from center of earth
m_1 = 5.972*10^24; % mass of earth (kg)
%% force acting on soyuz at any time = gravitational force,air drag, rocket thrust
%% F_total = rocket thrust-gravitational force - atmospheric pressure force
%% F_airdrag = r*A*v^2
F_thrust_boosters = 838500 + 792480; % newtons
% 423*10*10^3; % in newtons
F_thrust_core = 792480;
% 79*10*10^3;
F_thrust_m = 297930;
% 3.92*10^3;
% h(1) = 6371090; % meters of launch pad from center of earth
m_1 = 5.972*10^24; % mass of earth (kg)
% acceleration(k)= -k_s*pos(k-1);%+delta(k-1);
% velocity(k)= velocity(k-1)+TS*acceleration(k-1);
% pos(k)=pos(k-1)+TS*velocity(k-1);
for t = 2:1:305 % 5.08 minutes %% should maybe be t
if t<=118 % stage I(boosters and core)
m_1 = 4*(m_boosters_drymass + m_boosters_propel) + m_core_drymass + m_core_propel + m3_drymass + m3_propel + m_soyuz_drymass + m_soyuz_propel; % total mass
propel_dec = ((4*m_boosters_propel) + m_core_propel)/118; % decay of propellant mass
m_1new = m_1-t*propel_dec; % current total mass
F_G = G*m_1*m_1new/(h(t-1)^2); % current force due to gravity
if (h(t-1)-6371090)<= 11000 % troposphere
T_tropo = 15.04 - .00649 * (h(1)-6371090); % current temperature of troposphere
p_tropo = 101.29 * [(T_tropo + 273.1)/288.08]^5.256; % current pressure of troposphere
r_tropo = p_tropo / [.2869 * (T_tropo + 273.1)]; % current densityof troposphere
F_airdrag = r_tropo*A*v(t-1)^2; % current air reistnace of troposphere
F_current= F_thrust_boosters - F_G - F_airdrag; % current forces
acc(t) = (6*F_current/(2*m_1new));
v(t) = v(t-1) + (acc(t)-acc(t-1))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
if (h(t-1)-6371090)>11000 && (h(t-1)-6371090)<=25000 % lower stratosphere
T_lowstrato = -56.46;
p_lowstrato = 22.65 * exp(1.73 - .000157 * (h(1)-6371090));
r_lowstrato = p_lowstrato / [.2869 * (T_lowstrato + 273.1)];
F_airdrag2 = r_lowstrato*A*v(t-1)^2;
F_current2= F_thrust_boosters - F_G - F_airdrag;
acc(t) = (F_current2/m_2new);
v(t) = v(t-1) + (acc(t)-acc(t-1))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
end
if t>118 && t<=160 % stage II(core only)
m_2 = m_core_drymass + m_core_propel + m3_drymass + m3_propel + m_soyuz_drymass + m_soyuz_propel;
propel_dec2 = m_core_propel/(160-118);
m_2new = m_2-(t-117)*propel_dec2;
F_G2 = G*m_2*m_2new/(h(t-1)^2);
if (h(t-1)-6371090)>11000 && (h(t-1)-6371090)<=25000 % lower stratosphere
T_lowstrato = -56.46;
p_lowstrato = 22.65 * exp(1.73 - .000157 * (h(1)-6371090));
r_lowstrato = p_lowstrato / [.2869 * (T_lowstrato + 273.1)];
F_airdrag3 = (10e-6)*r_lowstrato*A*v(t-1)^2;
F_current3= F_thrust_core - F_G2 - F_airdrag3;
acc(t) = F_current3*10^3/m_2new;
v(t) = v(t-1) + (acc(t)-acc(t-1))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
if (h(t-1)-6371090)>25000 % upper stratosphere
T_upstrato = -131.21 + .00299 * (h(1)-6371090);
p_upstrato = 2.488 * [(T_upstrato + 273.1)/ 216.6]^-11.388;
r_upstrato = p_upstrato / [.2869 * (T_upstrato + 273.1)];
F_airdrag4 = (10e-6)*r_upstrato*A*v(t-1)^2;
F_current4= F_thrust_core - F_G2 - F_airdrag4;
acc(t) = F_current4*10^3/m_2new;
v(t) = v(t-1) + ((acc(t)-acc(t-1)))/2;
h(t) = h(t-1)+ (v(t)-v(t-1))/2;
end
end
if t>160 && t<=295 % stage III - all upper stratosphere
m_3 = m3_drymass + m3_propel + m_soyuz_drymass + m_soyuz_propel;
propel_dec3 = m3_propel/(160-118);
m_3new = m_3-(t-117)*propel_dec3;
F_G3 = G*m_3*m_3new/(h(t-1)^2);
T_upstrato = -131.21 + .00299 * (h(1)-6371090);
p_upstrato = 2.488 * [(T_upstrato + 273.1)/ 216.6]^-11.388;
r_upstrato = p_upstrato / [.2869 * (T_upstrato + 273.1)];
F_airdrag5 = r_upstrato*A*v(t-1)^2;
F_current5= F_thrust_m - F_G3 - F_airdrag5;
acc(t) = abs(F_current5/m_3new);
v(t) = v(t-1) + abs((acc(t)-acc(t-1)))/2;
h(t) = h(t-1)+ abs((v(t)-v(t-1)))/2;
end
if t>295 % spacecraft - all upper stratosphere
m_4 = m_soyuz_drymass + m_soyuz_propel;
propel_dec4 = m_soyuz_propel/(160-118);
m_4new = m_4-(t-117)*propel_dec4;
F_G4 = G*m_4*m_4new/(h(t-1)^2);
T_upstrato = -131.21 + .00299 * (h(1)-6371090);
p_upstrato = 2.488 * [(T_upstrato + 273.1)/ 216.6]^-11.388;
r_upstrato = p_upstrato / [.2869 * (T_upstrato + 273.1)];
F_airdrag6 = r_upstrato*A*v(t-1)^2;
F_current6= F_thrust_m - F_G4 - F_airdrag6;
acc(t) = F_current6/m_4new;
v(t) = v(t-1) + abs((acc(t)-acc(t-1)))/2;
h(t) = h(t-1)+ abs((v(t)-v(t-1)))/2;
end
end
figure(1)
plot(t,h)
figure(2)
plot(t,v)
figure(3)
plot(t,acc)