% 10.34 - Fall 2006
% HW Set#2, Problem #2a
% Rob Ashcraft - Sept. 11, 2006

% This function can do two things: (1) SOlve for the polynomial
% coefficients of an Nth-order polynomial given N+1 data points, (2) Find
% the best fit coefficients for an Nth-order polynomial given more than N+1
% data points.  
% f(x) = a0 + a1*x + a2*x^2 + ... + aN*x^N
%If you want to always use the largest polynomial possible for the data,
%N_order does no have to be passed to the function

function a_vec = calc_poly_coeff(x,f,N_order)
% Input: 
%   x = data for the independent variable 
%   f = data for dependent variable
%   N_order = order of the polynomial you want to fit
%       
% Output
%   a_vec = polynomial coefficients = [a0; a1; ... aN]
%       

% First, check to see if the number of arguments passed to the function is
% equal to 2 (i.e. the user only specified the x and f data).  If this is
% true, the user wants to largest polynomial possible and the coefficients
% can be found by solving a linear system: A*x = f. "nargin" in a built-in
% matlab command that tell you have many inputs were passed to a function

if(nargin == 2)
    
    % Generate the system matrix, with rows of the form: [1 x x^2 ... x^N]
    for i=1:length(x)
        for j=1:length(x)
            X_mat(j,i) = x(j)^(i-1);
        end
    end
    
    % determine and display the condition number of the system matrix
    disp(['The condition number for the polynomial is:  ',num2str(cond(X_mat),'%5.3e')]);
    disp('  ');
    
    % use the backslash command solve for the parameters
    a_vec = X_mat\f;

% If the user supplied the desired order, there are three possibilities:
% (1) you have too many parameters for the amount of data, (2) the number
% of parameters and data is equal and you can solve a linear system, or (3)
% the data set is larger than the number of parameters and a linear
% regression must be performed to find the parameter.  
elseif(nargin == 3)
    if(N_order+1 > length(f))
        disp(['You do not have enough data to solve for the coefficients of the order ',num2str(N_order),' polynomial.']);
        return;  % kicks out of the function... probably would cause errors.
        
    elseif(N_order+1 == length(f))
        for i=1:length(x)
            for j=1:length(x)
                X_mat(j,i) = x(j)^(i-1);
            end
        end

        a_vec = X_mat\f;       

    else
        
        for i=1:(N_order+1)
            for j=1:length(x)
                X_mat(j,i) = x(j)^(i-1);
            end
        end

        % in this case, a linear regression is performed to get the parameters
        a_vec = X_mat'*X_mat\X_mat'*f
        
    end
end