1 files changed, 160 insertions, 0 deletions
diff --git a/noao/onedspec/scombine/icstat.x b/noao/onedspec/scombine/icstat.x
new file mode 100644
index 00000000..3fce4165
--- /dev/null
+++ b/noao/onedspec/scombine/icstat.x
@@ -0,0 +1,160 @@
+include	<smw.h>
+include	"icombine.h"
+
+
+# IC_STATR -- Compute image statistics within spectrum.
+
+procedure ic_statr (sh, lflag, rg, domode, domedian, domean, mode, median, mean)
+
+pointer	sh			 # Spectrum structure
+int	lflag			 # Data flag
+pointer	rg			 # Wavelength ranges
+bool	domode, domedian, domean # Statistics to compute
+real	mode, median, mean	 # Statistics
+
+int	i, n, npts
+real	a, w
+pointer	sp, data, dp, lp, mp
+real	ic_moder(), asumr()
+bool	ic_wisinrange()
+double	shdr_lw()
+
+include	"icombine.com"
+
+begin
+	mp = SX(sh)
+	lp = SY(sh)
+	npts = SN(sh)
+
+	call smark (sp)
+	call salloc (data, npts, TY_REAL)
+
+	dp = data
+	if (lflag == D_ALL && rg == NULL) {
+	    if (dothresh) {
+		do i = 1, npts {
+		    a = Memr[lp]
+		    if (a >= lthresh && a <= hthresh) {
+			Memr[dp] = a
+			dp = dp + 1
+		    }
+		    lp = lp + 1
+		}
+	    } else {
+		do i = 1, npts {
+		    Memr[dp] = Memr[lp]
+		    dp = dp + 1
+		    lp = lp + 1
+		}
+	    }
+	} else if (lflag == D_MIX || rg != NULL) {
+	    if (dothresh) {
+		do i = 1, npts {
+		    if (Memr[mp] == 0) {
+			a = Memr[lp]
+			if (a >= lthresh && a <= hthresh) {
+			    w = shdr_lw (sh, double (i))
+			    if (ic_wisinrange (rg, w)) {
+				Memr[dp] = a
+				dp = dp + 1
+			    }
+			}
+		    }
+		    mp = mp + 1
+		    lp = lp + 1
+		}
+	    } else {
+		do i = 1, npts {
+		    if (Memr[mp] == 0) {
+			w = shdr_lw (sh, double (i))
+			if (ic_wisinrange (rg, w)) {
+			    Memr[dp] = Memr[lp]
+			    dp = dp + 1
+			}
+		    }
+		    mp = mp + 1
+		    lp = lp + 1
+		}
+	    }
+	}
+
+	n = dp - data
+	if (n > 0) {
+	    # Compute only statistics needed.
+	    if (domode || domedian) {
+		call asrtr (Memr[data], Memr[data], n)
+		mode = ic_moder (Memr[data], n)
+		median = Memr[data+n/2-1]
+	    }
+	    if (domean)
+		mean = asumr (Memr[data], n) / n
+	} else {
+	    mode = INDEF
+	    median = INDEF
+	    mean = INDEF
+	}
+
+	call sfree (sp)
+end
+
+
+define	NMIN	10	# Minimum number of pixels for mode calculation
+define	ZRANGE	0.8	# Fraction of pixels about median to use
+define	ZSTEP	0.01	# Step size for search for mode
+define	ZBIN	0.1	# Bin size for mode.
+
+# IC_MODE -- Compute mode of an array.  The mode is found by binning
+# with a bin size based on the data range over a fraction of the
+# pixels about the median and a bin step which may be smaller than the
+# bin size.  If there are too few points the median is returned.
+# The input array must be sorted.
+
+real procedure ic_moder (a, n)
+
+real	a[n]			# Data array
+int	n			# Number of points
+
+int	i, j, k, nmax
+real	z1, z2, zstep, zbin
+real	mode
+bool	fp_equalr()
+
+begin
+	if (n < NMIN)
+	    return (a[n/2])
+
+	# Compute the mode.  The array must be sorted.  Consider a
+	# range of values about the median point.  Use a bin size which
+	# is ZBIN of the range.  Step the bin limits in ZSTEP fraction of
+	# the bin size.
+
+	i = 1 + n * (1. - ZRANGE) / 2.
+	j = 1 + n * (1. + ZRANGE) / 2.
+	z1 = a[i]
+	z2 = a[j]
+	if (fp_equalr (z1, z2)) {
+	    mode = z1
+	    return (mode)
+	}
+
+	zstep = ZSTEP * (z2 - z1)
+	zbin = ZBIN * (z2 - z1)
+
+	z1 = z1 - zstep
+	k = i
+	nmax = 0
+	repeat {
+	    z1 = z1 + zstep
+	    z2 = z1 + zbin
+	    for (; i < j && a[i] < z1; i=i+1)
+		;
+	    for (; k < j && a[k] < z2; k=k+1)
+		;
+	    if (k - i > nmax) {
+	        nmax = k - i
+	        mode = a[(i+k)/2]
+	    }
+	} until (k >= j)
+
+	return (mode)
+end