% fitFall.m
% Fit 'experimental' data of a falling body in air, assuming a force
% F=m*g-c*v, with c a drag coefficient to be determined. We approximate
% v(t)=g*t-c*g/2/m*t^2 and x(t)=c2*t^2+c3*t^3, with c2=g/2 and c3=-c*g/6/m
% Jose M Soler. Nov.2015

% Set data (all in MKS system units)
clear all
m = 0.325;
X = [0.0  3.2  6.4  9.6  12.8 16.0]';  % 'measured' positions
T = [0.00 0.83 1.16 1.45 1.65 1.88]';  % 'measured' times

% Fit data, excluding one point each time
n = numel(X);
for i = 1:n
    t = [T(1:i-1);T(i+1:n)];
    x = [X(1:i-1);X(i+1:n)];
    M = [t.^2 t.^3];
    C = M\x;
    c2 = C(1);
    c3 = C(2);
    g(i) = 2*c2;
    c(i) = -6*m*c3/g(i);
end

% Print mean values of g,c and their standard deviations
dg = standDev(g);  % notice: this must be done before overwritting g
dc = standDev(c);
g = meanValue(g);  % overwrites vector g
c = meanValue(c);
xFit = (g/2)*T.^2 - (c*g/6/m)*T.^3;
xErr = standDev(X-xFit);
fprintf('g = %4.2f +/- %4.2f m/s^2 \n',g,dg)
fprintf('c = %5.3f +/- %5.3f N*s/m \n',c,dc)
fprintf('mean fit error = %5.2f m \n',xErr)

% Plot data and fit
t = 0:0.01:2;
x = (g/2)*t.^2 - (c*g/6/m)*t.^3;
plot(T,X,'o',t,x)
legend('data','fit')
xlabel('time (s)')
ylabel('position (m)')
