; the following is just a helper function; scroll down for the main body of fitplane

FUNCTION evalfitplane,P
	common _fitplane,img,psf,count,silentflag,x,y
	; P is [ratio, background]
	
	diff=img - (P[0]+ x*P[1] #y*P[2])
	
	score=total(abs(diff),/nan)
	count=count+1;
	if (count mod 100) eq 0 then  $
		if not(keyword_set(silentflag)) then print,count,score,P

;	message,/info,"Score is "+strc(score)
	return,score


end

;+
;
; FUNCTION: fitplane
; PURPOSE:
;
;		Fits a plane to an image
;	
; INPUTS:
; 	im1		2D images
; 	
; KEYWORDS:
; 	silent		supress printed output
; 	method		Registration method. See subreg.pro for documentation
; 	interp_type	interpolation method; see StarFinder image_shift.pro
; 	ftol		tolerance for POWELL optimization
; OUTPUTS:
; 	
; 	params		a float array containing the parameters of the subtraction fit.
;	scale		divide im2 by this for best subtraction.
;
; HISTORY:
; 	August 2001		Marshall Perrin
;
;-

PRO fitplane,im1o, output,silent=silent,ftol=ftol,scale=scale,$
	interp_type=interp_type,method=method,subtracted=sub
	common _fitplane,img,psf,count,silentflag,x,y

	if not(keyword_set(ftol)) then ftol=0.0001
	if keyword_set(silent) then silentflag=1 
	sz=size(im1o)

message,/info,"Finding optimum subtraction"

	sz = size(im1o)

	x = findgen(sz[1])
	y = findgen(sz[2])

count=0

	P=[1.0,1.0,1.0]

img = im1o

   xi = fltarr(3,3)
   xi[0,0]=.10 & xi[1,1]=.10 & xi[2,2]=.1; & xi[3,3]=1.


   POWELL, P, xi, ftol, fmin, 'evalfitplane',iter=iter,/double

   PRINT, 'Value at solution point: ', fmin
   print, "Coefficients: ",p
   scale=p

   plane = P[0]+P[1]*x # y*P[2]

	if arg_present(output) then begin
		sub=img - (plane)
		sz = size(sub)
		output = [[[im1o]],[[plane]],[[sub]]]
	endif

end