% 1.72 Problem Set #2
% Problem 1, Least-squares fit of falling permeameter data

% All data are in inches; t is in seconds
data = [0	5
5	6
15	6.875
30	8
45	9.25
60	10.375
75	11.5
90	12.625
120	14.375
300	20.5
540	23.375
780	25
1200	26.375
2400	27.5
3600	27.5];

t = data(:,1); h = data(:,2);
H0 = 32; h0 = h(1); L = 28.5;
D = 2.5; d = 0.75;

a = (H0 + h0*(d/D)^2)/(1+(d/D)^2);
b = -(H0 - h0)/(1 + (d/D)^2);
c = (1 + (D/d)^2)/L;

% Fit K using all data points
y1 = -log((h-a)/b)/c;
K1 = t\y1;			% least-squares solution in Matlab
yhat1 = t*K1;		% predicted values for y using estimated K

% Fit K without last 5 points
y2 = y1(1:end-5); t2 = t(1:end-5);
K2 = t2\y2;
yhat2 = t2*K2;

% Plot results
subplot(211)
plot(t,y1,'o',t,yhat1,'-')
legend('Measured','Predicted',4)
title(['Fit using all data: K = ',num2str(K1),' in/sec'])
xlabel('t [s]'); ylabel('-log((h-a)/b)/c');
subplot(212)
plot(t2,y2,'o',t2,yhat2,'-')
legend('Measured','Predicted',4)
title(['Fit for t<=300: K = ',num2str(K2),' in/sec']);
xlabel('t [s]'); ylabel('-log((h-a)/b)/c');