clear; close all
%
%  Rocket problem p. 204
%
%  Using vectors
%
m = 50;             %% g
m = m/1000;         %% kg
tEngine = 0.15;     %% engine burn time, sec
thrust = 16;        %% thrust, N
vChute = -20;       %% chute deployment velocity, m/sec (down)
g = 9.81;           %% grav, m/sec^2
dt = 0.01;
longtime = 500;
n = 1;
%
%  Powered flight
%
acc = (thrust - m*g)/m;
t1 = 0:dt:tEngine;
v1 = acc*t1;
h1 = 0.5*acc*t1.^2;
%
%  Free flight
%
t2 = (tEngine+dt):dt:longtime;
v2 = v1(end) - g*(t2-t1(end));
ind = find(v2 < vChute);
h2 = h1(end) + v1(end)*(t2-t1(end)) - 0.5*g*(t2 - t1(end)).^2;
v2(ind) = []; t2(ind) = []; h2(ind) = [];
%
%  Chute open
%
t3 = (t2(end)+dt):dt:longtime;
h3 = h2(end) + vChute*(t3-t2(end));
ind = find(h3 < 0);
t3(ind) = []; h3(ind) = [];
v3 = vChute*ones(size(t3));
t = [t1 t2 t3]; v = [v1 v2 v3]; h = [h1 h2 h3];
subplot(2,1,1); plot(t,h,t2(end),h2(end),'o');
xlabel('Time, sec'); ylabel('height, m'); grid
subplot(2,1,2); plot(t,v,t2(end),v2(end),'o');
xlabel('Time, sec'); ylabel('velocity, m/sec'); grid