Initial commit

14 files changed, 6648 insertions, 0 deletions

diff --git a/noao/imred/ccdred/src/combine/generic/icaclip.x b/noao/imred/ccdred/src/combine/generic/icaclip.x
new file mode 100644
index 00000000..1530145c
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icaclip.x
@@ -0,0 +1,1102 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+define	MINCLIP		3	# Minimum number of images for this algorithm
+
+
+# IC_AAVSIGCLIP -- Reject pixels using an average sigma about the average
+# The average sigma is normalized by the expected poisson sigma.
+
+procedure ic_aavsigclips (d, m, n, scales, zeros, nimages, npts, average)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+real	d1, low, high, sum, a, s, s1, r, one
+data	one /1.0/
+pointer	sp, sums, resid, dp1, dp2, mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (sums, npts, TY_REAL)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Since the unweighted average is computed here possibly skip combining
+	if (dowts || combine != AVERAGE)
+	    docombine = true
+	else
+	    docombine = false
+
+	# Compute the unweighted average with the high and low rejected and
+	# the poisson scaled average sigma.  There must be at least three
+	# pixels at each point to define the average and contributions to
+	# the mean sigma.  Corrections for differences in the image
+	# scale factors are selected by the doscale1 flag.
+
+	nin = n[1]
+	s = 0.
+	n2 = 0
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (n1 < 3)
+		next
+
+	    # Unweighted average with the high and low rejected
+	    low = Mems[d[1]+k]
+	    high = Mems[d[2]+k]
+	    if (low > high) {
+		d1 = low
+		low = high
+		high = d1
+	    }
+	    sum = 0.
+	    do j = 3, n1 {
+	       d1 = Mems[d[j]+k]
+	       if (d1 < low) {
+		   sum = sum + low
+		   low = d1
+		} else if (d1 > high) {
+		    sum = sum + high
+		    high = d1
+		} else
+		    sum = sum + d1
+	    }
+	    a = sum / (n1 - 2)
+	    sum = sum + low + high
+
+	    # Poisson scaled sigma accumulation
+	    if (doscale1) {
+		do j = 1, n1 {
+		    dp1 = d[j] + k
+		    mp1 = m[j] + k
+
+		    d1 = Mems[dp1]
+		    l = Memi[mp1]
+		    s1 = max (one, (a + zeros[l]) /  scales[l])
+		    s = s + (d1 - a) ** 2 / s1
+		}
+	    } else {
+		s1 = max (one, a)
+		do j = 1, n1
+		    s = s + (Mems[d[j]+k] - a) ** 2 / s1
+	    }
+	    n2 = n2 + n1
+
+	    # Save the average and sum for later.
+	    average[i] = a
+	    Memr[sums+k] = sum
+	}
+
+	# Here is the final sigma.
+	if (n2 > 1)
+	    s = sqrt (s / (n2 - 1))
+
+	# Reject pixels and compute the final average (if needed).
+	# There must be at least three pixels at each point for rejection.
+	# Iteratively scale the mean sigma and reject pixels
+	# Compact the data and keep track of the image IDs if needed.
+
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    if (n1 <= max (2, maxkeep)) {
+		if (!docombine) {
+		    if (n1 == 0)
+			average[i] = blank
+		    else {
+			sum = Mems[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Mems[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    a = average[i]
+	    sum = Memr[sums+k]
+
+	    repeat {
+		n2 = n1
+		if (s > 0.) {
+		    if (doscale1) {
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    mp1 = m[j] + k
+
+			    d1 = Mems[dp1]
+			    l = Memi[mp1]
+			    s1 = s * sqrt (max (one, (a+zeros[l]) /  scales[l]))
+			    r = (d1 - a) / s1
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs(r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    mp2 = m[n1] + k
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = l
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    } else {
+			s1 = s * sqrt (max (one, a))
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    d1 = Mems[dp1]
+			    r = (d1 - a) / s1
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs(r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[n1] + k
+					l = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = l
+				    }
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		}
+		if (n1 > 1)
+		    a = sum / n1
+	    } until (n1 == n2 || n1 <= max (2, maxkeep))
+
+	    # If too many are rejected add some back in.
+	    # Pixels with equal residuals are added together.
+	    if (n1 < maxkeep) {
+		nk = maxkeep
+		if (doscale1) {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Mems[dp1]
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Mems[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		} else {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Mems[dp1]
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Mems[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		}
+		if (n1 > 1)
+		    a = sum / n1
+	    }
+
+	    # Save the average if needed.
+	    n[i] = n1
+	    if (!docombine) {
+		if (n1 > 0)
+		    average[i] = a
+		else
+		    average[i] = blank
+	    }
+	}
+
+	# Check if the data flag has to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# IC_MAVSIGCLIP -- Reject pixels using an average sigma about the median
+# The average sigma is normalized by the expected poisson sigma.
+
+procedure ic_mavsigclips (d, m, n, scales, zeros, nimages, npts, median)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+pointer	sp, resid, mp1, mp2
+real	med, low, high, r, s, s1, one
+data	one /1.0/
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Compute the poisson scaled average sigma about the median.
+	# There must be at least three pixels at each point to define
+	# the mean sigma.  Corrections for differences in the image
+	# scale factors are selected by the doscale1 flag.
+
+	s = 0.
+	n2 = 0
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (n1 < 3) {
+		if (n1 == 0)
+		    median[i] = blank
+		else if (n1 == 1)
+		    median[i] = Mems[d[1]+k]
+		else {
+		    low = Mems[d[1]+k]
+		    high = Mems[d[2]+k]
+		    median[i] = (low + high) / 2.
+		}
+		next
+	    }
+
+	    # Median
+	    n3 = 1 + n1 / 2
+	    if (mod (n1, 2) == 0) {
+		low = Mems[d[n3-1]+k]
+		high = Mems[d[n3]+k]
+		med = (low + high) / 2.
+	    } else
+		med = Mems[d[n3]+k]
+
+	    # Poisson scaled sigma accumulation
+	    if (doscale1) {
+		do j = 1, n1 {
+		    l = Memi[m[j]+k]
+		    s1 = max (one, (med + zeros[l]) /  scales[l])
+		    s = s + (Mems[d[j]+k] - med) ** 2 / s1
+		}
+	    } else {
+		s1 = max (one, med)
+		do j = 1, n1
+		    s = s + (Mems[d[j]+k] - med) ** 2 / s1
+	    }
+	    n2 = n2 + n1
+
+	    # Save the median for later.
+	    median[i] = med
+	}
+
+	# Here is the final sigma.
+	if (n2 > 1)
+	    s = sqrt (s / (n2 - 1))
+	else
+	    return
+
+	# Compute individual sigmas and iteratively clip.
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    if (n1 < max (3, maxkeep+1))
+		next
+	    nl = 1
+	    nh = n1
+	    med = median[i]
+
+	    repeat {
+		n2 = n1
+		n3 = nl + n1 / 2
+
+		if (n1 >= max (MINCLIP, maxkeep+1) && s > 0.) {
+		    if (doscale1) {
+			for (; nl <= n2; nl = nl + 1) {
+			    l = Memi[m[nl]+k]
+			    s1 = s * sqrt (max (one, (med+zeros[l])/scales[l]))
+			    r = (med - Mems[d[nl]+k]) / s1
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    l = Memi[m[nh]+k]
+			    s1 = s * sqrt (max (one, (med+zeros[l])/scales[l]))
+			    r = (Mems[d[nh]+k] - med) / s1
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    } else {
+			s1 = s * sqrt (max (one, med))
+			for (; nl <= n2; nl = nl + 1) {
+			    r = (med - Mems[d[nl]+k]) / s1
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    r = (Mems[d[nh]+k] - med) / s1
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    }
+
+		    # Recompute median
+		    if (n1 < n2) {
+			if (n1 > 0) {
+			    n3 = nl + n1 / 2
+			    if (mod (n1, 2) == 0) {
+				low = Mems[d[n3-1]+k]
+				high = Mems[d[n3]+k]
+				med = (low + high) / 2.
+			    } else
+				med = Mems[d[n3]+k]
+			} else
+			    med = blank
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    # If too many are rejected add some back in.
+	    # Pixels with equal residuals are added together.
+	    while (n1 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n1 = nh - nl + 1
+
+		# Recompute median
+		if (n1 < n2) {
+		    if (n1 > 0) {
+			n3 = nl + n1 / 2
+			if (mod (n1, 2) == 0) {
+			    low = Mems[d[n3-1]+k]
+			    high = Mems[d[n3]+k]
+			    med = (low + high) / 2.
+			} else
+			    med = Mems[d[n3]+k]
+		    } else
+			med = blank
+		}
+	    }
+
+	    # Only set median and reorder if needed
+	    n[i] = n1
+	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Mems[d[l]+k] = Mems[d[j]+k]
+			if (grow > 0) {
+			    mp1 = m[l] + k
+			    mp2 = m[j] + k
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+k] = Memi[m[j]+k]
+			j = j + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Mems[d[l]+k] = Mems[d[j]+k]
+			j = j + 1
+		    }
+		}
+	    }
+
+	    if (combine == MEDIAN)
+		median[i] = med
+	}
+
+	# Check if data flag needs to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag that the median is computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
+
+# IC_AAVSIGCLIP -- Reject pixels using an average sigma about the average
+# The average sigma is normalized by the expected poisson sigma.
+
+procedure ic_aavsigclipr (d, m, n, scales, zeros, nimages, npts, average)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+real	d1, low, high, sum, a, s, s1, r, one
+data	one /1.0/
+pointer	sp, sums, resid, dp1, dp2, mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (sums, npts, TY_REAL)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Since the unweighted average is computed here possibly skip combining
+	if (dowts || combine != AVERAGE)
+	    docombine = true
+	else
+	    docombine = false
+
+	# Compute the unweighted average with the high and low rejected and
+	# the poisson scaled average sigma.  There must be at least three
+	# pixels at each point to define the average and contributions to
+	# the mean sigma.  Corrections for differences in the image
+	# scale factors are selected by the doscale1 flag.
+
+	nin = n[1]
+	s = 0.
+	n2 = 0
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (n1 < 3)
+		next
+
+	    # Unweighted average with the high and low rejected
+	    low = Memr[d[1]+k]
+	    high = Memr[d[2]+k]
+	    if (low > high) {
+		d1 = low
+		low = high
+		high = d1
+	    }
+	    sum = 0.
+	    do j = 3, n1 {
+	       d1 = Memr[d[j]+k]
+	       if (d1 < low) {
+		   sum = sum + low
+		   low = d1
+		} else if (d1 > high) {
+		    sum = sum + high
+		    high = d1
+		} else
+		    sum = sum + d1
+	    }
+	    a = sum / (n1 - 2)
+	    sum = sum + low + high
+
+	    # Poisson scaled sigma accumulation
+	    if (doscale1) {
+		do j = 1, n1 {
+		    dp1 = d[j] + k
+		    mp1 = m[j] + k
+
+		    d1 = Memr[dp1]
+		    l = Memi[mp1]
+		    s1 = max (one, (a + zeros[l]) /  scales[l])
+		    s = s + (d1 - a) ** 2 / s1
+		}
+	    } else {
+		s1 = max (one, a)
+		do j = 1, n1
+		    s = s + (Memr[d[j]+k] - a) ** 2 / s1
+	    }
+	    n2 = n2 + n1
+
+	    # Save the average and sum for later.
+	    average[i] = a
+	    Memr[sums+k] = sum
+	}
+
+	# Here is the final sigma.
+	if (n2 > 1)
+	    s = sqrt (s / (n2 - 1))
+
+	# Reject pixels and compute the final average (if needed).
+	# There must be at least three pixels at each point for rejection.
+	# Iteratively scale the mean sigma and reject pixels
+	# Compact the data and keep track of the image IDs if needed.
+
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    if (n1 <= max (2, maxkeep)) {
+		if (!docombine) {
+		    if (n1 == 0)
+			average[i] = blank
+		    else {
+			sum = Memr[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Memr[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    a = average[i]
+	    sum = Memr[sums+k]
+
+	    repeat {
+		n2 = n1
+		if (s > 0.) {
+		    if (doscale1) {
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    mp1 = m[j] + k
+
+			    d1 = Memr[dp1]
+			    l = Memi[mp1]
+			    s1 = s * sqrt (max (one, (a+zeros[l]) /  scales[l]))
+			    r = (d1 - a) / s1
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs(r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    mp2 = m[n1] + k
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = l
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    } else {
+			s1 = s * sqrt (max (one, a))
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    d1 = Memr[dp1]
+			    r = (d1 - a) / s1
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs(r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[n1] + k
+					l = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = l
+				    }
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		}
+		if (n1 > 1)
+		    a = sum / n1
+	    } until (n1 == n2 || n1 <= max (2, maxkeep))
+
+	    # If too many are rejected add some back in.
+	    # Pixels with equal residuals are added together.
+	    if (n1 < maxkeep) {
+		nk = maxkeep
+		if (doscale1) {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Memr[dp1]
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Memr[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		} else {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Memr[dp1]
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Memr[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		}
+		if (n1 > 1)
+		    a = sum / n1
+	    }
+
+	    # Save the average if needed.
+	    n[i] = n1
+	    if (!docombine) {
+		if (n1 > 0)
+		    average[i] = a
+		else
+		    average[i] = blank
+	    }
+	}
+
+	# Check if the data flag has to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# IC_MAVSIGCLIP -- Reject pixels using an average sigma about the median
+# The average sigma is normalized by the expected poisson sigma.
+
+procedure ic_mavsigclipr (d, m, n, scales, zeros, nimages, npts, median)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+pointer	sp, resid, mp1, mp2
+real	med, low, high, r, s, s1, one
+data	one /1.0/
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Compute the poisson scaled average sigma about the median.
+	# There must be at least three pixels at each point to define
+	# the mean sigma.  Corrections for differences in the image
+	# scale factors are selected by the doscale1 flag.
+
+	s = 0.
+	n2 = 0
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (n1 < 3) {
+		if (n1 == 0)
+		    median[i] = blank
+		else if (n1 == 1)
+		    median[i] = Memr[d[1]+k]
+		else {
+		    low = Memr[d[1]+k]
+		    high = Memr[d[2]+k]
+		    median[i] = (low + high) / 2.
+		}
+		next
+	    }
+
+	    # Median
+	    n3 = 1 + n1 / 2
+	    if (mod (n1, 2) == 0) {
+		low = Memr[d[n3-1]+k]
+		high = Memr[d[n3]+k]
+		med = (low + high) / 2.
+	    } else
+		med = Memr[d[n3]+k]
+
+	    # Poisson scaled sigma accumulation
+	    if (doscale1) {
+		do j = 1, n1 {
+		    l = Memi[m[j]+k]
+		    s1 = max (one, (med + zeros[l]) /  scales[l])
+		    s = s + (Memr[d[j]+k] - med) ** 2 / s1
+		}
+	    } else {
+		s1 = max (one, med)
+		do j = 1, n1
+		    s = s + (Memr[d[j]+k] - med) ** 2 / s1
+	    }
+	    n2 = n2 + n1
+
+	    # Save the median for later.
+	    median[i] = med
+	}
+
+	# Here is the final sigma.
+	if (n2 > 1)
+	    s = sqrt (s / (n2 - 1))
+	else
+	    return
+
+	# Compute individual sigmas and iteratively clip.
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    if (n1 < max (3, maxkeep+1))
+		next
+	    nl = 1
+	    nh = n1
+	    med = median[i]
+
+	    repeat {
+		n2 = n1
+		n3 = nl + n1 / 2
+
+		if (n1 >= max (MINCLIP, maxkeep+1) && s > 0.) {
+		    if (doscale1) {
+			for (; nl <= n2; nl = nl + 1) {
+			    l = Memi[m[nl]+k]
+			    s1 = s * sqrt (max (one, (med+zeros[l])/scales[l]))
+			    r = (med - Memr[d[nl]+k]) / s1
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    l = Memi[m[nh]+k]
+			    s1 = s * sqrt (max (one, (med+zeros[l])/scales[l]))
+			    r = (Memr[d[nh]+k] - med) / s1
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    } else {
+			s1 = s * sqrt (max (one, med))
+			for (; nl <= n2; nl = nl + 1) {
+			    r = (med - Memr[d[nl]+k]) / s1
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    r = (Memr[d[nh]+k] - med) / s1
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    }
+
+		    # Recompute median
+		    if (n1 < n2) {
+			if (n1 > 0) {
+			    n3 = nl + n1 / 2
+			    if (mod (n1, 2) == 0) {
+				low = Memr[d[n3-1]+k]
+				high = Memr[d[n3]+k]
+				med = (low + high) / 2.
+			    } else
+				med = Memr[d[n3]+k]
+			} else
+			    med = blank
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    # If too many are rejected add some back in.
+	    # Pixels with equal residuals are added together.
+	    while (n1 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n1 = nh - nl + 1
+
+		# Recompute median
+		if (n1 < n2) {
+		    if (n1 > 0) {
+			n3 = nl + n1 / 2
+			if (mod (n1, 2) == 0) {
+			    low = Memr[d[n3-1]+k]
+			    high = Memr[d[n3]+k]
+			    med = (low + high) / 2.
+			} else
+			    med = Memr[d[n3]+k]
+		    } else
+			med = blank
+		}
+	    }
+
+	    # Only set median and reorder if needed
+	    n[i] = n1
+	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Memr[d[l]+k] = Memr[d[j]+k]
+			if (grow > 0) {
+			    mp1 = m[l] + k
+			    mp2 = m[j] + k
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+k] = Memi[m[j]+k]
+			j = j + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Memr[d[l]+k] = Memr[d[j]+k]
+			j = j + 1
+		    }
+		}
+	    }
+
+	    if (combine == MEDIAN)
+		median[i] = med
+	}
+
+	# Check if data flag needs to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag that the median is computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icaverage.x b/noao/imred/ccdred/src/combine/generic/icaverage.x
new file mode 100644
index 00000000..3646b725
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icaverage.x
@@ -0,0 +1,163 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	"../icombine.h"
+
+
+# IC_AVERAGE -- Compute the average image line.
+# Options include a weight average.
+
+procedure ic_averages (d, m, n, wts, npts, average)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image ID pointers
+int	n[npts]			# Number of points
+real	wts[ARB]		# Weights
+int	npts			# Number of output points per line
+real	average[npts]		# Average (returned)
+
+int	i, j, k
+real	sumwt, wt
+real	sum
+
+include	"../icombine.com"
+
+begin
+	# If no data has been excluded do the average without checking the
+	# number of points and using the fact that the weights are normalized.
+	# If all the data has been excluded set the average to the blank value.
+
+	if (dflag == D_ALL) {
+	    if (dowts) {
+		do i = 1, npts {
+		    k = i - 1
+		    wt = wts[Memi[m[1]+k]]
+		    sum = Mems[d[1]+k] * wt
+		    do j = 2, n[i] {
+			wt = wts[Memi[m[j]+k]]
+			sum = sum + Mems[d[j]+k] * wt
+		    }
+		    average[i] = sum
+		}
+	    } else {
+		do i = 1, npts {
+		    k = i - 1
+		    sum = Mems[d[1]+k]
+		    do j = 2, n[i]
+			sum = sum + Mems[d[j]+k]
+		    average[i] = sum / n[i]
+		}
+	    }
+	} else if (dflag == D_NONE) {
+	    do i = 1, npts
+		average[i] = blank
+	} else {
+	    if (dowts) {
+		do i = 1, npts {
+		    if (n[i] > 0) {
+			k = i - 1
+			wt = wts[Memi[m[1]+k]]
+			sum = Mems[d[1]+k] * wt
+			sumwt = wt
+			do j = 2, n[i] {
+			    wt = wts[Memi[m[j]+k]]
+			    sum = sum + Mems[d[j]+k] * wt
+			    sumwt = sumwt + wt
+			}
+			average[i] = sum / sumwt
+		    } else
+			average[i] = blank
+		}
+	    } else {
+		do i = 1, npts {
+		    if (n[i] > 0) {
+			k = i - 1
+			sum = Mems[d[1]+k]
+			do j = 2, n[i]
+			    sum = sum + Mems[d[j]+k]
+			average[i] = sum / n[i]
+		    } else
+			average[i] = blank
+		}
+	    }
+	}
+end
+
+# IC_AVERAGE -- Compute the average image line.
+# Options include a weight average.
+
+procedure ic_averager (d, m, n, wts, npts, average)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image ID pointers
+int	n[npts]			# Number of points
+real	wts[ARB]		# Weights
+int	npts			# Number of output points per line
+real	average[npts]		# Average (returned)
+
+int	i, j, k
+real	sumwt, wt
+real	sum
+
+include	"../icombine.com"
+
+begin
+	# If no data has been excluded do the average without checking the
+	# number of points and using the fact that the weights are normalized.
+	# If all the data has been excluded set the average to the blank value.
+
+	if (dflag == D_ALL) {
+	    if (dowts) {
+		do i = 1, npts {
+		    k = i - 1
+		    wt = wts[Memi[m[1]+k]]
+		    sum = Memr[d[1]+k] * wt
+		    do j = 2, n[i] {
+			wt = wts[Memi[m[j]+k]]
+			sum = sum + Memr[d[j]+k] * wt
+		    }
+		    average[i] = sum
+		}
+	    } else {
+		do i = 1, npts {
+		    k = i - 1
+		    sum = Memr[d[1]+k]
+		    do j = 2, n[i]
+			sum = sum + Memr[d[j]+k]
+		    average[i] = sum / n[i]
+		}
+	    }
+	} else if (dflag == D_NONE) {
+	    do i = 1, npts
+		average[i] = blank
+	} else {
+	    if (dowts) {
+		do i = 1, npts {
+		    if (n[i] > 0) {
+			k = i - 1
+			wt = wts[Memi[m[1]+k]]
+			sum = Memr[d[1]+k] * wt
+			sumwt = wt
+			do j = 2, n[i] {
+			    wt = wts[Memi[m[j]+k]]
+			    sum = sum + Memr[d[j]+k] * wt
+			    sumwt = sumwt + wt
+			}
+			average[i] = sum / sumwt
+		    } else
+			average[i] = blank
+		}
+	    } else {
+		do i = 1, npts {
+		    if (n[i] > 0) {
+			k = i - 1
+			sum = Memr[d[1]+k]
+			do j = 2, n[i]
+			    sum = sum + Memr[d[j]+k]
+			average[i] = sum / n[i]
+		    } else
+			average[i] = blank
+		}
+	    }
+	}
+end
diff --git a/noao/imred/ccdred/src/combine/generic/iccclip.x b/noao/imred/ccdred/src/combine/generic/iccclip.x
new file mode 100644
index 00000000..57709064
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/iccclip.x
@@ -0,0 +1,898 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+define	MINCLIP		2	# Mininum number of images for algorithm
+
+
+# IC_ACCDCLIP -- Reject pixels using CCD noise parameters about the average
+
+procedure ic_accdclips (d, m, n, scales, zeros, nm, nimages, npts, average)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+real	nm[3,nimages]		# Noise model parameters
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+real	d1, low, high, sum, a, s, r, zero
+data	zero /0.0/
+pointer	sp, resid, dp1, dp2, mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	# If there are no pixels go on to the combining.  Since the unweighted
+	# average is computed here possibly skip the combining later.
+
+	# There must be at least max (1, nkeep) pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	} else if (dowts || combine != AVERAGE)
+	    docombine = true
+	else
+	    docombine = false
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# There must be at least two pixels for rejection.  The initial
+	# average is the low/high rejected average except in the case of
+	# just two pixels.  The rejections are iterated and the average
+	# is recomputed.  Corrections for scaling may be performed.
+	# Depending on other flags the image IDs may also need to be adjusted.
+
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    if (n1 <= max (MINCLIP-1, maxkeep)) {
+		if (!docombine) {
+		    if (n1 == 0)
+			average[i] = blank
+		    else {
+			sum = Mems[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Mems[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    repeat {
+		if (n1 == 2) {
+		    sum = Mems[d[1]+k]
+		    sum = sum + Mems[d[2]+k]
+		    a = sum / 2
+		} else {
+		    low = Mems[d[1]+k]
+		    high = Mems[d[2]+k]
+		    if (low > high) {
+			d1 = low
+			low = high
+			high = d1
+		    }
+		    sum = 0.
+		    do j = 3, n1 {
+		       d1 = Mems[d[j]+k]
+		       if (d1 < low) {
+			   sum = sum + low
+			   low = d1
+			} else if (d1 > high) {
+			    sum = sum + high
+			    high = d1
+			} else
+			    sum = sum + d1
+		    }
+		    a = sum / (n1 - 2)
+		    sum = sum + low + high
+		}
+		n2 = n1
+		if (doscale1) {
+		    for (j=1; j<=n1; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+
+			l = Memi[mp1]
+			s = scales[l]
+			d1 = max (zero, s * (a + zeros[l]))
+			s = sqrt (nm[1,l] + d1/nm[2,l] + (d1*nm[3,l])**2) / s
+
+			d1 = Mems[dp1]
+			r = (d1 - a) / s
+			if (r < -lsigma || r > hsigma) {
+			    Memr[resid+n1] = abs(r)
+			    if (j < n1) {
+				dp2 = d[n1] + k
+				Mems[dp1] = Mems[dp2]
+				Mems[dp2] = d1
+				mp2 = m[n1] + k
+				Memi[mp1] = Memi[mp2]
+				Memi[mp2] = l
+				j = j - 1
+			    }
+			    sum = sum - d1
+			    n1 = n1 - 1
+			}
+		    }
+		} else {
+		    if (!keepids) {
+			s = max (zero, a)
+			s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2)
+		    }
+		    for (j=1; j<=n1; j=j+1) {
+			if (keepids) {
+			    l = Memi[m[j]+k]
+			    s = max (zero, a)
+			    s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2)
+			}
+			dp1 = d[j] + k
+			d1 = Mems[dp1]
+			r = (d1 - a) / s
+			if (r < -lsigma || r > hsigma) {
+			    Memr[resid+n1] = abs(r)
+			    if (j < n1) {
+				dp2 = d[n1] + k
+				Mems[dp1] = Mems[dp2]
+				Mems[dp2] = d1
+				if (keepids) {
+				    mp1 = m[j] + k
+				    mp2 = m[n1] + k
+				    l = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = l
+				}
+				j = j - 1
+			    }
+			    sum = sum - d1
+			    n1 = n1 - 1
+			}
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    if (n1 < maxkeep) {
+		nk = maxkeep
+		if (doscale1) {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Mems[dp1]
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Mems[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		} else {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Mems[dp1]
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Mems[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		}
+	    }
+
+	    n[i] = n1
+	    if (!docombine)
+		if (n1 > 0)
+		    average[i] = sum / n1
+		else
+		    average[i] = blank
+	}
+
+	# Check if the data flag has to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# IC_CCDCLIP -- Reject pixels using CCD noise parameters about the median
+
+procedure ic_mccdclips (d, m, n, scales, zeros, nm, nimages, npts, median)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+real	nm[3,nimages]		# Noise model
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+real	r, s
+pointer	sp, resid, mp1, mp2
+real	med, zero
+data	zero /0.0/
+
+include	"../icombine.com"
+
+begin
+	# There must be at least max (MINCLIP, nkeep+1) pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Compute median and sigma and iteratively clip.
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    nl = 1
+	    nh = n1
+
+	    repeat {
+		n2 = n1
+		n3 = nl + n1 / 2
+
+		if (n1 == 0)
+		    med = blank
+		else if (mod (n1, 2) == 0) {
+		    med = Mems[d[n3-1]+k]
+		    med = (med + Mems[d[n3]+k]) / 2.
+		} else
+		    med = Mems[d[n3]+k]
+
+		if (n1 >= max (MINCLIP, maxkeep+1)) {
+		    if (doscale1) {
+			for (; nl <= n2; nl = nl + 1) {
+			    l = Memi[m[nl]+k]
+			    s = scales[l]
+			    r = max (zero, s * (med + zeros[l]))
+			    s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s
+			    r = (med - Mems[d[nl]+k]) / s
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    l = Memi[m[nh]+k]
+			    s = scales[l]
+			    r = max (zero, s * (med + zeros[l]))
+			    s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s
+			    r = (Mems[d[nh]+k] - med) / s
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    } else {
+			if (!keepids) {
+			    s = max (zero, med)
+			    s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2)
+			}
+			for (; nl <= n2; nl = nl + 1) {
+			    if (keepids) {
+				l = Memi[m[nl]+k]
+				s = max (zero, med)
+				s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2)
+			    }
+			    r = (med - Mems[d[nl]+k]) / s
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    if (keepids) {
+				l = Memi[m[nh]+k]
+				s = max (zero, med)
+				s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2)
+			    }
+			    r = (Mems[d[nh]+k] -  med) / s
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    while (n1 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n1 = nh - nl + 1
+	    }
+
+	    # Only set median and reorder if needed
+	    n[i] = n1
+	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Mems[d[l]+k] = Mems[d[j]+k]
+			if (grow > 0) {
+			    mp1 = m[l] + k
+			    mp2 = m[j] + k
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+k] = Memi[m[j]+k]
+			j = j + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Mems[d[l]+k] = Mems[d[j]+k]
+			j = j + 1
+		    }
+		}
+	    }
+
+	    if (combine == MEDIAN)
+		median[i] = med
+	}
+
+	# Check if data flag needs to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag that the median is computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
+
+# IC_ACCDCLIP -- Reject pixels using CCD noise parameters about the average
+
+procedure ic_accdclipr (d, m, n, scales, zeros, nm, nimages, npts, average)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+real	nm[3,nimages]		# Noise model parameters
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+real	d1, low, high, sum, a, s, r, zero
+data	zero /0.0/
+pointer	sp, resid, dp1, dp2, mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	# If there are no pixels go on to the combining.  Since the unweighted
+	# average is computed here possibly skip the combining later.
+
+	# There must be at least max (1, nkeep) pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	} else if (dowts || combine != AVERAGE)
+	    docombine = true
+	else
+	    docombine = false
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# There must be at least two pixels for rejection.  The initial
+	# average is the low/high rejected average except in the case of
+	# just two pixels.  The rejections are iterated and the average
+	# is recomputed.  Corrections for scaling may be performed.
+	# Depending on other flags the image IDs may also need to be adjusted.
+
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    if (n1 <= max (MINCLIP-1, maxkeep)) {
+		if (!docombine) {
+		    if (n1 == 0)
+			average[i] = blank
+		    else {
+			sum = Memr[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Memr[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    repeat {
+		if (n1 == 2) {
+		    sum = Memr[d[1]+k]
+		    sum = sum + Memr[d[2]+k]
+		    a = sum / 2
+		} else {
+		    low = Memr[d[1]+k]
+		    high = Memr[d[2]+k]
+		    if (low > high) {
+			d1 = low
+			low = high
+			high = d1
+		    }
+		    sum = 0.
+		    do j = 3, n1 {
+		       d1 = Memr[d[j]+k]
+		       if (d1 < low) {
+			   sum = sum + low
+			   low = d1
+			} else if (d1 > high) {
+			    sum = sum + high
+			    high = d1
+			} else
+			    sum = sum + d1
+		    }
+		    a = sum / (n1 - 2)
+		    sum = sum + low + high
+		}
+		n2 = n1
+		if (doscale1) {
+		    for (j=1; j<=n1; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+
+			l = Memi[mp1]
+			s = scales[l]
+			d1 = max (zero, s * (a + zeros[l]))
+			s = sqrt (nm[1,l] + d1/nm[2,l] + (d1*nm[3,l])**2) / s
+
+			d1 = Memr[dp1]
+			r = (d1 - a) / s
+			if (r < -lsigma || r > hsigma) {
+			    Memr[resid+n1] = abs(r)
+			    if (j < n1) {
+				dp2 = d[n1] + k
+				Memr[dp1] = Memr[dp2]
+				Memr[dp2] = d1
+				mp2 = m[n1] + k
+				Memi[mp1] = Memi[mp2]
+				Memi[mp2] = l
+				j = j - 1
+			    }
+			    sum = sum - d1
+			    n1 = n1 - 1
+			}
+		    }
+		} else {
+		    if (!keepids) {
+			s = max (zero, a)
+			s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2)
+		    }
+		    for (j=1; j<=n1; j=j+1) {
+			if (keepids) {
+			    l = Memi[m[j]+k]
+			    s = max (zero, a)
+			    s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2)
+			}
+			dp1 = d[j] + k
+			d1 = Memr[dp1]
+			r = (d1 - a) / s
+			if (r < -lsigma || r > hsigma) {
+			    Memr[resid+n1] = abs(r)
+			    if (j < n1) {
+				dp2 = d[n1] + k
+				Memr[dp1] = Memr[dp2]
+				Memr[dp2] = d1
+				if (keepids) {
+				    mp1 = m[j] + k
+				    mp2 = m[n1] + k
+				    l = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = l
+				}
+				j = j - 1
+			    }
+			    sum = sum - d1
+			    n1 = n1 - 1
+			}
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    if (n1 < maxkeep) {
+		nk = maxkeep
+		if (doscale1) {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Memr[dp1]
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Memr[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		} else {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Memr[dp1]
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Memr[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		}
+	    }
+
+	    n[i] = n1
+	    if (!docombine)
+		if (n1 > 0)
+		    average[i] = sum / n1
+		else
+		    average[i] = blank
+	}
+
+	# Check if the data flag has to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# IC_CCDCLIP -- Reject pixels using CCD noise parameters about the median
+
+procedure ic_mccdclipr (d, m, n, scales, zeros, nm, nimages, npts, median)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+real	nm[3,nimages]		# Noise model
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+real	r, s
+pointer	sp, resid, mp1, mp2
+real	med, zero
+data	zero /0.0/
+
+include	"../icombine.com"
+
+begin
+	# There must be at least max (MINCLIP, nkeep+1) pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Compute median and sigma and iteratively clip.
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    nl = 1
+	    nh = n1
+
+	    repeat {
+		n2 = n1
+		n3 = nl + n1 / 2
+
+		if (n1 == 0)
+		    med = blank
+		else if (mod (n1, 2) == 0) {
+		    med = Memr[d[n3-1]+k]
+		    med = (med + Memr[d[n3]+k]) / 2.
+		} else
+		    med = Memr[d[n3]+k]
+
+		if (n1 >= max (MINCLIP, maxkeep+1)) {
+		    if (doscale1) {
+			for (; nl <= n2; nl = nl + 1) {
+			    l = Memi[m[nl]+k]
+			    s = scales[l]
+			    r = max (zero, s * (med + zeros[l]))
+			    s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s
+			    r = (med - Memr[d[nl]+k]) / s
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    l = Memi[m[nh]+k]
+			    s = scales[l]
+			    r = max (zero, s * (med + zeros[l]))
+			    s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s
+			    r = (Memr[d[nh]+k] - med) / s
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    } else {
+			if (!keepids) {
+			    s = max (zero, med)
+			    s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2)
+			}
+			for (; nl <= n2; nl = nl + 1) {
+			    if (keepids) {
+				l = Memi[m[nl]+k]
+				s = max (zero, med)
+				s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2)
+			    }
+			    r = (med - Memr[d[nl]+k]) / s
+			    if (r <= lsigma)
+				break
+			    Memr[resid+nl] = r
+			    n1 = n1 - 1
+			}
+			for (; nh >= nl; nh = nh - 1) {
+			    if (keepids) {
+				l = Memi[m[nh]+k]
+				s = max (zero, med)
+				s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2)
+			    }
+			    r = (Memr[d[nh]+k] -  med) / s
+			    if (r <= hsigma)
+				break
+			    Memr[resid+nh] = r
+			    n1 = n1 - 1
+			}
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    while (n1 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n1 = nh - nl + 1
+	    }
+
+	    # Only set median and reorder if needed
+	    n[i] = n1
+	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Memr[d[l]+k] = Memr[d[j]+k]
+			if (grow > 0) {
+			    mp1 = m[l] + k
+			    mp2 = m[j] + k
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+k] = Memi[m[j]+k]
+			j = j + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Memr[d[l]+k] = Memr[d[j]+k]
+			j = j + 1
+		    }
+		}
+	    }
+
+	    if (combine == MEDIAN)
+		median[i] = med
+	}
+
+	# Check if data flag needs to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag that the median is computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icgdata.x b/noao/imred/ccdred/src/combine/generic/icgdata.x
new file mode 100644
index 00000000..5c6ac18c
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icgdata.x
@@ -0,0 +1,459 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<mach.h>
+include	"../icombine.h"
+
+
+# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
+# Entirely empty lines are excluded.  The data are compacted within the
+# input data buffers.  If it is required, the connection to the original
+# image index is keeped in the returned m data pointers.
+
+procedure ic_gdatas (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
+	zeros, nimages, npts, v1, v2)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output images
+pointer	dbuf[nimages]		# Data buffers for nonaligned images
+pointer	d[nimages]		# Data pointers
+pointer	id[nimages]		# ID pointers
+int	n[npts]			# Number of good pixels
+pointer	m[nimages]		# Mask pointers
+int	lflag[nimages]		# Empty mask flags
+int	offsets[nimages,ARB]	# Image offsets
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero offset factors
+int	nimages			# Number of input images
+int	npts			# NUmber of output points per line
+long	v1[ARB], v2[ARB]	# Line vectors
+
+int	i, j, k, l, ndim, nused
+real	a, b
+pointer	buf, dp, ip, mp, imgnls()
+
+include	"../icombine.com"
+
+begin
+	# Get masks and return if there is no data
+	call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
+	if (dflag == D_NONE)
+	    return
+
+	# Get data and fill data buffers.  Correct for offsets if needed.
+	ndim = IM_NDIM(out[1])
+	do i = 1, nimages {
+	    if (lflag[i] == D_NONE)
+		next
+	    if (aligned) {
+		call amovl (v1, v2, IM_MAXDIM)
+		if (project)
+		    v2[ndim+1] = i
+		j = imgnls (in[i], d[i], v2)
+	    } else {
+		v2[1] = v1[1]
+		do j = 2, ndim
+		    v2[j] = v1[j] - offsets[i,j]
+		if (project)
+		    v2[ndim+1] = i
+		j = imgnls (in[i], buf, v2)
+		call amovs (Mems[buf], Mems[dbuf[i]+offsets[i,1]],
+		    IM_LEN(in[i],1))
+		d[i] = dbuf[i]
+	    }
+	}
+
+	# Apply threshold if needed
+	if (dothresh) {
+	    do i = 1, nimages {
+		dp = d[i]
+		if (lflag[i] == D_ALL) {
+		    do j = 1, npts {
+			a = Mems[dp]
+			if (a < lthresh || a > hthresh) {
+			    Memi[m[i]+j-1] = 1
+			    lflag[i] = D_MIX
+			    dflag = D_MIX
+			}
+			dp = dp + 1
+		    }
+		} else if (lflag[i] == D_MIX) {
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    a = Mems[dp]
+			    if (a < lthresh || a > hthresh) {
+				Memi[m[i]+j-1] = 1
+				dflag = D_MIX
+			    }
+			}
+			dp = dp + 1
+			mp = mp + 1
+		    }
+		}
+
+		# Check for completely empty lines
+		if (lflag[i] == D_MIX) {
+		    lflag[i] = D_NONE
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    lflag[i] = D_MIX
+			    break
+			}
+			mp = mp + 1
+		    }
+		}
+	    }
+	}
+
+	# Apply scaling (avoiding masked pixels which might overflow?)
+	if (doscale) {
+	    if (dflag == D_ALL) {
+		do i = 1, nimages {
+		    dp = d[i]
+		    a = scales[i]
+		    b = -zeros[i]
+		    do j = 1, npts {
+			Mems[dp] = Mems[dp] / a + b
+			dp = dp + 1
+		    }
+		}
+	    } else if (dflag == D_MIX) {
+		do i = 1, nimages {
+		    dp = d[i]
+		    a = scales[i]
+		    b = -zeros[i]
+		    if (lflag[i] == D_ALL) {
+			do j = 1, npts {
+			    Mems[dp] = Mems[dp] / a + b
+			    dp = dp + 1
+			}
+		    } else if (lflag[i] == D_MIX) {
+			mp = m[i]
+			do j = 1, npts {
+			    if (Memi[mp] == 0)
+				Mems[dp] = Mems[dp] / a + b
+			    dp = dp + 1
+			    mp = mp + 1
+			}
+		    }
+		}
+	    }
+	}
+
+	# Sort pointers to exclude unused images.
+	# Use the lflag array to keep track of the image index.
+
+	if (dflag == D_ALL)
+	    nused = nimages
+	else {
+	    nused = 0
+	    do i = 1, nimages
+		if (lflag[i] != D_NONE) {
+		    nused = nused + 1
+		    d[nused] = d[i]
+		    m[nused] = m[i]
+		    lflag[nused] = i
+		}
+	    if (nused == 0)
+		dflag = D_NONE
+	}
+
+	# Compact data to remove bad pixels
+	# Keep track of the image indices if needed
+	# If growing mark the end of the included image indices with zero
+
+	if (dflag == D_ALL) {
+	    call amovki (nused, n, npts)
+	    if (keepids)
+		do i = 1, nimages
+		    call amovki (i, Memi[id[i]], npts)
+	} else if (dflag == D_NONE)
+	    call aclri (n, npts)
+	else {
+	    call aclri (n, npts)
+	    if (keepids) {
+		do i = 1, nused {
+		    l = lflag[i]
+		    dp = d[i]
+		    ip = id[i]
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    n[j] = n[j] + 1
+			    k = n[j]
+			    if (k < i) {
+				Mems[d[k]+j-1] = Mems[dp]
+				Memi[id[k]+j-1] = l
+			    } else
+				Memi[ip] = l
+			}
+			dp = dp + 1
+			ip = ip + 1
+			mp = mp + 1
+		    }
+		}
+		if (grow > 0) {
+		    do j = 1, npts {
+			do i = n[j]+1, nimages
+			    Memi[id[i]+j-1] = 0
+		    }
+		}
+	    } else {
+		do i = 1, nused {
+		    dp = d[i]
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    n[j] = n[j] + 1
+			    k = n[j]
+			    if (k < i)
+				Mems[d[k]+j-1] = Mems[dp]
+			}
+			dp = dp + 1
+			mp = mp + 1
+		    }
+		}
+	    }
+	}
+
+	# Sort the pixels and IDs if needed
+	if (mclip) {
+	    call malloc (dp, nimages, TY_SHORT)
+	    if (keepids) {
+		call malloc (ip, nimages, TY_INT)
+		call ic_2sorts (d, Mems[dp], id, Memi[ip], n, npts)
+		call mfree (ip, TY_INT)
+	    } else
+		call ic_sorts (d, Mems[dp], n, npts)
+	    call mfree (dp, TY_SHORT)
+	}
+end
+
+# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
+# Entirely empty lines are excluded.  The data are compacted within the
+# input data buffers.  If it is required, the connection to the original
+# image index is keeped in the returned m data pointers.
+
+procedure ic_gdatar (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
+	zeros, nimages, npts, v1, v2)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output images
+pointer	dbuf[nimages]		# Data buffers for nonaligned images
+pointer	d[nimages]		# Data pointers
+pointer	id[nimages]		# ID pointers
+int	n[npts]			# Number of good pixels
+pointer	m[nimages]		# Mask pointers
+int	lflag[nimages]		# Empty mask flags
+int	offsets[nimages,ARB]	# Image offsets
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero offset factors
+int	nimages			# Number of input images
+int	npts			# NUmber of output points per line
+long	v1[ARB], v2[ARB]	# Line vectors
+
+int	i, j, k, l, ndim, nused
+real	a, b
+pointer	buf, dp, ip, mp, imgnlr()
+
+include	"../icombine.com"
+
+begin
+	# Get masks and return if there is no data
+	call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
+	if (dflag == D_NONE)
+	    return
+
+	# Get data and fill data buffers.  Correct for offsets if needed.
+	ndim = IM_NDIM(out[1])
+	do i = 1, nimages {
+	    if (lflag[i] == D_NONE)
+		next
+	    if (aligned) {
+		call amovl (v1, v2, IM_MAXDIM)
+		if (project)
+		    v2[ndim+1] = i
+		j = imgnlr (in[i], d[i], v2)
+	    } else {
+		v2[1] = v1[1]
+		do j = 2, ndim
+		    v2[j] = v1[j] - offsets[i,j]
+		if (project)
+		    v2[ndim+1] = i
+		j = imgnlr (in[i], buf, v2)
+		call amovr (Memr[buf], Memr[dbuf[i]+offsets[i,1]],
+		    IM_LEN(in[i],1))
+		d[i] = dbuf[i]
+	    }
+	}
+
+	# Apply threshold if needed
+	if (dothresh) {
+	    do i = 1, nimages {
+		dp = d[i]
+		if (lflag[i] == D_ALL) {
+		    do j = 1, npts {
+			a = Memr[dp]
+			if (a < lthresh || a > hthresh) {
+			    Memi[m[i]+j-1] = 1
+			    lflag[i] = D_MIX
+			    dflag = D_MIX
+			}
+			dp = dp + 1
+		    }
+		} else if (lflag[i] == D_MIX) {
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    a = Memr[dp]
+			    if (a < lthresh || a > hthresh) {
+				Memi[m[i]+j-1] = 1
+				dflag = D_MIX
+			    }
+			}
+			dp = dp + 1
+			mp = mp + 1
+		    }
+		}
+
+		# Check for completely empty lines
+		if (lflag[i] == D_MIX) {
+		    lflag[i] = D_NONE
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    lflag[i] = D_MIX
+			    break
+			}
+			mp = mp + 1
+		    }
+		}
+	    }
+	}
+
+	# Apply scaling (avoiding masked pixels which might overflow?)
+	if (doscale) {
+	    if (dflag == D_ALL) {
+		do i = 1, nimages {
+		    dp = d[i]
+		    a = scales[i]
+		    b = -zeros[i]
+		    do j = 1, npts {
+			Memr[dp] = Memr[dp] / a + b
+			dp = dp + 1
+		    }
+		}
+	    } else if (dflag == D_MIX) {
+		do i = 1, nimages {
+		    dp = d[i]
+		    a = scales[i]
+		    b = -zeros[i]
+		    if (lflag[i] == D_ALL) {
+			do j = 1, npts {
+			    Memr[dp] = Memr[dp] / a + b
+			    dp = dp + 1
+			}
+		    } else if (lflag[i] == D_MIX) {
+			mp = m[i]
+			do j = 1, npts {
+			    if (Memi[mp] == 0)
+				Memr[dp] = Memr[dp] / a + b
+			    dp = dp + 1
+			    mp = mp + 1
+			}
+		    }
+		}
+	    }
+	}
+
+	# Sort pointers to exclude unused images.
+	# Use the lflag array to keep track of the image index.
+
+	if (dflag == D_ALL)
+	    nused = nimages
+	else {
+	    nused = 0
+	    do i = 1, nimages
+		if (lflag[i] != D_NONE) {
+		    nused = nused + 1
+		    d[nused] = d[i]
+		    m[nused] = m[i]
+		    lflag[nused] = i
+		}
+	    if (nused == 0)
+		dflag = D_NONE
+	}
+
+	# Compact data to remove bad pixels
+	# Keep track of the image indices if needed
+	# If growing mark the end of the included image indices with zero
+
+	if (dflag == D_ALL) {
+	    call amovki (nused, n, npts)
+	    if (keepids)
+		do i = 1, nimages
+		    call amovki (i, Memi[id[i]], npts)
+	} else if (dflag == D_NONE)
+	    call aclri (n, npts)
+	else {
+	    call aclri (n, npts)
+	    if (keepids) {
+		do i = 1, nused {
+		    l = lflag[i]
+		    dp = d[i]
+		    ip = id[i]
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    n[j] = n[j] + 1
+			    k = n[j]
+			    if (k < i) {
+				Memr[d[k]+j-1] = Memr[dp]
+				Memi[id[k]+j-1] = l
+			    } else
+				Memi[ip] = l
+			}
+			dp = dp + 1
+			ip = ip + 1
+			mp = mp + 1
+		    }
+		}
+		if (grow > 0) {
+		    do j = 1, npts {
+			do i = n[j]+1, nimages
+			    Memi[id[i]+j-1] = 0
+		    }
+		}
+	    } else {
+		do i = 1, nused {
+		    dp = d[i]
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memi[mp] == 0) {
+			    n[j] = n[j] + 1
+			    k = n[j]
+			    if (k < i)
+				Memr[d[k]+j-1] = Memr[dp]
+			}
+			dp = dp + 1
+			mp = mp + 1
+		    }
+		}
+	    }
+	}
+
+	# Sort the pixels and IDs if needed
+	if (mclip) {
+	    call malloc (dp, nimages, TY_REAL)
+	    if (keepids) {
+		call malloc (ip, nimages, TY_INT)
+		call ic_2sortr (d, Memr[dp], id, Memi[ip], n, npts)
+		call mfree (ip, TY_INT)
+	    } else
+		call ic_sortr (d, Memr[dp], n, npts)
+	    call mfree (dp, TY_REAL)
+	}
+end
+
diff --git a/noao/imred/ccdred/src/combine/generic/icgrow.x b/noao/imred/ccdred/src/combine/generic/icgrow.x
new file mode 100644
index 00000000..b94e1cbc
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icgrow.x
@@ -0,0 +1,148 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+
+# IC_GROW --  Reject neigbors of rejected pixels.
+# The rejected pixels are marked by having nonzero ids beyond the number
+# of included pixels.  The pixels rejected here are given zero ids
+# to avoid growing of the pixels rejected here.  The unweighted average
+# can be updated but any rejected pixels requires the median to be
+# recomputed.  When the number of pixels at a grow point reaches nkeep 
+# no further pixels are rejected.  Note that the rejection order is not
+# based on the magnitude of the residuals and so a grow from a weakly
+# rejected image pixel may take precedence over a grow from a strongly
+# rejected image pixel.
+
+procedure ic_grows (d, m, n, nimages, npts, average)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image id pointers
+int	n[npts]			# Number of good pixels
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i1, i2, j1, j2, k1, k2, l, is, ie, n2, maxkeep
+pointer	mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	if (dflag == D_NONE)
+	    return
+
+	do i1 = 1, npts {
+	    k1 = i1 - 1
+	    is = max (1, i1 - grow)
+	    ie = min (npts, i1 + grow)
+	    do j1 = n[i1]+1, nimages {
+		l = Memi[m[j1]+k1]
+		if (l == 0)
+		    next
+		if (combine == MEDIAN)
+		    docombine = true
+
+		do i2 = is, ie {
+		    if (i2 == i1)
+			next
+		    k2 = i2 - 1
+		    n2 = n[i2]
+		    if (nkeep < 0)
+			maxkeep = max (0, n2 + nkeep)
+		    else
+			maxkeep = min (n2, nkeep)
+		    if (n2 <= maxkeep)
+			next
+		    do j2 = 1, n2 {
+			mp1 = m[j2] + k2
+			if (Memi[mp1] == l) {
+			    if (!docombine && n2 > 1)
+				average[i2] =
+				    (n2*average[i2] - Mems[d[j2]+k2]) / (n2-1)
+			    mp2 = m[n2] + k2
+			    if (j2 < n2) {
+				Mems[d[j2]+k2] = Mems[d[n2]+k2]
+				Memi[mp1] = Memi[mp2]
+			    }
+			    Memi[mp2] = 0
+			    n[i2] = n2 - 1
+			    break
+			}
+		    }
+		}
+	    }
+	}
+end
+
+# IC_GROW --  Reject neigbors of rejected pixels.
+# The rejected pixels are marked by having nonzero ids beyond the number
+# of included pixels.  The pixels rejected here are given zero ids
+# to avoid growing of the pixels rejected here.  The unweighted average
+# can be updated but any rejected pixels requires the median to be
+# recomputed.  When the number of pixels at a grow point reaches nkeep 
+# no further pixels are rejected.  Note that the rejection order is not
+# based on the magnitude of the residuals and so a grow from a weakly
+# rejected image pixel may take precedence over a grow from a strongly
+# rejected image pixel.
+
+procedure ic_growr (d, m, n, nimages, npts, average)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image id pointers
+int	n[npts]			# Number of good pixels
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i1, i2, j1, j2, k1, k2, l, is, ie, n2, maxkeep
+pointer	mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	if (dflag == D_NONE)
+	    return
+
+	do i1 = 1, npts {
+	    k1 = i1 - 1
+	    is = max (1, i1 - grow)
+	    ie = min (npts, i1 + grow)
+	    do j1 = n[i1]+1, nimages {
+		l = Memi[m[j1]+k1]
+		if (l == 0)
+		    next
+		if (combine == MEDIAN)
+		    docombine = true
+
+		do i2 = is, ie {
+		    if (i2 == i1)
+			next
+		    k2 = i2 - 1
+		    n2 = n[i2]
+		    if (nkeep < 0)
+			maxkeep = max (0, n2 + nkeep)
+		    else
+			maxkeep = min (n2, nkeep)
+		    if (n2 <= maxkeep)
+			next
+		    do j2 = 1, n2 {
+			mp1 = m[j2] + k2
+			if (Memi[mp1] == l) {
+			    if (!docombine && n2 > 1)
+				average[i2] =
+				    (n2*average[i2] - Memr[d[j2]+k2]) / (n2-1)
+			    mp2 = m[n2] + k2
+			    if (j2 < n2) {
+				Memr[d[j2]+k2] = Memr[d[n2]+k2]
+				Memi[mp1] = Memi[mp2]
+			    }
+			    Memi[mp2] = 0
+			    n[i2] = n2 - 1
+			    break
+			}
+		    }
+		}
+	    }
+	}
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icmedian.x b/noao/imred/ccdred/src/combine/generic/icmedian.x
new file mode 100644
index 00000000..ec0166ba
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icmedian.x
@@ -0,0 +1,343 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+
+# IC_MEDIAN -- Median of lines
+
+procedure ic_medians (d, n, npts, median)
+
+pointer	d[ARB]			# Input data line pointers
+int	n[npts]			# Number of good pixels
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, j1, j2, n1, lo, up, lo1, up1
+bool	even
+real	val1, val2, val3
+short	temp, wtemp
+
+include	"../icombine.com"
+
+begin
+	# If no data return after possibly setting blank values.
+	if (dflag == D_NONE) {
+	    do i = 1, npts
+		median[i]= blank
+	    return
+	}
+
+	# If the data were previously sorted then directly compute the median.
+	if (mclip) {
+	    if (dflag == D_ALL) {
+		n1 = n[1]
+		even = (mod (n1, 2) == 0)
+		j1 = n1 / 2 + 1
+		j2 = n1 / 2
+		do i = 1, npts {
+		    k = i - 1
+		    if (even) {
+			val1 = Mems[d[j1]+k]
+			val2 = Mems[d[j2]+k]
+			median[i] = (val1 + val2) / 2.
+		    } else
+			median[i] = Mems[d[j1]+k]
+		}
+	    } else {
+		do i = 1, npts {
+		    k = i - 1
+		    n1 = n[i]
+		    if (n1 > 0) {
+			j1 = n1 / 2 + 1
+			if (mod (n1, 2) == 0) {
+			    j2 = n1 / 2
+			    val1 = Mems[d[j1]+k]
+			    val2 = Mems[d[j2]+k]
+			    median[i] = (val1 + val2) / 2.
+			} else
+			    median[i] = Mems[d[j1]+k]
+		    } else
+			median[i] = blank
+		}
+	    }
+	    return
+	}
+
+	# Compute the median.
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+
+	    # If there are more than 3 points use Wirth algorithm.  This
+	    # is the same as vops$amed.gx except for an even number of
+	    # points it selects the middle two and averages.
+	    if (n1 > 3) {
+		lo = 1
+		up = n1
+		j  = max (lo, min (up, (up+1)/2))
+
+		while (lo < up) {
+		    if (! (lo < up))
+			break
+
+		    temp = Mems[d[j]+k];  lo1 = lo;  up1 = up
+
+		    repeat {
+			while (Mems[d[lo1]+k] < temp)
+			    lo1 = lo1 + 1
+			while (temp < Mems[d[up1]+k])
+			    up1 = up1 - 1
+			if (lo1 <= up1) {
+			    wtemp = Mems[d[lo1]+k]
+			    Mems[d[lo1]+k] = Mems[d[up1]+k]
+			    Mems[d[up1]+k] = wtemp
+			    lo1 = lo1 + 1;  up1 = up1 - 1
+			}
+		    } until (lo1 > up1)
+
+		    if (up1 < j)
+			lo = lo1
+		    if (j < lo1)
+			up = up1
+		}
+
+		median[i] = Mems[d[j]+k]
+
+		if (mod (n1,2) == 0) {
+		    lo = 1
+		    up = n1
+		    j  = max (lo, min (up, (up+1)/2)+1)
+
+		    while (lo < up) {
+			if (! (lo < up))
+			    break
+
+			temp = Mems[d[j]+k];  lo1 = lo;  up1 = up
+
+			repeat {
+			    while (Mems[d[lo1]+k] < temp)
+				lo1 = lo1 + 1
+			    while (temp < Mems[d[up1]+k])
+				up1 = up1 - 1
+			    if (lo1 <= up1) {
+				wtemp = Mems[d[lo1]+k]
+				Mems[d[lo1]+k] = Mems[d[up1]+k]
+				Mems[d[up1]+k] = wtemp
+				lo1 = lo1 + 1;  up1 = up1 - 1
+			    }
+			} until (lo1 > up1)
+
+			if (up1 < j)
+			    lo = lo1
+			if (j < lo1)
+			    up = up1
+		    }
+		    median[i] = (median[i] + Mems[d[j]+k]) / 2
+		}
+
+	    # If 3 points find the median directly.
+	    } else if (n1 == 3) {
+	        val1 = Mems[d[1]+k]
+	        val2 = Mems[d[2]+k]
+	        val3 = Mems[d[3]+k]
+	        if (val1 < val2) {
+		    if (val2 < val3)		# abc
+		        median[i] = val2
+		    else if (val1 < val3)	# acb
+		        median[i] = val3
+		    else			# cab
+		        median[i] = val1
+	        } else {
+		    if (val2 > val3)		# cba
+		        median[i] = val2
+		    else if (val1 > val3)	# bca
+		        median[i] = val3
+		    else			# bac
+		        median[i] = val1
+	        }
+
+	    # If 2 points average.
+	    } else if (n1 == 2) {
+		val1 = Mems[d[1]+k]
+		val2 = Mems[d[2]+k]
+		median[i] = (val1 + val2) / 2
+
+	    # If 1 point return the value.
+	    } else if (n1 == 1)
+		median[i] = Mems[d[1]+k]
+
+	    # If no points return with a possibly blank value.
+	    else
+		median[i] = blank
+	}
+end
+
+# IC_MEDIAN -- Median of lines
+
+procedure ic_medianr (d, n, npts, median)
+
+pointer	d[ARB]			# Input data line pointers
+int	n[npts]			# Number of good pixels
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, j1, j2, n1, lo, up, lo1, up1
+bool	even
+real	val1, val2, val3
+real	temp, wtemp
+
+include	"../icombine.com"
+
+begin
+	# If no data return after possibly setting blank values.
+	if (dflag == D_NONE) {
+	    do i = 1, npts
+		median[i]= blank
+	    return
+	}
+
+	# If the data were previously sorted then directly compute the median.
+	if (mclip) {
+	    if (dflag == D_ALL) {
+		n1 = n[1]
+		even = (mod (n1, 2) == 0)
+		j1 = n1 / 2 + 1
+		j2 = n1 / 2
+		do i = 1, npts {
+		    k = i - 1
+		    if (even) {
+			val1 = Memr[d[j1]+k]
+			val2 = Memr[d[j2]+k]
+			median[i] = (val1 + val2) / 2.
+		    } else
+			median[i] = Memr[d[j1]+k]
+		}
+	    } else {
+		do i = 1, npts {
+		    k = i - 1
+		    n1 = n[i]
+		    if (n1 > 0) {
+			j1 = n1 / 2 + 1
+			if (mod (n1, 2) == 0) {
+			    j2 = n1 / 2
+			    val1 = Memr[d[j1]+k]
+			    val2 = Memr[d[j2]+k]
+			    median[i] = (val1 + val2) / 2.
+			} else
+			    median[i] = Memr[d[j1]+k]
+		    } else
+			median[i] = blank
+		}
+	    }
+	    return
+	}
+
+	# Compute the median.
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+
+	    # If there are more than 3 points use Wirth algorithm.  This
+	    # is the same as vops$amed.gx except for an even number of
+	    # points it selects the middle two and averages.
+	    if (n1 > 3) {
+		lo = 1
+		up = n1
+		j  = max (lo, min (up, (up+1)/2))
+
+		while (lo < up) {
+		    if (! (lo < up))
+			break
+
+		    temp = Memr[d[j]+k];  lo1 = lo;  up1 = up
+
+		    repeat {
+			while (Memr[d[lo1]+k] < temp)
+			    lo1 = lo1 + 1
+			while (temp < Memr[d[up1]+k])
+			    up1 = up1 - 1
+			if (lo1 <= up1) {
+			    wtemp = Memr[d[lo1]+k]
+			    Memr[d[lo1]+k] = Memr[d[up1]+k]
+			    Memr[d[up1]+k] = wtemp
+			    lo1 = lo1 + 1;  up1 = up1 - 1
+			}
+		    } until (lo1 > up1)
+
+		    if (up1 < j)
+			lo = lo1
+		    if (j < lo1)
+			up = up1
+		}
+
+		median[i] = Memr[d[j]+k]
+
+		if (mod (n1,2) == 0) {
+		    lo = 1
+		    up = n1
+		    j  = max (lo, min (up, (up+1)/2)+1)
+
+		    while (lo < up) {
+			if (! (lo < up))
+			    break
+
+			temp = Memr[d[j]+k];  lo1 = lo;  up1 = up
+
+			repeat {
+			    while (Memr[d[lo1]+k] < temp)
+				lo1 = lo1 + 1
+			    while (temp < Memr[d[up1]+k])
+				up1 = up1 - 1
+			    if (lo1 <= up1) {
+				wtemp = Memr[d[lo1]+k]
+				Memr[d[lo1]+k] = Memr[d[up1]+k]
+				Memr[d[up1]+k] = wtemp
+				lo1 = lo1 + 1;  up1 = up1 - 1
+			    }
+			} until (lo1 > up1)
+
+			if (up1 < j)
+			    lo = lo1
+			if (j < lo1)
+			    up = up1
+		    }
+		    median[i] = (median[i] + Memr[d[j]+k]) / 2
+		}
+
+	    # If 3 points find the median directly.
+	    } else if (n1 == 3) {
+	        val1 = Memr[d[1]+k]
+	        val2 = Memr[d[2]+k]
+	        val3 = Memr[d[3]+k]
+	        if (val1 < val2) {
+		    if (val2 < val3)		# abc
+		        median[i] = val2
+		    else if (val1 < val3)	# acb
+		        median[i] = val3
+		    else			# cab
+		        median[i] = val1
+	        } else {
+		    if (val2 > val3)		# cba
+		        median[i] = val2
+		    else if (val1 > val3)	# bca
+		        median[i] = val3
+		    else			# bac
+		        median[i] = val1
+	        }
+
+	    # If 2 points average.
+	    } else if (n1 == 2) {
+		val1 = Memr[d[1]+k]
+		val2 = Memr[d[2]+k]
+		median[i] = (val1 + val2) / 2
+
+	    # If 1 point return the value.
+	    } else if (n1 == 1)
+		median[i] = Memr[d[1]+k]
+
+	    # If no points return with a possibly blank value.
+	    else
+		median[i] = blank
+	}
+end
+
diff --git a/noao/imred/ccdred/src/combine/generic/icmm.x b/noao/imred/ccdred/src/combine/generic/icmm.x
new file mode 100644
index 00000000..259759bd
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icmm.x
@@ -0,0 +1,300 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+
+# IC_MM --  Reject a specified number of high and low pixels
+
+procedure ic_mms (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
+short	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
+		d1 = Mems[k]
+		dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k
+		do j = 2, n1 {
+		    d2 = d1
+		    k = d[j] + i1
+		    d1 = Mems[k]
+		    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) {
+			    Mems[kmax] = d2
+			    Memi[m[jmax]+i1] = Memi[m[j]+i1]
+			} else {
+			    Mems[kmax] = d1
+			    Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+			}
+		    }
+		    if (jmin < j) {
+			if (jmax != n1) {
+			    Mems[kmin] = d1
+			    Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+			} else {
+			    Mems[kmin] = d2
+			    Memi[m[jmin]+i1] = Memi[m[j]+i1]
+			}
+		    }
+		} else {
+		    if (jmax < j) {
+			if (jmin != j)
+			    Mems[kmax] = d2
+			else
+			    Mems[kmax] = d1
+		    }
+		    if (jmin < j) {
+			if (jmax != n1)
+			    Mems[kmin] = d1
+			else
+			    Mems[kmin] = d2
+		    }
+		}
+		n1 = n1 - 2
+	    }
+
+	    # Reject the excess low points.
+	    do np = 1, nlow {
+		k = d[1] + i1
+		d1 = Mems[k]
+		dmin = d1; jmin = 1; kmin = k
+		do j = 2, n1 {
+		    k = d[j] + i1
+		    d1 = Mems[k]
+		    if (d1 < dmin) {
+			dmin = d1; jmin = j; kmin = k
+		    }
+		}
+		if (keepids) {
+		    if (jmin < n1) {
+			Mems[kmin] = d1
+			Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+		    }
+		} else {
+		    if (jmin < n1)
+			Mems[kmin] = d1
+		}
+		n1 = n1 - 1 
+	    }
+
+	    # Reject the excess high points.
+	    do np = 1, nhigh {
+		k = d[1] + i1
+		d1 = Mems[k]
+		dmax = d1; jmax = 1; kmax = k
+		do j = 2, n1 {
+		    k = d[j] + i1
+		    d1 = Mems[k]
+		    if (d1 > dmax) {
+			dmax = d1; jmax = j; kmax = k
+		    }
+		}
+		if (keepids) {
+		    if (jmax < n1) {
+			Mems[kmax] = d1
+			Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+		    }
+		} else {
+		    if (jmax < n1)
+			Mems[kmax] = d1
+		}
+		n1 = n1 - 1 
+	    }
+	    n[i] = n1
+	}
+
+	if (dflag == D_ALL && npairs + nlow + nhigh > 0)
+		dflag = D_MIX
+end
+
+# IC_MM --  Reject a specified number of high and low pixels
+
+procedure ic_mmr (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
+real	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
+		d1 = Memr[k]
+		dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k
+		do j = 2, n1 {
+		    d2 = d1
+		    k = d[j] + i1
+		    d1 = Memr[k]
+		    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) {
+			    Memr[kmax] = d2
+			    Memi[m[jmax]+i1] = Memi[m[j]+i1]
+			} else {
+			    Memr[kmax] = d1
+			    Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+			}
+		    }
+		    if (jmin < j) {
+			if (jmax != n1) {
+			    Memr[kmin] = d1
+			    Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+			} else {
+			    Memr[kmin] = d2
+			    Memi[m[jmin]+i1] = Memi[m[j]+i1]
+			}
+		    }
+		} else {
+		    if (jmax < j) {
+			if (jmin != j)
+			    Memr[kmax] = d2
+			else
+			    Memr[kmax] = d1
+		    }
+		    if (jmin < j) {
+			if (jmax != n1)
+			    Memr[kmin] = d1
+			else
+			    Memr[kmin] = d2
+		    }
+		}
+		n1 = n1 - 2
+	    }
+
+	    # Reject the excess low points.
+	    do np = 1, nlow {
+		k = d[1] + i1
+		d1 = Memr[k]
+		dmin = d1; jmin = 1; kmin = k
+		do j = 2, n1 {
+		    k = d[j] + i1
+		    d1 = Memr[k]
+		    if (d1 < dmin) {
+			dmin = d1; jmin = j; kmin = k
+		    }
+		}
+		if (keepids) {
+		    if (jmin < n1) {
+			Memr[kmin] = d1
+			Memi[m[jmin]+i1] = Memi[m[n1]+i1]
+		    }
+		} else {
+		    if (jmin < n1)
+			Memr[kmin] = d1
+		}
+		n1 = n1 - 1 
+	    }
+
+	    # Reject the excess high points.
+	    do np = 1, nhigh {
+		k = d[1] + i1
+		d1 = Memr[k]
+		dmax = d1; jmax = 1; kmax = k
+		do j = 2, n1 {
+		    k = d[j] + i1
+		    d1 = Memr[k]
+		    if (d1 > dmax) {
+			dmax = d1; jmax = j; kmax = k
+		    }
+		}
+		if (keepids) {
+		    if (jmax < n1) {
+			Memr[kmax] = d1
+			Memi[m[jmax]+i1] = Memi[m[n1]+i1]
+		    }
+		} else {
+		    if (jmax < n1)
+			Memr[kmax] = d1
+		}
+		n1 = n1 - 1 
+	    }
+	    n[i] = n1
+	}
+
+	if (dflag == D_ALL && npairs + nlow + nhigh > 0)
+		dflag = D_MIX
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icombine.x b/noao/imred/ccdred/src/combine/generic/icombine.x
new file mode 100644
index 00000000..b4ff60be
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icombine.x
@@ -0,0 +1,607 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<imset.h>
+include	<error.h>
+include	<syserr.h>
+include	<mach.h>
+include	"../icombine.h"
+
+
+# ICOMBINE -- Combine images
+#
+# The memory and open file descriptor limits are checked and an attempt
+# to recover is made either by setting the image pixel files to be
+# closed after I/O or by notifying the calling program that memory
+# ran out and the IMIO buffer size should be reduced.  After the checks
+# a procedure for the selected combine option is called.
+# Because there may be several failure modes when reaching the file
+# limits we first assume an error is due to the file limit, except for
+# out of memory, and close some pixel files.  If the error then repeats
+# on accessing the pixels the error is passed back.
+
+
+procedure icombines (in, out, offsets, nimages, bufsize)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output images
+int	offsets[nimages,ARB]	# Input image offsets
+int	nimages			# Number of input images
+int	bufsize			# IMIO buffer size
+
+char	str[1]
+int	i, j, npts, fd, stropen(), errcode(), imstati()
+pointer	sp, d, id, n, m, lflag, scales, zeros, wts, dbuf
+pointer	buf, imgl1s(), impl1i()
+errchk	stropen, imgl1s, impl1i
+pointer	impl1r()
+errchk	impl1r
+
+include	"../icombine.com"
+
+begin
+	npts = IM_LEN(out[1],1)
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (d, nimages, TY_POINTER)
+	call salloc (id, nimages, TY_POINTER)
+	call salloc (n, npts, TY_INT)
+	call salloc (m, nimages, TY_POINTER)
+	call salloc (lflag, nimages, TY_INT)
+	call salloc (scales, nimages, TY_REAL)
+	call salloc (zeros, nimages, TY_REAL)
+	call salloc (wts, nimages, TY_REAL)
+	call amovki (D_ALL, Memi[lflag], nimages)
+
+	# If aligned use the IMIO buffer otherwise we need vectors of
+	# output length.
+
+	if (!aligned) {
+	    call salloc (dbuf, nimages, TY_POINTER)
+	    do i = 1, nimages
+		call salloc (Memi[dbuf+i-1], npts, TY_SHORT)
+	}
+
+	if (project) {
+	    call imseti (in[1], IM_NBUFS, nimages)
+	    call imseti (in[1], IM_BUFSIZE, bufsize)
+	    do i = 1, 3 {
+		if (out[i] != NULL)
+		    call imseti (out[i], IM_BUFSIZE, bufsize)
+	    }
+	} else {
+	    # Reserve FD for string operations.
+	    fd = stropen (str, 1, NEW_FILE)
+
+	    # Do I/O to the images.
+	    do i = 1, 3 {
+		if (out[i] != NULL)
+		    call imseti (out[i], IM_BUFSIZE, bufsize)
+	    }
+	    buf = impl1r (out[1])
+	    call aclrr (Memr[buf], npts)
+	    if (out[3] != NULL) {
+		buf = impl1r (out[3])
+		call aclrr (Memr[buf], npts)
+	    }
+	    if (out[2] != NULL) {
+		buf = impl1i (out[2])
+		call aclri (Memi[buf], npts)
+	    }
+
+	    do i = 1, nimages {
+		call imseti (in[i], IM_BUFSIZE, bufsize)
+		iferr (buf = imgl1s (in[i])) {
+		    switch (errcode()) {
+		    case SYS_MFULL:
+			call sfree (sp)
+			call strclose (fd)
+			call erract (EA_ERROR)
+		    case SYS_FTOOMANYFILES, SYS_IKIOPIX:
+			if (imstati (in[i], IM_CLOSEFD) == YES) {
+			    call sfree (sp)
+			    call strclose (fd)
+			    call erract (EA_ERROR)
+			}
+			do j = i-2, nimages
+			    call imseti (in[j], IM_CLOSEFD, YES)
+			buf = imgl1s (in[i])
+		    default:
+			call sfree (sp)
+			call strclose (fd)
+			call erract (EA_ERROR)
+		    }
+		}
+	    }
+
+	    call strclose (fd)
+	}
+
+	call ic_combines (in, out, Memi[dbuf], Memi[d], Memi[id], Memi[n],
+	    Memi[m], Memi[lflag], offsets, Memr[scales], Memr[zeros],
+	    Memr[wts], nimages, npts)
+end
+
+
+# IC_COMBINE -- Combine images.
+
+procedure ic_combines (in, out, dbuf, d, id, n, m, lflag, offsets,
+	scales, zeros, wts, nimages, npts)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output image
+pointer	dbuf[nimages]		# Data buffers for nonaligned images
+pointer	d[nimages]		# Data pointers
+pointer	id[nimages]		# Image index ID pointers
+int	n[npts]			# Number of good pixels
+pointer	m[nimages]		# Mask pointers
+int	lflag[nimages]		# Line flags
+int	offsets[nimages,ARB]	# Input image offsets
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero offset factors
+real	wts[nimages]		# Combining weights
+int	nimages			# Number of input images
+int	npts			# Number of points per output line
+
+int	i, ctor()
+real	r, imgetr()
+pointer	sp, v1, v2, v3, outdata, buf, nm, impnli()
+pointer	impnlr()
+errchk	ic_scale, imgetr
+
+include	"../icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (v3, IM_MAXDIM, TY_LONG)
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+	call ic_scale (in, out, offsets, scales, zeros, wts, nimages)
+
+	# Set combine parameters
+	switch (combine) {
+	case AVERAGE:
+	    if (dowts)
+		keepids = true
+	    else
+		keepids = false
+	case MEDIAN:
+	    dowts = false
+	    keepids = false
+	}
+	docombine = true
+
+	# Set rejection algorithm specific parameters
+	switch (reject) {
+	case CCDCLIP, CRREJECT:
+	    call salloc (nm, 3*nimages, TY_REAL)
+	    i = 1
+	    if (ctor (Memc[rdnoise], i, r) > 0) {
+		do i = 1, nimages
+		    Memr[nm+3*(i-1)] = r
+	    } else {
+		do i = 1, nimages
+		    Memr[nm+3*(i-1)] = imgetr (in[i], Memc[rdnoise])
+	    }
+	    i = 1
+	    if (ctor (Memc[gain], i, r) > 0) {
+		do i = 1, nimages {
+		    Memr[nm+3*(i-1)+1] = r
+		    Memr[nm+3*(i-1)] =
+			max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+		}
+	    } else {
+		do i = 1, nimages {
+		    r = imgetr (in[i], Memc[gain])
+		    Memr[nm+3*(i-1)+1] = r
+		    Memr[nm+3*(i-1)] =
+			max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+		}
+	    }
+	    i = 1
+	    if (ctor (Memc[snoise], i, r) > 0) {
+		do i = 1, nimages
+		    Memr[nm+3*(i-1)+2] = r
+	    } else {
+		do i = 1, nimages {
+		    r = imgetr (in[i], Memc[snoise])
+		    Memr[nm+3*(i-1)+2] = r
+		}
+	    }
+	    if (!keepids) {
+		if (doscale1 || grow > 0)
+		    keepids = true
+		else {
+		    do i = 2, nimages {
+			if (Memr[nm+3*(i-1)] != Memr[nm] ||
+			    Memr[nm+3*(i-1)+1] != Memr[nm+1] ||
+			    Memr[nm+3*(i-1)+2] != Memr[nm+2]) {
+			    keepids = true
+			    break
+			}
+		    }
+		}
+	    }
+	    if (reject == CRREJECT)
+		lsigma = MAX_REAL
+	case MINMAX:
+	    mclip = false
+	    if (grow > 0)
+		keepids = true
+	case PCLIP:
+	    mclip = true
+	    if (grow > 0)
+		keepids = true
+	case AVSIGCLIP, SIGCLIP:
+	    if (doscale1 || grow > 0)
+		keepids = true
+	case NONE:
+	    mclip = false
+	    grow = 0
+	}
+
+	if (keepids) {
+	    do i = 1, nimages
+		call salloc (id[i], npts, TY_INT)
+	}
+
+	while (impnlr (out[1], outdata, Meml[v1]) != EOF) {
+	    call ic_gdatas (in, out, dbuf, d, id, n, m, lflag, offsets,
+		scales, zeros, nimages, npts, Meml[v2], Meml[v3])
+
+	    switch (reject) {
+	    case CCDCLIP, CRREJECT:
+		if (mclip)
+		    call ic_mccdclips (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Memr[outdata])
+		else
+		    call ic_accdclips (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Memr[outdata])
+	    case MINMAX:
+		call ic_mms (d, id, n, npts)
+	    case PCLIP:
+		call ic_pclips (d, id, n, nimages, npts, Memr[outdata])
+	    case SIGCLIP:
+		if (mclip)
+		    call ic_msigclips (d, id, n, scales, zeros, nimages, npts,
+			Memr[outdata])
+		else
+		    call ic_asigclips (d, id, n, scales, zeros, nimages, npts,
+			Memr[outdata])
+	    case AVSIGCLIP:
+		if (mclip)
+		    call ic_mavsigclips (d, id, n, scales, zeros, nimages,
+			npts, Memr[outdata])
+		else
+		    call ic_aavsigclips (d, id, n, scales, zeros, nimages,
+			npts, Memr[outdata])
+	    }
+
+	    if (grow > 0)
+		call ic_grows (d, id, n, nimages, npts, Memr[outdata])
+
+	    if (docombine) {
+		switch (combine) {
+		case AVERAGE:
+		    call ic_averages (d, id, n, wts, npts, Memr[outdata])
+		case MEDIAN:
+		    call ic_medians (d, n, npts, Memr[outdata])
+		}
+	    }
+
+	    if (out[2] != NULL) {
+		call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		i = impnli (out[2], buf, Meml[v1])
+		call amovki (nimages, Memi[buf], npts)
+		call asubi (Memi[buf], n, Memi[buf], npts)
+	    }
+		    
+	    if (out[3] != NULL) {
+		call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		i = impnlr (out[3], buf, Meml[v1])
+		call ic_sigmas (d, id, n, wts, npts, Memr[outdata],
+		    Memr[buf])
+	    }
+	    call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+	}
+
+	call sfree (sp)
+end
+
+procedure icombiner (in, out, offsets, nimages, bufsize)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output images
+int	offsets[nimages,ARB]	# Input image offsets
+int	nimages			# Number of input images
+int	bufsize			# IMIO buffer size
+
+char	str[1]
+int	i, j, npts, fd, stropen(), errcode(), imstati()
+pointer	sp, d, id, n, m, lflag, scales, zeros, wts, dbuf
+pointer	buf, imgl1r(), impl1i()
+errchk	stropen, imgl1r, impl1i
+pointer	impl1r()
+errchk	impl1r
+
+include	"../icombine.com"
+
+begin
+	npts = IM_LEN(out[1],1)
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (d, nimages, TY_POINTER)
+	call salloc (id, nimages, TY_POINTER)
+	call salloc (n, npts, TY_INT)
+	call salloc (m, nimages, TY_POINTER)
+	call salloc (lflag, nimages, TY_INT)
+	call salloc (scales, nimages, TY_REAL)
+	call salloc (zeros, nimages, TY_REAL)
+	call salloc (wts, nimages, TY_REAL)
+	call amovki (D_ALL, Memi[lflag], nimages)
+
+	# If aligned use the IMIO buffer otherwise we need vectors of
+	# output length.
+
+	if (!aligned) {
+	    call salloc (dbuf, nimages, TY_POINTER)
+	    do i = 1, nimages
+		call salloc (Memi[dbuf+i-1], npts, TY_REAL)
+	}
+
+	if (project) {
+	    call imseti (in[1], IM_NBUFS, nimages)
+	    call imseti (in[1], IM_BUFSIZE, bufsize)
+	    do i = 1, 3 {
+		if (out[i] != NULL)
+		    call imseti (out[i], IM_BUFSIZE, bufsize)
+	    }
+	} else {
+	    # Reserve FD for string operations.
+	    fd = stropen (str, 1, NEW_FILE)
+
+	    # Do I/O to the images.
+	    do i = 1, 3 {
+		if (out[i] != NULL)
+		    call imseti (out[i], IM_BUFSIZE, bufsize)
+	    }
+	    buf = impl1r (out[1])
+	    call aclrr (Memr[buf], npts)
+	    if (out[3] != NULL) {
+		buf = impl1r (out[3])
+		call aclrr (Memr[buf], npts)
+	    }
+	    if (out[2] != NULL) {
+		buf = impl1i (out[2])
+		call aclri (Memi[buf], npts)
+	    }
+
+	    do i = 1, nimages {
+		call imseti (in[i], IM_BUFSIZE, bufsize)
+		iferr (buf = imgl1r (in[i])) {
+		    switch (errcode()) {
+		    case SYS_MFULL:
+			call sfree (sp)
+			call strclose (fd)
+			call erract (EA_ERROR)
+		    case SYS_FTOOMANYFILES, SYS_IKIOPIX:
+			if (imstati (in[i], IM_CLOSEFD) == YES) {
+			    call sfree (sp)
+			    call strclose (fd)
+			    call erract (EA_ERROR)
+			}
+			do j = i-2, nimages
+			    call imseti (in[j], IM_CLOSEFD, YES)
+			buf = imgl1r (in[i])
+		    default:
+			call sfree (sp)
+			call strclose (fd)
+			call erract (EA_ERROR)
+		    }
+		}
+	    }
+
+	    call strclose (fd)
+	}
+
+	call ic_combiner (in, out, Memi[dbuf], Memi[d], Memi[id], Memi[n],
+	    Memi[m], Memi[lflag], offsets, Memr[scales], Memr[zeros],
+	    Memr[wts], nimages, npts)
+end
+
+
+# IC_COMBINE -- Combine images.
+
+procedure ic_combiner (in, out, dbuf, d, id, n, m, lflag, offsets,
+	scales, zeros, wts, nimages, npts)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output image
+pointer	dbuf[nimages]		# Data buffers for nonaligned images
+pointer	d[nimages]		# Data pointers
+pointer	id[nimages]		# Image index ID pointers
+int	n[npts]			# Number of good pixels
+pointer	m[nimages]		# Mask pointers
+int	lflag[nimages]		# Line flags
+int	offsets[nimages,ARB]	# Input image offsets
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero offset factors
+real	wts[nimages]		# Combining weights
+int	nimages			# Number of input images
+int	npts			# Number of points per output line
+
+int	i, ctor()
+real	r, imgetr()
+pointer	sp, v1, v2, v3, outdata, buf, nm, impnli()
+pointer	impnlr()
+errchk	ic_scale, imgetr
+
+include	"../icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (v3, IM_MAXDIM, TY_LONG)
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+	call ic_scale (in, out, offsets, scales, zeros, wts, nimages)
+
+	# Set combine parameters
+	switch (combine) {
+	case AVERAGE:
+	    if (dowts)
+		keepids = true
+	    else
+		keepids = false
+	case MEDIAN:
+	    dowts = false
+	    keepids = false
+	}
+	docombine = true
+
+	# Set rejection algorithm specific parameters
+	switch (reject) {
+	case CCDCLIP, CRREJECT:
+	    call salloc (nm, 3*nimages, TY_REAL)
+	    i = 1
+	    if (ctor (Memc[rdnoise], i, r) > 0) {
+		do i = 1, nimages
+		    Memr[nm+3*(i-1)] = r
+	    } else {
+		do i = 1, nimages
+		    Memr[nm+3*(i-1)] = imgetr (in[i], Memc[rdnoise])
+	    }
+	    i = 1
+	    if (ctor (Memc[gain], i, r) > 0) {
+		do i = 1, nimages {
+		    Memr[nm+3*(i-1)+1] = r
+		    Memr[nm+3*(i-1)] =
+			max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+		}
+	    } else {
+		do i = 1, nimages {
+		    r = imgetr (in[i], Memc[gain])
+		    Memr[nm+3*(i-1)+1] = r
+		    Memr[nm+3*(i-1)] =
+			max ((Memr[nm+3*(i-1)] / r) ** 2, 1e4 / MAX_REAL)
+		}
+	    }
+	    i = 1
+	    if (ctor (Memc[snoise], i, r) > 0) {
+		do i = 1, nimages
+		    Memr[nm+3*(i-1)+2] = r
+	    } else {
+		do i = 1, nimages {
+		    r = imgetr (in[i], Memc[snoise])
+		    Memr[nm+3*(i-1)+2] = r
+		}
+	    }
+	    if (!keepids) {
+		if (doscale1 || grow > 0)
+		    keepids = true
+		else {
+		    do i = 2, nimages {
+			if (Memr[nm+3*(i-1)] != Memr[nm] ||
+			    Memr[nm+3*(i-1)+1] != Memr[nm+1] ||
+			    Memr[nm+3*(i-1)+2] != Memr[nm+2]) {
+			    keepids = true
+			    break
+			}
+		    }
+		}
+	    }
+	    if (reject == CRREJECT)
+		lsigma = MAX_REAL
+	case MINMAX:
+	    mclip = false
+	    if (grow > 0)
+		keepids = true
+	case PCLIP:
+	    mclip = true
+	    if (grow > 0)
+		keepids = true
+	case AVSIGCLIP, SIGCLIP:
+	    if (doscale1 || grow > 0)
+		keepids = true
+	case NONE:
+	    mclip = false
+	    grow = 0
+	}
+
+	if (keepids) {
+	    do i = 1, nimages
+		call salloc (id[i], npts, TY_INT)
+	}
+
+	while (impnlr (out[1], outdata, Meml[v1]) != EOF) {
+	    call ic_gdatar (in, out, dbuf, d, id, n, m, lflag, offsets,
+		scales, zeros, nimages, npts, Meml[v2], Meml[v3])
+
+	    switch (reject) {
+	    case CCDCLIP, CRREJECT:
+		if (mclip)
+		    call ic_mccdclipr (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Memr[outdata])
+		else
+		    call ic_accdclipr (d, id, n, scales, zeros, Memr[nm],
+			nimages, npts, Memr[outdata])
+	    case MINMAX:
+		call ic_mmr (d, id, n, npts)
+	    case PCLIP:
+		call ic_pclipr (d, id, n, nimages, npts, Memr[outdata])
+	    case SIGCLIP:
+		if (mclip)
+		    call ic_msigclipr (d, id, n, scales, zeros, nimages, npts,
+			Memr[outdata])
+		else
+		    call ic_asigclipr (d, id, n, scales, zeros, nimages, npts,
+			Memr[outdata])
+	    case AVSIGCLIP:
+		if (mclip)
+		    call ic_mavsigclipr (d, id, n, scales, zeros, nimages,
+			npts, Memr[outdata])
+		else
+		    call ic_aavsigclipr (d, id, n, scales, zeros, nimages,
+			npts, Memr[outdata])
+	    }
+
+	    if (grow > 0)
+		call ic_growr (d, id, n, nimages, npts, Memr[outdata])
+
+	    if (docombine) {
+		switch (combine) {
+		case AVERAGE:
+		    call ic_averager (d, id, n, wts, npts, Memr[outdata])
+		case MEDIAN:
+		    call ic_medianr (d, n, npts, Memr[outdata])
+		}
+	    }
+
+	    if (out[2] != NULL) {
+		call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		i = impnli (out[2], buf, Meml[v1])
+		call amovki (nimages, Memi[buf], npts)
+		call asubi (Memi[buf], n, Memi[buf], npts)
+	    }
+		    
+	    if (out[3] != NULL) {
+		call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		i = impnlr (out[3], buf, Meml[v1])
+		call ic_sigmar (d, id, n, wts, npts, Memr[outdata],
+		    Memr[buf])
+	    }
+	    call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+	}
+
+	call sfree (sp)
+end
+
diff --git a/noao/imred/ccdred/src/combine/generic/icpclip.x b/noao/imred/ccdred/src/combine/generic/icpclip.x
new file mode 100644
index 00000000..da09bb75
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icpclip.x
@@ -0,0 +1,442 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+define	MINCLIP		3	# Minimum number for clipping
+
+
+# IC_PCLIP -- Percentile clip
+#
+# 1) Find the median
+# 2) Find the pixel which is the specified order index away
+# 3) Use the data value difference as a sigma and apply clipping
+# 4) Since the median is known return it so it does not have to be recomputed
+
+procedure ic_pclips (d, m, n, nimages, npts, median)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image id pointers
+int	n[npts]			# Number of good pixels
+int	nimages			# Number of input images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, n4, n5, nl, nh, nin, maxkeep
+bool	even, fp_equalr()
+real	sigma, r, s, t
+pointer	sp, resid, mp1, mp2
+real	med
+
+include	"../icombine.com"
+
+begin
+	# There must be at least MINCLIP and more than nkeep pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Set sign of pclip parameter
+	if (pclip < 0)
+	    t = -1.
+	else
+	    t = 1.
+
+	# If there are no rejected pixels compute certain parameters once.
+	if (dflag == D_ALL) {
+	    n1 = n[1]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    n2 = 1 + n1 / 2
+	    even = (mod (n1, 2) == 0)
+	    if (pclip < 0.) {
+		if (even)
+		    n3 = max (1, nint (n2 - 1 + pclip))
+		else
+		    n3 = max (1, nint (n2 + pclip))
+	    } else
+		n3 = min (n1, nint (n2 + pclip))
+	    nin = n1
+	}
+
+	# Now apply clipping.
+	do i = 1, npts {
+	    # Compute median.
+	    if (dflag == D_MIX) {
+		n1 = n[i]
+		if (nkeep < 0)
+		    maxkeep = max (0, n1 + nkeep)
+		else
+		    maxkeep = min (n1, nkeep)
+		if (n1 == 0) {
+		    if (combine == MEDIAN)
+			median[i] = blank
+		    next
+		}
+		n2 = 1 + n1 / 2
+		even = (mod (n1, 2) == 0)
+		if (pclip < 0) {
+		    if (even)
+			n3 = max (1, nint (n2 - 1 + pclip))
+		    else
+			n3 = max (1, nint (n2 + pclip))
+		} else
+		    n3 = min (n1, nint (n2 + pclip))
+	    }
+
+	    j = i - 1 
+	    if (even) {
+		med = Mems[d[n2-1]+j]
+		med = (med + Mems[d[n2]+j]) / 2.
+	    } else
+		med = Mems[d[n2]+j]
+
+	    if (n1 < max (MINCLIP, maxkeep+1)) {
+		if (combine == MEDIAN)
+		    median[i] = med
+		next
+	    }
+
+	    # Define sigma for clipping
+	    sigma = t * (Mems[d[n3]+j] - med)
+	    if (fp_equalr (sigma, 0.)) {
+		if (combine == MEDIAN)
+		    median[i] = med
+		next
+	    }
+
+	    # Reject pixels and save residuals.
+	    # Check if any pixels are clipped.
+	    # If so recompute the median and reset the number of good pixels.
+	    # Only reorder if needed.
+
+	    for (nl=1; nl<=n1; nl=nl+1) {
+		r = (med - Mems[d[nl]+j]) / sigma
+		if (r < lsigma)
+		    break
+		Memr[resid+nl] = r
+	    }
+	    for (nh=n1; nh>=1; nh=nh-1) {
+		r = (Mems[d[nh]+j] -  med) / sigma
+		if (r < hsigma)
+		    break
+		Memr[resid+nh] = r
+	    }
+	    n4 = nh - nl + 1
+
+	    # If too many pixels are rejected add some back in.
+	    # All pixels with the same residual are added.
+	    while (n4 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n4 = nh - nl + 1
+	    }
+
+	    # If any pixels are rejected recompute the median.
+	    if (nl > 1 || nh < n1) {
+		n5 = nl + n4 / 2
+		if (mod (n4, 2) == 0) {
+		    med = Mems[d[n5-1]+j]
+		    med = (med + Mems[d[n5]+j]) / 2.
+		} else
+		    med = Mems[d[n5]+j]
+		n[i] = n4
+	    }
+	    if (combine == MEDIAN)
+		median[i] = med
+
+	    # Reorder if pixels only if necessary.
+	    if (nl > 1 && (combine != MEDIAN || grow > 0)) {
+		k = max (nl, n4 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Mems[d[l]+j] = Mems[d[k]+j]
+			if (grow > 0) {
+			    mp1 = m[l] + j
+			    mp2 = m[k] + j
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+j] = Memi[m[k]+j]
+			k = k + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Mems[d[l]+j] = Mems[d[k]+j]
+			k = k + 1
+		    }
+		}
+	    }
+	}
+
+	# Check if data flag needs to be reset for rejected pixels.
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag whether the median has been computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
+
+# IC_PCLIP -- Percentile clip
+#
+# 1) Find the median
+# 2) Find the pixel which is the specified order index away
+# 3) Use the data value difference as a sigma and apply clipping
+# 4) Since the median is known return it so it does not have to be recomputed
+
+procedure ic_pclipr (d, m, n, nimages, npts, median)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image id pointers
+int	n[npts]			# Number of good pixels
+int	nimages			# Number of input images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, n4, n5, nl, nh, nin, maxkeep
+bool	even, fp_equalr()
+real	sigma, r, s, t
+pointer	sp, resid, mp1, mp2
+real	med
+
+include	"../icombine.com"
+
+begin
+	# There must be at least MINCLIP and more than nkeep pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Set sign of pclip parameter
+	if (pclip < 0)
+	    t = -1.
+	else
+	    t = 1.
+
+	# If there are no rejected pixels compute certain parameters once.
+	if (dflag == D_ALL) {
+	    n1 = n[1]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    n2 = 1 + n1 / 2
+	    even = (mod (n1, 2) == 0)
+	    if (pclip < 0.) {
+		if (even)
+		    n3 = max (1, nint (n2 - 1 + pclip))
+		else
+		    n3 = max (1, nint (n2 + pclip))
+	    } else
+		n3 = min (n1, nint (n2 + pclip))
+	    nin = n1
+	}
+
+	# Now apply clipping.
+	do i = 1, npts {
+	    # Compute median.
+	    if (dflag == D_MIX) {
+		n1 = n[i]
+		if (nkeep < 0)
+		    maxkeep = max (0, n1 + nkeep)
+		else
+		    maxkeep = min (n1, nkeep)
+		if (n1 == 0) {
+		    if (combine == MEDIAN)
+			median[i] = blank
+		    next
+		}
+		n2 = 1 + n1 / 2
+		even = (mod (n1, 2) == 0)
+		if (pclip < 0) {
+		    if (even)
+			n3 = max (1, nint (n2 - 1 + pclip))
+		    else
+			n3 = max (1, nint (n2 + pclip))
+		} else
+		    n3 = min (n1, nint (n2 + pclip))
+	    }
+
+	    j = i - 1 
+	    if (even) {
+		med = Memr[d[n2-1]+j]
+		med = (med + Memr[d[n2]+j]) / 2.
+	    } else
+		med = Memr[d[n2]+j]
+
+	    if (n1 < max (MINCLIP, maxkeep+1)) {
+		if (combine == MEDIAN)
+		    median[i] = med
+		next
+	    }
+
+	    # Define sigma for clipping
+	    sigma = t * (Memr[d[n3]+j] - med)
+	    if (fp_equalr (sigma, 0.)) {
+		if (combine == MEDIAN)
+		    median[i] = med
+		next
+	    }
+
+	    # Reject pixels and save residuals.
+	    # Check if any pixels are clipped.
+	    # If so recompute the median and reset the number of good pixels.
+	    # Only reorder if needed.
+
+	    for (nl=1; nl<=n1; nl=nl+1) {
+		r = (med - Memr[d[nl]+j]) / sigma
+		if (r < lsigma)
+		    break
+		Memr[resid+nl] = r
+	    }
+	    for (nh=n1; nh>=1; nh=nh-1) {
+		r = (Memr[d[nh]+j] -  med) / sigma
+		if (r < hsigma)
+		    break
+		Memr[resid+nh] = r
+	    }
+	    n4 = nh - nl + 1
+
+	    # If too many pixels are rejected add some back in.
+	    # All pixels with the same residual are added.
+	    while (n4 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n4 = nh - nl + 1
+	    }
+
+	    # If any pixels are rejected recompute the median.
+	    if (nl > 1 || nh < n1) {
+		n5 = nl + n4 / 2
+		if (mod (n4, 2) == 0) {
+		    med = Memr[d[n5-1]+j]
+		    med = (med + Memr[d[n5]+j]) / 2.
+		} else
+		    med = Memr[d[n5]+j]
+		n[i] = n4
+	    }
+	    if (combine == MEDIAN)
+		median[i] = med
+
+	    # Reorder if pixels only if necessary.
+	    if (nl > 1 && (combine != MEDIAN || grow > 0)) {
+		k = max (nl, n4 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Memr[d[l]+j] = Memr[d[k]+j]
+			if (grow > 0) {
+			    mp1 = m[l] + j
+			    mp2 = m[k] + j
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+j] = Memi[m[k]+j]
+			k = k + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Memr[d[l]+j] = Memr[d[k]+j]
+			k = k + 1
+		    }
+		}
+	    }
+	}
+
+	# Check if data flag needs to be reset for rejected pixels.
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag whether the median has been computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icsclip.x b/noao/imred/ccdred/src/combine/generic/icsclip.x
new file mode 100644
index 00000000..d7ccfd84
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icsclip.x
@@ -0,0 +1,964 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+define	MINCLIP		3	# Mininum number of images for algorithm
+
+
+# IC_ASIGCLIP -- Reject pixels using sigma clipping about the average
+# The initial average rejects the high and low pixels.  A correction for
+# different scalings of the images may be made.  Weights are not used.
+
+procedure ic_asigclips (d, m, n, scales, zeros, nimages, npts, average)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+real	d1, low, high, sum, a, s, r, one
+data	one /1.0/
+pointer	sp, resid, w, wp, dp1, dp2, mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	} 
+	
+	# Flag whether returned average needs to be recomputed.
+	if (dowts || combine != AVERAGE)
+	    docombine = true
+	else
+	    docombine = false
+
+	# Save the residuals and the sigma scaling corrections if needed.
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+	if (doscale1)
+	    call salloc (w, nimages, TY_REAL)
+
+	# Do sigma clipping.
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+
+	    # If there are not enough pixels simply compute the average.
+	    if (n1 < max (3, maxkeep)) {
+		if (!docombine) {
+		    if (n1 == 0)
+			average[i] = blank
+		    else {
+			sum = Mems[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Mems[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    # Compute average with the high and low rejected.
+	    low = Mems[d[1]+k]
+	    high = Mems[d[2]+k]
+	    if (low > high) {
+		d1 = low
+		low = high
+		high = d1
+	    }
+	    sum = 0.
+	    do j = 3, n1 {
+	       d1 = Mems[d[j]+k]
+	       if (d1 < low) {
+		   sum = sum + low
+		   low = d1
+		} else if (d1 > high) {
+		    sum = sum + high
+		    high = d1
+		} else
+		    sum = sum + d1
+	    }
+	    a = sum / (n1 - 2)
+	    sum = sum + low + high
+
+	    # Iteratively reject pixels and compute the final average if needed.
+	    # Compact the data and keep track of the image IDs if needed.
+
+	    repeat {
+		n2 = n1
+		if (doscale1) {
+		    # Compute sigma corrected for scaling.
+		    s = 0.
+		    wp = w - 1
+		    do j = 1, n1 {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			wp = wp + 1
+
+			d1 = Mems[dp1]
+			l = Memi[mp1]
+			r = sqrt (max (one, (a + zeros[l]) / scales[l]))
+			s = s + ((d1 - a) / r) ** 2
+			Memr[wp] = r
+		    }
+		    s = sqrt (s / (n1 - 1))
+
+		    # Reject pixels.  Save the residuals and data values.
+		    wp = w - 1
+		    if (s > 0.) {
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    mp1 = m[j] + k
+			    wp = wp + 1
+
+			    d1 = Mems[dp1]
+			    r = (d1 - a) / (s * Memr[wp])
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs (r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    Memr[wp] = Memr[w+n1-1]
+				    mp2 = m[n1] + k
+				    l = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = l
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		} else {
+		    # Compute the sigma without scale correction.
+		    s = 0.
+		    do j = 1, n1
+			s = s + (Mems[d[j]+k] - a) ** 2
+		    s = sqrt (s / (n1 - 1))
+
+		    # Reject pixels.  Save the residuals and data values.
+		    if (s > 0.) {
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    d1 = Mems[dp1]
+			    r = (d1 - a) / s
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs (r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[n1] + k
+					l = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = l
+				    }
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		}
+
+		# Recompute the average.
+		if (n1 > 1)
+		    a = sum / n1
+	    } until (n1 == n2 || n1 <= max (2, maxkeep))
+
+	    # If too many pixels are rejected add some back.
+	    # All pixels with equal residuals are added back.
+	    if (n1 < maxkeep) {
+		nk = maxkeep
+		if (doscale1) {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Mems[dp1]
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Mems[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		} else {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Mems[dp1]
+				    Mems[dp1] = Mems[dp2]
+				    Mems[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Mems[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		}
+
+		# Recompute the average.
+		if (n1 > 1)
+		    a = sum / n1
+	    }
+
+	    # Save the average if needed.
+	    n[i] = n1
+	    if (!docombine) {
+		if (n1 > 0)
+		    average[i] = a
+		else
+		    average[i] = blank
+	    }
+	}
+
+	# Check if the data flag has to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# IC_MSIGCLIP -- Reject pixels using sigma clipping about the median
+
+procedure ic_msigclips (d, m, n, scales, zeros, nimages, npts, median)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+real	r, s
+pointer	sp, resid, w, mp1, mp2
+real	med, one
+data	one /1.0/
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	# Save the residuals and sigma scaling corrections if needed.
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+	if (doscale1)
+	    call salloc (w, nimages, TY_REAL)
+
+	# Compute median and sigma and iteratively clip.
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    nl = 1
+	    nh = n1
+
+	    repeat {
+		n2 = n1
+		n3 = nl + n1 / 2
+
+		if (n1 == 0)
+		    med = blank
+		else if (mod (n1, 2) == 0)
+		    med = (Mems[d[n3-1]+k] + Mems[d[n3]+k]) / 2.
+		else
+		    med = Mems[d[n3]+k]
+
+		if (n1 >= max (MINCLIP, maxkeep+1)) {
+		    if (doscale1) {
+			# Compute the sigma with scaling correction.
+			s = 0.
+			do j = nl, nh {
+			    l = Memi[m[j]+k]
+			    r = sqrt (max (one, (med + zeros[l]) / scales[l]))
+			    s = s + ((Mems[d[j]+k] - med) / r) ** 2
+			    Memr[w+j-1] = r
+			}
+			s = sqrt (s / (n1 - 1))
+
+			# Reject pixels and save the residuals.
+			if (s > 0.) {
+			    for (; nl <= n2; nl = nl + 1) {
+				r = (med - Mems[d[nl]+k]) / (s * Memr[w+nl-1])
+				if (r <= lsigma)
+				    break
+				Memr[resid+nl] = r
+				n1 = n1 - 1
+			    }
+			    for (; nh >= nl; nh = nh - 1) {
+				r = (Mems[d[nh]+k] - med) / (s * Memr[w+nh-1])
+				if (r <= hsigma)
+				    break
+				Memr[resid+nh] = r
+				n1 = n1 - 1
+			    }
+			}
+		    } else {
+			# Compute the sigma without scaling correction.
+			s = 0.
+			do j = nl, nh
+			    s = s + (Mems[d[j]+k] - med) ** 2
+			s = sqrt (s / (n1 - 1))
+
+			# Reject pixels and save the residuals.
+			if (s > 0.) {
+			    for (; nl <= n2; nl = nl + 1) {
+				r = (med - Mems[d[nl]+k]) / s
+				if (r <= lsigma)
+				    break
+				Memr[resid+nl] = r
+				n1 = n1 - 1
+			    }
+			    for (; nh >= nl; nh = nh - 1) {
+				r = (Mems[d[nh]+k] -  med) / s
+				if (r <= hsigma)
+				    break
+				Memr[resid+nh] = r
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    # If too many pixels are rejected add some back.
+	    # All pixels with equal residuals are added back.
+	    while (n1 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n1 = nh - nl + 1
+	    }
+
+	    # Only set median and reorder if needed
+	    n[i] = n1
+	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Mems[d[l]+k] = Mems[d[j]+k]
+			if (grow > 0) {
+			    mp1 = m[l] + k
+			    mp2 = m[j] + k
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+k] = Memi[m[j]+k]
+			j = j + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Mems[d[l]+k] = Mems[d[j]+k]
+			j = j + 1
+		    }
+		}
+	    }
+
+	    if (combine == MEDIAN)
+		median[i] = med
+	}
+
+	# Check if data flag needs to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag that the median has been computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
+
+# IC_ASIGCLIP -- Reject pixels using sigma clipping about the average
+# The initial average rejects the high and low pixels.  A correction for
+# different scalings of the images may be made.  Weights are not used.
+
+procedure ic_asigclipr (d, m, n, scales, zeros, nimages, npts, average)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+real	d1, low, high, sum, a, s, r, one
+data	one /1.0/
+pointer	sp, resid, w, wp, dp1, dp2, mp1, mp2
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	} 
+	
+	# Flag whether returned average needs to be recomputed.
+	if (dowts || combine != AVERAGE)
+	    docombine = true
+	else
+	    docombine = false
+
+	# Save the residuals and the sigma scaling corrections if needed.
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+	if (doscale1)
+	    call salloc (w, nimages, TY_REAL)
+
+	# Do sigma clipping.
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+
+	    # If there are not enough pixels simply compute the average.
+	    if (n1 < max (3, maxkeep)) {
+		if (!docombine) {
+		    if (n1 == 0)
+			average[i] = blank
+		    else {
+			sum = Memr[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Memr[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    # Compute average with the high and low rejected.
+	    low = Memr[d[1]+k]
+	    high = Memr[d[2]+k]
+	    if (low > high) {
+		d1 = low
+		low = high
+		high = d1
+	    }
+	    sum = 0.
+	    do j = 3, n1 {
+	       d1 = Memr[d[j]+k]
+	       if (d1 < low) {
+		   sum = sum + low
+		   low = d1
+		} else if (d1 > high) {
+		    sum = sum + high
+		    high = d1
+		} else
+		    sum = sum + d1
+	    }
+	    a = sum / (n1 - 2)
+	    sum = sum + low + high
+
+	    # Iteratively reject pixels and compute the final average if needed.
+	    # Compact the data and keep track of the image IDs if needed.
+
+	    repeat {
+		n2 = n1
+		if (doscale1) {
+		    # Compute sigma corrected for scaling.
+		    s = 0.
+		    wp = w - 1
+		    do j = 1, n1 {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			wp = wp + 1
+
+			d1 = Memr[dp1]
+			l = Memi[mp1]
+			r = sqrt (max (one, (a + zeros[l]) / scales[l]))
+			s = s + ((d1 - a) / r) ** 2
+			Memr[wp] = r
+		    }
+		    s = sqrt (s / (n1 - 1))
+
+		    # Reject pixels.  Save the residuals and data values.
+		    wp = w - 1
+		    if (s > 0.) {
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    mp1 = m[j] + k
+			    wp = wp + 1
+
+			    d1 = Memr[dp1]
+			    r = (d1 - a) / (s * Memr[wp])
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs (r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    Memr[wp] = Memr[w+n1-1]
+				    mp2 = m[n1] + k
+				    l = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = l
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		} else {
+		    # Compute the sigma without scale correction.
+		    s = 0.
+		    do j = 1, n1
+			s = s + (Memr[d[j]+k] - a) ** 2
+		    s = sqrt (s / (n1 - 1))
+
+		    # Reject pixels.  Save the residuals and data values.
+		    if (s > 0.) {
+			for (j=1; j<=n1; j=j+1) {
+			    dp1 = d[j] + k
+			    d1 = Memr[dp1]
+			    r = (d1 - a) / s
+			    if (r < -lsigma || r > hsigma) {
+				Memr[resid+n1] = abs (r)
+				if (j < n1) {
+				    dp2 = d[n1] + k
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[n1] + k
+					l = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = l
+				    }
+				    j = j - 1
+				}
+				sum = sum - d1
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		}
+
+		# Recompute the average.
+		if (n1 > 1)
+		    a = sum / n1
+	    } until (n1 == n2 || n1 <= max (2, maxkeep))
+
+	    # If too many pixels are rejected add some back.
+	    # All pixels with equal residuals are added back.
+	    if (n1 < maxkeep) {
+		nk = maxkeep
+		if (doscale1) {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			mp1 = m[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Memr[dp1]
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Memr[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		} else {
+		    for (j=n1+1; j<=nk; j=j+1) {
+			dp1 = d[j] + k
+			r = Memr[resid+j]
+			jj = 0
+			do l = j+1, n2 {
+			    s = Memr[resid+l]
+			    if (s < r + TOL) {
+				if (s > r - TOL)
+				    jj = jj + 1
+				else {
+				    jj = 0
+				    Memr[resid+l] = r
+				    r = s
+				    dp2 = d[l] + k
+				    d1 = Memr[dp1]
+				    Memr[dp1] = Memr[dp2]
+				    Memr[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Memr[dp1]
+			n1 = n1 + 1
+			nk = max (nk, j+jj)
+		    }
+		}
+
+		# Recompute the average.
+		if (n1 > 1)
+		    a = sum / n1
+	    }
+
+	    # Save the average if needed.
+	    n[i] = n1
+	    if (!docombine) {
+		if (n1 > 0)
+		    average[i] = a
+		else
+		    average[i] = blank
+	    }
+	}
+
+	# Check if the data flag has to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# IC_MSIGCLIP -- Reject pixels using sigma clipping about the median
+
+procedure ic_msigclipr (d, m, n, scales, zeros, nimages, npts, median)
+
+pointer	d[nimages]		# Data pointers
+pointer	m[nimages]		# Image id pointers
+int	n[npts]			# Number of good pixels
+real	scales[nimages]		# Scales
+real	zeros[nimages]		# Zeros
+int	nimages			# Number of images
+int	npts			# Number of output points per line
+real	median[npts]		# Median
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+real	r, s
+pointer	sp, resid, w, mp1, mp2
+real	med, one
+data	one /1.0/
+
+include	"../icombine.com"
+
+begin
+	# If there are insufficient pixels go on to the combining
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	# Save the residuals and sigma scaling corrections if needed.
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+	if (doscale1)
+	    call salloc (w, nimages, TY_REAL)
+
+	# Compute median and sigma and iteratively clip.
+	nin = n[1]
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    nl = 1
+	    nh = n1
+
+	    repeat {
+		n2 = n1
+		n3 = nl + n1 / 2
+
+		if (n1 == 0)
+		    med = blank
+		else if (mod (n1, 2) == 0)
+		    med = (Memr[d[n3-1]+k] + Memr[d[n3]+k]) / 2.
+		else
+		    med = Memr[d[n3]+k]
+
+		if (n1 >= max (MINCLIP, maxkeep+1)) {
+		    if (doscale1) {
+			# Compute the sigma with scaling correction.
+			s = 0.
+			do j = nl, nh {
+			    l = Memi[m[j]+k]
+			    r = sqrt (max (one, (med + zeros[l]) / scales[l]))
+			    s = s + ((Memr[d[j]+k] - med) / r) ** 2
+			    Memr[w+j-1] = r
+			}
+			s = sqrt (s / (n1 - 1))
+
+			# Reject pixels and save the residuals.
+			if (s > 0.) {
+			    for (; nl <= n2; nl = nl + 1) {
+				r = (med - Memr[d[nl]+k]) / (s * Memr[w+nl-1])
+				if (r <= lsigma)
+				    break
+				Memr[resid+nl] = r
+				n1 = n1 - 1
+			    }
+			    for (; nh >= nl; nh = nh - 1) {
+				r = (Memr[d[nh]+k] - med) / (s * Memr[w+nh-1])
+				if (r <= hsigma)
+				    break
+				Memr[resid+nh] = r
+				n1 = n1 - 1
+			    }
+			}
+		    } else {
+			# Compute the sigma without scaling correction.
+			s = 0.
+			do j = nl, nh
+			    s = s + (Memr[d[j]+k] - med) ** 2
+			s = sqrt (s / (n1 - 1))
+
+			# Reject pixels and save the residuals.
+			if (s > 0.) {
+			    for (; nl <= n2; nl = nl + 1) {
+				r = (med - Memr[d[nl]+k]) / s
+				if (r <= lsigma)
+				    break
+				Memr[resid+nl] = r
+				n1 = n1 - 1
+			    }
+			    for (; nh >= nl; nh = nh - 1) {
+				r = (Memr[d[nh]+k] -  med) / s
+				if (r <= hsigma)
+				    break
+				Memr[resid+nh] = r
+				n1 = n1 - 1
+			    }
+			}
+		    }
+		}
+	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))
+
+	    # If too many pixels are rejected add some back.
+	    # All pixels with equal residuals are added back.
+	    while (n1 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n1 = nh - nl + 1
+	    }
+
+	    # Only set median and reorder if needed
+	    n[i] = n1
+	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Memr[d[l]+k] = Memr[d[j]+k]
+			if (grow > 0) {
+			    mp1 = m[l] + k
+			    mp2 = m[j] + k
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+k] = Memi[m[j]+k]
+			j = j + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Memr[d[l]+k] = Memr[d[j]+k]
+			j = j + 1
+		    }
+		}
+	    }
+
+	    if (combine == MEDIAN)
+		median[i] = med
+	}
+
+	# Check if data flag needs to be reset for rejected pixels
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag that the median has been computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icsigma.x b/noao/imred/ccdred/src/combine/generic/icsigma.x
new file mode 100644
index 00000000..bc0d9788
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icsigma.x
@@ -0,0 +1,205 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	"../icombine.h"
+
+
+# IC_SIGMA -- Compute the sigma image line.
+# The estimated sigma includes a correction for the finite population.
+# Weights are used if desired.
+
+procedure ic_sigmas (d, m, n, wts, npts, average, sigma)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image ID pointers
+int	n[npts]			# Number of points
+real	wts[ARB]		# Weights
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+real	sigma[npts]		# Sigma line (returned)
+
+int	i, j, k, n1
+real	wt, sigcor, sumwt
+real	a, sum
+
+include	"../icombine.com"
+
+begin
+	if (dflag == D_ALL) {
+	    n1 = n[1]
+	    if (dowts) {
+		if (n1 > 1)
+		    sigcor = real (n1) / real (n1 - 1)
+		else
+		    sigcor = 1.
+		do i = 1, npts {
+		    k = i - 1
+		    a = average[i]
+		    wt = wts[Memi[m[1]+k]]
+		    sum = (Mems[d[1]+k] - a) ** 2 * wt
+		    do j = 2, n1 {
+			wt = wts[Memi[m[j]+k]]
+			sum = sum + (Mems[d[j]+k] - a) ** 2 * wt
+		    }
+		    sigma[i] = sqrt (sum * sigcor)
+		}
+	    } else {
+		if (n1 > 1)
+		    sigcor = 1. / real (n1 - 1)
+		else
+		    sigcor = 1.
+		do i = 1, npts {
+		    k = i - 1
+		    a = average[i]
+		    sum = (Mems[d[1]+k] - a) ** 2
+		    do j = 2, n1
+			sum = sum + (Mems[d[j]+k] - a) ** 2
+		    sigma[i] = sqrt (sum * sigcor)
+		}
+	    }
+	} else if (dflag == D_NONE) {
+	    do i = 1, npts
+		sigma[i] = blank
+	} else {
+	    if (dowts) {
+		do i = 1, npts {
+		    n1 = n[i]
+		    if (n1 > 0) {
+			k = i - 1
+			if (n1 > 1)
+			    sigcor = real (n1) / real (n1 -1)
+			else
+			    sigcor = 1
+			a = average[i]
+			wt = wts[Memi[m[1]+k]]
+			sum = (Mems[d[1]+k] - a) ** 2 * wt
+			sumwt = wt
+			do j = 2, n1 {
+			    wt = wts[Memi[m[j]+k]]
+			    sum = sum + (Mems[d[j]+k] - a) ** 2 * wt
+			    sumwt = sumwt + wt
+			}
+			sigma[i] = sqrt (sum / sumwt * sigcor)
+		    } else
+			sigma[i] = blank
+		}
+	    } else {
+		do i = 1, npts {
+		    n1 = n[i]
+		    if (n1 > 0) {
+			k = i - 1
+			if (n1 > 1)
+			    sigcor = 1. / real (n1 - 1)
+			else
+			    sigcor = 1.
+			a = average[i]
+			sum = (Mems[d[1]+k] - a) ** 2
+			do j = 2, n1
+			    sum = sum + (Mems[d[j]+k] - a) ** 2
+			sigma[i] = sqrt (sum * sigcor)
+		    } else
+			sigma[i] = blank
+		}
+	    }
+	}
+end
+
+# IC_SIGMA -- Compute the sigma image line.
+# The estimated sigma includes a correction for the finite population.
+# Weights are used if desired.
+
+procedure ic_sigmar (d, m, n, wts, npts, average, sigma)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image ID pointers
+int	n[npts]			# Number of points
+real	wts[ARB]		# Weights
+int	npts			# Number of output points per line
+real	average[npts]		# Average
+real	sigma[npts]		# Sigma line (returned)
+
+int	i, j, k, n1
+real	wt, sigcor, sumwt
+real	a, sum
+
+include	"../icombine.com"
+
+begin
+	if (dflag == D_ALL) {
+	    n1 = n[1]
+	    if (dowts) {
+		if (n1 > 1)
+		    sigcor = real (n1) / real (n1 - 1)
+		else
+		    sigcor = 1.
+		do i = 1, npts {
+		    k = i - 1
+		    a = average[i]
+		    wt = wts[Memi[m[1]+k]]
+		    sum = (Memr[d[1]+k] - a) ** 2 * wt
+		    do j = 2, n1 {
+			wt = wts[Memi[m[j]+k]]
+			sum = sum + (Memr[d[j]+k] - a) ** 2 * wt
+		    }
+		    sigma[i] = sqrt (sum * sigcor)
+		}
+	    } else {
+		if (n1 > 1)
+		    sigcor = 1. / real (n1 - 1)
+		else
+		    sigcor = 1.
+		do i = 1, npts {
+		    k = i - 1
+		    a = average[i]
+		    sum = (Memr[d[1]+k] - a) ** 2
+		    do j = 2, n1
+			sum = sum + (Memr[d[j]+k] - a) ** 2
+		    sigma[i] = sqrt (sum * sigcor)
+		}
+	    }
+	} else if (dflag == D_NONE) {
+	    do i = 1, npts
+		sigma[i] = blank
+	} else {
+	    if (dowts) {
+		do i = 1, npts {
+		    n1 = n[i]
+		    if (n1 > 0) {
+			k = i - 1
+			if (n1 > 1)
+			    sigcor = real (n1) / real (n1 -1)
+			else
+			    sigcor = 1
+			a = average[i]
+			wt = wts[Memi[m[1]+k]]
+			sum = (Memr[d[1]+k] - a) ** 2 * wt
+			sumwt = wt
+			do j = 2, n1 {
+			    wt = wts[Memi[m[j]+k]]
+			    sum = sum + (Memr[d[j]+k] - a) ** 2 * wt
+			    sumwt = sumwt + wt
+			}
+			sigma[i] = sqrt (sum / sumwt * sigcor)
+		    } else
+			sigma[i] = blank
+		}
+	    } else {
+		do i = 1, npts {
+		    n1 = n[i]
+		    if (n1 > 0) {
+			k = i - 1
+			if (n1 > 1)
+			    sigcor = 1. / real (n1 - 1)
+			else
+			    sigcor = 1.
+			a = average[i]
+			sum = (Memr[d[1]+k] - a) ** 2
+			do j = 2, n1
+			    sum = sum + (Memr[d[j]+k] - a) ** 2
+			sigma[i] = sqrt (sum * sigcor)
+		    } else
+			sigma[i] = blank
+		}
+	    }
+	}
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icsort.x b/noao/imred/ccdred/src/combine/generic/icsort.x
new file mode 100644
index 00000000..a39b68e2
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icsort.x
@@ -0,0 +1,550 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+define	LOGPTR	32			# log2(maxpts) (4e9)
+
+
+# IC_SORT -- Quicksort.  This is based on the VOPS asrt except that
+# the input is an array of pointers to image lines and the sort is done
+# across the image lines at each point along the lines.  The number of
+# valid pixels at each point is allowed to vary.  The cases of 1, 2, and 3
+# pixels per point are treated specially.
+
+procedure ic_sorts (a, b, nvecs, npts)
+
+pointer	a[ARB]			# pointer to input vectors
+short	b[ARB]			# work array
+int	nvecs[npts]		# number of vectors
+int	npts			# number of points in vectors
+
+short	pivot, temp, temp3
+int	i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR]
+define	swap {temp=$1;$1=$2;$2=temp}
+define	copy_	10
+
+begin
+	do l = 0, npts-1 {
+	    npix = nvecs[l+1]
+	    if (npix <= 1)
+		next
+
+	    do i = 1, npix
+		b[i] = Mems[a[i]+l]
+
+	    # Special cases
+	    if (npix <= 3) {
+		pivot = b[1]
+		temp = b[2]
+		if (npix == 2) {
+		    if (temp < pivot) {
+			b[1] = temp
+			b[2] = pivot
+		    } else
+			next
+		} else {
+		    temp3 = b[3]
+		    if (temp < pivot) {				# bac|bca|cba
+			if (temp < temp3) {			# bac|bca
+			    b[1] = temp
+			    if (pivot < temp3)			# bac
+				b[2] = pivot
+			    else {				# bca
+				b[2] = temp3
+				b[3] = pivot
+			    }
+			} else {				# cba
+			    b[1] = temp3
+			    b[3] = pivot
+			}
+		    } else if (temp3 < temp) {			# acb|cab
+			b[3] = temp
+			if (pivot < temp3)			# acb
+			    b[2] = temp3
+			else {					# cab
+			    b[1] = temp3
+			    b[2] = pivot
+			}
+		    } else
+			next
+		}
+		goto copy_
+	    }
+
+	    # General case
+	    do i = 1, npix
+		b[i] = Mems[a[i]+l]
+
+	    lv[1] = 1
+	    uv[1] = npix
+	    p = 1
+
+	    while (p > 0) {
+		if (lv[p] >= uv[p])		# only one elem in this subset
+		    p = p - 1			# pop stack
+		else {
+		    # Dummy do loop to trigger the Fortran optimizer.
+		    do p = p, ARB {
+			i = lv[p] - 1
+			j = uv[p]
+
+			# Select as the pivot the element at the center of the
+			# array, to avoid quadratic behavior on an already
+			# sorted array.
+
+			k = (lv[p] + uv[p]) / 2
+			swap (b[j], b[k])
+			pivot = b[j]		   # pivot line
+
+			while (i < j) {
+			    for (i=i+1;  b[i] < pivot;  i=i+1)
+				;
+			    for (j=j-1;  j > i;  j=j-1)
+				if (b[j] <= pivot)
+				    break
+			    if (i < j)		   # out of order pair
+				swap (b[i], b[j])  # interchange elements
+			}
+
+			j = uv[p]		   # move pivot to position i
+			swap (b[i], b[j])	   # interchange elements
+
+			if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+			    lv[p+1] = lv[p]
+			    uv[p+1] = i - 1
+			    lv[p] = i + 1
+			} else {
+			    lv[p+1] = i + 1
+			    uv[p+1] = uv[p]
+			    uv[p] = i - 1
+			}
+
+			break
+		    }
+		    p = p + 1			   # push onto stack
+		}
+	    }
+
+copy_
+	    do i = 1, npix
+		Mems[a[i]+l] = b[i]
+	}
+end
+
+
+# IC_2SORT -- Quicksort.  This is based on the VOPS asrt except that
+# the input is an array of pointers to image lines and the sort is done
+# across the image lines at each point along the lines.  The number of
+# valid pixels at each point is allowed to vary.  The cases of 1, 2, and 3
+# pixels per point are treated specially.  A second integer set of
+# vectors is sorted.
+
+procedure ic_2sorts (a, b, c, d, nvecs, npts)
+
+pointer	a[ARB]			# pointer to input vectors
+short	b[ARB]			# work array
+pointer	c[ARB]			# pointer to associated integer vectors
+int	d[ARB]			# work array
+int	nvecs[npts]		# number of vectors
+int	npts			# number of points in vectors
+
+short	pivot, temp, temp3
+int	i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR], itemp
+define	swap {temp=$1;$1=$2;$2=temp}
+define	iswap {itemp=$1;$1=$2;$2=itemp}
+define	copy_	10
+
+begin
+	do l = 0, npts-1 {
+	    npix = nvecs[l+1]
+	    if (npix <= 1)
+		next
+
+	    do i = 1, npix {
+		b[i] = Mems[a[i]+l]
+		d[i] = Memi[c[i]+l]
+	    }
+
+	    # Special cases
+	    if (npix <= 3) {
+		pivot = b[1]
+		temp = b[2]
+		if (npix == 2) {
+		    if (temp < pivot) {
+			b[1] = temp
+			b[2] = pivot
+			iswap (d[1], d[2])
+		    } else
+			next
+		} else {
+		    temp3 = b[3]
+		    if (temp < pivot) {				# bac|bca|cba
+			if (temp < temp3) {			# bac|bca
+			    b[1] = temp
+			    if (pivot < temp3) {		# bac
+				b[2] = pivot
+				iswap (d[1], d[2])
+			    } else {				# bca
+				b[2] = temp3
+				b[3] = pivot
+				itemp = d[2]
+				d[2] = d[3]
+				d[3] = d[1]
+				d[1] = itemp
+			    }
+			} else {				# cba
+			    b[1] = temp3
+			    b[3] = pivot
+			    iswap (d[1], d[3])
+			}
+		    } else if (temp3 < temp) {			# acb|cab
+			b[3] = temp
+			if (pivot < temp3) {			# acb
+			    b[2] = temp3
+			    iswap (d[2], d[3])
+			} else {				# cab
+			    b[1] = temp3
+			    b[2] = pivot
+			    itemp = d[2]
+			    d[2] = d[1]
+			    d[1] = d[3]
+			    d[3] = itemp
+			}
+		    } else
+			next
+		}
+		goto copy_
+	    }
+
+	    # General case
+	    lv[1] = 1
+	    uv[1] = npix
+	    p = 1
+
+	    while (p > 0) {
+		if (lv[p] >= uv[p])		# only one elem in this subset
+		    p = p - 1			# pop stack
+		else {
+		    # Dummy do loop to trigger the Fortran optimizer.
+		    do p = p, ARB {
+			i = lv[p] - 1
+			j = uv[p]
+
+			# Select as the pivot the element at the center of the
+			# array, to avoid quadratic behavior on an already
+			# sorted array.
+
+			k = (lv[p] + uv[p]) / 2
+			swap (b[j], b[k]); swap (d[j], d[k])
+			pivot = b[j]		   # pivot line
+
+			while (i < j) {
+			    for (i=i+1;  b[i] < pivot;  i=i+1)
+				;
+			    for (j=j-1;  j > i;  j=j-1)
+				if (b[j] <= pivot)
+				    break
+			    if (i < j) {	   # out of order pair
+				swap (b[i], b[j])  # interchange elements
+				swap (d[i], d[j])
+			    }
+			}
+
+			j = uv[p]		   # move pivot to position i
+			swap (b[i], b[j])	   # interchange elements
+			swap (d[i], d[j])
+
+			if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+			    lv[p+1] = lv[p]
+			    uv[p+1] = i - 1
+			    lv[p] = i + 1
+			} else {
+			    lv[p+1] = i + 1
+			    uv[p+1] = uv[p]
+			    uv[p] = i - 1
+			}
+
+			break
+		    }
+		    p = p + 1			   # push onto stack
+		}
+	    }
+
+copy_	   
+	    do i = 1, npix {
+		Mems[a[i]+l] = b[i]
+		Memi[c[i]+l] = d[i]
+	    }
+	}
+end
+
+# IC_SORT -- Quicksort.  This is based on the VOPS asrt except that
+# the input is an array of pointers to image lines and the sort is done
+# across the image lines at each point along the lines.  The number of
+# valid pixels at each point is allowed to vary.  The cases of 1, 2, and 3
+# pixels per point are treated specially.
+
+procedure ic_sortr (a, b, nvecs, npts)
+
+pointer	a[ARB]			# pointer to input vectors
+real	b[ARB]			# work array
+int	nvecs[npts]		# number of vectors
+int	npts			# number of points in vectors
+
+real	pivot, temp, temp3
+int	i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR]
+define	swap {temp=$1;$1=$2;$2=temp}
+define	copy_	10
+
+begin
+	do l = 0, npts-1 {
+	    npix = nvecs[l+1]
+	    if (npix <= 1)
+		next
+
+	    do i = 1, npix
+		b[i] = Memr[a[i]+l]
+
+	    # Special cases
+	    if (npix <= 3) {
+		pivot = b[1]
+		temp = b[2]
+		if (npix == 2) {
+		    if (temp < pivot) {
+			b[1] = temp
+			b[2] = pivot
+		    } else
+			next
+		} else {
+		    temp3 = b[3]
+		    if (temp < pivot) {				# bac|bca|cba
+			if (temp < temp3) {			# bac|bca
+			    b[1] = temp
+			    if (pivot < temp3)			# bac
+				b[2] = pivot
+			    else {				# bca
+				b[2] = temp3
+				b[3] = pivot
+			    }
+			} else {				# cba
+			    b[1] = temp3
+			    b[3] = pivot
+			}
+		    } else if (temp3 < temp) {			# acb|cab
+			b[3] = temp
+			if (pivot < temp3)			# acb
+			    b[2] = temp3
+			else {					# cab
+			    b[1] = temp3
+			    b[2] = pivot
+			}
+		    } else
+			next
+		}
+		goto copy_
+	    }
+
+	    # General case
+	    do i = 1, npix
+		b[i] = Memr[a[i]+l]
+
+	    lv[1] = 1
+	    uv[1] = npix
+	    p = 1
+
+	    while (p > 0) {
+		if (lv[p] >= uv[p])		# only one elem in this subset
+		    p = p - 1			# pop stack
+		else {
+		    # Dummy do loop to trigger the Fortran optimizer.
+		    do p = p, ARB {
+			i = lv[p] - 1
+			j = uv[p]
+
+			# Select as the pivot the element at the center of the
+			# array, to avoid quadratic behavior on an already
+			# sorted array.
+
+			k = (lv[p] + uv[p]) / 2
+			swap (b[j], b[k])
+			pivot = b[j]		   # pivot line
+
+			while (i < j) {
+			    for (i=i+1;  b[i] < pivot;  i=i+1)
+				;
+			    for (j=j-1;  j > i;  j=j-1)
+				if (b[j] <= pivot)
+				    break
+			    if (i < j)		   # out of order pair
+				swap (b[i], b[j])  # interchange elements
+			}
+
+			j = uv[p]		   # move pivot to position i
+			swap (b[i], b[j])	   # interchange elements
+
+			if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+			    lv[p+1] = lv[p]
+			    uv[p+1] = i - 1
+			    lv[p] = i + 1
+			} else {
+			    lv[p+1] = i + 1
+			    uv[p+1] = uv[p]
+			    uv[p] = i - 1
+			}
+
+			break
+		    }
+		    p = p + 1			   # push onto stack
+		}
+	    }
+
+copy_
+	    do i = 1, npix
+		Memr[a[i]+l] = b[i]
+	}
+end
+
+
+# IC_2SORT -- Quicksort.  This is based on the VOPS asrt except that
+# the input is an array of pointers to image lines and the sort is done
+# across the image lines at each point along the lines.  The number of
+# valid pixels at each point is allowed to vary.  The cases of 1, 2, and 3
+# pixels per point are treated specially.  A second integer set of
+# vectors is sorted.
+
+procedure ic_2sortr (a, b, c, d, nvecs, npts)
+
+pointer	a[ARB]			# pointer to input vectors
+real	b[ARB]			# work array
+pointer	c[ARB]			# pointer to associated integer vectors
+int	d[ARB]			# work array
+int	nvecs[npts]		# number of vectors
+int	npts			# number of points in vectors
+
+real	pivot, temp, temp3
+int	i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR], itemp
+define	swap {temp=$1;$1=$2;$2=temp}
+define	iswap {itemp=$1;$1=$2;$2=itemp}
+define	copy_	10
+
+begin
+	do l = 0, npts-1 {
+	    npix = nvecs[l+1]
+	    if (npix <= 1)
+		next
+
+	    do i = 1, npix {
+		b[i] = Memr[a[i]+l]
+		d[i] = Memi[c[i]+l]
+	    }
+
+	    # Special cases
+	    if (npix <= 3) {
+		pivot = b[1]
+		temp = b[2]
+		if (npix == 2) {
+		    if (temp < pivot) {
+			b[1] = temp
+			b[2] = pivot
+			iswap (d[1], d[2])
+		    } else
+			next
+		} else {
+		    temp3 = b[3]
+		    if (temp < pivot) {				# bac|bca|cba
+			if (temp < temp3) {			# bac|bca
+			    b[1] = temp
+			    if (pivot < temp3) {		# bac
+				b[2] = pivot
+				iswap (d[1], d[2])
+			    } else {				# bca
+				b[2] = temp3
+				b[3] = pivot
+				itemp = d[2]
+				d[2] = d[3]
+				d[3] = d[1]
+				d[1] = itemp
+			    }
+			} else {				# cba
+			    b[1] = temp3
+			    b[3] = pivot
+			    iswap (d[1], d[3])
+			}
+		    } else if (temp3 < temp) {			# acb|cab
+			b[3] = temp
+			if (pivot < temp3) {			# acb
+			    b[2] = temp3
+			    iswap (d[2], d[3])
+			} else {				# cab
+			    b[1] = temp3
+			    b[2] = pivot
+			    itemp = d[2]
+			    d[2] = d[1]
+			    d[1] = d[3]
+			    d[3] = itemp
+			}
+		    } else
+			next
+		}
+		goto copy_
+	    }
+
+	    # General case
+	    lv[1] = 1
+	    uv[1] = npix
+	    p = 1
+
+	    while (p > 0) {
+		if (lv[p] >= uv[p])		# only one elem in this subset
+		    p = p - 1			# pop stack
+		else {
+		    # Dummy do loop to trigger the Fortran optimizer.
+		    do p = p, ARB {
+			i = lv[p] - 1
+			j = uv[p]
+
+			# Select as the pivot the element at the center of the
+			# array, to avoid quadratic behavior on an already
+			# sorted array.
+
+			k = (lv[p] + uv[p]) / 2
+			swap (b[j], b[k]); swap (d[j], d[k])
+			pivot = b[j]		   # pivot line
+
+			while (i < j) {
+			    for (i=i+1;  b[i] < pivot;  i=i+1)
+				;
+			    for (j=j-1;  j > i;  j=j-1)
+				if (b[j] <= pivot)
+				    break
+			    if (i < j) {	   # out of order pair
+				swap (b[i], b[j])  # interchange elements
+				swap (d[i], d[j])
+			    }
+			}
+
+			j = uv[p]		   # move pivot to position i
+			swap (b[i], b[j])	   # interchange elements
+			swap (d[i], d[j])
+
+			if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+			    lv[p+1] = lv[p]
+			    uv[p+1] = i - 1
+			    lv[p] = i + 1
+			} else {
+			    lv[p+1] = i + 1
+			    uv[p+1] = uv[p]
+			    uv[p] = i - 1
+			}
+
+			break
+		    }
+		    p = p + 1			   # push onto stack
+		}
+	    }
+
+copy_	   
+	    do i = 1, npix {
+		Memr[a[i]+l] = b[i]
+		Memi[c[i]+l] = d[i]
+	    }
+	}
+end
diff --git a/noao/imred/ccdred/src/combine/generic/icstat.x b/noao/imred/ccdred/src/combine/generic/icstat.x
new file mode 100644
index 00000000..41512ccb
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/icstat.x
@@ -0,0 +1,444 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	"../icombine.h"
+
+define	NMAX	10000	# Maximum number of pixels to sample
+
+
+# IC_STAT -- Compute image statistics within specified section.
+# The image section is relative to a reference image which may be
+# different than the input image and may have an offset.  Only a
+# subsample of pixels is used.  Masked and thresholded pixels are
+# ignored.  Only the desired statistics are computed to increase
+# efficiency.
+
+procedure ic_stats (im, imref, section, offsets, image, nimages,
+	domode, domedian, domean, mode, median, mean)
+
+pointer	im			 # Data image
+pointer	imref			 # Reference image for image section
+char	section[ARB]		 # Image section
+int	offsets[nimages,ARB]	 # Image section offset from data to reference
+int	image			 # Image index (for mask I/O)
+int	nimages			 # Number of images in  offsets.
+bool	domode, domedian, domean # Statistics to compute
+real	mode, median, mean	 # Statistics
+
+int	i, j, ndim, n, nv
+real	a
+pointer	sp, v1, v2, dv, va, vb
+pointer	data, mask, dp, lp, mp, imgnls()
+short	ic_modes()
+real    asums()
+
+
+include	"../icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (dv, IM_MAXDIM, TY_LONG)
+	call salloc (va, IM_MAXDIM, TY_LONG)
+	call salloc (vb, IM_MAXDIM, TY_LONG)
+
+	# Determine the image section parameters.  This must be in terms of
+	# the data image pixel coordinates though the section may be specified
+	# in terms of the reference image coordinates.  Limit the number of
+	# pixels in each dimension to a maximum.
+
+	ndim = IM_NDIM(im)
+	if (project)
+	    ndim = ndim - 1
+	call amovki (1, Memi[v1], IM_MAXDIM)
+	call amovki (1, Memi[va], IM_MAXDIM)
+	call amovki (1, Memi[dv], IM_MAXDIM)
+	call amovi (IM_LEN(imref,1), Memi[vb], ndim)
+	call ic_section (section, Memi[va], Memi[vb], Memi[dv], ndim)
+	if (im != imref)
+	    do i = 1, ndim {
+		Memi[va+i-1] = Memi[va+i-1] - offsets[image,i]
+		Memi[vb+i-1] = Memi[vb+i-1] - offsets[image,i]
+	    }
+
+	do j = 1, 10 {
+	    n = 1
+	    do i = 0, ndim-1 {
+		Memi[v1+i] = max (1, min (Memi[va+i], Memi[vb+i]))
+		Memi[v2+i] = min (IM_LEN(im,i+1), max (Memi[va+i], Memi[vb+i]))
+		Memi[dv+i] = j
+		nv = max (1, (Memi[v2+i] - Memi[v1+i]) / Memi[dv+i] + 1)
+		Memi[v2+i] = Memi[v1+i] + (nv - 1) * Memi[dv+i]
+		n = n * nv
+	    }
+	    if (n < NMAX)
+		break
+	}
+
+	call amovl (Memi[v1], Memi[va], IM_MAXDIM)
+	Memi[va] = 1
+	if (project)
+	   Memi[va+ndim] = image
+	call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+
+	# Accumulate the pixel values within the section.  Masked pixels and
+	# thresholded pixels are ignored.
+
+	call salloc (data, n, TY_SHORT)
+	dp = data
+	while (imgnls (im, lp, Memi[vb]) != EOF) {
+	    call ic_mget1 (im, image, offsets[image,1], Memi[va], mask)
+	    lp = lp + Memi[v1] - 1
+	    if (dflag == D_ALL) {
+		if (dothresh) {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			a = Mems[lp]
+			if (a >= lthresh && a <= hthresh) {
+			    Mems[dp] = a
+			    dp = dp + 1
+			}
+			lp = lp + Memi[dv]
+		    }
+		} else {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			Mems[dp] = Mems[lp]
+			dp = dp + 1
+			lp = lp + Memi[dv]
+		    }
+		}
+	    } else if (dflag == D_MIX) {
+		mp = mask + Memi[v1] - 1
+		if (dothresh) {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			if (Memi[mp] == 0) {
+			    a = Mems[lp]
+			    if (a >= lthresh && a <= hthresh) {
+				Mems[dp] = a
+				dp = dp + 1
+			    }
+			}
+			mp = mp + Memi[dv]
+			lp = lp + Memi[dv]
+		    }
+		} else {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			if (Memi[mp] == 0) {
+			    Mems[dp] = Mems[lp]
+			    dp = dp + 1
+			}
+			mp = mp + Memi[dv]
+			lp = lp + Memi[dv]
+		    }
+		}
+	    }
+	    for (i=2; i<=ndim; i=i+1) {
+		Memi[va+i-1] = Memi[va+i-1] + Memi[dv+i-1]
+		if (Memi[va+i-1] <= Memi[v2+i-1])
+		    break
+		Memi[va+i-1] = Memi[v1+i-1]
+	    }
+	    if (i > ndim)
+		break
+	    call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+	}
+
+	n = dp - data
+	if (n < 1) {
+	    call sfree (sp)
+	    call error (1, "Image section contains no pixels")
+	}
+
+	# Compute only statistics needed.
+	if (domode || domedian) {
+	    call asrts (Mems[data], Mems[data], n)
+	    mode = ic_modes (Mems[data], n)
+	    median = Mems[data+n/2-1]
+	}
+	if (domean)
+	    mean = asums (Mems[data], n) / n
+
+	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.
+
+short procedure ic_modes (a, n)
+
+short	a[n]			# Data array
+int	n			# Number of points
+
+int	i, j, k, nmax
+real	z1, z2, zstep, zbin
+short	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)
+	zstep = max (1., zstep)
+	zbin = max (1., zbin)
+
+	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
+
+# IC_STAT -- Compute image statistics within specified section.
+# The image section is relative to a reference image which may be
+# different than the input image and may have an offset.  Only a
+# subsample of pixels is used.  Masked and thresholded pixels are
+# ignored.  Only the desired statistics are computed to increase
+# efficiency.
+
+procedure ic_statr (im, imref, section, offsets, image, nimages,
+	domode, domedian, domean, mode, median, mean)
+
+pointer	im			 # Data image
+pointer	imref			 # Reference image for image section
+char	section[ARB]		 # Image section
+int	offsets[nimages,ARB]	 # Image section offset from data to reference
+int	image			 # Image index (for mask I/O)
+int	nimages			 # Number of images in  offsets.
+bool	domode, domedian, domean # Statistics to compute
+real	mode, median, mean	 # Statistics
+
+int	i, j, ndim, n, nv
+real	a
+pointer	sp, v1, v2, dv, va, vb
+pointer	data, mask, dp, lp, mp, imgnlr()
+real	ic_moder()
+real    asumr()
+
+
+include	"../icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (dv, IM_MAXDIM, TY_LONG)
+	call salloc (va, IM_MAXDIM, TY_LONG)
+	call salloc (vb, IM_MAXDIM, TY_LONG)
+
+	# Determine the image section parameters.  This must be in terms of
+	# the data image pixel coordinates though the section may be specified
+	# in terms of the reference image coordinates.  Limit the number of
+	# pixels in each dimension to a maximum.
+
+	ndim = IM_NDIM(im)
+	if (project)
+	    ndim = ndim - 1
+	call amovki (1, Memi[v1], IM_MAXDIM)
+	call amovki (1, Memi[va], IM_MAXDIM)
+	call amovki (1, Memi[dv], IM_MAXDIM)
+	call amovi (IM_LEN(imref,1), Memi[vb], ndim)
+	call ic_section (section, Memi[va], Memi[vb], Memi[dv], ndim)
+	if (im != imref)
+	    do i = 1, ndim {
+		Memi[va+i-1] = Memi[va+i-1] - offsets[image,i]
+		Memi[vb+i-1] = Memi[vb+i-1] - offsets[image,i]
+	    }
+
+	do j = 1, 10 {
+	    n = 1
+	    do i = 0, ndim-1 {
+		Memi[v1+i] = max (1, min (Memi[va+i], Memi[vb+i]))
+		Memi[v2+i] = min (IM_LEN(im,i+1), max (Memi[va+i], Memi[vb+i]))
+		Memi[dv+i] = j
+		nv = max (1, (Memi[v2+i] - Memi[v1+i]) / Memi[dv+i] + 1)
+		Memi[v2+i] = Memi[v1+i] + (nv - 1) * Memi[dv+i]
+		n = n * nv
+	    }
+	    if (n < NMAX)
+		break
+	}
+
+	call amovl (Memi[v1], Memi[va], IM_MAXDIM)
+	Memi[va] = 1
+	if (project)
+	   Memi[va+ndim] = image
+	call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+
+	# Accumulate the pixel values within the section.  Masked pixels and
+	# thresholded pixels are ignored.
+
+	call salloc (data, n, TY_REAL)
+	dp = data
+	while (imgnlr (im, lp, Memi[vb]) != EOF) {
+	    call ic_mget1 (im, image, offsets[image,1], Memi[va], mask)
+	    lp = lp + Memi[v1] - 1
+	    if (dflag == D_ALL) {
+		if (dothresh) {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			a = Memr[lp]
+			if (a >= lthresh && a <= hthresh) {
+			    Memr[dp] = a
+			    dp = dp + 1
+			}
+			lp = lp + Memi[dv]
+		    }
+		} else {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			Memr[dp] = Memr[lp]
+			dp = dp + 1
+			lp = lp + Memi[dv]
+		    }
+		}
+	    } else if (dflag == D_MIX) {
+		mp = mask + Memi[v1] - 1
+		if (dothresh) {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			if (Memi[mp] == 0) {
+			    a = Memr[lp]
+			    if (a >= lthresh && a <= hthresh) {
+				Memr[dp] = a
+				dp = dp + 1
+			    }
+			}
+			mp = mp + Memi[dv]
+			lp = lp + Memi[dv]
+		    }
+		} else {
+		    do i = Memi[v1], Memi[v2], Memi[dv] {
+			if (Memi[mp] == 0) {
+			    Memr[dp] = Memr[lp]
+			    dp = dp + 1
+			}
+			mp = mp + Memi[dv]
+			lp = lp + Memi[dv]
+		    }
+		}
+	    }
+	    for (i=2; i<=ndim; i=i+1) {
+		Memi[va+i-1] = Memi[va+i-1] + Memi[dv+i-1]
+		if (Memi[va+i-1] <= Memi[v2+i-1])
+		    break
+		Memi[va+i-1] = Memi[v1+i-1]
+	    }
+	    if (i > ndim)
+		break
+	    call amovl (Memi[va], Memi[vb], IM_MAXDIM)
+	}
+
+	n = dp - data
+	if (n < 1) {
+	    call sfree (sp)
+	    call error (1, "Image section contains no pixels")
+	}
+
+	# 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
+
+	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
+
diff --git a/noao/imred/ccdred/src/combine/generic/mkpkg b/noao/imred/ccdred/src/combine/generic/mkpkg
new file mode 100644
index 00000000..63695459
--- /dev/null
+++ b/noao/imred/ccdred/src/combine/generic/mkpkg
@@ -0,0 +1,23 @@
+# Make CCDRED Package.
+
+$checkout libpkg.a ../../..
+$update   libpkg.a
+$checkin  libpkg.a ../../..
+$exit
+
+libpkg.a:
+	icaclip.x	../icombine.com ../icombine.h
+	icaverage.x	../icombine.com ../icombine.h <imhdr.h>
+	iccclip.x	../icombine.com ../icombine.h
+	icgdata.x	../icombine.com ../icombine.h <imhdr.h> <mach.h>
+	icgrow.x	../icombine.com ../icombine.h
+	icmedian.x	../icombine.com ../icombine.h
+	icmm.x		../icombine.com ../icombine.h
+	icombine.x	../icombine.com ../icombine.h <error.h> <syserr.h>\
+			<imhdr.h> <imset.h> <mach.h>
+	icpclip.x	../icombine.com ../icombine.h
+	icsclip.x	../icombine.com ../icombine.h
+	icsigma.x	../icombine.com ../icombine.h <imhdr.h>
+	icsort.x	
+	icstat.x	../icombine.com ../icombine.h <imhdr.h>
+	;