;+
; NAME:
;	psf_merge
; PURPOSE:
;	Merge two psf matrices together by smoothly splicing/glueing.
; CALLING EXAMPLE:
;	psf_combo = psf_merge( psf_inner, psf_outer, RADIUS=7 )
; INPUTS:
;	psf_inner = will be used inside (at the peak).
;	psf_outer = will be spliced to outside (this one must be centered)
; KEYWORDS:
;	LEVEL_SPLICE = magnitude at which psf_outer is spliced with psf_inner,
;			overrides and determines RADIUS_SPLICE.
;	RADIUS_SPLICE = radius at which splicing occurs (default = 3*FWHM).
;	FACTOR_SPLICE = control steepness of glueing function (default = 2).
; OUTPUT:
;	Function returns matrix of new spliced PSF.
; EXTERNAL CALLS:
;	function dist
;	function FullWid_HalfMax (if RADIUS_SPLICE not specified)
; PROCEDURE:
;	Use the function 1/(1+exp(x)) to smoothly glue psf_outer to psf_inner.
; HISTORY:
;	Written, Frank Varosi NASA/GSFC 1993.
;-

function psf_merge, psf_inner, psf_outer, FACTOR_SPLICE=splice_Fac, $
					  RADIUS_SPLICE=splice_Rad, $
					  LEVEL_SPLICE=splice_Lev
	si = size( psf_inner )
	so = size( psf_outer )
	cm = so(1:2)/2
	dim = shift( dist( so(1),so(2) ), cm(0),cm(1) )

	psfex = fltarr( so(1), so(2) )
	Loc = ( so/2 - si/2 ) > 0	;center PSF in new array,
	s = (si/2 - so/2) > 0	   ;handle all cases: smaller or bigger
	L = (s + so-1) < (si-1)
	psfex( Loc(1), Loc(2) ) = psf_inner( s(1):L(1), s(2):L(2) )

	if N_elements( splice_Lev ) EQ 1 then begin
		message,"using splice Level = " + strtrim(splice_Lev,2),/INFO
		w = where( psfex LE splice_Lev, nw )
		if (nw GT 9) then begin
			splice_Rad = fix( min( dim(w) ) )
		    message,"found splice Radius = "+strtrim(splice_Rad,2),/INFO
		  endif else message,"cannot find splice Radius",/INFO
	   endif

	if N_elements( splice_Fac ) NE 1 then splice_Fac=2
	splice_Fac = ( splice_Fac > 0.1 ) < 10

	if N_elements( splice_Rad ) NE 1 then $
			splice_Rad = 3 * FullWid_HalfMax( psf_outer,/AVERAGE )
	splice_Rad = ( splice_Rad > 2 ) < ( min( cm )-1 )

	dim = (dim - splice_Rad) * splice_Fac
	alpha = fltarr( so(1), so(2) )
	c = [ cm-splice_Rad, cm+splice_Rad ]
	alpha( c(0):c(2), c(1):c(3) ) = 1	
	w = where( abs( dim ) LE 21 )
	alpha(w) = 1/( 1 + exp( dim(w) ) )

	if !DEBUG GT 2 then stop

return, alpha*psfex + (1-alpha)*psf_outer
end