function [a,b,da,db]=fit_linear(x,y,varargin)
%fit_linear: 1st order linear regression.

%This function fits a dataset of x and y values to the linear equation
%y=ax+b. it is optional for the user to insert also y error values for the
%estimation of a and b errors.
%

%Yair Judkovsky, 21.8.2013

if nargin==2     %case of no error input - assuming that all errors are equal
 dy = ones(size(y));
else
 dy = varargin{1};
end

dy2=dy.^2;
x2=x.^2;
xy=x.*y;

S=sum(1./dy2);
Sx=sum(x./dy2);
Sy=sum(y./dy2);
Sxx=sum(x2./dy2);
Sxy=sum(xy./dy2);
D=S*Sxx-Sx^2;

a=(S*Sxy-Sx*Sy)/D;
b=(Sxx*Sy-Sxy*Sx)/D;

if nargin==2    %no error given - need to estimate it
    Ytheo=a*x+b;
    dy2=mean((Ytheo-y).^2)*ones(size(y));    
    S=sum(1./dy2);
    Sx=sum(x./dy2);
    Sxx=sum(x2./dy2);
    D=S*Sxx-Sx^2;
    da=sqrt(S/D);
    db=sqrt(Sxx/D);
else
    da=sqrt(S/D);
    db=sqrt(Sxx/D);
end

return