Initial commit

1 files changed, 177 insertions, 0 deletions

diff --git a/noao/imred/ccdred/src/combine/icmm.gx b/noao/imred/ccdred/src/combine/icmm.gx
new file mode 100644
index 00000000..90837ae5
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/icmm.gx
@@ -0,0 +1,177 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+$for (sr)
+# IC_MM --  Reject a specified number of high and low pixels
+
+procedure ic_mm$t (d, m, n, npts)
+
+pointer	d[ARB]		# Data pointers
+pointer	m[ARB]		# Image ID pointers
+int	n[npts]			# Number of good pixels
+int	npts			# Number of output points per line
+
+int	n1, ncombine, npairs, nlow, nhigh, np
+int	i, i1, j, jmax, jmin
+pointer	k, kmax, kmin
+PIXEL	d1, d2, dmin, dmax
+
+include	"../icombine.com"
+
+begin
+	if (dflag == D_NONE)
+	    return
+
+	if (dflag == D_ALL) {
+	    n1 = n[1]
+	    nlow = flow * n1 + 0.001
+	    nhigh = fhigh * n1 + 0.001
+	    ncombine = n1 - nlow -  nhigh
+	    npairs = min (nlow, nhigh)
+	    nlow = nlow - npairs
+	    nhigh = nhigh - npairs
+	}
+
+	do i = 1, npts {
+	    i1 = i - 1
+	    n1 = n[i]
+	    if (dflag == D_MIX) {
+		nlow = flow * n1 + 0.001
+		nhigh = fhigh * n1 + 0.001
+		ncombine = max (ncombine, n1 - nlow - nhigh)
+		npairs = min (nlow, nhigh)
+		nlow = nlow - npairs
+		nhigh = nhigh - npairs
+	    }
+
+	    # Reject the npairs low and high points.
+	    do np = 1, npairs {
+		k = d[1] + i1
+		$if (datatype == x)
+		d1 = abs (Mem$t[k])
+		$else
+		d1 = Mem$t[k]
+		$endif
+		dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k
+		do j = 2, n1 {
+		    d2 = d1
+		    k = d[j] + i1
+		    $if (datatype == x)
+		    d1 = abs (Mem$t[k])
+		    $else
+		    d1 = Mem$t[k]
+		    $endif
+		    if (d1 > dmax) {
+			dmax = d1; jmax = j; kmax = k
+		    } else if (d1 < dmin) {
+			dmin = d1; jmin = j; kmin = k
+		    }
+		}
+		j = n1 - 1
+		if (keepids) {
+		    if (jmax < j) {
+			if (jmin != j) {
+			    Mem$t[kmax] = d2
+			    Memi[m[jmax]+i1] = Memi[m[j]+i1]
+			} else {
+			    Mem$t[kmax] = d1
+			    Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+			}
+		    }
+		    if (jmin < j) {
+			if (jmax != n1) {
+			    Mem$t[kmin] = d1
+			    Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+			} else {
+			    Mem$t[kmin] = d2
+			    Memi[m[jmin]+i1] = Memi[m[j]+i1]
+			}
+		    }
+		} else {
+		    if (jmax < j) {
+			if (jmin != j)
+			    Mem$t[kmax] = d2
+			else
+			    Mem$t[kmax] = d1
+		    }
+		    if (jmin < j) {
+			if (jmax != n1)
+			    Mem$t[kmin] = d1
+			else
+			    Mem$t[kmin] = d2
+		    }
+		}
+		n1 = n1 - 2
+	    }
+
+	    # Reject the excess low points.
+	    do np = 1, nlow {
+		k = d[1] + i1
+		$if (datatype == x)
+		d1 = abs (Mem$t[k])
+		$else
+		d1 = Mem$t[k]
+		$endif
+		dmin = d1; jmin = 1; kmin = k
+		do j = 2, n1 {
+		    k = d[j] + i1
+		    $if (datatype == x)
+		    d1 = abs (Mem$t[k])
+		    $else
+		    d1 = Mem$t[k]
+		    $endif
+		    if (d1 < dmin) {
+			dmin = d1; jmin = j; kmin = k
+		    }
+		}
+		if (keepids) {
+		    if (jmin < n1) {
+			Mem$t[kmin] = d1
+			Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+		    }
+		} else {
+		    if (jmin < n1)
+			Mem$t[kmin] = d1
+		}
+		n1 = n1 - 1 
+	    }
+
+	    # Reject the excess high points.
+	    do np = 1, nhigh {
+		k = d[1] + i1
+		$if (datatype == x)
+		d1 = abs (Mem$t[k])
+		$else
+		d1 = Mem$t[k]
+		$endif
+		dmax = d1; jmax = 1; kmax = k
+		do j = 2, n1 {
+		    k = d[j] + i1
+		    $if (datatype == x)
+		    d1 = abs (Mem$t[k])
+		    $else
+		    d1 = Mem$t[k]
+		    $endif
+		    if (d1 > dmax) {
+			dmax = d1; jmax = j; kmax = k
+		    }
+		}
+		if (keepids) {
+		    if (jmax < n1) {
+			Mem$t[kmax] = d1
+			Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+		    }
+		} else {
+		    if (jmax < n1)
+			Mem$t[kmax] = d1
+		}
+		n1 = n1 - 1 
+	    }
+	    n[i] = n1
+	}
+
+	if (dflag == D_ALL && npairs + nlow + nhigh > 0)
+		dflag = D_MIX
+end
+$endfor