% Example 2.10
% Program Example2_10.m , for solution of Example Problem 2.10, cooling billet.
clear all; format short g
n=17; %Measured temperature in Celsius. Time step is 96 sec.
y=[137.55 129.37 121.96 115.17 108.99 102.91 97.7 92.98 88.69 84.69 ...
    80.91 77.45 73.91 71.05 68.44 65.99 63.71]';
R=y'*y; t=0:1:n-1; t=t'; %Only for compactness is more than one statement on line.
figure(1); clf;
plot(t*96, y)
hold on
for p=1:6
  X=[ones(n,1)];
  for jj=1:p-1
      tj=t.^jj; X=[X tj];
  end
  b=X\y; e=y-X*b; sum=e'*e; DR=R-sum; 
  s2=sum/(n-p); Con=cond(X); [S]=svd(X);
  F_test = finv( 0.95, 1, n-p );
  mm=n-p; F=DR/s2; B(p,:)=[p mm sum s2 DR F F_test Con]; R=sum;
  if( F/F_test > 1.0 )
      e_opt = e;
      b_opt = b;
  end
end
figure(1);
plot( t*96, 100*e_opt, 'o', 'MarkerEdgeColor', [ 0.1 0.5 .1], 'linewidth', 1.0 )
plot( [ 0 1600], [0 0], 'r-', 'linewidth', 1.0 )
xlabel('time, seconds')
ylabel('Temperature, Celsius')
text(400,20,'100 times residuals')
text(500,110,'billet temperature history')
fprintf('  p   dof.    R_LS   s-squared      Delta-R       F      F_test Cond(X)\n'); 
fprintf(' %2d   %2d  %8.2f   %8.4f   %12.4f   %7.3f   %6.2f  %4.2f\n',B')
fprintf('\nb_opt: '); fprintf(' %9.6f', b_opt); fprintf('\n');
% print the 600 dpi figure to tiff format
fig=figure(1);
print(fig,fullfile(pwd, 'Figure2_4.tif'),'-dtiff','-r600')