function  [coef, RSquared]=plotfit(xi,yi,polyorder)
% plotfit accepts data in the form of a single vector, a pair of vectors
% "x" and "y", or a 2xn or nx2 matrix with x in first row or column, 
% optionally fits the data to a polynomial of order "polyorder", plots 
% the data in red dots and the fit as a greeen line, and displays the 
% fit coefficients and R-squared in the upper left corner of the 
% graph. Polyorder=1 for straight line, =2 for quadratic (parabola) etc.
% Tom O'Haver, 2008. Version 2, October 2011 
%
% Examples
% x=[1 2 3 4 5];y=[ 0 2 3 2 0];
% plotfit(y); % plot y only vs index number, no fit.
% plotfit(x,y); % plot x vs y only, data in separate vectors
% plotfit([x;y]); % plot data only, data in matrix with x in first row
% plotfit([x;y]'); % plot data only, data in matrix with x in first column
% plotfit(y,2); % plot y vs index and fit to second order polynomial
% plotfit(x,y,3); % plot x vs y and fit to third order polynomial
% plotfit([x;y],3);  % plot and fit data in matrix
% [coef, RSquared]=plotfit([x;y],2) % return fit coefficients and r-squared

% process input arguments
if nargin==1, % Single argument, might be vector or matrix
    polyorder=0;
    datasize=size(xi);
    if isvector(xi), % x is vector, not matrix
        y=xi;
        x=1:length(y); % Use this only to create an x vector if needed
    else % x is a matrix, split it into x and y vectors
        if datasize(1)<datasize(2),xi=xi';end
        x=xi(:,1);
        y=xi(:,2);
        polyorder=0;
    end
end

if nargin==2, % Two arguments, might be x,y vectors or matrix and polyorder
    datasize=size(xi);
    if isscalar(yi),   %  matrix and polyorder
        if datasize(1)<datasize(2),xi=xi';end
        if isvector(xi),
            x=1:length(xi);
            y=xi';
            polyorder=yi;
        else
            x=xi(:,1);
            y=xi(:,2);
            polyorder=yi;
        end
    else %  x,y vectors
        x=xi;
        y=yi;
        polyorder=0;
    end
end

if nargin==3, % Tree arguments, must be x,y vectors and polyorder
    x=xi;
    y=yi;
end

if length(x)==length(y),
    if polyorder>0, % Plot and fit the data
        % Compute the fit
        coef=polyfit(x,y,polyorder);
        
        % Plot the data and the polynomial fit
        xx=linspace(min(x),max(x));
        yhat=polyval(coef,xx);
        plot(x,y,'.r',xx,yhat,'-g')
        axis([min(x) max(x) min(y) max(y)]);
        xlabel('x');ylabel('y')
        title(['Number of data points= ' num2str(length(x)) ] )
        % Compute the correlation coefficient and R-Squared
        cc=corrcoef(polyval(coef,x),y);
        RSquared=cc(2).^2;
        
        % Label the graph with the fit information
        text(min(x),max(y)-.05.*(max(y)-min(y)),['   Polynomial Order of fit = ' num2str(polyorder)] );
        text(min(x),max(y)-.1.*(max(y)-min(y)),['   Fit coefficients = ' num2str(coef)] );
        text(min(x),max(y)-.15.*(max(y)-min(y)),['   R-Squared = ' num2str(RSquared)] );
    else % Just plot the data, no fit.
        plot(x,y,'.r')
        axis([min(x) max(x) min(y) max(y)]);
        xlabel('x');ylabel('y')
        title(['Number of data points= ' num2str(length(x)) ] )
        coef=[];
        RSquared=[];
    end
else
    disp('Error: x and y must be the same size.')
    sizex=size(x)
    sizey=size(y)
end