function [a,b,da,db]=linreg(x,y,dy)
%  linear regression 
% usage:[a,b,da,db]=linreg(x,y,dy)
% where a is the intercept and b is the slope and da and db are 
% one sigma uncertainty estimates.

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


  sigsq=dy.^(2);
  Sy=sum(y./sigsq);
		Sx=sum(x./sigsq);
		Sxx=sum((x.^2)./sigsq);
		Sxy=sum((x.*y)./sigsq);
  S=sum(sigsq.^(-1));
  del=S*Sxx-Sx^2;

		a=(Sxx*Sy-Sx*Sxy)/del;
		b=(S*Sxy-Sx*Sy)/del;
		
		da=sqrt(Sxx/del);
		db=sqrt(S/del);

	
