clear
load thrust
load velocity

mass = 55;
t(1) = 0;
a(1) = thrust(1);
v(1) = 0;
d(1) = 0;
r = 32.174*.00273;
A = 2;
color = 'rgbyk';
disp('Drag      T/O Velocity (ft/s)')
for drag = .01:.01:.05
for ii = 2:10000
    t(ii) = (ii-1)/100;
    a(ii) = interp1(velocity,thrust,v(ii-1))/mass*32.174 - (drag*A*r*v(ii-1)^2/2 + .02*mass);
    v(ii) = v(ii-1) + (a(ii-1)+a(ii))/2*(t(ii)-t(ii-1));
    d(ii) = d(ii-1) + (v(ii-1)+v(ii))/2*(t(ii)-t(ii-1));
    if d(ii) > 170,break,end
end

% figure(1)
% plot(t,[a',v',d'])
% legend('Acceleration','Velocity','Distance')
% xlabel('Time (s)')
% ylabel('^f^t/_s2    |    ^f^t/_s    |    ft')
figure(2)
plot(d,v,color(drag*100))
hold on
xlabel('Groundroll Distance (ft)')
ylabel('Velocity (^f^t/_s)')
disp(sprintf('%g            %g',drag,max(v(ii))));
end
hold off