; linear fit a*x+b to data cube
; vectors for x and w can be provided if possible (to save memory)

; $id: lincfit.pro
function lincfit,x,y,w=w
; cube or vector of x values
; cube of y values
; cube or vector of weights
; returns cube slope, constant term, slope error, constant error, chi square
; written by Bringfried Stecklum

xs=size(x)
ys=size(y)
ws=size(w)
n=ys(3)

; vectors provided?
if (xs(0) eq 1) then begin
	; weights provided?
	if not (ws(0) eq 1) then begin
		w=x
		w(*)=1.
	endif
	; make arrays
	p=fltarr(ys(1),ys(2))
	px=p
	px2=p
	px(*)=total(w*x)
	px2(*)=total(w*x^2)
	py=p
	pxy=p
	p(*)=total(w)
	for i=0,n-1 do begin
		py=temporary(py)+w(i)*y(*,*,i)
		pxy=temporary(pxy)+w(i)*x(i)*y(*,*,i)			
	endfor
endif else begin
	; test if weights defined
	if not (ws(0) eq 3) then begin
		w=x
		w(*,*,*)=1.
	endif
	; generate moments
	p=total(w,3)
	px=total(w*x,3)
	py=total(w*y,3)
	pxy=total(w*x*y,3)
	px2=total(w*x*x,3)
end

; compute coefficients
b=(pxy*p-px*py)/(px2*p-px^2)
a=py/p-b*px/p

; get chisquare
if (xs(0) eq 1) then begin
	chi2=fltarr(ys(1),ys(2))
	for i=0,n-1 do chi2=temporary(chi2)+w(i)*(y(*,*,i)-b*x(i)-a)^2
endif else chi2=total(w*(y-ima2cube(b,n)*x-ima2cube(a,n))^2,3)

; mean error
m=sqrt(chi2/(n-2.0))

; error of coefficients
ma=m*sqrt(px2/(p*px2-px^2))
mb=m*sqrt(p/(p*px2-px^2))

; return cube
return,[[[b]],[[a]],[[mb]],[[ma]],[[chi2]]]
end