% Unit 4 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%===============================================================
% free-fall in air (F = mg - cv) determining g and c
%===============================================================
clear all, clc

% experimental data
m =   0.325;                            % kg
x = -[0.0  3.2  6.4  9.6  12.8 16.0];   % m
t =  [0.00 0.83 1.16 1.45 1.65 1.88];   % sec

figure(2)
plot(t,x,'o','MarkerSize',10,'LineWidth',3), hold on
xlabel('t (s)'), ylabel('x (m)')


% Fit data to analytical function x(t) = a1 t^2 + a2 t^3 (with a1=g/2 and a2=-cg/6m)

% 1. use all points
i = [1:length(x)];
M = [t(i).^2; t(i).^3]';
a = inv(transpose(M)*M)*transpose(M)*x(i)';
g = 2*a(1)
c = -6*m*a(2)/g


time = linspace(min(t),max(t),100);
plot(time,a(1)*time.^2+a(2)*time.^3,'-r')

% 2. use sub-set of points
isubset = [1, 2, 6];
M       = [t(isubset).^2; t(isubset).^3]';
a       = inv(transpose(M)*M)*transpose(M)*x(isubset)';
g       = 2*a(1)
c       = -6*m*a(2)/g
time    = linspace(min(t),max(t),100);
time    = linspace(min(t),2.5,100);
plot(time,a(1)*time.^2+a(2)*time.^3,'-g')

legend('experimental data','best-fit (all points)','best-fit (subset: [1,2,6])')

% root-finding for ground level (-20m)
g = @(t)(a(1).*t.^2+a(2).*t.^3-(-20));

a=1.88;
b=2.5;

while b-a>1e-3
   t0 = (a+b)/2;
   if g(a)*g(t0) > 0
       a = t0;
   else
       b = t0;
   end
end
t0



