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

FUNCTION evaloptsubwithoffset,P
	common _optsub,img,psf,count,silentflag
	; P is [ratio, background]
	
	diff=img - (psf/p[0]) - P[1] 
	
	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: optsubwithoffset
; PURPOSE:
; 	PSF Subtraction optimizer. Given two images already aligned, 
; 	what scale factor optimizes their subtraction?
;
; 	This version also solves for a constant shift between the two images.
;	
; 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
; 	/show		display on screen the two images
; 	/all		use the whole image, not a subimage, for the 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.
;	offset=		Add this to im2 for best subtraction.
;
; HISTORY:
; 	August 2001		Marshall Perrin
;
;	2006-05-25		Added /show and /all
;-

PRO optsubwithoffset,im1o,im2o,output,silent=silent,ftol=ftol,scale=scale,$
	interp_type=interp_type,boxsize=boxsize,method=method,offset=offset,$
	badmask = badmask,show=show,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)

	; no registration done here! Just the scaling.

if ~(keyword_set(silent)) then message,/info,"Finding optimum subtraction, using boxsize="+strc(boxsize)
	; do aperture photometry to get rough guesses at the flux ratio
	im1 = im1o
	im2 = im2o

	; remove NANs before APER, since it chokes on them.
	wnan1 = where(~ finite(im1),nanct1)
	if nanct1 gt 0 then im1[wnan1] = median(im1) 
	wnan2 = where(~ finite(im2),nanct2)
	if nanct2 gt 0 then im2[wnan2] = median(im2) 
	
	findmaxstar,im1,x1,y1,silent=silent
	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
	if ~(keyword_set(silent)) then 	message,/info,"Aperture photometry gave NaNs; using initial guess 1.0 instead."
		P = [1.0,0.0]
	endif
	if ~(keyword_set(silent)) then print,"Initial Fluxes are: TARGET "+strc(flux_img)+", PSF "+strc(flux_psf)+"  Ratio="+strc(flux_psf[0]/flux_img[0])

	; now put NANs back in to ignore those pixels.
	if nanct1 gt 0 then im1[wnan1] =!values.f_nan
	if nanct2 gt 0 then im2[wnan2] =!values.f_nan

if keyword_set(all) then begin
	; now using entire image??!
	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

	if keyword_set(show) then begin
		pm = !p.multi
		!p.multi=[0,2,1]
		imdisp,alogscale(img),title="Image 1",/axis
		imdisp,alogscale(psf),title="Image 2",/axis
		!p.multi=pm
	endif

	if keyword_set(badmask) then begin
		wbad = where(badmask eq 0)
		img[wbad] = !values.f_nan
		psf[wbad] = !values.f_nan
	endif

	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, 'evaloptsubwithoffset',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:
	if ~(keyword_set(silent)) then  begin
	   PRINT, 'Value at solution point: ', fmin
	   print, "Scale factor: ",p[0]
	   print, "Offset :       ",p[1]
   endif
   scale=p[0]
   offset=P[1]

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

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

end