procedure getsky (image)

string	image	=""	{prompt="image name"}
real	smean	=0.	{prompt="average sky background level"}
real	smode	=0.	{prompt="modal sky background level"}
real	smidpt	=0.	{prompt="median sky background level"}
real	srms	=0.	{prompt="rms (sigma) of sky background"}
int	niter	=3	{prompt="number of rejection iterations"}
bool	verbose	=no	{prompt="print results to screen?"}

# @(#) task getsky                   Author: R.A. Jansen -- Jan 29, 2002
# @(#)
# @(#) task to measure the aproximate sky background level in an image,
# @(#) assuming that the background is fairly flat (has no large scale
# @(#) gradient or prominent higher order structure).  Uses: 'imgstat'.

begin
	string	img
	real	zavg, zrms, zmin, zmax, z1, z2
	int	nn, i

	img	= image
	nn	= niter

#       Check existence of input image...
	unlearn chkimg; chkimg (img, "access", verbose=yes)
	if (!chkimg.ok) { return }

#	Opening message...
        if (verbose) { print ("GETSKY: "//img) }

#	Measure raw image statistics...
	imgstat (img, domode=no, verbose=no)
	zavg = imgstat.mean
	zrms = imgstat.stddev
	z1   = zavg - 4*zrms
	z2   = zavg + 4*zrms

#	Measure image statistics after niter rejections iterations...
	for ( i = 1 ; i <= nn ; i = i + 1 ) {
	    imgstat (img, lower=z1, upper=z2, domode=yes, binwidth=0.1, verbose=no)
	    zavg = imgstat.mean
	    zrms = imgstat.stddev
	    zmin = imgstat.min
	    zmax = imgstat.max
	    z1   = zavg - (3./sqrt(i-0.5))*zrms
	    z2   = zavg + (3./sqrt(i-0.5))*zrms

	    # Print results...
	    if (verbose) {
		printf ("Iter= %2g, Mean = %.4f, Median= %.4f, Mode = %.4f, RMS = %.4f\n",
			i, zavg, imgstat.midpt, imgstat.rmode, zrms)
	    }
	}

#       Return task parameters...
	smean	= imgstat.mean
	smidpt	= imgstat.midpt
	smode	= imgstat.rmode
	srms	= imgstat.stddev

end