clear all; close all; clc;

g = -9.81;
mass = 1.2;
mdot = 0.0182;
dt = 0.001;

vprev = 0;
hprev = 0;
i = 1;

for t = 0 : dt : 3.45
    m = mass - mdot*t;
    if t < 0.148
        thrust = 0;
    elseif (t > 0.147) && (t < 0.419)
        thrust = 65.165.*t - 2.3921;
    elseif (t > 0.420) && (t < 3.382)
        thrust = 0.8932*t^6 - 11.609*t^5 + 60.739*t^4 - 162.99*t^3 ...
            + 235.6*t^2 - 174.43*t + 67.17;
    elseif (t > 3.381) && (t < 3.451)
        thrust = -195.78*t + 675.11;
    elseif t > 3.45
        thrust = 0;
    end
   
    acceleration(i) = (thrust / mass) + g;
    velocity(i) = acceleration(i) * dt + vprev;
    height(i) = velocity(i) * dt + hprev;
    
    vprev = velocity(i);
    hprev = height(i);
    i = i + 1;
end
figure(1)
hold on
t = 0 : dt : 3.45;
plot(t + 0.7712, velocity - max(velocity))

t1 = 0 : 0.001 : 0.7712;
v = -9.81*t1;
plot(t1, v)

figure(2)
h = 0.5*-9.81*t1.^2
plot(t1, h)