% Himmelblau & Riggs, 7th ed
% Workbook problem 23.2
% given T (degC) and Cp (J/gmol/degC)
% find coefficients of best-fit 3rd-order (degree) polynomial
% Cp = coef(1)*T^3 + coef(2)*T^2 + coef(3)*T^1 + coef(4) 

format short g

% enter input data
T = [0 18 25 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500]';
Cp = [29.062 29.075 29.075 29.142 29.292 29.514 29.782 30.083 30.401 30.711 ...
    31.020 31.317 31.585 31.865 32.108 32.338 32.556 32.761]';

% fitting a polynomial function to data is a linear algebra
% problem because the system is linear in the unknown coefficients
% of the polynomial

% it is such a common problem that Matlab has a standard function
% to do this - polyfit - at command line, enter >> help polyfit

fprintf('\n3rd-order polynomial coefficients using polyfit \n')
coef1 = polyfit(T,Cp,3)'

% to plot the data and the curve you should compute the "model" values
% at more points that the data points in order to make sure the function
% doesn't "oscillate" between the data points

% define a new array Tp for T prime with more points than data T
Tp = 0:1500;

% use std function polyval to compute Cpmodel at points Tp
Cpmodel = polyval(coef1,Tp);

plot(T,Cp,'o',Tp,Cpmodel)
title('H&R Workbook, 23.2, data fit to 3rd-order polynomial with polyfit')
ylabel('Cp (J/gmol/degC)')
xlabel('T (degC)')

% in addition, solve using matrix algebra to see if we get same results

X = [T.^3 T.^2 T.^1 T.^0]; % note use of dot .^ - see Matlab help > operators

fprintf('3rd-order polynomial coefficients using matrix algebra \n\n')
coef2 = inv(X'*X)*X'*Cp

fprintf('3rd-order polynomial coefficients using matrix and \\ operator \n')
coef3 = X\Cp