call linfit(x,y,sigmay,npts,mode,int,sigmaa,slope,sigmab,r)
	write(6,50)int,sigmaa
50	format(1x,'y intercept =',f10.4,'     stand. dev. =',f10.4)
	write(6,60)slope,sigmab
60	format(1x,'slope =',f10.4,'    stand. dev. =',f10.4)
	write(6,70)r
70	format(1x,'corr. coef. = ',f10.4)
	subroutine linfit(x,y,sigmay,npts,mode,int,sigmaa,slope,sigmab,r)
c   least squares straight line fitting subroutine
c   from bevington
	real int
	dimension x(1),y(1),sigmay(1)
c	accumulate weighted sums
11	sum=0.
	sumx=0.
	sumy=0.
	sumx2=0.
	sumxy=0.
	sumy2=0.
21	do 50 i=1,npts
	x1=x(i)
	y1=y(i)
	if (mode.eq.0) goto 36
	if (mode.gt.0) goto 38
	if (y1.lt.0.0)goto 34
	if (y1.eq.0.0)goto 36
	weight=1./y1
	goto 41
34	weight=1./(-y1)
	goto 41
36	weight=1.
	goto 41
38	weight=1./sigmay(i)**2.
41	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
50	continue
c	calculate coefficients
51	delta=sum*sumx2-sumx*sumx
	int=(sumx2*sumy-sumx*sumxy)/delta
	slope=(sumxy*sum-sumx*sumy)/delta
	if (mode.eq.0) goto 64
	varnce=1.
	goto 67
64	c=npts-2
	varnce=(sumy2+int*int*sum+slope*slope*sumx2-2.*(int*sumy+slope*sumxy
     1   -int*slope*sumx))/c
67	sigmaa=(varnce*sumx2/delta)**.5
	sigmab=(varnce*sum/delta)**.5
	r=(sum*sumxy-sumx*sumy)/(delta*(sum*sumy2-sumy*sumy))**.5
	return
	end