1 files changed, 120 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitsky/apmean.x b/noao/digiphot/apphot/fitsky/apmean.x
new file mode 100644
index 00000000..f134f02c
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apmean.x
@@ -0,0 +1,120 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+# AP_MEAN -- Procedure to calculate the mean of the sky pixel array.
+
+int procedure ap_mean (skypix, coords, wgt, index, nskypix, snx, sny, losigma,
+	hisigma, rgrow, maxiter, sky_mean, sky_sigma, sky_skew, nsky,
+	nsky_reject)
+
+real	skypix[ARB]		# unsorted array of skypixels
+int	coords[ARB]		# coordinate array for region growing
+real	wgt[ARB]		# the weight array for rejection
+int	index[ARB]		# sorted array of indices
+int	nskypix			# total number of sky pixels
+int	snx, sny		# dimensions of the sky subraster
+real	losigma, hisigma	# number of sky_sigma for rejection
+real	rgrow			# radius of region growing
+int	maxiter			# maximum number of cycles of rejection
+real	sky_mean		# computed sky value
+real	sky_sigma		# the computed sky sigma
+real	sky_skew		# skewness of sky distribution
+int	nsky			# the number of sky pixels used
+int	nsky_reject		# the number of sky pixels rejected
+
+double  dsky, sumpx, sumsqpx, sumcbpx
+int	i, j, nreject
+real	sky_zero, dmin, dmax, locut, hicut
+int	ap_grow_regions()
+real	ap_asumr()
+
+begin
+	# Intialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mean = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Compute the mean, sigma and skew.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+
+	# Decide whether to do the rejection cycle.
+	if (maxiter < 1 || (IS_INDEFR(losigma) && IS_INDEFR(hisigma)) ||
+	    sky_sigma <= 0.0)
+	    return (AP_OK)
+
+	# Reject points within k1 * sky_sigma of the median.
+	do i = 1, maxiter {
+
+	    # Compute the new rejection limits.
+	    if (IS_INDEFR(losigma))
+		locut = max (-MAX_REAL, dmin)
+	    else if (i == 1)
+		locut = sky_mean -  min (sky_mean - dmin, dmax - sky_mean,
+		    losigma * sky_sigma)
+	    else
+	        locut = sky_mean - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = min (MAX_REAL, dmax)
+	    else if (i == 1)
+		hicut = sky_mean + min (dmax - sky_mean, sky_mean - dmin,
+		    hisigma * sky_sigma)
+	    else
+	        hicut = sky_mean + hisigma * sky_sigma
+
+	    nreject = 0
+	    do j = 1, nskypix {
+	        if (wgt[index[j]] <= 0.0)
+		    next
+	        if (skypix[index[j]] < locut || skypix[index[j]] > hicut) {
+		    if (rgrow > 0.0) {
+		        nreject = ap_grow_regions (skypix, coords, wgt,
+			    nskypix, sky_zero, index[j], snx, sny, rgrow,
+			    sumpx, sumsqpx, sumcbpx)
+		    } else {
+			dsky = (skypix[index[j]] - sky_zero)
+		        sumpx = sumpx - dsky
+		        sumsqpx = sumsqpx - dsky ** 2
+		        sumcbpx = sumcbpx - dsky ** 3
+		        wgt[index[j]] = 0.0 
+		        nreject = nreject + 1
+		    }
+	        }
+	    }
+
+	    # Test the number of rejected pixels.
+	    if (nreject <= 0)
+		break
+	    nsky_reject = nsky_reject + nreject
+
+	    # Test that some pixels actually remain.
+	    nsky = nskypix - nsky_reject
+	    if (nsky <= 0)
+		break
+
+	    # Recompute the mean, sigma and skew.
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	    if (sky_sigma <= 0.0)
+		break
+	    sky_mean = max (dmin, min (sky_mean, dmax))
+	}
+
+	# Return an appropriate error code.
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mean = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end