%---------------------------------------------------------------
%					linfit				(L. Braile 1/20/95)
%---------------------------------------------------------------
function [a,sigmaa,b,sigmab,r]=linfit(x,y,sigmay,npts,mode)
% function to calculate least squares straight line
% y=a + b*x and std dev of coefficients.
%---------------------------------------------------------------
sum=0;
sumx=0;
sumy=0;
sumx2=0;
sumxy=0;
sumy2=0;
for i=1:npts
	x1=x(i);
	y1=y(i);
	if mode<0
	if(y1)<0; weight=1/(-y1); end
	if(y1)==0; weight=1; end
	if(y1)>0; weight=1/y1; end
	end
	if mode==0; weight=1;end
	if mode>0;weight=1/sigmay(i)^2; end
	sum=sum + weight;
	sumx=sumx + weight*x1;
	sumy=sumy + weight*y1;
	sumx2=sumx2 + weight*x1*x1;
	sumxy=sumxy + weight*x1*y1;
	sumy2=sumy2 + weight*y1*y1;
end
%---------------------------------------------------------------
% Calculate coefficients and standard deviations
delta=sum*sumx2 - sumx*sumx;
a=(sumx2*sumy - sumx*sumxy)/delta;
b=(sumxy*sum - sumx*sumy)/delta;
if mode<0; varnce=1; end
if mode>0; varnce=1; end
if mode==0
	c=npts-2;
	varnce=(sumy2 + a*a*sum + b*b*sumx2 - 2*(a*sumy + b*sumxy - a*b*sumx))/c;
end
sigmaa=sqrt(varnce*sumx2/delta);
sigmab=sqrt(varnce*sum/delta);
r=(sum*sumxy - sumx*sumy)/sqrt(delta*(sum*sumy2 - sumy*sumy));
%---------------------------------------------------------------