% M3eP182.m
% Program on page 182.
% Script file rocket1.m
% Computes flight duration as a function of burn time.
% Basic data values.
m_e = 100; q = 1; u = 8000; g = 32.2;
dt = 0.1; h_desired = 50000;
for b = 1:100 % Loop over burn time.
burn_time(b) = b;
% The following lines implement the formulas in the text.
m_0 = m_e + q*b; v_b = u*log(m_0/m_e) - g*b;
h_b = ((u*m_e)/q)*log(m_e/(m_e+q*b))+u*b - 0.5*g*b^2;
h_p = h_b + v_b^2/(2*g);
if h_p >= h_desired
% Calculate only if peak height > desired height.
t_p = b + v_b/g; % Compute peak time.
t_hit = t_p + sqrt(2*h_p/g); % Compute time to hit.
for p = 0:t_hit/dt
% Use a loop to compute the height vector.
k = p + 1; t = p*dt; time(k) = t;
if t <= b
% Burnout has not yet occurred.
h(k) = (u/q)*(m_0 - q*t)*log(m_0 - q*t)...
+ u*(log(m_0) + 1)*t - 0.5*g*t^2 ...
- (m_0*u/q)*log(m_0);
else
% Burnout has occurred.
h(k) = h_b - 0.5*g*(t - b)^2 + v_b*(t - b);
end
end
% Compute the duration.
duration(b) = length(find(h>=h_desired))*dt;
else
% Rocket did not reach the desired height.
duration(b) = 0;
end
end % Plot the results.
plot(burn_time,duration),xlabel('Burn Time (sec)'),...
ylabel('Duration (sec)'),title('Duration Above 50 000 Feet')