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

FUNCTION evalcenpsfsub,P
	common _cenpsfsub,img,psf,count,silentflag
	; P is [ratio, background]
	
	diff=img - (psf/p[0]) - P[1]
	
	score=total(abs(diff))
	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: cenpsfsub
; PURPOSE:
;	centroiding psf subtraction routine. 
; NOTES:
;	Registers images using the
;	subreg procedure, then uses Powell minimization to find the bias offset
;	and scale factors which result in the best PSF subtraction, as measured
;	by minimizing the sum of the absolute value of the difference image.
;	
; INPUTS:
; 	im1,im2		2D images
; 	
; KEYWORDS:
; 	silent		supress printed output
; 	method		Registration method. See subreg.pro for documentation
; 	interp_type	interpolation method; see StarFinder image_shift.pro
; 	boxsize		radius of square box around star to be used
; 	ftol		tolerance for POWELL optimization
; OUTPUTS:
; 	output		subtracted overlap region and originals
; 		output[*,*,0] = subtracted (  = im1 - scale&shift(im2))
; 		output[*,*,1] = image 1
; 		output[*,*,2] = scaled and shifted image 2
; 	
; 	params		a float array containing the parameters of the subtraction fit.
;
; HISTORY:
; 	August 2001		Marshall Perrin
;
;-

PRO cenpsfsub,im1,im2,output,silent=silent,ftol=ftol,params=params,$
	interp_type=interp_type,boxsize=boxsize,method=method
	common _cenpsfsub,img,psf,count,silentflag

	if not(keyword_set(ftol)) then ftol=0.001
	if not(keyword_set(boxsize)) then boxsize=60
	if keyword_set(silent) then silentflag=1 
	sz=size(im1)

	subreg,im1,im2,dxy,method=method ; get initial offsets
	c_iter=0
	repeat begin
		im2s=image_shift(im2,dxy[0],dxy[1],interp_type=interp_type,data)
		subreg,im1,im2s,shifts,/silent,method=method,boxsize=boxsize; calculate the new relative shifts
		dxy=dxy+shifts  		; new overall shifts for next round
		error=sqrt(total(shifts^2))
		c_iter=c_iter+1
		print,"Iter "+strc(c_iter)+", shifts="+strc(shifts[0])+", "+strc(shifts[1])+"      Error: "+strc(error)
	endrep until (error lt ftol)

	message,/info, "Image registration complete"

	print,"Image shift found: "+printcoo(dxy[0],dxy[1])

	; now that we have the images aligned, use some form of minimization 
	; to get the ratio and offset right.
	;
	; do aperture photometry to get rough guesses at the flux ratio
	findmaxstar,im1,x1,y1
	aper, im1, [x1], [y1], flux_img,errap,sky,skyerr,1.0,[10],[25,35],[0,0],/flux,/silent
	aper, im2s, [x1], [y1], flux_psf,errap,sky,skyerr,1.0,[10],[25,35],[0,0],/flux,/silent
	P=[flux_psf[0]/flux_img[0],0]
	print,"Initial Fluxes are: TARGET "+strc(flux_img)+", PSF "+strc(flux_psf)+"  Ratio="+strc(flux_psf[0]/flux_img[0])

	img=im1[0>x1-boxsize:x1+boxsize< sz[1],0 >y1-boxsize:y1+boxsize < sz[2]]
	psf=im2s[0>x1-boxsize:x1+boxsize< sz[1],0 >y1-boxsize:y1+boxsize < sz[2]]
	count=0
	
if not(keyword_set(ftol)) then begin
   ftol = 1.0e-12
endif

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


   POWELL, P, xi, ftol, fmin, 'evalcenpsfsub',iter=iter,/double
   PRINT, 'Solution point, found after '+strc(iter)+" iterations: "
	print,"		: ",P[0],P[1],dxy[0],dxy[1]
   ; Print the value at the solution point:
   PRINT, 'Value at solution point: ', fmin
	params=[p,dxy]

	if arg_present(output) then begin
		sub=img - (psf/p[0] + P[1])
		output=fltarr(2*boxsize+1,2*boxsize+1,3)

		output(*,*,0)=sub
		output(*,*,1)=img
		output(*,*,2)=psf/p[0]+P[1]
	endif

end