Initial commit

25 files changed, 4676 insertions, 0 deletions

diff --git a/noao/onedspec/scombine/README b/noao/onedspec/scombine/README
new file mode 100644
index 00000000..e4031e6a
--- /dev/null
+++ b/noao/onedspec/scombine/README
@@ -0,0 +1,17 @@
+SCOMBINE -- Combine spectra
+
+This routine is based in large part on IMCOMBINE.  The routines in the
+generic directory are identical to those in that task except that they
+only contain routines for real data.  The ic routines in this directory
+are similar though modified from IMCOMBINE.
+
+The iscombine files are for an interactive combine task based on work
+by CTIO.  Because it is limited currently to linear spectra and is not
+organized to take advantage of the IMCOMBINE options it is not installed.
+A version of this may someday be added based on the current software.
+
+
+=======
+
+This version was renamed to OSCOMBINE.  It is obsolete and may be removed
+at some future time.  (4/14/04, Valdes)
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	
+	;
diff --git a/noao/onedspec/scombine/icgdata.x b/noao/onedspec/scombine/icgdata.x
new file mode 100644
index 00000000..907adc5e
--- /dev/null
+++ b/noao/onedspec/scombine/icgdata.x
@@ -0,0 +1,199 @@
+include	<smw.h>
+include "icombine.h"
+
+
+# IC_GDATAR - Apply threshold, scaling, and masking
+
+procedure ic_gdatar (sh, d, id, n, m, lflag, scales, zeros, nimages, npts)
+
+pointer	sh[nimages]		# Input spectra structures
+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
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero offset factors
+int	nimages			# Number of spectra
+int	npts			# NUmber of output points
+
+int	i, j, k, l, nused
+real	a, b
+pointer	dp, ip, mp
+
+include	"icombine.com"
+	
+begin
+	# Set data vectors
+	do i = 1, nimages {
+	    d[i] = SY(sh[i])
+	    m[i] = SX(sh[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) {
+			    Memr[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 (Memr[mp] == 0) {
+			    a = Memr[dp]
+			    if (a < lthresh || a > hthresh) {
+				Memr[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 (Memr[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 (Memr[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 (Memr[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 = 0, npts-1 {
+			do i = n[i]+1, nimages
+			    Memi[id[i]+j] = 0
+		    }
+		}
+	    } else {
+		do i = 1, nused {
+		    dp = d[i]
+		    mp = m[i]
+		    do j = 1, npts {
+			if (Memr[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/onedspec/scombine/iclog.x b/noao/onedspec/scombine/iclog.x
new file mode 100644
index 00000000..29002c0f
--- /dev/null
+++ b/noao/onedspec/scombine/iclog.x
@@ -0,0 +1,301 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<smw.h>
+include	"icombine.h"
+
+# IC_LOG -- Output log information is a log file has been specfied.
+
+procedure ic_log (sh, shout, ncombine, exptime, sname, zname, wname,
+	mode, median, mean, scales, zeros, wts, nimages,
+	dozero, nout, expname, exposure)
+
+pointer	sh[nimages]		# Input spectra
+pointer	shout			# Output spectrum
+int	ncombine[nimages]	# Number of previous combined images
+real	exptime[nimages]	# Exposure times
+char	sname[ARB]              # Scale name
+char	zname[ARB]              # Zero name
+char	wname[ARB]              # Weight name
+real	mode[nimages]		# Modes
+real	median[nimages]		# Medians
+real	mean[nimages]		# Means
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero or sky levels
+real	wts[nimages]		# Weights
+int	nimages			# Number of images
+bool	dozero			# Zero flag
+int	nout			# Number of images combined in output
+char	expname[ARB]		# Exposure name
+real	exposure		# Output exposure
+
+int	i, j, ctor()
+real	rval
+long	clktime()
+bool	prncombine, prexptime, prmode, prmedian, prmean
+bool	prrdn, prgain, prsn, prscale, przero, prwts, strne()
+pointer	sp, fname
+
+include	"icombine.com"
+
+begin
+	if (logfd == NULL)
+	    return
+
+	call smark (sp)
+	call salloc (fname, SZ_LINE, TY_CHAR)
+
+	# Time stamp the log and print parameter information.
+
+	call cnvdate (clktime(0), Memc[fname], SZ_LINE)
+	call fprintf (logfd, "\n%s: SCOMBINE\n")
+	    call pargstr (Memc[fname])
+	switch (combine) {
+	case AVERAGE:
+	    call fprintf (logfd, "  combine = average,")
+	case MEDIAN:
+	    call fprintf (logfd, "  combine = median,")
+	case SUM:
+	    call fprintf (logfd, "  combine = sum\n")
+	}
+	if (combine != SUM) {
+	    call fprintf (logfd, " scale = %s, zero = %s, weight = %s\n")
+		call pargstr (sname)
+		call pargstr (zname)
+		call pargstr (wname)
+	}
+
+	switch (reject) {
+	case MINMAX:
+	    call fprintf (logfd, "  reject = minmax, nlow = %d, nhigh = %d\n")
+		call pargi (nint (flow * nimages))
+		call pargi (nint (fhigh * nimages))
+	case CCDCLIP:
+	    call fprintf (logfd, "  reject = ccdclip, mclip = %b, nkeep = %d\n")
+		call pargb (mclip)
+		call pargi (nkeep)
+	    call fprintf (logfd,
+	    "  rdnoise = %s, gain = %s, snoise = %s, sigma = %g, hsigma = %g\n")
+		call pargstr (Memc[rdnoise])
+		call pargstr (Memc[gain])
+		call pargstr (Memc[snoise])
+		call pargr (lsigma)
+		call pargr (hsigma)
+	case CRREJECT:
+	    call fprintf (logfd,
+		"  reject = crreject, mclip = %b, nkeep = %d\n")
+		call pargb (mclip)
+		call pargi (nkeep)
+	    call fprintf (logfd,
+		"  rdnoise = %s, gain = %s, snoise = %s, hsigma = %g\n")
+		call pargstr (Memc[rdnoise])
+		call pargstr (Memc[gain])
+		call pargstr (Memc[snoise])
+		call pargr (hsigma)
+	case PCLIP:
+	    call fprintf (logfd, "  reject = pclip, nkeep = %d\n")
+		call pargi (nkeep)
+	    call fprintf (logfd, "  pclip = %g, lsigma = %g, hsigma = %g\n")
+		call pargr (pclip)
+		call pargr (lsigma)
+		call pargr (hsigma)
+	case SIGCLIP:
+	    call fprintf (logfd, "  reject = sigclip, mclip = %b, nkeep = %d\n")
+		call pargb (mclip)
+		call pargi (nkeep)
+	    call fprintf (logfd, "  lsigma = %g, hsigma = %g\n")
+		call pargr (lsigma)
+		call pargr (hsigma)
+	case AVSIGCLIP:
+	    call fprintf (logfd,
+		"  reject = avsigclip, mclip = %b, nkeep = %d\n")
+		call pargb (mclip)
+		call pargi (nkeep)
+	    call fprintf (logfd, "  lsigma = %g, hsigma = %g\n")
+		call pargr (lsigma)
+		call pargr (hsigma)
+	}
+	if (reject != NONE && grow > 0) {
+	    call fprintf (logfd, "  grow = %d\n")
+		call pargi (grow)
+	}
+	if (dothresh) {
+	    if (lthresh > -MAX_REAL && hthresh < MAX_REAL) {
+		call fprintf (logfd, "  lthreshold = %g, hthreshold = %g\n")
+		    call pargr (lthresh)
+		    call pargr (hthresh)
+	    } else if (lthresh > -MAX_REAL) {
+		call fprintf (logfd, "  lthreshold = %g\n")
+		    call pargr (lthresh)
+	    } else {
+		call fprintf (logfd, "  hthreshold = %g\n")
+		    call pargr (hthresh)
+	    }
+	}
+	call fprintf (logfd, "  blank = %g\n")
+	    call pargr (blank)
+	call clgstr ("sample", Memc[fname], SZ_LINE)
+	if (Memc[fname] != EOS) {
+	    call fprintf (logfd, "  sample = %s\n")
+		call pargstr (Memc[fname])
+	}
+
+	# Print information pertaining to individual images as a set of
+	# columns with the image name being the first column.  Determine
+	# what information is relevant and print the appropriate header.
+
+	prncombine = false
+	prexptime = (expname[1] != EOS)
+	prscale = (doscale || strne (sname, "none"))
+	przero = (dozero || strne (zname, "none"))
+	prwts = (dowts || strne (wname, "none"))
+	prmode = false
+	prmedian = false
+	prmean = false
+	prrdn = false
+	prgain = false
+	prsn = false
+	do i = 1, nimages {
+	    if (ncombine[i] != ncombine[1])
+		prncombine = true
+	    if (exptime[i] != exptime[1])
+		prexptime = true
+	    if (mode[i] != mode[1])
+		prmode = true
+	    if (median[i] != median[1])
+		prmedian = true
+	    if (mean[i] != mean[1])
+		prmean = true
+	    if (reject == CCDCLIP || reject == CRREJECT) {
+		j = 1
+		if (ctor (Memc[rdnoise], j, rval) == 0)
+		    prrdn = true
+		j = 1
+		if (ctor (Memc[gain], j, rval) == 0)
+		    prgain = true
+		j = 1
+		if (ctor (Memc[snoise], j, rval) == 0)
+		    prsn = true
+	    }
+	}
+
+	call fprintf (logfd, "  %20s ")
+	    call pargstr ("Images")
+	if (prncombine) {
+	    call fprintf (logfd, " %6s")
+		call pargstr ("N")
+	}
+	if (prexptime) {
+	    call fprintf (logfd, " %6s")
+		call pargstr ("Exp")
+	}
+	if (prmode) {
+	    call fprintf (logfd, " %7s")
+		call pargstr ("Mode")
+	}
+	if (prmedian) {
+	    call fprintf (logfd, " %7s")
+		call pargstr ("Median")
+	}
+	if (prmean) {
+	    call fprintf (logfd, " %7s")
+		call pargstr ("Mean")
+	}
+	if (prrdn) {
+	    call fprintf (logfd, " %7s")
+		call pargstr ("Rdnoise")
+	}
+	if (prgain) {
+	    call fprintf (logfd, " %6s")
+		call pargstr ("Gain")
+	}
+	if (prsn) {
+	    call fprintf (logfd, " %6s")
+		call pargstr ("Snoise")
+	}
+	if (prscale) {
+	    call fprintf (logfd, " %6s")
+		call pargstr ("Scale")
+	}
+	if (przero) {
+	    call fprintf (logfd, " %7s")
+		call pargstr ("Zero")
+	}
+	if (prwts) {
+	    call fprintf (logfd, " %6s")
+		call pargstr ("Weight")
+	}
+	call fprintf (logfd, "\n")
+
+	do i = 1, nimages {
+	    call fprintf (logfd, "  %16s[%3d]")
+		call pargstr (IMNAME(sh[i]))
+		call pargi (AP(sh[i]))
+	    if (prncombine) {
+		call fprintf (logfd, " %6d")
+		    call pargi (ncombine[i])
+	    }
+	    if (prexptime) {
+		call fprintf (logfd, " %6.1f")
+		    call pargr (exptime[i])
+	    }
+	    if (prmode) {
+		call fprintf (logfd, " %7.5g")
+		    call pargr (mode[i])
+	    }
+	    if (prmedian) {
+		call fprintf (logfd, " %7.5g")
+		    call pargr (median[i])
+	    }
+	    if (prmean) {
+		call fprintf (logfd, " %7.5g")
+		    call pargr (mean[i])
+	    }
+	    if (prrdn) {
+		call fprintf (logfd, " %7g")
+		    call pargr (RA(sh[i]))
+	    }
+	    if (prgain) {
+		call fprintf (logfd, " %6g")
+		    call pargr (DEC(sh[i]))
+	    }
+	    if (prsn) {
+		call fprintf (logfd, " %6g")
+		    call pargr (UT(sh[i]))
+	    }
+	    if (prscale) {
+		call fprintf (logfd, " %6.3f")
+		    call pargr (1./scales[i])
+	    }
+	    if (przero) {
+		call fprintf (logfd, " %7.5g")
+		    call pargr (-zeros[i])
+	    }
+	    if (prwts) {
+		call fprintf (logfd, " %6.3f")
+		    call pargr (wts[i])
+	    }
+	    call fprintf (logfd, "\n")
+	}
+
+	# Log information about the output images.
+	call fprintf (logfd, "\n  Output image = %s, ncombine = %d")
+	    call pargstr (IMNAME(shout))
+	    call pargi (nout)
+	if (expname[1] != EOS) {
+	    call fprintf (logfd, ", %s = %g")
+		call pargstr (expname)
+		call pargr (exposure)
+	}
+	call fprintf (logfd, "\n")
+	call fprintf (logfd,
+	    "  w1 = %g, w2 = %g, dw = %g, nw = %g, dtype = %d\n")
+	    call pargr (W0(shout))
+	    call pargr (W1(shout))
+	    call pargr (WP(shout))
+	    call pargi (SN(shout))
+	    call pargi (DC(shout))
+
+	call flush (logfd)
+	call sfree (sp)
+end
diff --git a/noao/onedspec/scombine/icombine.com b/noao/onedspec/scombine/icombine.com
new file mode 100644
index 00000000..771ada77
--- /dev/null
+++ b/noao/onedspec/scombine/icombine.com
@@ -0,0 +1,36 @@
+# SCOMBINE Common
+
+int	combine			# Combine algorithm
+int	reject			# Rejection algorithm
+real	blank			# Blank value
+pointer	rdnoise			# CCD read noise
+pointer	gain			# CCD gain
+pointer	snoise			# CCD sensitivity noise
+real	lthresh			# Low threshold
+real	hthresh			# High threshold
+int	nkeep			# Minimum to keep
+real	lsigma			# Low sigma cutoff
+real	hsigma			# High sigma cutoff
+real	pclip			# Number or fraction of pixels from median
+real	flow			# Fraction of low pixels to reject
+real	fhigh			# Fraction of high pixels to reject
+int	grow			# Grow radius
+bool	mclip			# Use median in sigma clipping?
+real	sigscale		# Sigma scaling tolerance
+int	logfd			# Log file descriptor
+
+# These flags allow special conditions to be optimized.
+
+int	dflag			# Data flag (D_ALL, D_NONE, D_MIX)
+bool	doscale			# Do the images have to be scaled?
+bool	doscale1		# Do the sigma calculations have to be scaled?
+bool	dothresh		# Check pixels outside specified thresholds?
+bool	dowts			# Does the final average have to be weighted?
+bool	keepids			# Keep track of the image indices?
+bool	docombine		# Call the combine procedure?
+bool	sort			# Sort data?
+
+common	/scbcom/ combine, reject, blank, rdnoise, gain, snoise, lsigma, hsigma,
+		 lthresh, hthresh, nkeep, pclip, flow, fhigh, grow, logfd,
+		 dflag, sigscale, mclip, doscale, doscale1,
+		 dothresh, dowts, keepids, docombine, sort
diff --git a/noao/onedspec/scombine/icombine.h b/noao/onedspec/scombine/icombine.h
new file mode 100644
index 00000000..8a45a673
--- /dev/null
+++ b/noao/onedspec/scombine/icombine.h
@@ -0,0 +1,74 @@
+# SCOMBINE Definitions
+
+# Grouping options
+define	GROUP	"|all|images|apertures|"
+define	GRP_ALL		1
+define	GRP_IMAGES	2
+define	GRP_APERTURES	3
+
+# Sorting options
+define	SORT		"|none|increasing|decreasing|"
+define	SORT_NONE	1
+define	SORT_INC	2
+define	SORT_DEC	3
+
+# Combining modes in interactive mode
+define	CMB_AGAIN	0
+define	CMB_ALL		1
+define	CMB_FIRST	2
+define	CMB_NEXT	3
+define	CMB_SKIP	4
+
+# Rejection options:
+define	REJECT	"|none|ccdclip|crreject|minmax|pclip|sigclip|avsigclip|"
+define	NONE		1	# No rejection algorithm
+define	CCDCLIP		2	# CCD noise function clipping
+define	CRREJECT	3	# CCD noise function clipping
+define	MINMAX		4	# Minmax rejection
+define	PCLIP		5	# Percentile clip
+define	SIGCLIP		6	# Sigma clip
+define	AVSIGCLIP	7	# Sigma clip with average poisson sigma
+
+# Combine options:
+define	COMBINE	"|average|median|sum|"
+define	AVERAGE		1
+define	MEDIAN		2
+define	SUM		3
+
+# Scaling options:
+define	STYPES		"|none|mode|median|mean|exposure|"
+define	ZTYPES		"|none|mode|median|mean|"
+define	WTYPES		"|none|mode|median|mean|exposure|"
+define	S_NONE		1
+define	S_MODE		2
+define	S_MEDIAN	3
+define	S_MEAN		4
+define	S_EXPOSURE	5
+define	S_FILE		6
+define	S_KEYWORD	7
+define	S_SECTION	"|input|output|overlap|"
+define	S_INPUT		1
+define	S_OUTPUT	2
+define	S_OVERLAP	3
+
+# Data flag
+define	D_ALL		0	# All pixels are good
+define	D_NONE		1	# All pixels are bad or rejected
+define	D_MIX		2	# Mixture of good and bad pixels
+
+define	TOL		0.001	# Tolerance for equal residuals
+
+# Spectrum data structure
+define	NS	Memi[$1+$2-1]			# Number of spec of given ap
+define	SH	Memi[Memi[$1+$2-1]+$3-1]	# Spectrum header structure
+
+# Combining options
+#define	COMBINE		"|average|sum|"
+#define	CMB_AVERAGE	1
+#define	CMB_SUM		2
+
+# Weighting options
+#define	WT_TYPE		"|none|expo|user|"
+#define	WT_NONE		1
+#define	WT_EXPO		2
+#define	WT_USER		3
diff --git a/noao/onedspec/scombine/icombine.x b/noao/onedspec/scombine/icombine.x
new file mode 100644
index 00000000..5650d3ab
--- /dev/null
+++ b/noao/onedspec/scombine/icombine.x
@@ -0,0 +1,174 @@
+include	<mach.h>
+include	<smw.h>
+include "icombine.h"
+
+
+# IC_COMBINE -- Combine images.
+
+procedure ic_combiner (sh, shout, d, id, n, m, lflag, scales, zeros, wts,
+	nimages, npts)
+
+pointer	sh[nimages]		# Input spectra
+pointer	shout			# Output spectrum
+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
+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
+pointer	sp, nm
+errchk	ic_scale
+
+include	"icombine.com"
+
+begin
+	call smark (sp)
+
+	# Rebin spectra and set mask arrays
+	call scb_rebin (sh, shout, lflag, nimages, npts)
+
+	# Set scale and weights and log
+	call ic_scale (sh, shout, lflag, scales, zeros, wts, nimages)
+
+	# Set combine parameters
+	switch (combine) {
+	case AVERAGE:
+	    if (dowts)
+		keepids = true
+	    else
+		keepids = false
+	case MEDIAN:
+	    dowts = false
+	    keepids = false
+	case SUM:
+	    keepids = false
+	    reject = NONE
+	    grow = 0
+	}
+	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)] = RA(sh[i])
+	    }
+	    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)] = (Memr[nm+3*(i-1)] / r) ** 2
+		}
+	    } else {
+		do i = 1, nimages {
+		    r = DEC(sh[i])
+		    Memr[nm+3*(i-1)+1] = r
+		    Memr[nm+3*(i-1)] = (Memr[nm+3*(i-1)] / r) ** 2
+		}
+	    }
+	    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 = UT(sh[i])
+		    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)
+	}
+
+	call ic_gdatar (sh, d, id, n, m, lflag, scales, zeros, nimages, npts)
+
+	switch (reject) {
+	case CCDCLIP, CRREJECT:
+	    if (mclip)
+		call ic_mccdclipr (d, id, n, scales, zeros, Memr[nm],
+		    nimages, npts, Memr[SY(shout)])
+	    else
+		call ic_accdclipr (d, id, n, scales, zeros, Memr[nm],
+		    nimages, npts, Memr[SY(shout)])
+	case MINMAX:
+	    call ic_mmr (d, id, n, npts)
+	case PCLIP:
+	    call ic_pclipr (d, id, n, nimages, npts, Memr[SY(shout)])
+	case SIGCLIP:
+	    if (mclip)
+		call ic_msigclipr (d, id, n, scales, zeros, nimages, npts,
+		    Memr[SY(shout)])
+	    else
+		call ic_asigclipr (d, id, n, scales, zeros, nimages, npts,
+		    Memr[SY(shout)])
+	case AVSIGCLIP:
+	    if (mclip)
+		call ic_mavsigclipr (d, id, n, scales, zeros, nimages,
+		    npts, Memr[SY(shout)])
+	    else
+		call ic_aavsigclipr (d, id, n, scales, zeros, nimages,
+		    npts, Memr[SY(shout)])
+	}
+
+	if (grow > 0)
+	    call ic_growr (d, id, n, nimages, npts, Memr[SY(shout)])
+
+	if (docombine) {
+	    switch (combine) {
+	    case AVERAGE:
+		call ic_averager (d, id, n, wts, npts, Memr[SY(shout)])
+	    case MEDIAN:
+		call ic_medianr (d, n, npts, Memr[SY(shout)])
+	    case SUM:
+		call ic_sumr (d, n, npts, Memr[SY(shout)])
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/scombine/icscale.x b/noao/onedspec/scombine/icscale.x
new file mode 100644
index 00000000..009b30c3
--- /dev/null
+++ b/noao/onedspec/scombine/icscale.x
@@ -0,0 +1,463 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<imset.h>
+include	<error.h>
+include	<ctype.h>
+include	<smw.h>
+include	"icombine.h"
+
+# IC_SCALE -- Get the scale factors for the spectra.
+# 1. This procedure does CLIO to determine the type of scaling desired.
+# 2. The output header parameters for exposure time and NCOMBINE are set.
+
+procedure ic_scale (sh, shout, lflags, scales, zeros, wts, nimages)
+
+pointer	sh[nimages]		# Input spectra
+pointer	shout			# Output spectrum
+int	lflags[nimages]		# Data flags
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero or sky levels
+real	wts[nimages]		# Weights
+int	nimages			# Number of images
+
+int	stype, ztype, wtype
+int	i, j, nout
+real	mode, median, mean, exposure, zmean
+pointer	sp, ncombine, exptime, modes, medians, means, expname
+pointer	str, sname, zname, wname, rg
+bool	domode, domedian, domean, dozero
+
+int	ic_gscale()
+real	asumr(), asumi()
+pointer	ic_wranges()
+errchk	ic_gscale, ic_statr
+
+include	"icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (ncombine, nimages, TY_INT)
+	call salloc (exptime, nimages, TY_REAL)
+	call salloc (modes, nimages, TY_REAL)
+	call salloc (medians, nimages, TY_REAL)
+	call salloc (means, nimages, TY_REAL)
+	call salloc (expname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (sname, SZ_FNAME, TY_CHAR)
+	call salloc (zname, SZ_FNAME, TY_CHAR)
+	call salloc (wname, SZ_FNAME, TY_CHAR)
+
+	# Set the defaults.
+	call amovki (1, Memi[ncombine], nimages)
+	call amovkr (0., Memr[exptime], nimages)
+	call amovkr (INDEF, Memr[modes], nimages)
+	call amovkr (INDEF, Memr[medians], nimages)
+	call amovkr (INDEF, Memr[means], nimages)
+	call amovkr (1., scales, nimages)
+	call amovkr (0., zeros, nimages)
+	call amovkr (1., wts, nimages)
+
+	# Set scaling factors.
+	if (combine == SUM) {
+	    stype = S_NONE
+	    ztype = S_NONE
+	    wtype = S_NONE
+	    do i = 1, nimages
+		Memr[exptime+i-1] = IT(sh[i])
+	} else {
+	    stype = ic_gscale ("scale", Memc[sname], STYPES, sh, Memr[exptime],
+		scales, nimages)
+	    ztype = ic_gscale ("zero", Memc[zname], ZTYPES, sh, Memr[exptime],
+		zeros, nimages)
+	    wtype = ic_gscale ("weight", Memc[wname], WTYPES, sh, Memr[exptime],
+		wts, nimages)
+	}
+
+	Memc[expname] = EOS
+	if (combine == SUM || stype == S_EXPOSURE || wtype == S_EXPOSURE) {
+	    call strcpy ("exptime", Memc[expname], SZ_FNAME)
+	    do i = 1, nimages
+		if (IS_INDEFR(Memr[exptime+i-1]))
+		    Memc[expname] = EOS
+	}
+
+	# Get image statistics only if needed.
+	domode = ((stype==S_MODE)||(ztype==S_MODE)||(wtype==S_MODE))
+	domedian = ((stype==S_MEDIAN)||(ztype==S_MEDIAN)||(wtype==S_MEDIAN))
+	domean = ((stype==S_MEAN)||(ztype==S_MEAN)||(wtype==S_MEAN))
+	if (domode || domedian || domean) {
+	    call clgstr ("sample", Memc[str], SZ_LINE)
+	    rg = ic_wranges (Memc[str])
+	    do i = 1, nimages {
+		call ic_statr (sh[i], lflags[i], rg, domode, domedian, domean,
+		    mode, median, mean)
+		if (domode) {
+		    Memr[modes+i-1] = mode
+		    if (stype == S_MODE)
+			scales[i] = mode
+		    if (ztype == S_MODE)
+			zeros[i] = mode
+		    if (wtype == S_MODE)
+			wts[i] = mode
+		}
+		if (domedian) {
+		    Memr[medians+i-1] = median
+		    if (stype == S_MEDIAN)
+			scales[i] = median
+		    if (ztype == S_MEDIAN)
+			zeros[i] = median
+		    if (wtype == S_MEDIAN)
+			wts[i] = median
+		}
+		if (domean) {
+		    Memr[means+i-1] = mean
+		    if (stype == S_MEAN)
+			scales[i] = mean
+		    if (ztype == S_MEAN)
+			zeros[i] = mean
+		    if (wtype == S_MEAN)
+			wts[i] = mean
+		}
+	    }
+	    call mfree (rg, TY_REAL)
+	}
+
+	do i = 1, nimages
+	    if (scales[i] <= 0.) {
+		call eprintf ("WARNING: Negative scale factors")
+		call eprintf (" -- ignoring scaling\n")
+		call amovkr (1., scales, nimages)
+		break
+	    }
+
+	# Convert to relative factors.
+	mean = asumr (scales, nimages) / nimages
+	call adivkr (scales, mean, scales, nimages)
+	call adivr (zeros, scales, zeros, nimages)
+	zmean = asumr (zeros, nimages) / nimages
+
+	if (wtype != S_NONE) {
+	    do i = 1, nimages {
+		if (wts[i] <= 0.) {
+		    call eprintf ("WARNING: Negative weights")
+		    call eprintf (" -- using only NCOMBINE weights\n")
+		    do j = 1, nimages
+			wts[j] = Memi[ncombine+j-1]
+		    break
+		}
+		if (ztype == S_NONE)
+	            wts[i] = Memi[ncombine+i-1] * wts[i]
+		else {
+		    if (zeros[i] <= 0.) {
+			call eprintf ("WARNING: Negative zero offsets")
+			call eprintf (" -- ignoring zero weight adjustments\n")
+			do j = 1, nimages
+			    wts[j] = Memi[ncombine+j-1] * wts[j]
+			break
+		    }
+		    wts[i] = Memi[ncombine+i-1] * wts[i] * zmean / zeros[i]
+		}
+	    }
+	}
+
+	call asubkr (zeros, zmean, zeros, nimages)
+	mean = asumr (wts, nimages)
+	call adivkr (wts, mean, wts, nimages)
+
+	# Because of finite arithmetic it is possible for the zero offsets to
+	# be nonzero even when they are all equal.  Just for the sake of
+	# a nice log set the zero offsets in this case.
+
+	for (i=2; (i<=nimages)&&(zeros[i]==zeros[1]); i=i+1)
+	    ;
+	if (i > nimages)
+	    call aclrr (zeros, nimages)
+
+	# Set flags for scaling, zero offsets, sigma scaling, weights.
+	# Sigma scaling may be suppressed if the scales or zeros are
+	# different by a specified tolerance.
+
+	doscale = false
+	dozero = false
+	doscale1 = false
+	dowts = false
+	do i = 2, nimages {
+	    if (scales[i] != scales[1])
+		doscale = true
+	    if (zeros[i] != zeros[1])
+		dozero = true
+	    if (wts[i] != wts[1])
+		dowts = true
+	}
+	if (doscale && sigscale != 0.) {
+	    do i = 1, nimages {
+		if (abs (scales[i] - 1) > sigscale) {
+		    doscale1 = true
+		    break
+		}
+	    }
+	    if (!doscale1 && zmean > 0.) {
+		do i = 1, nimages {
+		    if (abs (zeros[i] / zmean) > sigscale) {
+			doscale1 = true
+			break
+		    }
+		}
+	    }
+	}
+		    
+	# Set the output header parameters.
+	nout = asumi (Memi[ncombine], nimages)
+	call imaddi (IM(shout), "ncombine", nout)
+	if (Memc[expname] != EOS) {
+	    exposure = 0.
+	    if (combine == SUM) {
+		do i = 1, nimages
+		    exposure = exposure + Memr[exptime+i-1]
+	    } else {
+		do i = 1, nimages
+		    exposure = exposure + wts[i] * Memr[exptime+i-1] / scales[i]
+	    }
+	    call imaddr (IM(shout), Memc[expname], exposure)
+	} else
+	    exposure = INDEF
+
+	# Start the log here since much of the info is only available here.
+	call ic_log (sh, shout, Memi[ncombine], Memr[exptime], Memc[sname],
+	    Memc[zname], Memc[wname], Memr[modes], Memr[medians], Memr[means],
+	    scales, zeros, wts, nimages, dozero, nout, Memc[expname], exposure)
+
+	doscale = (doscale || dozero)
+
+	call sfree (sp)
+end
+
+
+# IC_GSCALE -- Get scale values as directed by CL parameter
+# The values can be one of those in the dictionary, from a file specified
+# with a @ prefix, or from an image header keyword specified by a ! prefix.
+
+int procedure ic_gscale (param, name, dic, sh, exptime, values, nimages)
+
+char	param[ARB]		#I CL parameter name
+char	name[SZ_FNAME]		#O Parameter value
+char	dic[ARB]		#I Dictionary string
+pointer	sh[nimages]		#I SHDR pointers
+real	exptime[nimages]	#I Exposure times
+real	values[nimages]		#O Values
+int	nimages			#I Number of images
+
+int	type			#O Type of value
+
+int	fd, i, nowhite(), open(), fscan(), nscan(), strdic()
+real	rval
+pointer	errstr
+errchk	open
+
+include	"icombine.com"
+
+begin
+	call clgstr (param, name, SZ_FNAME)
+	if (nowhite (name, name, SZ_FNAME) == 0)
+	    type = S_NONE
+	else if (name[1] == '@') {
+	    type = S_FILE
+	    fd = open (name[2], READ_ONLY, TEXT_FILE)
+	    i = 0
+	    while (fscan (fd) != EOF) {
+		call gargr (rval)
+		if (nscan() != 1)
+		    next
+		if (i == nimages) {
+		   call eprintf (
+		       "Warning: Ignoring additional %s values in %s\n")
+		       call pargstr (param)
+		       call pargstr (name[2])
+		   break
+		}
+		i = i + 1
+		values[i] = rval
+	    }
+	    call close (fd)
+	    if (i < nimages) {
+		call salloc (errstr, SZ_LINE, TY_CHAR)
+		call sprintf (Memc[errstr], SZ_FNAME,
+		    "Insufficient %s values in %s")
+		    call pargstr (param)
+		    call pargstr (name[2])
+		call error (1, Memc[errstr])
+	    }
+	} else if (name[1] == '!') {
+	    type = S_KEYWORD
+	    do i = 1, nimages {
+		switch (param[1]) {
+		case 's':
+		    values[i] = ST(sh[i])
+		case 'z':
+		    values[i] = HA(sh[i])
+		case 'w':
+		    values[i] = AM(sh[i])
+		}
+	    }
+	} else {
+	    type = strdic (name, name, SZ_FNAME, dic)
+	    if (type == 0)
+		call error (1, "Unknown scale, zero, or weight type")
+	    if (type==S_EXPOSURE) {
+		do i = 1, nimages {
+		    if (IS_INDEF(IT(sh[i])))
+			call error (1, "Exposure time not defined")
+		    exptime[i] = IT(sh[i])
+		    values[i] = max (0.001, exptime[i])
+		}
+	    }
+	}
+
+	return (type)
+end
+
+
+# IC_WRANGES -- Parse wavelength range string.
+# A wavelength range string  consists of colon delimited ranges with
+# multiple ranges separated by comma and/or whitespace.
+
+pointer procedure ic_wranges (rstr)
+
+char	rstr[ARB]		# Range string
+pointer	rg			# Range pointer
+
+int	i, fd, strlen(), open(), getline()
+pointer	sp, str, ptr
+errchk	open, ic_wadd
+
+begin
+	call smark (sp)
+	call salloc (str, max (strlen (rstr), SZ_LINE), TY_CHAR)
+	call calloc (rg, 1, TY_REAL)
+
+	i = 1
+	while (rstr[i] != EOS) {
+
+	    # Find beginning and end of a range and copy it to the work string
+	    while (IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n')
+	        i = i + 1
+	    if (rstr[i] == EOS)
+		break
+
+	    ptr = str
+	    while (!(IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n' ||
+		rstr[i]==EOS)) {
+		Memc[ptr] = rstr[i]
+	        i = i + 1
+		ptr = ptr + 1
+	    }
+	    Memc[ptr] = EOS
+
+	    # Add range(s)
+	    iferr {
+		if (Memc[str] == '@') {
+		    fd = open (Memc[str+1], READ_ONLY, TEXT_FILE)
+		    while (getline (fd, Memc[str]) != EOF) {
+			iferr (call ic_wadd (rg, Memc[str]))
+			    call erract (EA_WARN)
+		    }
+		    call close (fd)
+		} else
+		    call ic_wadd (rg, Memc[str])
+	    } then
+		call erract (EA_WARN)
+	}
+
+	call sfree (sp)
+
+	# Set final structure
+	i = Memr[rg]
+	if (i == 0)
+	    call mfree (rg, TY_REAL)
+	else
+	    call realloc (rg, 1 + 2 * i, TY_REAL)
+	return (rg)
+end
+
+
+# IC_WADD -- Add a range
+
+procedure ic_wadd (rg, rstr)
+
+pointer	rg			# Range descriptor
+char	rstr[ARB]		# Range string
+
+int	i, j, n, strlen(), ctor()
+real	w1, w2
+pointer	sp, str, ptr
+
+begin
+	call smark (sp)
+	call salloc (str, strlen (rstr), TY_CHAR)
+
+	i = 1
+	n = Memr[rg]
+	while (rstr[i] != EOS) {
+
+	    # Find beginning and end of a range and copy it to the work string
+	    while (IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n')
+	        i = i + 1
+	    if (rstr[i] == EOS)
+		break
+
+	    ptr = str
+	    while (!(IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n' ||
+		rstr[i]==EOS)) {
+		if (rstr[i] == ':')
+		    Memc[ptr] = ' '
+		else
+		    Memc[ptr] = rstr[i]
+	        i = i + 1
+		ptr = ptr + 1
+	    }
+	    Memc[ptr] = EOS
+
+	    # Parse range
+	    if (Memc[str] == '@')
+		call error (1, "Cannot nest @files")
+	    else {
+		# Get range
+		j = 1
+		if (ctor (Memc[str], j, w1) == 0)
+		    call error (1, "Range syntax error")
+		if (ctor (Memc[str], j, w2) == 0)
+		    call error (1, "Range syntax error")
+	    }
+
+	    if (mod (n, 10) == 0)
+		call realloc (rg, 1+2*(n+10), TY_REAL)
+	    n = n + 1
+	    Memr[rg+2*n-1] = min (w1, w2)
+	    Memr[rg+2*n] = max (w1, w2)
+	}
+	Memr[rg] = n
+
+	call sfree (sp)
+end
+
+
+# IC_WISINRANGE -- Is wavelength in range?
+
+bool procedure ic_wisinrange (rg, w)
+
+pointer	rg		# Wavelength range array
+real	w		# Wavelength
+
+int	i, n
+
+begin
+	if (rg == NULL)
+	    return (true)
+
+	n = nint (Memr[rg])
+	do i = 1, 2*n, 2
+	    if (w >= Memr[rg+i] && w <= Memr[rg+i+1])
+		return (true)
+	return (false)
+end
diff --git a/noao/onedspec/scombine/icstat.x b/noao/onedspec/scombine/icstat.x
new file mode 100644
index 00000000..3fce4165
--- /dev/null
+++ b/noao/onedspec/scombine/icstat.x
@@ -0,0 +1,160 @@
+include	<smw.h>
+include	"icombine.h"
+
+
+# IC_STATR -- Compute image statistics within spectrum.
+
+procedure ic_statr (sh, lflag, rg, domode, domedian, domean, mode, median, mean)
+
+pointer	sh			 # Spectrum structure
+int	lflag			 # Data flag
+pointer	rg			 # Wavelength ranges
+bool	domode, domedian, domean # Statistics to compute
+real	mode, median, mean	 # Statistics
+
+int	i, n, npts
+real	a, w
+pointer	sp, data, dp, lp, mp
+real	ic_moder(), asumr()
+bool	ic_wisinrange()
+double	shdr_lw()
+
+include	"icombine.com"
+
+begin
+	mp = SX(sh)
+	lp = SY(sh)
+	npts = SN(sh)
+
+	call smark (sp)
+	call salloc (data, npts, TY_REAL)
+
+	dp = data
+	if (lflag == D_ALL && rg == NULL) {
+	    if (dothresh) {
+		do i = 1, npts {
+		    a = Memr[lp]
+		    if (a >= lthresh && a <= hthresh) {
+			Memr[dp] = a
+			dp = dp + 1
+		    }
+		    lp = lp + 1
+		}
+	    } else {
+		do i = 1, npts {
+		    Memr[dp] = Memr[lp]
+		    dp = dp + 1
+		    lp = lp + 1
+		}
+	    }
+	} else if (lflag == D_MIX || rg != NULL) {
+	    if (dothresh) {
+		do i = 1, npts {
+		    if (Memr[mp] == 0) {
+			a = Memr[lp]
+			if (a >= lthresh && a <= hthresh) {
+			    w = shdr_lw (sh, double (i))
+			    if (ic_wisinrange (rg, w)) {
+				Memr[dp] = a
+				dp = dp + 1
+			    }
+			}
+		    }
+		    mp = mp + 1
+		    lp = lp + 1
+		}
+	    } else {
+		do i = 1, npts {
+		    if (Memr[mp] == 0) {
+			w = shdr_lw (sh, double (i))
+			if (ic_wisinrange (rg, w)) {
+			    Memr[dp] = Memr[lp]
+			    dp = dp + 1
+			}
+		    }
+		    mp = mp + 1
+		    lp = lp + 1
+		}
+	    }
+	}
+
+	n = dp - data
+	if (n > 0) {
+	    # Compute only statistics needed.
+	    if (domode || domedian) {
+		call asrtr (Memr[data], Memr[data], n)
+		mode = ic_moder (Memr[data], n)
+		median = Memr[data+n/2-1]
+	    }
+	    if (domean)
+		mean = asumr (Memr[data], n) / n
+	} else {
+	    mode = INDEF
+	    median = INDEF
+	    mean = INDEF
+	}
+
+	call sfree (sp)
+end
+
+
+define	NMIN	10	# Minimum number of pixels for mode calculation
+define	ZRANGE	0.8	# Fraction of pixels about median to use
+define	ZSTEP	0.01	# Step size for search for mode
+define	ZBIN	0.1	# Bin size for mode.
+
+# IC_MODE -- Compute mode of an array.  The mode is found by binning
+# with a bin size based on the data range over a fraction of the
+# pixels about the median and a bin step which may be smaller than the
+# bin size.  If there are too few points the median is returned.
+# The input array must be sorted.
+
+real procedure ic_moder (a, n)
+
+real	a[n]			# Data array
+int	n			# Number of points
+
+int	i, j, k, nmax
+real	z1, z2, zstep, zbin
+real	mode
+bool	fp_equalr()
+
+begin
+	if (n < NMIN)
+	    return (a[n/2])
+
+	# Compute the mode.  The array must be sorted.  Consider a
+	# range of values about the median point.  Use a bin size which
+	# is ZBIN of the range.  Step the bin limits in ZSTEP fraction of
+	# the bin size.
+
+	i = 1 + n * (1. - ZRANGE) / 2.
+	j = 1 + n * (1. + ZRANGE) / 2.
+	z1 = a[i]
+	z2 = a[j]
+	if (fp_equalr (z1, z2)) {
+	    mode = z1
+	    return (mode)
+	}
+
+	zstep = ZSTEP * (z2 - z1)
+	zbin = ZBIN * (z2 - z1)
+
+	z1 = z1 - zstep
+	k = i
+	nmax = 0
+	repeat {
+	    z1 = z1 + zstep
+	    z2 = z1 + zbin
+	    for (; i < j && a[i] < z1; i=i+1)
+		;
+	    for (; k < j && a[k] < z2; k=k+1)
+		;
+	    if (k - i > nmax) {
+	        nmax = k - i
+	        mode = a[(i+k)/2]
+	    }
+	} until (k >= j)
+
+	return (mode)
+end
diff --git a/noao/onedspec/scombine/icsum.x b/noao/onedspec/scombine/icsum.x
new file mode 100644
index 00000000..f038b37b
--- /dev/null
+++ b/noao/onedspec/scombine/icsum.x
@@ -0,0 +1,48 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"icombine.h"
+
+
+# IC_SUM -- Compute the summed image line.
+
+procedure ic_sumr (d, n, npts, sum)
+
+pointer	d[ARB]			# Data pointers
+int	n[npts]			# Number of points
+int	npts			# Number of output points per line
+real	sum[npts]		# Average (returned)
+
+int	i, j, k
+real	s
+
+include	"icombine.com"
+
+begin
+	# If no data has been excluded do the sum without checking the
+	# number of points.  If all the data has been excluded set the
+	# sum to the blank value.
+
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		k = i - 1
+		s = Memr[d[1]+k]
+		do j = 2, n[i]
+		    s = s + Memr[d[j]+k]
+		sum[i] = s
+	    }
+	} else if (dflag == D_NONE) {
+	    do i = 1, npts
+		sum[i] = blank
+	} else {
+	    do i = 1, npts {
+		if (n[i] > 0) {
+		    k = i - 1
+		    s = Memr[d[1]+k]
+		    do j = 2, n[i]
+			s = s + Memr[d[j]+k]
+		    sum[i] = s
+		} else
+		    sum[i] = blank
+	    }
+	}
+end
diff --git a/noao/onedspec/scombine/iscombine.key b/noao/onedspec/scombine/iscombine.key
new file mode 100644
index 00000000..91d7876b
--- /dev/null
+++ b/noao/onedspec/scombine/iscombine.key
@@ -0,0 +1,23 @@
+		SCOMBINE CURSOR KEYS
+
+a - Mark scaling ranges in overlap region for 'v' key
+b - Cancel scaling ranges
+c - Print cursor position
+d - Replace range of pixels by straight cursor line
+e - Replace range of pixels by linear interpolation from the endpoint pixels
+f - Start over from the first spectrum
+j - Fudge a point to vertical cursor value
+n - Go on to next spectrum
+o - Reset data for current spectrum to initial values
+p - Don't include current spectrum in combined image and go on to next spectrum
+q - Quit and combine remaining spectra noninteractively (no|yes|YES)
+s - Mark accumulation ranges in current spectrum
+t - Cancel accumulation ranges
+v - Shift overlap average of spectrum vertically to accumulated spectrum
+w - Window the graph using gtools commands
+x - Shift spectrum horizontally to cursor position
+y - Shift spectrum vertically to cursor position
+z - Shift spectrum vertically to accumulated spectrum
++ - Set additive scaling for 'v' key
+* - Set multiplicative scaling 
+? - This help page
diff --git a/noao/onedspec/scombine/iscombine.par b/noao/onedspec/scombine/iscombine.par
new file mode 100644
index 00000000..a9d8846e
--- /dev/null
+++ b/noao/onedspec/scombine/iscombine.par
@@ -0,0 +1,18 @@
+input,s,a,"",,,List of input spectra
+output,s,a,"",,,List of output spectra
+woutput,s,h,"",,,List of output weight spectra
+apertures,s,h,"",,,Apertures to combine
+group,s,h,"apertures","all|images|apertures",,Grouping option
+combine,s,h,"average","average|sum",,Combining option
+scale,s,h,"",,,Header keyword for scaling
+weight,s,h,"","",,"Header keyword for weighting
+"
+w1,r,h,INDEF,,,Starting wavelength of output spectra
+w2,r,h,INDEF,,,Ending wavelength of output spectra
+dw,r,h,INDEF,,,Wavelength increment of output spectra
+nw,i,h,INDEF,,,Length of output spectra
+log,b,h,no,,,"Logarithmic increments?
+"
+interactive,b,h,no,,,Adjust spectra interactively?
+sort,s,h,"none","none|increasing|decreasing",,Interactive combining order
+cursor,*gcur,h,"",,,Graphics cursor input
diff --git a/noao/onedspec/scombine/mkpkg b/noao/onedspec/scombine/mkpkg
new file mode 100644
index 00000000..ab60e45b
--- /dev/null
+++ b/noao/onedspec/scombine/mkpkg
@@ -0,0 +1,35 @@
+# SCOMBINE
+
+$call	relink
+$exit
+
+update:
+	$call	relink
+	$call	install
+	;
+
+relink:
+	$update	libpkg.a
+	$omake	x_scombine.x
+	$link	x_scombine.o libpkg.a -lsmw -lxtools -liminterp \
+	    -o xx_scombine.e
+	;
+
+install:
+	$move	xx_scombine.e noaobin$x_scombine.e
+	;
+
+
+libpkg.a:
+	@generic
+
+	icgdata.x	<smw.h> icombine.com icombine.h
+	iclog.x		<smw.h> icombine.com icombine.h <mach.h>
+	icombine.x	<smw.h> icombine.com <mach.h> icombine.h
+	icscale.x	<smw.h> icombine.com icombine.h <ctype.h> <error.h>\
+			<imhdr.h> <imset.h>
+	icstat.x	<smw.h> icombine.com icombine.h
+	icsum.x		icombine.com icombine.h
+	t_scombine.x	<smw.h> icombine.h icombine.com <error.h> <imhdr.h>\
+			<mach.h>
+	;
diff --git a/noao/onedspec/scombine/scombine.par b/noao/onedspec/scombine/scombine.par
new file mode 100644
index 00000000..932e6e31
--- /dev/null
+++ b/noao/onedspec/scombine/scombine.par
@@ -0,0 +1,37 @@
+input,s,a,"",,,List of input spectra
+output,s,a,"",,,List of output spectra
+noutput,s,h,"",,,List of output number combined spectra
+logfile,s,h,"STDOUT",,,"Log file
+"
+apertures,s,h,"",,,Apertures to combine
+group,s,h,"apertures","all|images|apertures",,"Grouping option"
+combine,s,h,"average","average|median|sum",,Type of combine operation
+reject,s,h,"none","none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip",,"Type of rejection
+"
+first,b,h,no,,,Use first spectrum for dispersion?
+w1,r,h,INDEF,,,Starting wavelength of output spectra
+w2,r,h,INDEF,,,Ending wavelength of output spectra
+dw,r,h,INDEF,,,Wavelength increment of output spectra
+nw,i,h,INDEF,,,Length of output spectra
+log,b,h,no,,,"Logarithmic increments?
+"
+scale,s,h,"none",,,Image scaling
+zero,s,h,"none",,,Image zero point offset
+weight,s,h,"none",,,Image weights
+sample,s,h,"",,,"Wavelength sample regions for statistics
+"
+lthreshold,r,h,INDEF,,,Lower threshold
+hthreshold,r,h,INDEF,,,Upper threshold
+nlow,i,h,1,0,,minmax: Number of low pixels to reject
+nhigh,i,h,1,0,,minmax: Number of high pixels to reject
+nkeep,i,h,1,,,Minimum to keep (pos) or maximum to reject (neg)
+mclip,b,h,yes,,,Use median in sigma clipping algorithms?
+lsigma,r,h,3.,0.,,Lower sigma clipping factor
+hsigma,r,h,3.,0.,,Upper sigma clipping factor
+rdnoise,s,h,"0.",,,ccdclip: CCD readout noise (electrons)
+gain,s,h,"1.",,,ccdclip: CCD gain (electrons/DN)
+snoise,s,h,"0.",,,ccdclip: Sensitivity noise (fraction)
+sigscale,r,h,0.1,0.,,Tolerance for sigma clipping scaling corrections
+pclip,r,h,-0.5,,,pclip: Percentile clipping parameter
+grow,i,h,0,,,Radius (pixels) for 1D neighbor rejection
+blank,r,h,0.,,,Value if there are no pixels
diff --git a/noao/onedspec/scombine/t_scombine.x b/noao/onedspec/scombine/t_scombine.x
new file mode 100644
index 00000000..774a5f87
--- /dev/null
+++ b/noao/onedspec/scombine/t_scombine.x
@@ -0,0 +1,630 @@
+include	<imhdr.h>
+include	<error.h>
+include	<mach.h>
+include	<smw.h>
+include "icombine.h"
+
+
+# T_SCOMBINE - Combine spectra
+# The input spectra are combined by medianing, averaging or summing
+# with optional rejection, scaling and weighting.  The input may be
+# grouped by aperture or by image.  The combining algorithms are
+# similar to those in IMCOMBINE.
+
+procedure t_scombine()
+
+int	ilist			# list of input images
+int	olist			# list of output images
+pointer	nlist			# image name for number combined
+pointer	aps			# aperture ranges
+int	group			# grouping option
+
+int	reject1
+real	flow1, fhigh1, pclip1, nkeep1
+
+real	rval
+bool	grdn, ggain, gsn
+int	i, j, k, l, n, naps, npts
+pointer	im, mw, nout, refim, shin, shout
+pointer	sp, input, output, noutput, scale, zero, weight, str, logfile, sh, ns
+pointer	sp1, d, id, nc, m, lflag, scales, zeros, wts
+
+real	clgetr(), imgetr()
+bool	clgetb(), rng_elementi()
+int	clgeti(), clgwrd(), ctor()
+int	imtopenp(), imtgetim(), open(), nowhite()
+pointer	rng_open(), immap(), smw_openim(), impl2i(), impl2r()
+errchk	open, immap, smw_openim, shdr_open, imgetr
+errchk	scb_output, scb_combine, ic_combiner
+
+include	"icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (noutput, SZ_FNAME, TY_CHAR)
+	call salloc (scale, SZ_FNAME, TY_CHAR)
+	call salloc (zero, SZ_FNAME, TY_CHAR)
+	call salloc (weight, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (gain, SZ_FNAME, TY_CHAR)
+	call salloc (snoise, SZ_FNAME, TY_CHAR)
+	call salloc (rdnoise, SZ_FNAME, TY_CHAR)
+	call salloc (logfile, SZ_FNAME, TY_CHAR)
+
+	# Get parameters
+	ilist = imtopenp ("input")
+	olist = imtopenp ("output")
+	nlist = imtopenp ("noutput")
+	call clgstr ("apertures", Memc[str], SZ_LINE)
+	group = clgwrd ("group", Memc[input], SZ_FNAME, GROUP)
+
+	# IMCOMBINE parameters
+	call clgstr ("logfile", Memc[logfile], SZ_FNAME)
+	combine = clgwrd ("combine", Memc[input], SZ_FNAME, COMBINE)
+	reject1 = clgwrd ("reject", Memc[input], SZ_FNAME, REJECT)
+	blank = clgetr ("blank")
+	call clgstr ("scale", Memc[scale], SZ_FNAME)
+	call clgstr ("zero", Memc[zero], SZ_FNAME)
+	call clgstr ("weight", Memc[weight], SZ_FNAME)
+	call clgstr ("gain", Memc[gain], SZ_FNAME)
+	call clgstr ("rdnoise", Memc[rdnoise], SZ_FNAME)
+	call clgstr ("snoise", Memc[snoise], SZ_FNAME)
+	lthresh = clgetr ("lthreshold")
+	hthresh = clgetr ("hthreshold")
+	lsigma = clgetr ("lsigma")
+	hsigma = clgetr ("hsigma")
+	pclip1 = clgetr ("pclip")
+	flow1 = clgetr ("nlow")
+	fhigh1 = clgetr ("nhigh")
+	nkeep1 = clgeti ("nkeep")
+	grow = clgeti ("grow")
+	mclip = clgetb ("mclip")
+	sigscale = clgetr ("sigscale")
+
+	i = nowhite (Memc[scale], Memc[scale], SZ_FNAME)
+	i = nowhite (Memc[zero], Memc[zero], SZ_FNAME)
+	i = nowhite (Memc[weight], Memc[weight], SZ_FNAME)
+
+	# Check parameters, map INDEFs, and set threshold flag
+	if (combine == SUM)
+	    reject1 = NONE
+	if (pclip1 == 0. && reject1 == PCLIP)
+	    call error (1, "Pclip parameter may not be zero")
+	if (IS_INDEFR (blank))
+	    blank = 0.
+	if (IS_INDEFR (lsigma))
+	    lsigma = MAX_REAL
+	if (IS_INDEFR (hsigma))
+	    hsigma = MAX_REAL
+	if (IS_INDEFR (pclip1))
+	    pclip1 = -0.5
+	if (IS_INDEFI (nkeep1))
+	    nkeep1 = 0
+	if (IS_INDEFR (flow1))
+	    flow1 = 0
+	if (IS_INDEFR (fhigh))
+	    fhigh = 0
+	if (IS_INDEFI (grow))
+	    grow = 0
+	if (IS_INDEF (sigscale))
+	    sigscale = 0.
+
+	if (IS_INDEF(lthresh) && IS_INDEF(hthresh))
+	    dothresh = false
+	else {
+	    dothresh = true
+	    if (IS_INDEF(lthresh))
+		lthresh = -MAX_REAL
+	    if (IS_INDEF(hthresh))
+		hthresh = MAX_REAL
+	}
+
+	# Get read noise and gain?
+	grdn = false
+	ggain = false
+	gsn = false
+	if (reject1 == CCDCLIP || reject1 == CRREJECT) {
+	    i = 1
+	    if (ctor (Memc[rdnoise], i, rval) == 0)
+		grdn = true
+	    i = 1
+	    if (ctor (Memc[gain], i, rval) == 0)
+		ggain = true
+	    i = 1
+	    if (ctor (Memc[snoise], i, rval) == 0)
+		gsn = true
+	}
+
+	# Open the log file.
+	logfd = NULL
+	if (Memc[logfile] != EOS) {
+	    iferr (logfd = open (Memc[logfile], APPEND, TEXT_FILE)) {
+	        logfd = NULL
+	        call erract (EA_WARN)
+	    }
+	}
+
+	iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF))
+	    call error (1, "Error in aperture list")
+
+	# Loop through input images.
+	while (imtgetim (ilist, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (olist, Memc[output], SZ_FNAME) == EOF) {
+		call eprintf ("No output image\n")
+		break
+	    }
+	    if (imtgetim (nlist, Memc[noutput], SZ_FNAME) == EOF)
+		Memc[noutput] = EOS
+
+	    # Get spectra to combine.
+	    # Because the input images are unmapped we must get all the
+	    # data we need for combining into the spectrum data structures.
+	    # In particular any header keyword parameters that will be
+	    # used.  We save the header values in unused elements of
+	    # the spectrum data structure.
+
+	    naps = 0
+	    repeat {
+		iferr (im = immap (Memc[input], READ_ONLY, 0)) {
+		    if (group == GRP_IMAGES) {
+			call erract (EA_WARN)
+			next
+		    } else {
+			call erract (EA_ERROR)
+		    }
+		}
+		mw = smw_openim (im)
+		shin = NULL
+
+		do i = 1, SMW_NSPEC(mw) {
+		    call shdr_open (im, mw, i, 1, INDEFI, SHDATA, shin)
+		    if (Memc[scale] == '!')
+			ST(shin) = imgetr (im, Memc[scale+1])
+		    if (Memc[zero] == '!')
+			HA(shin) = imgetr (im, Memc[zero+1])
+		    if (Memc[weight] == '!')
+			AM(shin) = imgetr (im, Memc[weight+1])
+		    if (grdn)
+			RA(shin) = imgetr (im, Memc[rdnoise])
+		    if (ggain)
+			DEC(shin) = imgetr (im, Memc[gain])
+		    if (gsn)
+			UT(shin) = imgetr (im, Memc[snoise])
+		    if (!rng_elementi (aps, AP(shin)))
+			next
+		    if (group == GRP_APERTURES) {
+			for (j=1; j<=naps; j=j+1)
+			    if (AP(shin) == AP(SH(sh,j,1)))
+				break
+			n = 10
+		    } else {
+			j = 1
+			n = 1
+		    }
+
+		    if (naps == 0) {
+			call calloc (sh, n, TY_POINTER)
+			call calloc (ns, n, TY_INT)
+		    } else if (j > naps && mod (naps, n) == 0) {
+			call realloc (sh, naps+n, TY_POINTER)
+			call realloc (ns, naps+n, TY_INT)
+			call aclri (Memi[sh+naps], n)
+			call aclri (Memi[ns+naps], n)
+		    }
+		    if (j > naps)
+			naps = naps + 1
+		    n = NS(ns,j)
+		    if (n == 0)
+			call malloc (Memi[sh+j-1], 10, TY_POINTER)
+		    else if (mod (n, 10) == 0)
+			call realloc (Memi[sh+j-1], n+10, TY_POINTER)
+
+		    n = n + 1
+		    SH(sh,j,n) = NULL
+		    NS(ns,j) = n
+		    call shdr_copy (shin, SH(sh,j,n), NO)
+		}
+
+		call imunmap (IM(shin))
+		MW(shin) = NULL
+		call shdr_close (shin)
+
+		if (group == GRP_IMAGES)
+		    break
+	    } until (imtgetim (ilist, Memc[input], SZ_FNAME) == EOF)
+
+	    if (naps < 1) {
+		call eprintf ("No input spectra to combine\n")
+		next
+	    }
+
+	    # Set the output and combine the spectra.
+	    call scb_output (sh, ns, naps, Memc[output], Memc[noutput],
+		im, mw, nout, refim)
+
+	    do j = 1, naps {
+		call shdr_open (im, mw, j, 1, INDEFI, SHHDR, shout)
+		npts = SN(shout)
+		n = NS(ns,j)
+
+		# Allocate additional memory
+		call smark (sp1)
+		call salloc (d, n, TY_POINTER)
+		call salloc (id, n, TY_POINTER)
+		call salloc (nc, npts, TY_INT)
+		call salloc (m, n, TY_POINTER)
+		call salloc (lflag, n, TY_INT)
+		call salloc (scales, n, TY_REAL)
+		call salloc (zeros, n, TY_REAL)
+		call salloc (wts, n, TY_REAL)
+		call calloc (SX(shout), npts, TY_REAL)
+		call calloc (SY(shout), npts, TY_REAL)
+		call amovki (D_ALL, Memi[lflag], n)
+
+		# Convert the pclip parameter to a number of pixels rather than
+		# a fraction.  This number stays constant even if pixels are
+		# rejected.  The number of low and high pixel rejected, however,
+		# are converted to a fraction of the valid pixels.
+
+		reject = reject1
+		nkeep = nkeep1
+		if (nkeep < 0)
+		    nkeep = n + nkeep
+		if (reject == PCLIP) {
+		    pclip = pclip1
+		    i = (n - 1) / 2.
+		    if (abs (pclip) < 1.)
+			pclip = pclip * i
+		    if (pclip < 0.)
+			pclip = min (-1, max (-i, int (pclip)))
+		    else
+			pclip = max (1, min (i, int (pclip)))
+		}
+		if (reject == MINMAX) {
+		    flow = flow1
+		    fhigh = fhigh1
+		    if (flow >= 1)
+			flow = flow / n
+		    if (fhigh >= 1)
+			fhigh = fhigh / n
+		    i = flow * n + fhigh * n
+		    if (i == 0)
+			reject = NONE
+		    else if (i >= n) {
+			call eprintf ("Bad minmax rejection parameters\n")
+			call eprintf ("Using no rejection\n")
+			reject = NONE
+		    }
+		}
+
+		# Combine spectra
+		call ic_combiner (SH(sh,j,1), shout, Memi[d], Memi[id],
+		    Memi[nc], Memi[m], Memi[lflag], Memr[scales], Memr[zeros],
+		    Memr[wts], n, npts)
+
+		# Write the results
+		call amovr (Memr[SY(shout)], Memr[impl2r(im,j)], npts)
+		if (nout != NULL)
+		    call amovi (Memi[nc], Memi[impl2i(nout,j)], npts)
+		call sfree (sp1)
+	    }
+
+	    # Finish up
+	    call shdr_close (shout)
+	    call smw_close (mw)
+	    call imunmap (im)
+	    call imunmap (refim)
+	    if (nout != NULL)
+		call imunmap (nout)
+
+	    # Find all the distinct SMW pointers and free them.
+	    do j = 1, naps {
+		do i = 1, NS(ns,j) {
+		    mw = MW(SH(sh,j,i))
+		    if (mw != NULL) {
+			do k = 1, naps  {
+			    do l = 1, NS(ns,k) {
+				shin = SH(sh,k,l)
+				if (MW(shin) == mw)
+				    MW(shin) = NULL
+			    }
+			}
+			call smw_close (mw)
+		    }
+		}
+	    }
+	    do j = 1, naps {
+		do i = 1, NS(ns,j)
+		    call shdr_close (SH(sh,j,i))
+		call mfree (Memi[sh+j-1], TY_POINTER)
+	    }
+	    call mfree (sh, TY_POINTER)
+	    call mfree (ns, TY_INT)
+	}
+
+	call rng_close (aps)
+	call imtclose (ilist)
+	call imtclose (olist)
+	call imtclose (nlist)
+
+	call sfree (sp)
+end
+
+
+# SCB_REBIN - Rebin input spectra to output dispersion
+# Use the SX array as mask.  If less than 1% of an input
+# pixel contributes to an output pixel then flag it as missing data.
+
+procedure scb_rebin (sh, shout, lflag, ns, npts)
+
+pointer	sh[ns]			# Input spectra structures
+pointer	shout			# Output spectrum structure
+int	lflag[ns]		# Empty mask flags
+int	ns			# Number of spectra
+int	npts			# NUmber of output points
+
+int	i, j
+real	a, b, c
+pointer	shin
+double	shdr_wl(), shdr_lw()
+
+include	"icombine.com"
+	
+begin
+	# Rebin to common dispersion
+	# Determine overlap with output and set mask arrays
+
+	do i = 1, ns {
+	    shin = sh[i]
+	    c = shdr_wl (shout, shdr_lw (shin, double(0.5)))
+	    b = shdr_wl (shout, shdr_lw (shin, double(SN(shin)+0.5)))
+	    a = max (1, nint (min (b, c) + 0.01))
+	    b = min (npts, nint (max (b, c) - 0.01))
+	    j = b - a + 1
+	    if (j < 1) {
+		lflag[i] = D_NONE
+		next
+	    }
+	    else if (j < npts)
+		lflag[i] = D_MIX
+	    else
+		lflag[i] = D_ALL
+
+	    call shdr_rebin (shin, shout)
+	    call aclrr (Memr[SX(shin)], SN(shin))
+	    j = a - 1
+	    if (j > 0)
+		call amovkr (1.0, Memr[SX(shin)], j)
+	    j = SN(shin) - b
+	    if (j > 0)
+		call amovkr (1.0, Memr[SX(shin)+SN(shin)-j], j)
+	}
+
+	dflag = lflag[1]
+	do i = 2, ns {
+	    if (dflag != lflag[i]) {
+		dflag = D_MIX
+		break
+	    }
+	}
+end
+
+
+# SCB_OUTPUT - Set the output spectrum
+
+procedure scb_output (sh, ns, naps, output, noutput, im, mw, nout, refim)
+
+pointer	sh			# spectra structures
+int	ns			# number of spectra
+int	naps			# number of apertures
+char	output[SZ_FNAME]	# output spectrum name
+char	noutput[SZ_FNAME]	# output number combined image name
+pointer	im			# output IMIO pointer
+pointer	mw			# output MWCS pointer
+pointer	nout			# output number combined IMIO pointer
+pointer	refim			# reference image for output image
+
+int	i, ap, beam, dtype, nw, nmax, axis[2]
+double	w1, dw, z
+real	aplow[2], aphigh[2]
+pointer	sp, key, coeff, sh1
+pointer	immap(), mw_open(), smw_openim()
+errchk	immap, smw_openim
+data	axis/1,2/
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	coeff = NULL
+
+	# Create output image using the first input image as a reference
+	refim = immap (IMNAME(SH(sh,1,1)), READ_ONLY, 0)
+	im = immap (output, NEW_COPY, refim)
+
+	# Use smw_openim to clean up old keywords(?).
+	mw = smw_openim (im)
+	call smw_close (mw)
+
+	if (naps == 1)
+	    IM_NDIM(im) = 1
+	else
+	    IM_NDIM(im) = 2
+	call imaddi (im, "SMW_NDIM", IM_NDIM(im))
+	IM_LEN(im,2) = naps
+	if (IM_PIXTYPE(im) != TY_DOUBLE)
+	    IM_PIXTYPE(im) = TY_REAL
+
+	# Set new header.
+	mw = mw_open (NULL, 2)
+	call mw_newsystem (mw, "multispec", 2)
+	call mw_swtype (mw, axis, 2, "multispec",
+	    "label=Wavelength units=Angstroms")
+	call smw_open (mw, NULL, im)
+
+	nmax = 0
+	do i = 1, naps {
+	    sh1 = SH(sh,i,1)
+	    call smw_gwattrs (MW(sh1), APINDEX(sh1), 1, ap, beam, dtype,
+		w1, dw, nw, z, aplow, aphigh, coeff)
+	    call scb_default (SH(sh,i,1), NS(ns,i),
+		dtype, w1, dw, nw, z, Memc[coeff])
+	    call smw_swattrs (mw, i, 1, ap, beam, dtype,
+		w1, dw, nw, z, aplow, aphigh, Memc[coeff])
+	    call smw_sapid (mw, i, 1, TITLE(sh1))
+	    nmax = max (nmax, nw)
+	}
+
+	IM_LEN(im,1) = nmax
+
+	# Set MWCS header.
+	call smw_saveim (mw, im)
+	call smw_close (mw)
+	mw = smw_openim (im)
+
+	# Create number combined image
+	if (noutput[1] != EOS) {
+	    nout = immap (noutput, NEW_COPY, im)
+	    IM_PIXTYPE(nout) = TY_INT
+	    call sprintf (IM_TITLE(nout), SZ_LINE, "Number combined for %s")
+		call pargstr (output)
+	}
+
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end
+
+
+# SCB_DEFAULT - Set default values for the starting wavelength, ending
+# wavelength, wavelength increment and spectrum length for the output
+# spectrum.
+
+procedure scb_default (shdr, ns, dtype, w1, dw, nw, z, coeff)
+
+pointer	shdr[ARB]		# spectra structures
+int	ns			# number of spectra
+int	dtype			# dispersion type
+double	w1			# starting wavelength
+double	dw			# wavelength increment
+int	nw			# spectrum length
+double	z			# redshift
+char	coeff[ARB]		# nonlinear coefficient array
+
+bool	clgetb()
+int	i, nwa, clgeti()
+double	w2, aux, w1a, w2a, dwa, clgetd()
+pointer	sh
+
+begin
+	if (clgetb ("first"))
+	    return
+
+	w1a = clgetd ("w1")
+	w2a = clgetd ("w2")
+	dwa = clgetd ("dw")
+	nwa = clgeti ("nw")
+	if (clgetb ("log"))
+	    dtype = DCLOG
+	else
+	    dtype = DCLINEAR
+	z = 0.
+	coeff[1] = EOS
+
+	# Dispersion type
+	if (dtype == DCLINEAR) {
+	    do i = 1, ns {
+		if (DC(shdr[i]) == DCNO) {
+		    dtype = DCNO
+		    break
+		}
+	    }
+	}
+
+	w1 = w1a
+	w2 = w2a
+	dw = dwa
+	nw = nwa
+
+	# Starting wavelength
+	if (IS_INDEFD (w1)) {
+	    if (IS_INDEFD (dw) || dw > 0.) {
+		w1 = MAX_REAL
+		do i = 1, ns {
+		    sh = shdr[i]
+		    if (WP(sh) > 0.)
+			aux = W0(sh)
+		    else
+			aux = W1(sh)
+		    if (aux < w1)
+			w1 = aux
+		}
+	    } else {
+		w1 = -MAX_REAL
+		do i = 1, ns {
+		    sh = shdr[i]
+		    if (WP(sh) > 0.)
+			aux = W1(sh)
+		    else
+			aux = W0(sh)
+		    if (aux > w1)
+			w1 = aux
+		}
+	    }
+	}
+
+	# Ending wavelength
+	if (IS_INDEFD (w2)) {
+	    if (IS_INDEFD (dw) || dw > 0.) {
+		w2 = -MAX_REAL
+		do i = 1, ns {
+		    sh = shdr[i]
+		    if (WP(sh) > 0.)
+			aux = W1(sh)
+		    else
+			aux = W0(sh)
+		    if (aux > w2)
+			w2 = aux
+		}
+	    } else {
+		w2 = MAX_REAL
+		do i = 1, ns {
+		    sh = shdr[i]
+		    if (WP(sh) > 0.)
+			aux = W0(sh)
+		    else
+			aux = W1(sh)
+		    if (aux < w2)
+			w2 = aux
+		}
+	    }
+	}
+
+	# Wavelength increment
+	if (IS_INDEFD (dw)) {
+	    dw = MAX_REAL
+	    do i = 1, ns {
+		aux = abs (WP(shdr[i]))
+		if (aux < dw)
+		    dw = aux
+	    }
+	}
+	if ((w2 - w1) / dw < 0.)
+	    dw = -dw
+
+	# Spectrum length
+	if (IS_INDEFI (nw))
+	    nw = int ((w2 - w1) / dw + 0.5) + 1
+
+	# Adjust the values.
+	if (IS_INDEFD (dwa))
+	    dw = (w2 - w1) / (nw - 1)
+	else if (IS_INDEFD (w2a))
+	    w2 = w1 + (nw - 1) * dw
+	else if (IS_INDEFD (w1a))
+	    w1 = w2 - (nw - 1) * dw
+	else {
+	    nw = int ((w2 - w1) / dw + 0.5) + 1
+	    w2 = w1 + (nw - 1) * dw
+	}
+end
diff --git a/noao/onedspec/scombine/x_scombine.x b/noao/onedspec/scombine/x_scombine.x
new file mode 100644
index 00000000..33943271
--- /dev/null
+++ b/noao/onedspec/scombine/x_scombine.x
@@ -0,0 +1 @@
+task	scombine = t_scombine