% regressFall.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]';
t = [0.00 0.83 1.16 1.45 1.65 1.88]';

% Fit data by linear regression
T = [t.^2 t.^3];
[C,Cint] = regress(x,T);   % substitutes C = T\x;
c2 = C(1);
c3 = C(2);
c2int = Cint(1,:);
c3int = Cint(2,:);
g = 2*c2;
gint = 2*c2int;
dg = (gint(2)-gint(1))/3.92;
c = -6*m*c3/g;
cint = -6*m*c3int(2:-1:1)/g;
dc = (cint(2)-cint(1))/3.92;   % 3.92*sigma is the 95% confidence interval

% Find fit error
xFit = (g/2)*t.^2 - (c*g/6/m)*t.^3;
dx = standDev(x-xFit);

% Print results
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('g_interval = [%4.2f,%4.2f] m/s^2 \n',gint)
fprintf('c_interval = [%5.3f,%5.3f] N*s/m \n',cint)
fprintf('mean fit error = %5.2f m \n',dx)

% Plot data and fit
time = [0:0.1:2];
plot(time,c2*time.^2+c3*time.^3,'-', ...
     time,c2int(1)*time.^2+c3int(1)*time.^3, ...
     time,c2int(2)*time.^2+c3int(2)*time.^3, t,x,'ok')
xlabel('t (s)'), ylabel('x (m)')
