function [ par dpar ] = linfit( x, y, dy )
%LINFIT Linear fit
%   [par dpar] = LINFIT(x, y, dy) calculates the linear fit of the 
%   form f(x) = par(1) + par(2)*x according to eqns 8.37-8.39 in Taylor.
%   The 1 sigma uncertainties on the fit parameters are also calculated
%   (see problem 8.19 in Taylor) and are returned as dpar.
%   If the errors dy are not specified, unweighted data is assumed
%
if length(x) ~= length(y)
    disp('Length mismatch in data vectors!'); 
    return; 
end
if nargin<3 | isempty(dy)
    dy=ones(1,length(x)); 
end
if length(dy) ~= length(y)
    disp('Length mismatch in data vectors!'); 
    return; 
end

% Weights
wt = 1./dy.^2;

sumwt = sum(wt);
sumwtx = sum(wt.*x);
sumwtxx = sum(wt.*x.*x);
sumwtxy = sum(wt.*x.*y);
sumwty = sum(wt.*y);

delta = sumwt*sumwtxx - sumwtx^2;

par(1) = (sumwtxx*sumwty - sumwtx*sumwtxy)/delta;
par(2) = (sumwt*sumwtxy - sumwtx*sumwty)/delta;

dpar(1) = sqrt(sumwtxx/delta);
dpar(2) = sqrt(sumwt/delta);

end

