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

FUNCTION evaloptsub,P
	common _optsub,img,psf,count,silentflag
	; P is [ratio, background]
	
	diff=img - (psf/p[0]) 
	
	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: optsub
; PURPOSE:
;
; 	PSF Subtraction optimizer. Given two images already aligned, 
; 	what scale factor optimizes their subtraction?
;	
; 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.
;	scale		divide im2 by this for best subtraction.
;
; HISTORY:
; 	August 2001		Marshall Perrin
;
;-

PRO optsub,im1o,im2o,output,silent=silent,ftol=ftol,scale=scale,$
	interp_type=interp_type,boxsize=boxsize,method=method,all=all
	common _optsub,img,psf,count,silentflag

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

message,/info,"Finding optimum subtraction, using boxsize="+strc(boxsize)

	; no registration done here! Just the scaling.

	; do aperture photometry to get rough guesses at the flux ratio
	im1 = im1o
	im2 = im2o
	wnan = where(~ finite(im1),nanct)
	if nanct gt 0 then im1[wnan] = median(im1) 
	wnan = where(~ finite(im2),nanct)
	if nanct gt 0 then im2[wnan] = median(im2) 
	
	findmaxstar,im1,x1,y1
	aper, im1, [x1], [y1], flux_img,errap,sky,skyerr,1.0,[10],[25,35],[0,0],/flux,/silent
	aper, im2, [x1], [y1], flux_psf,errap,sky,skyerr,1.0,[10],[25,35],[0,0],/flux,/silent
	P=[flux_psf[0]/flux_img[0],0]
	if total(finite(P)) ne 2 then begin
		message,/info,"Aperture photometry gave NaNs; using initial guess 1.0 instead."
		P = [1.0,0.0]
	endif
	print,"Initial Fluxes are: TARGET "+strc(flux_img)+", PSF "+strc(flux_psf)+"  Ratio="+strc(flux_psf[0]/flux_img[0])

	if (keyword_set(all)) then begin
		img = im1
		psf=im2
	endif else begin
		img=im1[0>(x1-boxsize):(x1+boxsize)< (sz[1]-1),0 >(y1-boxsize):(y1+boxsize) < (sz[2]-1)]
		psf=im2[0>(x1-boxsize):(x1+boxsize)< (sz[1]-1),0 >(y1-boxsize):(y1+boxsize) < (sz[2]-1)]
	endelse
	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.


	; here we ignore the second parameter to P, but Powell seems to want at
	; least two options available. It doesn't do anything though...

   POWELL, P, xi, ftol, fmin, 'evaloptsub',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
   print, "Scale factor: ",p[0]
   scale=p[0]

	if arg_present(output) then begin
		sub=img - (psf/p[0] )
		sz = size(sub)
		output=fltarr(sz[1],sz[2],3)

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

end