;+ optpsfsub
;
;PRO optpsfsub,filename1,filename2,unsharp=unsharp,cbox=cbox,$
;	output=output,nosave=nosave,savecube=savecube,$
;	oversample=oversample
; PURPOSE:
;  Optimizing PSF subtraction routine.
;
; INPUTS: filename1,filename2	filenames of the two fits files to 
; 	subtract. Assumes suffix ".ks_mos.fits"
; OUTPUTS: 
; 	output		This keyword returns the subtracted image, or
; 				an image cube if /savecube is set.
; KEYWORDS:
; 	unsharp=	do unsharp masking with this radius
; 	cbox=		size of box around star to evaluate
; 	oversample=	do oversampling of this amount
; 	/nosave		don't save to disk, just keep in memory.
; 	/qzap		use qzap to kill CRs
; Interpolation:
; 	/fft 		use fft interpolation
; 	/cubic		use cubic splines for subpixel interpolation
; 	/interp		use bilinear interpolation
; 	/sinc		use sinc interpolation
; 	
; 	If multiple interpolation flags are passed, the highest one in the
; 	list above is honored.
; 	
; HISTORY: Began by MDP on 2001-07-09 as a wrapper for psfpowell
; 2001-08-22	Subpixel interpolation options added.
;
; 
;-


PRO optpsfsub,filename1,filename2,unsharp=unsharp,cbox=cbox,$
	output=output,nosave=nosave,$
	oversample=oversample,qzap=qzap,sinc=sinc,fft=fft,cubic=cubic,interp=interp

if n_params() lt 2 then begin
	print,"ERROR: incorrect number of arguments specified."
	print,"	optpsfsub,filename1,filename2 "
endif

img=readfits(filename1,head0)
psf=readfits(filename2,psfh)
;FIXME: add error checking here for nonexistant files.

if keyword_set(qzap) then begin
	print,"Removing cosmic rays with QZAP..."
	qzap,img,img2
	qzap,psf,psf2
	img=temporary(img2)
	psf=temporary(psf2)
endif

; always save cube if output to variable desired

print,"Now optimizing subtraction..."
psfpowell,img,psf,output,params=params,unsharp=unsharp,cbox=cbox,$
	oversample=oversample,sinc=sinc,fft=fft,interp=interp,cubic=cubic

paramtxt=["Amplitude Ratio","Background level","Dx","Dy","Iterations","Boxsize","Unsharp masking","Oversampling","Floating-point tolerance"]

sxaddhist,"<< OPTPSFSUB.PRO: "+systime()+ userid()+">>",head0
sxaddhist,"Performedoptimizing PSF subtraction on "+filename1,head0
sxaddhist,"  Using PSF star "+filename2,head0
sxaddhist,"  Subtraction Parameters: ",head0
sxaddhist,"    "+paramtxt+": "+strc(params),head0
if keyword_set(sint) then  begin
	sxaddhist,"    Sub-pixel sinc interpolation used.",head0
endif else if keyword_set(fft) then begin
	sxaddhist,"    sub-pixel fft interpolation used.",head0
endif else sxaddhist,"    No sub-pixel interp; ignore fractional part of dx & dy"

outfile="sub_"+filename1

if (not (keyword_set(nosave))) then begin
	print,"Now writing data to disk."
	if (filecheck(outfile)) then begin
			sxaddhist," Image Layers as follows:",head0
			sxaddhist,"   0=sub  1=original 2=psf",head0
			writefits,outfile,output,head0
	endif
endif

end