;
; fit a plane to xyz data
; A*x + B*y + c*z + D = 0
; or z=-(A/C *x + B/C * y + D/C)
; so params 
; A*xy[0,*] + B*x[1,*] + C*x[2,*] +D = 0
; a[0]=A/C
; a[1]=B/C
; a[2]=D/C
; the partial derivatives are:
; -[x[0,*],x[1,*],dblarr(npnts)+1d]
;
pro  fitplanefunc,x,a,f,pderr

;	print,a
	npnts=n_elements(x)/2l
	xy=reform(x,2,npnts)
	f=-(reform(a[0]*xy[0,*] + a[1]*xy[1,*]) + a[2])
	if n_params() ge 4 then begin
		pderr= -[[reform(xy[0,*])],[reform(xy[1,*])],[dblarr(npnts)+1d]]
	endif
	return
end