Initial commit

10 files changed, 2460 insertions, 0 deletions

diff --git a/noao/onedspec/scombine/generic/icaclip.x b/noao/onedspec/scombine/generic/icaclip.x
new file mode 100644
index 00000000..41432dd7
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icaclip.x
@@ -0,0 +1,555 @@
+# 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_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/onedspec/scombine/generic/icaverage.x b/noao/onedspec/scombine/generic/icaverage.x
new file mode 100644
index 00000000..6c5c870b
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icaverage.x
@@ -0,0 +1,84 @@
+# 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_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/onedspec/scombine/generic/iccclip.x b/noao/onedspec/scombine/generic/iccclip.x
new file mode 100644
index 00000000..26b17ba2
--- /dev/null
+++ b/noao/onedspec/scombine/generic/iccclip.x
@@ -0,0 +1,453 @@
+# 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_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/onedspec/scombine/generic/icgrow.x b/noao/onedspec/scombine/generic/icgrow.x
new file mode 100644
index 00000000..074bd8c3
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icgrow.x
@@ -0,0 +1,76 @@
+# 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_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/onedspec/scombine/generic/icmedian.x b/noao/onedspec/scombine/generic/icmedian.x
new file mode 100644
index 00000000..e7607340
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icmedian.x
@@ -0,0 +1,139 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+
+# 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, j1, j2, j3, k, n1
+bool	even
+real	val1, val2, val3
+
+include	"../icombine.com"
+
+begin
+	if (dflag == D_NONE) {
+	    do i = 1, npts
+		median[i]= blank
+	    return
+	}
+
+	# Check for previous sorting
+	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
+	}
+
+	# Repeatedly exchange the extreme values until there are three
+	# or fewer pixels.
+
+	do i = 1, npts {
+	    k = i - 1
+	    n1 = n[i]
+	    while (n1 > 3) {
+		j1 = 1
+		j2 = 1
+	        val1 = Memr[d[j1]+k]
+	        val2 = val1
+	        do j3 = 2, n1 {
+		    val3 = Memr[d[j3]+k]
+	            if (val3 > val1) {
+		        j1 = j3
+		        val1 = val3
+		    } else if (val3 < val2) {
+		        j2 = j3
+		        val2 = val3
+		    }
+	        }
+		j3 = n1 - 1
+		if (j1 < j3 && j2 < j3) {
+		    Memr[d[j1]+k] = val3
+		    Memr[d[j2]+k] = Memr[d[j3]+k]
+		    Memr[d[j3]+k] = val1
+		    Memr[d[n1]+k] = val2
+		} else if (j1 < j3) {
+		    if (j2 == j3) {
+			Memr[d[j1]+k] = val3
+			Memr[d[n1]+k] = val1
+		    } else {
+			Memr[d[j1]+k] = Memr[d[j3]+k]
+			Memr[d[j3]+k] = val1
+		    }
+		} else if (j2 < j3) {
+		    if (j1 == j3) {
+			Memr[d[j2]+k] = val3
+			Memr[d[n1]+k] = val2
+		    } else {
+			Memr[d[j2]+k] = Memr[d[j3]+k]
+			Memr[d[j3]+k] = val2
+		    }
+		}
+		n1 = n1 - 2
+	    }
+
+	    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
+	        }
+	    } else if (n1 == 2) {
+		val1 = Memr[d[1]+k]
+		val2 = Memr[d[2]+k]
+		median[i] = (val1 + val2) / 2
+	    } else if (n1 == 1)
+		median[i] = Memr[d[1]+k]
+	    else
+		median[i] = blank
+	}
+end
diff --git a/noao/onedspec/scombine/generic/icmm.x b/noao/onedspec/scombine/generic/icmm.x
new file mode 100644
index 00000000..1c314241
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icmm.x
@@ -0,0 +1,152 @@
+# 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_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/onedspec/scombine/generic/icpclip.x b/noao/onedspec/scombine/generic/icpclip.x
new file mode 100644
index 00000000..d9028a93
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icpclip.x
@@ -0,0 +1,224 @@
+# 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_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/onedspec/scombine/generic/icsclip.x b/noao/onedspec/scombine/generic/icsclip.x
new file mode 100644
index 00000000..e38f7935
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icsclip.x
@@ -0,0 +1,486 @@
+# 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_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/onedspec/scombine/generic/icsort.x b/noao/onedspec/scombine/generic/icsort.x
new file mode 100644
index 00000000..f3d2fb21
--- /dev/null
+++ b/noao/onedspec/scombine/generic/icsort.x
@@ -0,0 +1,275 @@
+# 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_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]
+			    }
+			} 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/onedspec/scombine/generic/mkpkg b/noao/onedspec/scombine/generic/mkpkg
new file mode 100644
index 00000000..4d371363
--- /dev/null
+++ b/noao/onedspec/scombine/generic/mkpkg
@@ -0,0 +1,16 @@
+$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
+	icgrow.x	../icombine.com ../icombine.h
+	icmedian.x	../icombine.com ../icombine.h
+	icmm.x		../icombine.com ../icombine.h
+	icpclip.x	../icombine.com ../icombine.h
+	icsclip.x	../icombine.com ../icombine.h
+	icsort.x	
+	;