;
; use sfit to compute the plane directly
;
; sigma=fitplanefast(xyz,coef,zfit=zfit,usemedian=usemedian)
;
;ARGS:
;  xyz[3,n]: double   points to fit
;KEYWORDS:
;usemed  :           if set then use median point to compute D
;                    default is to use the average
;zfit[n]:  double    return the z value from the fit
;RETURNS:
; sigma : double    sigma of  dif=xyz - yfit
; coef[4]: double   A,B,C,D  from fit: Ax  + By + Cz + D  = 0
;                   [coef[0:2] is normal to plane
;                    coef[3]   is distance plane to origin 
;KEYWORDS:
;zfit[n] : double    the fitted value for z

function fitplanefast,xyz,coef,zfit=zfit,usemedian=usemedian
;  

	deg=1
	zfit=sfit(xyz,deg,/irreg,kx=kx,/max_degree)
;
; kx = [D/C , B/C, A/C]
; to get C assume mean value of  xyz gives distance D to origin.
;
	a=size(xyz)
	npnts=a[2]
	pavg=(keyword_set(usemedian))?median(xyz,dimension=2):total(xyz,2)/npnts
	D=sqrt(total(pavg*pavg))
    if (kx[0] lt 0) then D*=-1d  ; this may not always  be correct if C < 0.
	C=1d/kx[0] * D
	coef=[kx[2],kx[1],1d,D/C]*C
	a=rms(zfit-reform(xyz[2,*]),/quiet)
	sigma=a[1]
	return,sigma
end