%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                  %  
%  MATLAB Fitting Template Script  %
%          MIT Junior Lab          %  
%       Spring 2005 Ver:1.0        %
%                                  %
%   Junior Lab Technical Staff     %
%                                  %   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%
%%% LOADING YOUR DATA %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

clear;   % Clears all variables from memory

load gauss3.dat;  % Loads sample data set 'gauss3.dat' into memory
                  % This script assumes that the data is in a text file, 
                  % divided into 2 or 3 columns (depending on whether there is 
		          % a column for the error)

x = gauss3(:,1);   % Assigns the first column of mydata to a vector called 'x'
y = gauss3(:,2);   % Assigns the second column of mydata to a vector called 'y'



% *** Note: There are several different ways you might assign your error vector
% ***       You might want to use one of these three methods below.

%sig = mydata(:,3); % Assigns the third column to a vector called 'sig'


sigma = 2.5;  % Constant error value, variance of gauss3 data set is 6.25.
sig = ones(size(x))*sigma;  % creates a constant vector the same size as 'x' 
                            % with the value sigma


%sig = sqrt(y);  % sets sigma to the square root of 'y' for Poisson statistics

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CREATING THE FIT MODEL %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% For non-linear fits, you will need to create a Matlab function for the 
% functional form you want to fit.  This function will need to be in its 
% own M-file.  For example, if you are trying to fit a sine wave, your 
% function M-file might look like this:
%       function f = sinfit(x,a)
%       f = a(1)*sin(a(2)*x+a(3));
%
% Save the M-file with the same name as the function.  In this case it
% would be "sinfit.m"
%
% For the gauss3 dataset, you could have the following M-file:
%       function f = g3function(x,a)
%       f = a(1)*exp(-a(2)*x)+a(3)*exp(-((x-a(4))/a(5)).^2)+a(6)*exp(-((x-a(7))/a(8)).^2);
%
% saved in a file called "g3function.m"

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% STARTING PARAMETERS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

a0 = [100 .01 100 100 20 40 150 20];  
      % This is a vector of the starting values of fitting parameters.  Assign 
      % them in the order that they appear in your function.

%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FITTING YOUR DATA %%%
%%%%%%%%%%%%%%%%%%%%%%%%% 

% Junior Lab supports two Matlab fitting scripts which you may find useful:
% 'fitnonlin' for non-linear fits and 'fitlin' for linear fits.  Both of
% these algorithms are based on the methods found in Bevington and
% Robinson.  These scripts are well commented and we encourage you to
% inspect them to strengthen your understanding of and intuition for
% fitting data.
% 

% ** NON-LINEAR FITS **

[a,aerr,chisq,yfit] = fitnonlin(x,y,sig,'g3function',a0);
   % Here the data in 'x' and 'y' is fit to the 
   % function represented in the file g3function.m.  Note the way g3function.m
   % invoked: in single quotes without the '.m'
   %    Inputs:  x -- the x data to fit
   %             y -- the y data to fit
   %             sig -- the uncertainties on the data points
   %             fitfun -- the name of the function to fit to
   %             a0 -- the initial guess at the parameters 
   %
   %    Outputs: a -- the best fit parameters
   %             aerr -- the errors on these parameters
   %             chisq -- the final value of chi-squared
   %             yfit -- the value of the fitted function
   %                     at the points in x
   %
   % * NOTE * - since fitnonlin is an iterative function, you can vary
   % the step size and tolerance used in the fitting process.  These are
   % specified within the 'fitnonlin.m' function
   
   
% ** LINEAR FITS **

% If you wish to fit your data to a straight line, you can use the function
% 'fitlin.m'  This function is also taken directly from the pages of Bevington 
% and Robinson Ch. 6 (p. 114) and it is very easy to follow.  The usage of
% fitlin is as follows:
%
% [a,aerr,chisq,yfit] = fitlin(x,y,sig)
%
%    Inputs:  x -- the x data to fit
%             y -- the y data to fit
%             sig -- the uncertainties on the data points
%
%    Outputs: a -- the best fit parameters
%             aerr -- the errors on these parameters
%             chisq -- the value of chi-squared
%             yfit -- the value of the fitted function
%                     at the points in x


[a', aerr']

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% EVALUATING THE GOODNESS OF YOUR FIT %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% It's useful to examine the Reduced Chi-Square of your fit.  Recall from
% Bevington that this is simply the Chi-Square divided by the number of
% degrees of freedom (ie the number of data points less the number of
% fitting parameters).

RChiSquare = chisq/(length(x)-length(a))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PLOTTING YOUR DATA AND FIT %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

figure (1);  % opens figure 1 and makes it active
clf;         % clears the active figure
hold on;     % Makes matlab plot without clearing the graph

errorbar (x,y,sig,'b.');  % Plots 'y' vs. 'x' with errobars of plus/minus 'sig'
plot(x,yfit,'r-');        % Plots  the fitted curve

axis([x(1),x(length(x)),min(y)-2*max(sig),max(y)+2*max(sig)]);  
    % sets visible range of the plot

title('Matlab Fit of NIST Generated Data Set','fontsize',16);
       % Places the title on the graph
xlabel('wavelength (nm)','fontsize',14);  % Labels the 'x' axis
ylabel('y axis label (units)','fontsize',14);  % Labels the 'y' axis

str1=num2str(RChiSquare,2);
text(5,120,['\chi^2_{\nu-1} = ' str1]); % Plots the reduced chi-square