pro sigma_clean, psf, im_resid, im_clean, SIGMA_IMAGE=im_sigma,  $
					  SKY_NOISE=sky_noise, $
				DET_GAIN=dgain, CFACTOR=factor, SMOOTH=smf, $
				POST_PROCESS=post, FMIN=fmin, FLUX=flux
;Frank Varosi NASA/GSFC 1992.
;(still in code development stage...)
; using W.C.Keel, Proc.Astro.Soc.Pacific, July 1991.

  common sigma_clean, sigma, px,py

	sd = size( im_resid )
	sp = size( psf )

	if N_elements( im_clean ) NE sd(sd(0)+2) then begin

		im_clean = fltarr( sd(1), sd(2) )
		if N_elements( sky_noise ) NE 1 then sky_noise=1

		if N_elements( im_sigma ) EQ sd(sd(0)+2) then sigma = im_sigma $
		 else if keyword_set( dgain ) then $
			sigma = sqrt( sky_noise^2 + (im_resid>0)/dgain ) $
		  else	sigma = replicate( sky_noise, sd(1), sd(2) )
	   endif

	max_snr = max( filter_image( im_resid/sigma, SM=smf,/ALL ), Lmax )
	px = Lmax MOD sd(1)
	py = Lmax / sd(1)

	if N_elements( factor ) NE 1 then factor=0.01
	if N_elements( fmin ) NE 1 then fmin=1.e27

	flux = im_resid(px,py) * factor
	fmin = fmin < abs( flux )

	im_clean(px,py) = im_clean(px,py) + flux

	sph = sp(1:2)/2
	b = [px,py] - sph
	L = [px,py] + sph

	im_resid(b(0),b(1)) = im_resid(b(0):L(0),b(1):L(1)) - flux * psf
return
end