Repatch (from linux) of OSX IRAF

1 files changed, 1071 insertions, 0 deletions

diff --git a/noao/onedspec/odcombine/t_odcombine.x b/noao/onedspec/odcombine/t_odcombine.x
new file mode 100644
index 00000000..9f5a58d8
--- /dev/null
+++ b/noao/onedspec/odcombine/t_odcombine.x
@@ -0,0 +1,1071 @@
+include	<imhdr.h>
+include	<error.h>
+include	<mach.h>
+include	<mwset.h>
+include	<smw.h>
+include	"src/icombine.h"
+
+# Grouping options
+define	GROUP	"|all|images|apertures|"
+define	GRP_ALL		1
+define	GRP_IMAGES	2
+define	GRP_APERTURES	3
+
+# Mask formats
+define	MASKFORMATS	"|bpmpixel|bpmspectrum|"
+define	BPMPIX		1
+define	BPMSPEC		2
+
+# 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
+
+
+# T_ODCOMBINE - Combine spectra matched in world coordinates.
+# 
+# The input spectra are combined by medianing, averaging or summing with
+# optional rejection, scaling and weighting.  The combining algorithms and
+# other features are the same as those in IMCOMBINE.
+# 
+# The main difference with IMCOMBINE is that the spectra are first resampled
+# to a common grid of pixels in dispersion.  To do this each spectrum
+# is resampled to a temporary file which are then combined and deleted.
+# When bad pixels are used they are also resampled in the same way.
+# 
+# Since there can be multiple spectra per file (each with different
+# aperture numbers) there are three ways to group the spectra for combining.
+# One is by image where all spectra in the file are combined.  The second
+# is by aperture where the same aperture across multple files are combine.
+# The third is to combine all spectra independent of aperture or file.
+# 
+# The structure of the program is to first internally collect all the
+# spectra from each input file.  When combining by image this is done file
+# by file otherwise all the files are collected together.  The reason for
+# this is to avoid opening, search, reading, and closing the same file
+# for each aperture Then the spectra for one output are rebinned to a
+# common dispersion and written to temporary files.  The same is done for
+# bad pixel masks if used.  The spectra in the files are then combined.
+# Finally the temporary files are deleted.  The rebinning and combine are
+# repeated for each output spectrum.
+
+procedure t_odcombine()
+
+pointer	aps			# aperture ranges
+int	group			# grouping option
+
+int	mformat, mtype, mvalue
+int	i, j, index, naps
+pointer	im, mw, refim, shout
+pointer	sp, input, output, headers, bmask, rmask, nrmask, emask, sigma, logfile
+pointer	tmp, str, s, b, ns
+int	ilist1, ilist2, ilist, olist, hlist, blist, rlist, slist, nrlist, elist
+
+bool	clgetb()
+int	clgeti(), clgwrd()
+int	imtopen(), imtopenp(), imtgetim(), imtlen()
+real	clgetr()
+pointer	rng_open()
+errchk	shdr_open, odc_gspec, odc_rebin, odc_output, odc_combine
+
+include	"src/icombine.com"
+
+begin
+	# Allocate stack memory.  Note some of the variables are declared in
+	# the icombine common block but still need to be allocated here.
+
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (headers, SZ_FNAME, TY_CHAR)
+	call salloc (bmask, SZ_FNAME, TY_CHAR)
+	call salloc (rmask, SZ_FNAME, TY_CHAR)
+	call salloc (nrmask, SZ_FNAME, TY_CHAR)
+	call salloc (emask, SZ_FNAME, TY_CHAR)
+	call salloc (sigma, SZ_FNAME, TY_CHAR)
+	call salloc (logfile, SZ_FNAME, TY_CHAR)
+	call salloc (tmp, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (ictask, SZ_FNAME, TY_CHAR)
+	call salloc (expkeyword, SZ_FNAME, TY_CHAR)
+	call salloc (statsec, SZ_FNAME, TY_CHAR)
+	call salloc (gain, SZ_FNAME, TY_CHAR)
+	call salloc (rdnoise, SZ_FNAME, TY_CHAR)
+	call salloc (snoise, SZ_FNAME, TY_CHAR)
+
+	# Set the IMCOMBINE parameters.
+
+	call strcpy ("ODCOMBINE", Memc[ictask], SZ_FNAME)
+	ilist = imtopenp ("input")
+	olist = imtopenp ("output")
+	hlist = imtopenp ("headers")
+	blist = imtopenp ("bpmasks")
+	rlist = imtopenp ("rejmasks")
+	nrlist = imtopenp ("nrejmasks")
+	elist = imtopenp ("expmasks")
+	slist = imtopenp ("sigmas")
+	call clgstr ("logfile", Memc[logfile], SZ_FNAME)
+
+	project = false
+	combine = clgwrd ("combine", Memc[str], SZ_LINE, COMBINE)
+	reject = clgwrd ("reject", Memc[str], SZ_LINE, REJECT)
+	blank = clgetr ("blank")
+	call clgstr ("expname", Memc[expkeyword], SZ_FNAME)
+	call clgstr ("statsec", Memc[statsec], 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")
+	pclip = clgetr ("pclip")
+	flow = clgetr ("nlow")
+	fhigh = clgetr ("nhigh")
+	nkeep = clgeti ("nkeep")
+	grow = clgetr ("grow")
+	mclip = clgetb ("mclip")
+	sigscale = clgetr ("sigscale")
+
+	# Check parameters, map INDEFs, and set threshold flag
+
+	if (pclip == 0. && reject == 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 (pclip))
+	    pclip = -0.5
+	if (IS_INDEFR (flow))
+	    flow = 0
+	if (IS_INDEFR (fhigh))
+	    fhigh = 0
+	if (IS_INDEFR (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 ODCOMBINE specific parameters.
+
+	call clgstr ("apertures", Memc[str], SZ_LINE)
+	group = clgwrd ("group", Memc[input], SZ_FNAME, GROUP)
+
+	# Expand aperture list.
+	iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF))
+	    call error (1, "Error in aperture list")
+
+	# We need to know about the mask in order to resample them.
+	# This does not support specifying a mask by name or keyword.
+
+	mformat = clgwrd ("smaskformat", Memc[str], SZ_LINE, MASKFORMATS)
+	mtype = clgwrd ("smasktype", Memc[str], SZ_LINE, MASKTYPES)
+	if (mtype == 0)
+	    call error (1, "Unsupported masktype")
+	mvalue = clgeti ("smaskvalue")
+	if (mtype == M_BADBITS && mvalue == 0) 
+	    mtype = M_NONE
+	if (mtype == M_NONE)
+	    call clpstr ("masktype", "none")
+	else
+	    call clpstr ("masktype", "goodvalue")
+	call clputi ("maskvalue", 0)
+
+	# Check lists.
+	i = imtlen (ilist)
+	if (i == 0)
+	    call error (1, "No input images to combine")
+	switch (group) {
+	case GRP_ALL, GRP_APERTURES:
+	    if (imtlen (olist) != 1)
+		call error (1, "Wrong number of output images")
+	    if (imtlen (hlist) > 1)
+		call error (1, "Wrong number of header files")
+	    if (imtlen (blist) > 1)
+		call error (1, "Wrong number of bad pixel masks")
+	    if (imtlen (rlist) > 1)
+		call error (1, "Wrong number of rejection masks")
+	    if (imtlen (nrlist) > 1)
+		call error (1, "Wrong number of number rejected masks")
+	    if (imtlen (elist) > 1)
+		call error (1, "Wrong number of exposure masks")
+	    if (imtlen (slist) > 1)
+		call error (1, "Wrong number of sigma images")
+	case GRP_IMAGES:
+	    if (imtlen (olist) != i)
+		call error (1, "Wrong number of output images")
+	    if (imtlen (hlist) > 0 && imtlen (hlist) != i)
+		call error (1, "Wrong number of header files")
+	    if (imtlen (blist) > 0 && imtlen (blist) != i)
+		call error (1, "Wrong number of bad pixel masks")
+	    if (imtlen (rlist) > 0 && imtlen (rlist) != i)
+		call error (1, "Wrong number of rejection masks")
+	    if (imtlen (nrlist) > 0 && imtlen (nrlist) != i)
+		call error (1, "Wrong number of number rejected masks")
+	    if (imtlen (elist) > 0 && imtlen (elist) != i)
+		call error (1, "Wrong number of exposure masks")
+	    if (imtlen (slist) > 1 && imtlen (slist) != i)
+		call error (1, "Wrong number of sigma images")
+	}
+
+	# Set temporary output rootname.
+	call mktemp ("tmp", Memc[tmp], SZ_FNAME)
+
+	# Loop through input images.
+	index = 0
+	while (imtgetim (ilist, Memc[input], SZ_FNAME) != EOF) {
+
+	    # Get all requested apertures from an image.  When not grouping
+	    # by image go through all images and exhaust the input list.
+
+	    naps = 0
+	    repeat {
+		iferr (call odc_gspec (Memc[input], aps, group, mtype, mformat,
+		    s, b, ns, naps)) {
+		    if (group == GRP_IMAGES) {
+			call erract (EA_WARN)
+			next
+		    } else {
+			call erract (EA_ERROR)
+		    }
+		}
+		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
+	    }
+
+	    # Create each output spectrum.  This involves rebinning to
+	    # temporary files, combining, and cleaning up.  The files are
+	    # deleted in the odc_combine routine.
+
+	    do i = 1, naps {
+
+		# Set the output dispersion in a temporary template image.
+		call odc_output (SH(s,i,1), NS(ns,i), Memc[tmp], im, mw, refim)
+		call shdr_open (im, mw, i, 1, INDEFI, SHDATA, shout)
+
+		# Rebin the spectra.
+		call odc_rebin (im, shout, SH(s,i,1), SH(b,i,1), NS(ns,i),
+		    mformat, mtype, mvalue, Memc[tmp])
+
+		# Close and delete the template image.
+		call shdr_close (shout)
+		call smw_close (mw)
+		call imunmap (im)
+		call imunmap (refim)
+		iferr (call imdelete (Memc[tmp]))
+		    ;
+
+		# Set lists to be combined.
+		call sprintf (Memc[str], SZ_LINE, "%s.*\\[^x]")
+		    call pargstr (Memc[tmp])
+		ilist1 = imtopen (Memc[str])
+		if (mtype != NONE) {
+		    call sprintf (Memc[str], SZ_LINE, "%sbpm.*\\[^x]")
+			call pargstr (Memc[tmp])
+		    ilist2 = imtopen (Memc[str])
+		} else
+		    ilist2 = imtopen ("")
+
+		# Set output names.
+		switch (group) {
+		case GRP_ALL:
+		    index = 1
+		    j = INDEFI
+		case GRP_IMAGES:
+		    index = index + 1
+		    j = INDEFI
+		case GRP_APERTURES:
+		    index = 1
+		    j = AP(SH(s,i,1))
+		}
+		call odc_imtgetim (olist, index, j, Memc[output], SZ_FNAME)
+		call odc_imtgetim (hlist, index, j, Memc[headers], SZ_FNAME)
+		call odc_imtgetim (blist, index, j, Memc[bmask], SZ_FNAME)
+		call odc_imtgetim (rlist, index, j, Memc[rmask], SZ_FNAME)
+		call odc_imtgetim (nrlist, index, j, Memc[nrmask], SZ_FNAME)
+		call odc_imtgetim (elist, index, j, Memc[emask], SZ_FNAME)
+		call odc_imtgetim (slist, index, j, Memc[sigma], SZ_FNAME)
+
+		# Combine and delete the lists.
+		iferr (call odc_combine (ilist1, ilist2, Memc[output],
+		    Memc[headers], Memc[bmask], Memc[rmask], Memc[nrmask],
+		    Memc[emask], Memc[sigma], Memc[logfile], YES))
+		    call erract (EA_WARN)
+
+		call imtclose (ilist1)
+		call imtclose (ilist2)
+	    }
+
+	    # Free all the spectrum data structures.
+	    call odc_fspec (s, b, ns, naps)
+	}
+
+	# Finish up.
+	call rng_close (aps)
+	call imtclose (ilist)
+	call imtclose (olist)
+	call imtclose (hlist)
+	call imtclose (blist)
+	call imtclose (rlist)
+	call imtclose (nrlist)
+	call imtclose (elist)
+	call imtclose (slist)
+	call sfree (sp)
+end
+
+
+# ODC_GSPEC -- Get spectra from an input image.
+#
+# This allocates and sets arrays of spectrum structures.  There is an array
+# for each "output" aperture and the number of elements is given by another
+# array (ns).  The number of output apertures is given by naps.  Note that each
+# call to this accumulates new spectra.
+
+procedure odc_gspec (input, aps, group, mtype, mformat, s, b, ns, naps)
+
+char	input[ARB]		#I Input spectrum file
+pointer	aps			#I Apertures to select
+int	group			#I Grouping for combining
+int	mtype			#I Mask type
+int	mformat			#I Mask format
+pointer	s			#U Spectra data structure
+pointer	b			#U Spectra data structure for pixel masks
+int	ns			#U Number of spectra per group
+int	naps			#U Number of output apertures
+
+int	i, j, k, n
+pointer	im, mw, sh, sh1, bpm, err
+
+bool	rng_elementi()
+pointer	immap(), smw_openim()
+errchk	immap, smw_openim, shdr_open
+
+begin
+	# Map the input spectrum file. Check format.
+	im = immap (input, READ_ONLY, 0)
+	mw = smw_openim (im)
+	if (SMW_FORMAT(mw) != SMW_ES && SMW_FORMAT(mw) != SMW_MS) {
+	    call smw_close (mw)
+	    call imunmap (im)
+	    call salloc (err, SZ_LINE, TY_CHAR)
+	    call sprintf (Memc[err], SZ_LINE,
+	       "Unsupported spectral format (%s)")
+	       call pargstr (input)
+	    call error (1, Memc[err])
+	}
+	sh = NULL
+
+	# Get the associated mask if requested.  It is not an error if there
+	# is mask.
+
+	if (mtype == M_NONE)
+	    bpm = NULL
+	else {
+	    switch (mformat) {
+	    case BPMPIX, BPMSPEC:
+		call malloc (bpm, SZ_FNAME, TY_CHAR)
+		iferr (call imgstr (im, "BPM", Memc[bpm], SZ_FNAME))
+		    call mfree (bpm, TY_CHAR)
+	    default:
+	        bpm = NULL
+	    }
+	}
+
+	# Select the requested apertures and group by output aperture.
+	do i = 1, SMW_NSPEC(mw) {
+	    call shdr_open (im, mw, i, 1, INDEFI, SHDATA, sh)
+	    if (!rng_elementi (aps, AP(sh)))
+		next
+
+	    if (group == GRP_APERTURES) {
+		for (j=1; j<=naps; j=j+1)
+		    if (AP(sh) == AP(SH(s,j,1)))
+			break
+		n = 10
+	    } else {
+		j = 1
+		n = 1
+	    }
+
+	    if (naps == 0) {
+		call calloc (s, n, TY_POINTER)
+		call calloc (b, n, TY_POINTER)
+		call calloc (ns, n, TY_INT)
+	    } else if (j > naps && mod (naps, n) == 0) {
+		call realloc (s, naps+n, TY_POINTER)
+		call realloc (b, naps+n, TY_POINTER)
+		call realloc (ns, naps+n, TY_INT)
+		call aclri (Memi[s+naps], n)
+		call aclri (Memi[b+naps], n)
+		call aclri (Memi[ns+naps], n)
+	    }
+	    if (j > naps)
+		naps = naps + 1
+	    n = NS(ns,j)
+	    if (n == 0) {
+		call malloc (Memi[s+j-1], 10, TY_POINTER)
+		call malloc (Memi[b+j-1], 10, TY_POINTER)
+	    } else if (mod (n, 10) == 0) {
+		call realloc (Memi[s+j-1], n+10, TY_POINTER)
+		call realloc (Memi[b+j-1], n+10, TY_POINTER)
+	    }
+
+	    n = n + 1
+	    SH(s,j,n) = NULL
+	    SH(b,j,n) = NULL
+	    call shdr_copy (sh, SH(s,j,n), NO)
+	    NS(ns,j) = n
+	}
+
+	call imunmap (IM(sh))
+	MW(sh) = NULL
+	call shdr_close (sh)
+
+	# Get BPMs if defined.
+	if (bpm != NULL) {
+	    im = immap (Memc[bpm], READ_ONLY, 0)
+	    mw = smw_openim (im)
+	    sh = NULL
+
+	    switch (mformat) {
+	    case BPMPIX:
+		do j = 1, naps {
+		    n = NS(ns,j)
+		    sh1 = SH(s,j,n)
+		    if (sh1 == NULL)
+			next
+		    k = LINDEX(sh1,1)
+		    call shdr_open (im, mw, k, 1, INDEFI, SHDATA, sh)
+		    if (LINDEX(sh,1) != k)
+			next
+		    call shdr_copy (sh1, SH(b,j,n), YES)
+		    sh1 = SH(b,j,n)
+		    call strcpy (IMNAME(sh), IMNAME(sh1), LEN_SHDRS)
+		    call strcpy (IMSEC(sh), IMSEC(sh1), LEN_SHDRS)
+		    call strcpy (TITLE(sh), TITLE(sh1), LEN_SHDRS)
+		    call amovi (LINDEX(sh,1), LINDEX(sh1,1), 2)
+		    call amovi (PINDEX(sh,1), PINDEX(sh1,1), 2)
+		    APINDEX(sh) = APINDEX(sh1)
+		    call amovr (Memr[SY(sh)], Memr[SY(sh1)],
+			min (SN(sh), SN(sh1)))
+		}
+		call smw_close (mw)
+	    case BPMSPEC:
+		do j = 1, naps {
+		    n = NS(ns,j)
+		    sh1 = SH(s,j,n)
+		    if (sh1 == NULL)
+			next
+		    k = AP(sh1)
+		    call shdr_open (im, mw, 1, 1, k, SHDATA, sh)
+		    if (AP(sh) != k)
+			next
+		    call shdr_copy (sh, SH(b,j,n), NO)
+		}
+	    }
+
+	    call imunmap (IM(sh))
+	    MW(sh) = NULL
+	    call shdr_close (sh)
+	    call mfree (bpm, TY_CHAR)
+	}
+end
+
+
+
+# ODC_FSPEC -- Free spectrum data structures.
+
+procedure odc_fspec (s, b, ns, naps)
+
+pointer	s		#U Spectrum data structures
+pointer	b		#U BPM data structures
+pointer	ns		#U Number of spectra per output aperture
+int	naps		#I Number of output apertures
+
+int	i, j, k, l
+pointer	sh, mw
+
+begin
+	# Find all the distinct SMW pointers and free them.
+	# Then free all the spectrum data pointers.
+
+	do j = 1, naps {
+	    do i = 1, NS(ns,j) {
+		sh = SH(s,j,i)
+		if (sh == NULL)
+		    next
+		mw = MW(sh)
+		if (mw != NULL) {
+		    do k = 1, naps  {
+			do l = 1, NS(ns,k) {
+			    sh = SH(s,k,l)
+			    if (sh == NULL)
+				next
+			    if (MW(sh) == mw)
+				MW(sh) = NULL
+			}
+		    }
+		    call smw_close (mw)
+		}
+	    }
+	}
+	do j = 1, naps {
+	    do i = 1, NS(ns,j) {
+		sh = SH(s,j,i)
+		if (sh == NULL)
+		    next
+		call shdr_close (sh)
+	    }
+	    call mfree (Memi[s+j-1], TY_POINTER)
+	}
+	call mfree (s, TY_POINTER)
+
+	do j = 1, naps {
+	    do i = 1, NS(ns,j) {
+		sh = SH(b,j,i)
+		if (sh == NULL)
+		    next
+		mw = MW(sh)
+		if (mw != NULL) {
+		    do k = 1, naps  {
+			do l = 1, NS(ns,k) {
+			    sh = SH(b,k,l)
+			    if (sh == NULL)
+				next
+			    if (MW(sh) == mw)
+				MW(sh) = NULL
+			}
+		    }
+		    call smw_close (mw)
+		}
+	    }
+	}
+	do j = 1, naps {
+	    do i = 1, NS(ns,j) {
+		sh = SH(b,j,i)
+		if (sh == NULL)
+		    next
+		call shdr_close (sh)
+	    }
+	    call mfree (Memi[b+j-1], TY_POINTER)
+	}
+	call mfree (b, TY_POINTER)
+
+	call mfree (ns, TY_INT)
+end
+
+
+# ODC_REBIN -- Rebin spectra and masks.
+
+procedure odc_rebin (refim, shout, s, b, n, mformat, mtype, mvalue, output)
+
+pointer	refim			#I Output reference image
+pointer	shout			#I Output spectrum structure
+pointer	s[ARB]			#I Array of spectrum structures
+pointer	b[ARB]			#I Array of BPM spectrum structures
+int	n			#I Number of spectra
+int	mformat			#I Mask format
+int	mtype			#I Mask type
+int	mvalue			#I Mask value
+char	output[ARB]		#I Output rootname
+
+int	i, j, k, p1, p2, npts
+double	c[3], d[3,3]
+pointer	sh, bpm, im, mw
+pointer	sp, str
+
+int	mw_stati()
+double	shdr_lw(), shdr_wl()
+pointer	immap(), mw_openim(), impl1r()
+errchk	immap, mw_openim, impl1r
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	j = 0
+	do i = 1, n {
+	    sh = s[i]
+	    bpm = b[i]
+
+	    # Determine limits of input spectrum relative to the output
+	    # spectrum.
+
+	    c[1] = shdr_wl (shout, shdr_lw (sh, double(0.5)))
+	    c[2] = shdr_wl (shout, shdr_lw (sh, double(SN(sh)+0.5)))
+	    p1 = max (1, nint (min (c[1], c[2]) + 0.01))
+	    p2 = min (SN(shout), nint (max (c[1], c[2]) - 0.01))
+	    npts = p2 - p1 + 1
+	    if (npts < 1)
+		next
+	    p1 = 1 - p1
+
+	    # Rebin the spectra and masks.
+	    call shdr_rebin (sh, shout)
+	    call odc_bpm (bpm, shout, mtype, mvalue)
+
+	    # Write the results.  We only write the part of the output
+	    # contained by the input spectrum and then let the combining deal
+	    # with the origin offsets.  This is done by setting the physical
+	    # pixel coordinate system to match the desired output system.
+	    # The main reason for this it to make the output of bounds
+	    # pixel implicitly bad or excluded.
+
+	    j = j + 1
+	    call sprintf (Memc[str], SZ_LINE, "%s.%04d")
+		call pargstr (output)
+		call pargi (j)
+	    im = immap (Memc[str], NEW_COPY, refim)
+	    call sprintf (Memc[str], SZ_LINE, "%s%s(%s)")
+		call pargstr (IMNAME(sh))
+		call pargstr (IMSEC(sh))
+		call pargi (AP(sh))
+	    call imastr (im, "ICFNAME", Memc[str])
+	    IM_LEN(im,1) = npts
+	    if (p1 != 0) {
+		mw = mw_openim (im)
+		k = mw_stati (mw, MW_NPHYSDIM)
+		call mw_gltermd (mw, d, c, k)
+		c[1] = c[1] + p1
+		call mw_sltermd (mw, d, c, k)
+		call mw_saveim (mw, im)
+	    }
+	    call amovr (Memr[SY(sh)-p1], Memr[impl1r(im)], npts)
+	    if (bpm != NULL) {
+		switch (mformat) {
+		case BPMPIX:
+		    call sprintf (Memc[str], SZ_LINE, "%s%s")
+			call pargstr (IMNAME(bpm))
+			call pargstr (IMSEC(bpm))
+		case BPMSPEC:
+		    call sprintf (Memc[str], SZ_LINE, "%s%s(%s)")
+			call pargstr (IMNAME(bpm))
+			call pargstr (IMSEC(bpm))
+			call pargi (AP(bpm))
+		}
+		call imastr (im, "ICBPM", Memc[str])
+		call sprintf (Memc[str], SZ_LINE, "%sbpm.%04d")
+		    call pargstr (output)
+		    call pargi (j)
+		call imastr (im, "BPM", Memc[str])
+	    } else iferr (call imdelf (im, "BPM"))
+		;
+	    call imunmap (im)
+
+	    if (bpm == NULL)
+		next
+
+	    im = immap (Memc[str], NEW_COPY, refim)
+	    IM_PIXTYPE(im) = TY_INT
+	    IM_LEN(im,1) = npts
+	    if (p1 != 0) {
+		mw = mw_openim (im)
+		k = mw_stati (mw, MW_NPHYSDIM)
+		call mw_gltermd (mw, d, c, k)
+		c[1] = c[1] + p1
+		call mw_sltermd (mw, d, c, k)
+		call mw_saveim (mw, im)
+	    }
+	    call amovr (Memr[SY(bpm)-p1], Memr[impl1r(im)], npts)
+	    iferr (call imdelf (im, "BPM"))
+	         ;
+	    call imunmap (im)
+	}
+
+	call sfree (sp)
+end
+
+
+# ODC_BPM -- Rebin the bad pixel masks.
+#
+# Even though the input mask can be specified by good or bad values or bits
+# the rebinned mask is created as a boolean mask.  Note that the rebinning
+# is done by setting a large mask value and then values computed from good
+# and bad pixels will have some intermediate value which we then threshold
+# to define good and bad.
+
+procedure odc_bpm (sh, shout, mtype, mvalue)
+
+pointer	sh		#I SHDR pointer for mask spectrum
+pointer	shout		#I SHDR pointer for template output spectrum
+int	mtype		#I Mask type
+int	mvalue		#I Mask value
+
+int	i, n, val, and()
+pointer	ptr
+
+begin
+	if (sh == NULL)
+	    return
+
+	n = SN(sh)
+	ptr = SY(sh)
+	switch (mtype) {
+	case M_GOODVAL:
+	    do i = 1, n {
+	        val = nint (Memr[ptr])
+	        if (val == mvalue)
+		    Memr[ptr] = 0
+		else
+		    Memr[ptr] = 1000
+		ptr = ptr + 1
+	    }
+	case M_BADVAL:
+	    do i = 1, n {
+	        val = nint (Memr[ptr])
+	        if (val != mvalue)
+		    Memr[ptr] = 0
+		else
+		    Memr[ptr] = 1000
+		ptr = ptr + 1
+	    }
+	case M_GOODBITS:
+	    do i = 1, n {
+	        val = nint (Memr[ptr])
+	        if (and (val, mvalue) != 0)
+		    Memr[ptr] = 0
+		else
+		    Memr[ptr] = 1000
+		ptr = ptr + 1
+	    }
+	case M_BADBITS:
+	    do i = 1, n {
+	        val = nint (Memr[ptr])
+	        if (and (val, mvalue) == 0)
+		    Memr[ptr] = 0
+		else
+		    Memr[ptr] = 1000
+		ptr = ptr + 1
+	    }
+	}
+
+	call shdr_rebin (sh, shout)
+
+	n = SN(sh)
+	ptr = SY(sh)
+	do i = 1, n {
+	    val = nint (Memr[ptr])
+	    if (val < 10)
+		Memr[ptr] = 0
+	    else
+		Memr[ptr] = 1
+	    ptr = ptr + 1
+	}
+end
+
+
+# ODC_OUTPUT - Set the output spectrum.
+
+procedure odc_output (sh, ns, output, im, mw, refim)
+
+pointer	sh[ARB]			# spectra structures
+int	ns			# number of spectra
+char	output[SZ_FNAME]	# output spectrum name
+pointer	im			# output IMIO pointer
+pointer	mw			# output MWCS pointer
+pointer	refim			# reference image for output image
+
+int	ap, beam, dtype, nw, axis[2]
+double	w1, dw, z
+real	aplow[2], aphigh[2]
+pointer	coeff
+pointer	immap(), mw_open(), smw_openim()
+errchk	immap, smw_openim
+data	axis/1,2/
+
+begin
+	coeff = NULL
+
+	# Create output image using the first input image as a reference
+	refim = immap (IMNAME(sh[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)
+
+	IM_NDIM(im) = 1
+	call imaddi (im, "SMW_NDIM", IM_NDIM(im))
+	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)
+
+	call smw_gwattrs (MW(sh[1]), APINDEX(sh[1]), 1, ap, beam, dtype,
+	    w1, dw, nw, z, aplow, aphigh, coeff)
+	call odc_default (sh, ns, dtype, w1, dw, nw, z, Memc[coeff])
+	call smw_swattrs (mw, 1, 1, ap, beam, dtype,
+	    w1, dw, nw, z, aplow, aphigh, Memc[coeff])
+	call smw_sapid (mw, 1, 1, TITLE(sh[1]))
+
+	IM_LEN(im,1) = nw
+
+	# Set MWCS header.
+	call smw_saveim (mw, im)
+	call smw_close (mw)
+	mw = smw_openim (im)
+
+	call mfree (coeff, TY_CHAR)
+end
+
+
+# ODC_DEFAULT - Set default values for the starting wavelength, ending
+# wavelength, wavelength increment and spectrum length for the output
+# spectrum.
+
+procedure odc_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")) {
+	    # For now we don't allow non-linear dispersions because the
+	    # generic combine routines don't understand multispec.
+	    if (dtype == DCFUNC) {
+	        dtype = DCLINEAR
+		coeff[1] = EOS
+		z = 0.
+	    }
+
+	    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
+
+
+# ODC_IMTGETIM -- Set output image from an list of root names.
+
+procedure odc_imtgetim (list, index, aperture, image, maxch)
+
+int	list			#I List of images
+int	index			#I List index
+int	aperture		#I Aperture
+char	image[maxch]		#O Image name
+int	maxch			#I Maximum character for image
+
+pointer	sp, root, extn
+
+int	imtrgetim()
+
+begin
+	if (imtrgetim (list, index, image, maxch) == EOF) {
+	    image[1] = EOS
+	    return
+	}
+
+	if (aperture == INDEFI)
+	    return
+
+	call smark (sp)
+	call salloc (root, SZ_FNAME, TY_CHAR)
+	call salloc (extn, SZ_FNAME, TY_CHAR)
+
+	call iki_init()
+	call iki_parse (image, Memc[root], Memc[extn])
+	if (Memc[extn] == EOS) {
+	    call sprintf (image, maxch, "%s.%04d")
+		call pargstr (Memc[root])
+		call pargi (aperture)
+	} else {
+	    call sprintf (image, maxch, "%s.%04d.%s")
+		call pargstr (Memc[root])
+		call pargi (aperture)
+		call pargstr (Memc[extn])
+	}
+
+	call sfree (sp)
+end
+
+
+# ODC_COMBINE -- Combine the spectra by calling the IMCOMBINE source.
+
+procedure odc_combine (slist, blist, output, headers, bmask, rmask, nrmask,
+	emask, sigma, logfile, delete)
+
+int	slist			#I List of 1D spectra to combine
+int	blist			#I List of 1D bad pixel spectra
+char	output[ARB]		#I Output combined spectrum
+char	headers[ARB]		#I Output headers
+char	bmask[ARB]		#I Output bad pixel mask
+char	rmask[ARB]		#I Output rejection mask
+char	nrmask[ARB]		#I Output number rejected mask
+char	emask[ARB]		#I Ouput exposure time mask
+char	sigma[ARB]		#I Output sigma
+char	logfile[ARB]		#I Logfile
+int	delete			#I Delete input spectra?
+
+int	n
+pointer	sp, fname, scales, zeros, wts
+
+int	imtlen(), imtgetim()
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Allocate and initialize scaling factors.
+	n = imtlen (slist)
+	call salloc (scales, 3*n, TY_REAL)
+	zeros = scales + n
+	wts = scales + 2 * n
+	call amovkr (INDEFR, Memr[scales], 3*n)
+
+	# Combine.
+	iferr (call icombine (slist, output, headers, bmask, rmask,
+	    nrmask, emask, sigma, logfile, Memr[scales], Memr[zeros],
+	    Memr[wts], NO, NO))
+	    call erract (EA_WARN)
+
+	# Delete the files.
+	if (delete == YES) {
+	    call imtrew (slist)
+	    while (imtgetim (slist, Memc[fname], SZ_FNAME) != EOF)
+		call imdelete (Memc[fname])
+	    call imtrew (blist)
+	    while (imtgetim (blist, Memc[fname], SZ_FNAME) != EOF)
+		call imdelete (Memc[fname])
+	}
+
+	call sfree (sp)
+end