%
%   AL28_HALFLIFE.M   Solution to curve fitting example problem 
%                Experimental data for Al28 activity vs time should allow
%                  estimation of the halflife of Al28
%
%   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 
%   Al28 activity vs time. Two approaches are compared:
%
%     A. Linear fit to a transformed exponential decay model
%             A(t) = Ao*exp(-lam*t)   and Thalf = ln(2)/lam
%        the nonlinear model can be linearized to give
%             y(t) = a + b*t    with  y = ln(A),  Ao = exp(a), and lam = -b
%
%     B. Use nonlinear fitting routine directly (nlinfit) to base nonlinear model
%             A(t) = Ao*exp(-lam*t)   and Thalf = ln(2)/lam
%
%   Both qualitative comparisons via graphical means and quantitative 
%   analysis by statisical measures are performed.
%
%   We will use MATLAB's polyfit function to do the linear curve fitting and its
%   nlinfit function to do the nonlinear fit.
%
%   Program written by J. R. White, UMass-Lowell (updated Dec. 2015)
%

      clear all,   close all,   nfig = 0;
%
%   write header for summary results
      fid = 1;
      fprintf(fid,'\n         Summary Results from the AL28_HalfLife.M Program \n');
%
%   base data
      N = 21;
      t = [  0.00  30.00  60.00  90.00 120.00 150.00 180.00 ...
           210.00 240.00 270.00 300.00 330.00 360.00 390.00 ...
           420.00 450.00 480.00 510.00 540.00 570.00 600.00]';  % time (seconds)
      A = [1.27E+04 8.26E+03 8.26E+03 6.75E+03 6.73E+03 5.50E+03 4.19E+03 ...
           2.95E+03 3.54E+03 2.62E+03 2.40E+03 2.20E+03 1.97E+03 1.58E+03 ...
           1.08E+03 1.02E+03 1.07E+03 9.08E+02 6.43E+02 6.61E+02 4.03E+02]'; % activity
      fprintf(fid,'\n    Experimental Al28 Activities (counts per second) vs Time Data  \n');
      fprintf(fid,'    t (sec)    = ');  fprintf(fid,'%5.0f ',t);  fprintf(fid,'\n');
      fprintf(fid,'    A(t) (cps) = ');  fprintf(fid,'%5.0f ',A);  fprintf(fid,'\n');
%
%   Model A  Linearized Exponential Decay Law (use polyfit)
%
%   this model must first be converted to linear form
      y = log(A);
%
%   now fit linear curve to the [t,y] data and convert back to original model
%   note: Matlab stores the coeffs of the resultant polynomial in decreasing order
      cc1 = polyfit(t,y,1);
      AoA = exp(cc1(2));     lamA = -cc1(1);     ThalfA = log(2)/lamA/60;
      fprintf(fid,'\n    Coeffs for Linearized Exponential Decay Law Model (Model A) \n');
      fprintf(fid,'    Ao = %7.3e   lam = %7.3e   Thalf = %7.3f min  \n',AoA,lamA,ThalfA);
%
%   do error analysis for model A (use the nonlinear math model)
      Am = sum(A)/N;           % mean value of activities
      St = sum((A-Am).^2);     % sum of squares of deviation (original data)
      fA = AoA*exp(-lamA*t);   % evaluate math model A at expt points
      SrA = sum((fA-A).^2);    % sum of squares of residuals (model - org data)
      r2A = (St - SrA)/St;     % coeff of determination (approx unity for good model)
% 
%   Model B  Nonliner Exponential Decay Law (use nlinfit)
%
      paramg = [A(1) log(2)/120];  % initial guess for model parameters
      param = nlinfit(t,A,@al28_nlmodel,paramg);
      AoB = param(1);  lamB = param(2);   ThalfB = log(2)/lamB/60;
      fprintf(fid,'\n    Coeffs for Nonliner Exponential Decay Law (Model B) \n');
      fprintf(fid,'    Ao = %7.3e   lam = %7.3e   Thalf = %7.3f min\n',AoB,lamB,ThalfB);
%
%   do error analysis for model B
      fB =  AoB*exp(-lamB*t);  % evaluate math model at expt points
      SrB = sum((fB-A).^2);    % sum of squares of residuals (model - org data)
      r2B = (St - SrB)/St;     % coeff of determination (approx unity for good model)
%
%   edit quantitative statisical info
      fprintf(fid,'\n        Summary Data        Linearized Model A   Nonlinear Model B \n');
      fprintf(fid,'    sum of squares of errors    %7.3e        %7.3e   \n',SrA,SrB);
      fprintf(fid,'    coeff of determination      %7.3e        %7.3e   \n',r2A,r2B);
%
%   plot data
      tt = linspace(t(1),t(N),50);
      AA = AoA*exp(-lamA*tt);
      AB = AoB*exp(-lamB*tt);
      nfig = nfig+1;   figure(nfig)
      plot(t,A,'ro',tt,AA,'b-',tt,AB,'g--','LineWidth',2),grid
      title('Al28\_HalfLife:  Model Prediction Versus Experimental Data')
      xlabel('Time (sec)'),ylabel('Al28 Activity (cps)')
      legend('Experimental Data','Model A (linearized)','Model B (nonlinear)')
%
%   end of problem 
