%CheE 216 - Polynomial Fitting example: using the normal equation method

% first, start with the input data
x = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50];
y = [17, 24, 31, 33, 37, 37, 40, 40, 42, 41];

% plot P vs V, with red squares
figure(1);
plot(x,y,'sr');

% The model is PV/RT = 1 + A_1 V^-1 + A_2 V^-2.
% To linearize the model, we change to PV/RT-1

Pfit = x;
Vfit = y;

% Now, construct the Z matrix of function values

Z = [Vfit' (Vfit.^2)'];

% use general linear fit theory to get the value for the parameters A
% Z'*Z A = Z' Y' 
%M = inv(Z'*Z)
%A = M*Z'*Pfit';
% X = A\B is the solution to the equation Ax = B (See 'help slash')
A = Z'*Z \ Z'*Pfit';
%compute the variance of the pressures

St = var(Pfit);
disp(['St =' num2str(St)]);

%Note -- can we use polyfit here?

%Compute the residuals
Pmod = A(1)*Vfit + A(2)*Vfit.^2;

Sr = sum((Pmod - Pfit).^(2));
disp(['Sr =' num2str(Sr)]);

% compute the error in the fit
numdat = size(P);
s = sqrt(Sr/(numdat(2)-2));

% compute the r^2 
r2 = (St-Sr)/St;
disp(['r^2 =' num2str(r2)]);

% compute the uncertiainty in the parameters
M = inv(Z'*Z)';
ErrA(1) = sqrt(M(1,1)*s^2);
ErrA(2) = sqrt(M(2,2)*s^2);

disp(['A1 = ' num2str(A(1)) '+/-' num2str(ErrA(1))])
disp(['A2 = ' num2str(A(2)) '+/-' num2str(ErrA(2))])

% plot the data with the fit

Vline = linspace(0.035,0.08,1000);
Vline2 = Vline.^(-1);

Pline = A(1)*Vline + A(2)*Vline.^2;

figure(2)
hold on;
plot(Pline,Vline)
plot(Pfit,Vfit,'sr');
hold off;

% plot the Pressure vs V with the fit. 

Pline2 = ((Pline+1)*R*T)./Vline2;

figure(3);
hold on;
plot(Pline2,Vline2);
plot(P,V,'sr');
hold off;