function [q,dq,xf,yf]=plot_polyfit_res(x,y,dy,n,xtxt,ytxt,xunt,yunt,axisv,ll)

%   PLOT_POLYFIT
%
%   [q,dq,xf,yf]=plot_polyfit(x,y,dy,n,xtxt,ytxt[,xunt,yunt,axisv,ll,res]);
%
%   Function to plot x data vs y data, and fit to a polynomial of order "n"
%   (with error bars if given), and order skipping available (see below).
%   Will plot AND return fit coefficients and their errors if output variables
%   are provided.  Otherwise only the plot is generated.
%
%   INPUTS:
%
%      x     = the X-axis data vector (column or row vector)
%      y     = the Y-axis data vector (column or row vector)
%      dy    = the Y-axis data vector error (or =1 for no errors given; in
%               this case the fit coefficient errors are rescaled such that
%               CHISQ/NDF = 1)
%      n     = polynomial order to fit; e.g. n=4 implies 0th, 1st, 2nd, 3rd,
%               & 4th order, while n=[1 0 1 1 0] implies only 0th, 2nd, & 3rd
%               order fit
%      xtxt  = (optional) the text string describing the X-axis data
%      ytxt  = (optional) the text string describing the Y-axis data
%      xunt  = (optional) the text string of X-axis units
%      yunt  = (optional) the text string of Y-axis units
%      axisv = (optional) 4-element vector of forced plot scale
%               [xmin xmax ymin ymax], or not used if 1-element
%      ll    = (optional) specifies location of fit coefficients on plot:
%               "ll"==0 -> upper right corner (the default)
%               "ll"==1 -> lower left  corner (historical)
%               "ll"==2 -> upper left  corner
%               "ll"==3 -> lower right corner
%      
%
%   OUTPUTS:
%
%      q     = (optional) if function is written as "q=polyfit(..." then
%               "q" will be a column vector of fit coefficients with row-1
%                as lowest order; if written as "plot_polyfit(...", then no
%                output is echoed (plot only)
%      dq    = (optional) error on "q" if function is written
%               "[q,dq]=polyfit(..."
%      xf    = (optional) more densely sampled abscissa if written as
%               "[q,dq,xf,yf]=polyfit(..."
%      yf    = (optional) fitted function if function is written as
%               "[q,dq,xf,yf]=polyfit(..."
%
%   P.Emma 11/13/91

%==========================================================================

% make column vectors out of everything

x=x(:);
y=y(:);
dy=dy(:);

% check input args for validity

nx=length(x);
ny=length(y);
if (length(dy)==1)
   bars=0;
else
   bars=1;
   if (ny~=length(dy))
      error('   Y-data and Y-error-data are unequal length vectors')
   end
end
if (nx~=ny)
   error('   X-data and Y-data are unequal length vectors')
end
if (ny<n+1)
   error('   Not enough data points to fit to this order')
end
if (~exist('xtxt'))
   xtxt=' ';
end
if (~exist('ytxt'))
   ytxt=' ';
end
if (~exist('xunt'))
   xunt=' ';
end
if (~exist('yunt'))
   yunt=' ';
end
if (~exist('ll'))
   ll=0;
end


% set up the fit

nn=length(n);
if (nn==1)
   mm=n+1;
   m=[0:1:mm];
else
   i=find(n);
   m=i-1;
   mm=length(m);
end
Q=[];
f=mean(abs(x));
x1=x/f;
for j=1:mm
   Q=[Q x1.^m(j)];
end

% do the fit

if (bars)
   [yfit,dyfit,p1,dp1,chisq]=fit(Q,y,dy);
else
   [yfit,dyfit,p1,dp1]=fit(Q,y);
end
for j=1:mm
   p(j)=p1(j)/(f^m(j));
   dp(j)=dp1(j)/(f^m(j));
end
if (mm==2)
   if ((m(1)==0)&(m(2)==1))
      r=p(2)*std(x)/std(y);  % linear correlation coefficient
   end
end
[Qr,Qc]=size(Q);
NDF=Qr-Qc
difs=(yfit-y).^2;
difbar=sqrt(mean(difs))
xsig=std(x);
ysig=std(y);

% do the plot

figure
subplot(2,1,1)

if (bars)
    plot_bars(x,y,dy,'*r')
else
   plot(x,y,'*r')
end
if (exist('axisv'))
   if (length(axisv)==4)
      axis(axisv)
   end
end
   lims=get(gca,'XLim');
   hold on
   xp=[lims(1):(lims(2)-lims(1))/250:lims(2)];
   Q=[];
   for j=1:mm
      Q=[Q xp(:).^m(j)];
   end
   p=p(:);
   yp=Q*p;
   plot(xp(:),yp(:),'-')
   hold off

title([ytxt ' VS ' xtxt])
xlabel([xtxt ' /' xunt])
ylabel([ytxt ' /' yunt])

% put polynomial coefficients on the plot

if (ll==1)
  xpn=0.20;  % lower left
  ypn=0.41;
elseif (ll==2)
  xpn=0.20;  % upper left
  ypn=0.91;
elseif (ll==3)
  xpn=0.60;  % lower right
  ypn=0.41;
else
  xpn=0.60;  % upper right
  ypn=0.91;
end
for j=1:mm
   text(scx(xpn),scy(ypn-j*0.025), ...
      [sprintf('p%1.0f=%9.5g+-%9.5g',m(j),p(j),dp(j))])
end
text(scx(0.20),scy(0.91),[sprintf('RMS=%5.3g ',difbar) yunt])
text(scx(0.20),scy(0.86),[sprintf('NDF=%5.0f ',NDF)])
if (bars)
   text(scx(0.2),scy(0.81),[sprintf('chisq/N = %5.3f',chisq)])
end
if (exist('r'))
   text(scx(0.2),scy(0.75),[sprintf('rho     = %6.3f',r)])
end

% download the results

if (nargout==1)
   q=p;
% Reverse order of coeff
q = q(end:-1:1);;
end
if (nargout==2)
   q=p;
   dq=dp;
   % Reverse order of coeff
   q = q(end:-1:1);;
   dq = dq(end:-1:1);
end
if (nargout>=3)
   q=p;
   dq=dp;
   xf=xp;
   yf=yp;
      % Reverse order of coeff
   q = q(end:-1:1);;
   dq = dq(end:-1:1);
   xf = xf(end:-1:1);;
   yf = yf(end:-1:1);
end


%  Plot residuals

subplot(2,1,2)

y_poly = polyval(q,x);
residuals = y-y_poly;
plot(x,residuals);

ytxt=[ytxt ' (FIT RESIDUALS)'];
title([ytxt ' VS ' xtxt])
xlabel([xtxt ' /' xunt])
ylabel([ytxt ' /' yunt])