%
%   CFIT_UO2K_1.M   Solution to curve fitting example problem 
%                Experimental data for UO2 thermal conductivity
%
%   The goal of this program is to determine the best mathematical model that fits
%   the set of experimental data.  The experimental data are associated with the 
%   temperature dependence of the thermal conductivity of uranium dioxide UO2. The 
%   two models to be compared are:
%
%     A. power law model       y(x) = a1*x^a2
%
%     B. quadratic polynomial  y(x) = c1 + c2*x + c3*x^2 
%
%        where x -> temperature (F) 
%              y -> thermal conductivity (Btu/hr-ft-F)
%
%   Both qualitative comparisons via graphical means and quantitative 
%   analysis by statisical measures are performed.
%
%   We will use MATLAB's polyfit function to do the curve fitting.
%
%   Program written by J. R. White, UMass-Lowell (last update: Nov. 2015)
%

      clear all,   close all,   nfig = 0;
%
%   write header for summary results
      fid = 1;
      fprintf(fid,'\n         Summary Results from the CFIT_UO2K_1.M Program \n');
%
%   base data
      N = 10;
      x = [200 400 600 800 1000 1200 1400 1600 1800 2000]';   % temp (F)
      y = [4.5 3.5 2.8 2.5  2.2  2.0  1.6  1.5  1.4  1.3]';   % k (Btu/hr-ft-F)
      fprintf(fid,'\n    Experimental UO2 Temperature and Thermal Conductivity Data  \n');
      fprintf(fid,'    T (F)           = ');  fprintf(fid,'%5.0f ',x);  fprintf(fid,'\n');
      fprintf(fid,'    k (Btu/hr-ft-F) = ');  fprintf(fid,'%5.1f ',y);  fprintf(fid,'\n');
%
%   Model A  Power Law
%
%   this model must first be converted to linear form
      t = log(x);    s = log(y);
%
%   now fit linear curve to the [t,s] data and convert back to original model
%   note that Matlab stores the coeffs of the resultant polynomial in decreasing order
      cc1 = polyfit(t,s,1);
      a1 = exp(cc1(2));     a2 = cc1(1);
      fprintf(fid,'\n    Coeffs for Power Law Model (Model A) \n');
      fprintf(fid,'    a1 = %7.3e   a2 = %7.3e  \n',a1,a2);
%
%   do error analysis for model A
      ym = sum(y)/N;            % mean value of y
      St = sum((y-ym).^2);      % sum of squares of deviation (original data)
      f1 = a1*x.^a2;            % evaluate math model at expt points
      Sr1 = sum((f1-y).^2);     % sum of squares of residuals (model - org data)
      r1 = (St - Sr1)/St;       % coeff of determination (approx unity for good model)
%
%   Model B  Quadratic polynomial
%
%   fit quadratic polynomial curve to the [x,y] data
      cc2 = polyfit(x,y,2);
      b3 = cc2(1);  b2 = cc2(2);  b1 = cc2(3);
      fprintf(fid,'\n    Coeffs for Quadratic Poly Model (Model B) \n');
      fprintf(fid,'    b1 = %7.3e   b2 = %7.3e   b3 = %7.3e\n',b1,b2,b3);
%
%   do error analysis for model B
      f2 =  b1 + b2*x + b3*x.^2;% evaluate math model at expt points
      Sr2 = sum((f2-y).^2);     % sum of squares of residuals (model - org data)
      r2 = (St - Sr2)/St;       % coeff of determination (approx unity for good model)
%
%   edit quantitative statisical info
      fprintf(fid,'\n        Summary Data            Power Law   Quadratic Fit   \n');
      fprintf(fid,'    sum of squares of errors    %7.3e    %7.3e   \n',Sr1,Sr2);
      fprintf(fid,'    coeff of determination      %7.3e    %7.3e   \n',r1,r2);
%
%   plot data (use more points for smooth curve)
      xx = linspace(x(1),x(N),50);
      y1 = a1*xx.^a2;
      y2 = polyval(cc2,xx);
      nfig = nfig+1;   figure(nfig)
      plot(x,y,'ro',xx,y1,'b-',xx,y2,'g--','LineWidth',2),grid
      title('CFIT\_UO2K\_1:  Model Prediction Versus Experimental Data')
      xlabel('Temperature (^oF)'),ylabel('UO_2 Thermal Conductivity (Btu/hr-ft-F)')
      legend('Experimental Data', ...
            ['Power Law (r^2 = ',num2str(r1,'%5.3f'),')'], ...
            ['Quadratic Poly (r^2 = ',num2str(r2,'%5.3f'),')'])
%
%   end of problem 
