Initial commit

22 files changed, 4907 insertions, 0 deletions

diff --git a/noao/onedspec/dispcor/dcio.x b/noao/onedspec/dispcor/dcio.x
new file mode 100644
index 00000000..b700da6a
--- /dev/null
+++ b/noao/onedspec/dispcor/dcio.x
@@ -0,0 +1,1155 @@
+include	<error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<pkg/dttext.h>
+include	<smw.h>
+include	<units.h>
+include	"dispcor.h"
+
+# Symbol table structure for the dispersion solutions.
+define	LEN_DC		11		# Length of dispersion solution struct.
+define	DC_FORMAT	Memi[$1]	# Type of dispersion
+define	DC_PAPS		Memi[$1+1]	# Pointer to aperture numbers
+define	DC_PAPCEN	Memi[$1+2]	# Pointer to aperture centers
+define	DC_PUN		Memi[$1+3]	# Pointer to units
+define	DC_PSHIFT	Memi[$1+4]	# Pointer to shifts
+define	DC_PCOEFF	Memi[$1+5]	# Pointer to coefficients
+define	DC_NAPS		Memi[$1+6]	# Number of apertures
+define  DC_OFFSET       Memi[$1+7]      # Aperture to order offset
+define  DC_SLOPE        Memi[$1+8]      # Aperture to order slope
+define  DC_COEFFS       Memi[$1+9]      # Dispersion coefficients
+define  DC_SHIFT        Memr[P2R($1+10)]# Dispersion function shift
+
+
+# DC_OPEN    -- Initialize the dispersion data structures
+# DC_CLOSE   -- Close the dispersion data structures
+# DC_GMS     -- Get a multispec spectrum
+# DC_GMSDB   -- Get a multispec dispersion database entry
+# DC_REFSHFT -- Get a reference shift
+# DC_GEC     -- Get an echelle spectrum
+# DC_GECDB   -- Get an echelle dispersion database entry
+# DC_ECMS    -- Convert echelle dispersion coeffs to multispec coeffs
+
+
+# DC_OPEN -- Initialize the dispersion routines.  This consists
+# of opening a symbol table for the dispersion solution functions.  A
+# symbol table is used since the same dispersion reference (arc image)
+# may be be used multiple times and the database access is slow.
+
+procedure dc_open (stp, db)
+
+pointer	stp			# Symbol table pointer
+char	db[SZ_FNAME]		# Database name
+
+pointer	sym, stopen(), stenter(), stpstr()
+
+begin
+	stp = stopen ("disp", 10, 10, 10*SZ_FNAME)
+	sym = stenter (stp, "database", 1)
+	Memi[sym] = stpstr (stp, db, 0)
+end
+
+
+# DC_CLOSE -- Close the dispersion data structures.
+
+procedure dc_close (stp)
+
+int	i
+pointer stp, sym, sthead, stnext
+
+begin
+	# Close each dispersion function and then the symbol table.
+	for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {
+	    if (DC_FORMAT(sym) == 1) {
+		do i = 1, DC_NAPS(sym) {
+		    call un_close (Memi[DC_PUN(sym)+i-1])
+		    call mfree (Memi[DC_PCOEFF(sym)+i-1], TY_DOUBLE)
+		}
+		call mfree (DC_PAPS(sym), TY_INT)
+		call mfree (DC_PAPCEN(sym), TY_REAL)
+		call mfree (DC_PUN(sym), TY_POINTER)
+		call mfree (DC_PSHIFT(sym), TY_DOUBLE)
+		call mfree (DC_PCOEFF(sym), TY_POINTER)
+	    } else if (DC_FORMAT(sym) == 2) {
+		call un_close (DC_PUN(sym))
+		call mfree (DC_COEFFS(sym), TY_DOUBLE)
+	    }
+	}
+	call stclose (stp)
+end
+
+
+# DC_GMS -- Get a multispec spectrum.  This consists of mapping the image
+# and setting a MWCS coordinate transformation.  If not dispersion corrected
+# the dispersion function is found in the database for the reference
+# spectra and set in the SMW.
+
+procedure dc_gms (spec, im, smw, stp, ignoreaps, ap, fd1, fd2)
+
+char	spec[ARB]	#I Spectrum name
+pointer	im		#I IMIO pointer
+pointer	smw		#I SMW pointer
+pointer	stp		#I Dispersion symbol table
+int	ignoreaps	#I Ignore aperture numbers?
+pointer	ap		#O Aperture data structure
+int	fd1		#I Logfile descriptor
+int	fd2		#I Logfile descriptor
+
+double	wt1, wt2, dval
+int	i, j, k, k1, k2, l, dc, sfd, naps, naps1, naps2, ncoeffs
+pointer	sp, str1, str2, papcen, pshift, coeffs, ct1, ct2, un, un1, un2
+pointer	paps1, paps2, punits1, punits2, pshift1, pshift2, pcoeff1, pcoeff2
+
+bool	un_compare()
+double	smw_c1trand()
+int	imaccf(), nscan(), stropen()
+pointer	smw_sctran(), un_open()
+errchk	dc_gmsdb, dc_refshft, imgstr, smw_sctran, un_open
+
+define	done_	90
+
+begin
+	call smark (sp)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (str2, SZ_LINE, TY_CHAR)
+
+	# Set WCS attributes
+	naps = IM_LEN(im,2)
+	call calloc (ap, LEN_AP(naps), TY_STRUCT)
+	do i = 1, naps {
+	    DC_PL(ap,i) = i
+	    DC_CF(ap,i) = NULL
+	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), DC_CF(ap,i))
+	    if (i == 1) {
+		iferr (call mw_gwattrs (SMW_MW(smw,0), 1, "units", Memc[str1],
+		    SZ_LINE))
+		    Memc[str1] = EOS
+		DC_UN(ap,i) = un_open (Memc[str1])
+	    }
+	    dc = DC_DT(ap,i)
+	}
+
+	# Check if the spectra have been dispersion corrected
+	# by an earlier version of DISPCOR.  If so then don't allow
+	# another database dispersion correction.  This assumes all
+	# spectra have the same dispersion type.  Check for a
+	# reference spectrum.
+
+	#if ((imaccf (im, "REFSPEC1") == NO) ||
+	#    (dc > -1 && imaccf (im, "DCLOG1") == NO)) {
+	if (imaccf (im, "REFSPEC1") == NO) {
+	    if (fd1 != NULL) {
+		call fprintf (fd1,
+		    "%s: Resampling using current coordinate system\n")
+		    call pargstr (spec)
+	    }
+	    if (fd2 != NULL) {
+		call fprintf (fd2,
+		    "%s: Resampling using current coordinate system\n")
+		    call pargstr (spec)
+	    }
+	    goto done_
+	}
+	    
+	# Get the reference spectra dispersion function from the database
+	# and determine a reference shift.
+
+	iferr {
+	    call imgstr (im, "REFSPEC1", Memc[str1], SZ_LINE)
+	    call sscan (Memc[str1])
+	    call gargwrd (Memc[str1], SZ_LINE)
+	    call gargd (wt1)
+	    if (nscan() == 1)
+		wt1 = 1.
+	} then {
+	    call strcpy (spec, Memc[str1], SZ_FNAME)
+	    wt1 = 1.
+	}
+	iferr (call dc_gmsdb (Memc[str1], stp, paps1, papcen, punits1, pshift,
+	    pcoeff1, naps1)) {
+	    call sfree (sp)
+	    call erract (EA_ERROR)
+	}
+	call salloc (pshift1, naps1, TY_DOUBLE)
+	call amovd (Memd[pshift], Memd[pshift1], naps1)
+	if (fd1 != NULL) {
+	    call fprintf (fd1, "%s: REFSPEC1 = '%s %.8g'\n")
+		call pargstr (spec)
+		call pargstr (Memc[str1])
+		call pargd (wt1)
+	}
+	if (fd2 != NULL) {
+	    call fprintf (fd2, "%s: REFSPEC1 = '%s %.8g'\n")
+		call pargstr (spec)
+		call pargstr (Memc[str1])
+		call pargd (wt1)
+	}
+
+	iferr (call  dc_refshft (spec, stp, Memc[str1], "REFSHFT1", im,
+	    Memi[paps1], Memr[papcen], Memd[pshift1], naps1, fd1, fd2))
+	    ;
+
+	iferr {
+            call imgstr (im, "REFSPEC2", Memc[str1], SZ_LINE)
+            call sscan (Memc[str1])
+            call gargwrd (Memc[str1], SZ_LINE)
+            call gargd (wt2)
+            if (nscan() == 1)
+                wt2 = 1.
+            call dc_gmsdb (Memc[str1], stp, paps2, papcen, punits2, pshift,
+		pcoeff2, naps2)
+            call salloc (pshift2, naps2, TY_DOUBLE)
+            call amovd (Memd[pshift], Memd[pshift2], naps2)
+	    if (fd1 != NULL) {
+		call fprintf (fd1, "%s: REFSPEC2 = '%s %.8g'\n")
+		    call pargstr (spec)
+		    call pargstr (Memc[str1])
+		    call pargd (wt2)
+	    }
+	    if (fd2 != NULL) {
+		call fprintf (fd2, "%s: REFSPEC2 = '%s %.8g'\n")
+		    call pargstr (spec)
+		    call pargstr (Memc[str1])
+		    call pargd (wt2)
+	    }
+	    iferr (call  dc_refshft (spec, stp, Memc[str1],
+		"REFSHFT2", im, Memi[paps2], Memr[papcen], Memd[pshift2],
+		naps2, fd1, fd2))
+	        ;
+	} then
+	    wt2 = 0.
+
+	# Adjust weights to unit sum.
+	dval = wt1 + wt2
+	wt1 = wt1 / dval
+	wt2 = wt2 / dval
+
+	# Enter dispersion function in the MWCS.
+	do i = 1, naps {
+	    j = DC_AP(ap,i)
+	    for (k1=0; k1<naps1 && Memi[paps1+k1]!=j; k1=k1+1)
+		;
+	    if (k1 == naps1)
+		for (k1=0; k1<naps1 && !IS_INDEFI(Memi[paps1+k1]); k1=k1+1)
+		    ;
+	    if (k1 == naps1) {
+		if (ignoreaps == YES)
+		    k1 = 0
+		else {
+		    call sprintf (Memc[str1], SZ_LINE,
+			"%s - Missing reference for aperture %d")
+			call pargstr (spec)
+			call pargi (j)
+		    call fatal (1, Memc[str1])
+		}
+	    }
+	    un1 = Memi[punits1+k1]
+
+	    # The following assumes some knowledge of the data structure in
+	    # order to shortten the the attribute string.
+	    coeffs = Memi[pcoeff1+k1]
+	    if (coeffs == NULL) {
+		if (DC_DT(ap,i) == 2) {
+		    sfd = NULL
+		    if (wt2 <= 0.)
+			call sshift1 (Memd[pshift1+k1], DC_CF(ap,i))
+		} else {
+		    ncoeffs = 6
+		    l = 20 * (ncoeffs + 2)
+		    if (wt2 > 0.)
+			l = 2 * l
+		    call realloc (DC_CF(ap,i), l, TY_CHAR)
+		    call aclrc (Memc[DC_CF(ap,i)], l)
+		    sfd = stropen (Memc[DC_CF(ap,i)], l, NEW_FILE)
+		    call fprintf (sfd, "%.8g %g")
+			call pargd (wt1)
+			call pargd (Memd[pshift1+k1])
+		    dval = DC_DW(ap,i) * (DC_NW(ap,i) - 1) / 2.
+		    call fprintf (sfd, " 1 2 1 %d %g %g")
+			call pargi (DC_NW(ap,i))
+			call pargd (DC_W1(ap,i) + dval)
+			call pargd (dval)
+		}
+	    } else {
+		ncoeffs = nint (Memd[coeffs])
+		l = 20 * (ncoeffs + 2)
+		if (wt2 > 0.)
+		    l = 2 * l
+		call realloc (DC_CF(ap,i), l, TY_CHAR)
+		call aclrc (Memc[DC_CF(ap,i)], l)
+		sfd = stropen (Memc[DC_CF(ap,i)], l, NEW_FILE)
+		call fprintf (sfd, "%.8g %g %d %d")
+		    call pargd (wt1)
+		    call pargd (Memd[pshift1+k1])
+		    call pargi (nint (Memd[coeffs+1]))
+		    call pargi (nint (Memd[coeffs+2]))
+		do k = 3, ncoeffs {
+		    call fprintf (sfd, " %.15g")
+			call pargd (Memd[coeffs+k])
+		}
+	    }
+
+	    if (wt2 > 0.) {
+		for (k2=0; k2<naps2 && Memi[paps2+k2]!=j; k2=k2+1)
+		    ;
+		if (k2 == naps2)
+		    for (k2=0; k2<naps2 && !IS_INDEFI(Memi[paps2+k2]); k2=k2+1)
+			;
+		if (k2 == naps2) {
+		    if (ignoreaps == YES)
+			k2 = 0
+		    else {
+			call sprintf (Memc[str1], SZ_LINE,
+			    "%s - Missing reference for aperture %d")
+			    call pargstr (spec)
+			    call pargi (j)
+			if (sfd != NULL)
+			    call strclose (sfd)
+			call sfree (sp)
+			call fatal (1, Memc[str1])
+		    }
+		}
+		un2 = Memi[punits2+k2]
+		if (!un_compare (un1, un2)) {
+		    call sfree (sp)
+		    call error (2,
+			"Can't combine references with different units")
+		}
+		if (DC_DT(ap,i)==2 && !(coeffs==NULL&&Memi[pcoeff2+k2]==NULL)) {
+		    call sfree (sp)
+		    call error (2,
+			"Can't combine references with non-linear dispersions")
+		}
+		coeffs = Memi[pcoeff2+k2]
+		if (coeffs == NULL) {
+		    if (DC_DT(ap,i) == 2) {
+			dval = (wt1*Memd[pshift1+k1] + wt2*Memd[pshift2+k2]) /
+			    (wt1 + wt2)
+			call sshift1 (dval, DC_CF(ap,i))
+		    } else {
+			call fprintf (sfd, " %.8g %g")
+			    call pargd (wt2)
+			    call pargd (Memd[pshift2+k2])
+			dval = DC_DW(ap,i) * (DC_NW(ap,i) - 1) / 2.
+			call fprintf (sfd, " 1 2 1 %d %g %g")
+			    call pargi (DC_NW(ap,i))
+			    call pargd (DC_W1(ap,i) + dval)
+			    call pargd (dval)
+		    }
+		} else {
+		    call fprintf (sfd, " %.8g %g %d %d")
+			call pargd (wt2)
+			call pargd (Memd[pshift2+k2])
+			call pargi (nint (Memd[coeffs+1]))
+			call pargi (nint (Memd[coeffs+2]))
+		    ncoeffs = nint (Memd[coeffs])
+		    do k = 3, ncoeffs {
+			call fprintf (sfd, " %.15g")
+			    call pargd (Memd[coeffs+k])
+		    }
+		}
+	    }
+
+	    if (i == 1) {
+		un = un1
+		if (UN_LABEL(un) != EOS)
+		    call mw_swattrs (SMW_MW(smw,0), 1, "label", UN_LABEL(un))
+		if (UN_UNITS(un) != EOS)
+		    call mw_swattrs (SMW_MW(smw,0), 1, "units", UN_UNITS(un))
+		call un_close (DC_UN(ap,i))
+		DC_UN(ap,i) = un
+	    } else if (!un_compare (un, un1)) {
+		call sfree (sp)
+		call error (3, "Units must be the same for all apertures")
+	    }
+	    DC_DT(ap,i) = 2
+	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), Memc[DC_CF(ap,i)])
+	    if (sfd != NULL)
+		call strclose (sfd)
+	}
+
+	# Update the linear part of WCS.
+	ct1 = smw_sctran (smw, "logical", "physical", 2)
+	ct2 = smw_sctran (smw, "physical", "world", 3)
+	do i = 1, naps {
+	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), DC_CF(ap,i))
+	    wt1 = nint (smw_c1trand (ct1, double(i)))
+	    call smw_c2trand (ct2, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
+	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
+	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), Memc[DC_CF(ap,i)])
+	}
+	call smw_ctfree (ct1)
+	call smw_ctfree (ct2)
+
+done_	# Set aperture parameters in terms of logical image.
+	ct1 = smw_sctran (smw, "physical", "logical", 1)
+	j = nint (smw_c1trand (ct1, 1D0))
+	do i = 1, naps {
+	    k = nint (smw_c1trand (ct1, double(DC_NW(ap,i))))
+	    DC_NW(ap,i) = min (IM_LEN(im,1), max (j, k))
+	}
+	call smw_ctfree (ct1)
+
+	ct1 = smw_sctran (smw, "logical", "world", 3)
+	do i = 1, naps {
+	    wt1 = i
+	    call smw_c2trand (ct1, 1D0, wt1, DC_W1(ap,i), wt2)
+	    call smw_c2trand (ct1, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
+	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
+	}
+	call smw_ctfree (ct1)
+
+	do i = 1, naps
+	    call mfree (DC_CF(ap,i), TY_CHAR)
+	call sfree (sp)
+end
+
+
+# DC_GMSDB -- Get a dispersion database entry.
+# The database entry is read only once from the database and stored in a
+# symbol table keyed by the spectrum name.  Subsequent requests for the
+# reference spectrum returns the data from the symbol table.
+
+procedure dc_gmsdb (spec, stp, paps, papcen, punits, pshift, pcoeff, naps)
+
+char	spec[ARB]	# Spectrum image name
+pointer	stp		# Symbol table pointer
+pointer	paps		# Pointer to aperture numbers
+pointer	papcen		# Pointer to aperture centers
+pointer	punits		# Pointer to units
+pointer	pshift		# Pointer to shifts
+pointer	pcoeff		# Pointer to coefficients
+int	naps		# Number of apertures
+
+double	dval
+int	i, n, dtgeti(), getline(), ctod()
+real	low, high, dtgetr()
+pointer	sp, str, coeffs, sym, db, dt, dt1
+pointer	stfind(), stenter(), strefsbuf(), dtmap1(), un_open()
+errchk	dtmap1, dtgeti, dtgad, un_open
+
+begin
+	# Check if dispersion solution is in the symbol table from a previous
+	# call.  If not in the symbol table get it from the database and
+	# store it in the symbol table.
+
+	sym = stfind (stp, spec)
+	if (sym == NULL) {
+	    call smark (sp)
+	    call salloc (str, SZ_LINE, TY_CHAR)
+	    call strcpy ("id", Memc[str], SZ_LINE)
+	    call imgcluster (spec, Memc[str+2], SZ_LINE-2)
+	    call xt_imroot (Memc[str+2], Memc[str+2], SZ_LINE-2)
+	    db = strefsbuf (stp, Memi[stfind (stp, "database")])
+	    dt = dtmap1 (Memc[db], Memc[str], READ_ONLY)
+	    call strcpy ("ec", Memc[str], SZ_LINE)
+	    call imgcluster (spec, Memc[str+2], SZ_LINE-2)
+	    call xt_imroot (Memc[str+2], Memc[str+2], SZ_LINE-2)
+	    ifnoerr (dt1 = dtmap1 (Memc[db], Memc[str], READ_ONLY)) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "Ambiguous database files: %s/%s and %s/%s")
+		    call pargstr (DT_DNAME(dt))
+		    call pargstr (DT_FNAME(dt))
+		    call pargstr (DT_DNAME(dt1))
+		    call pargstr (DT_FNAME(dt1))
+		call dtunmap (dt)
+		call dtunmap (dt1)
+		call fatal (3, Memc[str])
+	    }
+
+	    naps = max (1, DT_NRECS(dt))
+	    call calloc (paps, naps, TY_INT)
+	    call calloc (papcen, naps, TY_REAL)
+	    call calloc (punits, naps, TY_POINTER)
+	    call calloc (pshift, naps, TY_DOUBLE)
+	    call calloc (pcoeff, naps, TY_POINTER)
+	    if (DT_NRECS(dt) > 0) {
+		for (i = 1; i <= naps; i = i + 1) {
+		    iferr (Memi[paps+i-1] = dtgeti (dt, i, "aperture"))
+			Memi[paps+i-1] = INDEFI
+		    iferr (low = dtgetr (dt, i, "aplow"))
+			low = INDEF
+		    iferr (high = dtgetr (dt, i, "aphigh"))
+			high = INDEF
+		    if (IS_INDEF(low) || IS_INDEF(high))
+			Memr[papcen+i-1] = 0.
+		    else
+			Memr[papcen+i-1] = (low + high) / 2.
+		    iferr (call dtgstr (dt, i, "units", Memc[str], SZ_LINE))
+			call strcpy ("Angstroms", Memc[str], SZ_LINE)
+		    Memi[punits+i-1] = un_open (Memc[str])
+		    iferr (Memd[pshift+i-1] = dtgetr (dt, i, "shift"))
+			Memd[pshift+i-1] = 0.
+		    iferr {
+			n = dtgeti (dt, i, "coefficients")
+			call malloc (coeffs, 1+n, TY_DOUBLE)
+			Memd[coeffs] = n
+			call dtgad (dt, i, "coefficients", Memd[coeffs+1], n, n)
+			Memi[pcoeff+i-1] = coeffs
+		    } then
+			Memi[pcoeff+i-1] = NULL
+		}
+	    } else {
+		Memi[paps] = INDEFI
+		Memr[papcen] = INDEFR
+		Memi[punits] = un_open ("")
+		Memd[pshift] = 0.
+		call malloc (coeffs, 100, TY_DOUBLE)
+		n = 3
+		call seek (Memi[dt], BOF)
+		while (getline (Memi[dt], Memc[str]) != EOF) {
+		    i = 1
+		    if (ctod (Memc[str], i, dval) == 0)
+			next
+		    if (mod (n, 100) == 0)
+			call realloc (coeffs, n+100, TY_DOUBLE)
+		    Memd[coeffs+n] = dval
+		    n = n + 1
+		}
+		Memd[coeffs] = n - 1
+		Memd[coeffs+1] = 5
+		Memd[coeffs+2] = n - 3
+		Memi[pcoeff] = coeffs
+	    }
+
+	    call dtunmap (dt)
+	    call sfree (sp)
+
+	    sym = stenter (stp, spec, LEN_DC)
+	    DC_FORMAT(sym) = 1
+	    DC_PAPS(sym) = paps
+	    DC_PAPCEN(sym) = papcen
+	    DC_PUN(sym) = punits
+	    DC_PSHIFT(sym) = pshift
+	    DC_PCOEFF(sym) = pcoeff
+	    DC_NAPS(sym) = naps
+	} else {
+	    if (DC_FORMAT(sym) != 1)
+		call error (1, "Not a multispec dispersion function")
+	    paps = DC_PAPS(sym)
+	    papcen = DC_PAPCEN(sym)
+	    punits = DC_PUN(sym)
+	    pshift = DC_PSHIFT(sym)
+	    pcoeff = DC_PCOEFF(sym)
+	    naps = DC_NAPS(sym)
+	}
+end
+
+
+# DC_REFSHFT -- Compute dispersion shift.
+
+procedure dc_refshft (spec, stp, refspec, keywrd, im, aps, apcens, shifts,
+	naps, fd1, fd2)
+
+char	spec[ARB]		# Spectrum to be corrected
+pointer	stp			# Symbol table pointer
+char	refspec[ARB]		# Reference spectrum
+char	keywrd[ARB]		# Header keyword (for log only)
+pointer	im			# IMIO pointer to spectrum to be corrected
+int	aps[naps]		# Reference apertures
+real	apcens[naps]		# Reference aperture centers
+double	shifts[naps]		# Reference aperture shifts (to be modified)
+int	naps			# Number of refernce apertures
+int	fd1			# Logfile descriptor
+int	fd2			# Logfile descriptor
+
+int	i, j, k, pnaps
+double	apcen, shift, sumx, sumy, sumxx, sumyy, sumxy, a, b
+pointer	sp, refshft, option, paps, papcen, punits, pshift, pcoeff
+bool	streq()
+errchk	imgstr, dc_gmsdb
+
+begin
+	call smark (sp)
+	call salloc (refshft, SZ_FNAME, TY_CHAR)
+	call salloc (option, SZ_FNAME, TY_CHAR)
+
+	# Parse header parameter.
+	call imgstr (im, keywrd, Memc[refshft], SZ_FNAME)
+	call sscan (Memc[refshft])
+	call gargwrd (Memc[refshft], SZ_FNAME)
+	if (streq (Memc[refshft], refspec)) {
+	    call sfree (sp)
+	    return
+	}
+	call gargwrd (Memc[option], SZ_FNAME)
+
+	# Get reference shift apertures.
+	call dc_gmsdb (Memc[refshft], stp, paps, papcen, punits, pshift,
+	    pcoeff, pnaps)
+	if (pnaps == 0) {
+	    call sfree (sp)
+	    return
+	}
+	
+	# Compute mean shift and RMS.
+	sumy = 0.
+	sumyy = 0.
+	do i = 1, pnaps {
+	    sumy = sumy + Memd[pshift+i-1]
+	    sumyy = sumyy + Memd[pshift+i-1] ** 2
+	}
+	sumy = sumy / pnaps
+	sumyy = sqrt (max (0.D0, sumyy / pnaps - sumy ** 2))
+
+	# Print.
+	if (fd1 != NULL) {
+	    call fprintf (fd1, "%s: %s = '%s %s', shift = %.6g, rms = %.6g\n")
+		call pargstr (spec)
+		call pargstr (keywrd)
+		call pargstr (Memc[refshft])
+		call pargstr (Memc[option])
+		call pargd (sumy)
+		call pargd (sumyy)
+	}
+	if (fd2 != NULL) {
+	    call fprintf (fd2, "%s: %s = '%s %s', shift = %.6g, rms = %.6g\n")
+	        call pargstr (spec)
+		call pargstr (keywrd)
+	        call pargstr (Memc[refshft])
+	        call pargstr (Memc[option])
+	        call pargd (sumy)
+	        call pargd (sumyy)
+	}
+
+	if (streq (Memc[option], "interp")) {
+	    if (pnaps > 1) {
+	        sumx = 0.
+	        sumy = 0.
+	        sumxx = 0.
+	        sumyy = 0.
+	        sumxy = 0.
+	        do i = 0, pnaps-1 {
+		    apcen = Memr[papcen+i]
+		    shift = Memd[pshift+i]
+		    sumx = sumx + apcen
+		    sumy = sumy + shift
+		    sumxx = sumxx + apcen * apcen
+		    sumyy = sumyy + shift * shift
+		    sumxy = sumxy + apcen * shift
+	        }
+	        b = pnaps * sumxx - sumx * sumx
+	        a = (sumy * sumxx - sumx * sumxy) / b
+	        b = (pnaps * sumxy - sumx * sumy) / b
+	    } else {
+	        a = sumy
+	        b = 0.
+	    }
+	    do i = 1, naps
+	        shifts[i] = shifts[i] + a + b * apcens[i]
+	    if (fd1 != NULL) {
+	        call fprintf (fd1, "\tintercept = %.6g, slope = %.6g\n")
+		    call pargd (a)
+		    call pargd (b)
+	    }
+	    if (fd2 != NULL) {
+	        call fprintf (fd2, "\tintercept = %.6g, slope = %.6g\n")
+		    call pargd (a)
+		    call pargd (b)
+	    }
+	} else if (streq (Memc[option], "nearest")) {
+	    do i = 1, naps {
+	        k = 0
+	        sumy = abs (apcens[i] - Memr[papcen])
+	        for (j = 1; j < pnaps; j = j + 1)
+		    if (abs (apcens[i] - Memr[papcen+j]) < sumy) {
+		        k = j
+		        sumy = abs (apcens[i] - Memr[papcen+k])
+		    }
+	        shifts[i] = shifts[i] + Memd[pshift+k]
+	        if (fd1 != NULL) {
+		    call fprintf (fd1, "\t%4d %7.2f %4d %7.2f %.6g\n")
+		        call pargi (aps[i])
+		        call pargr (apcens[i])
+		        call pargi (Memi[paps+k])
+		        call pargr (Memr[papcen+k])
+		        call pargd (Memd[pshift+k])
+	        }
+	        if (fd2 != NULL) {
+		    call fprintf (fd2, "\t%4d %7.2f %4d %7.2f %.6g\n")
+		        call pargi (aps[i])
+		        call pargr (apcens[i])
+		        call pargi (Memi[paps+k])
+		        call pargr (Memr[papcen+k])
+		        call pargd (Memd[pshift+k])
+	        }
+	    }
+	} else
+	    call aaddkd (shifts, sumy, shifts, naps)
+
+	call sfree (sp)
+end
+
+
+# DC_GEC -- Get an echelle spectrum.  This consists of mapping the image
+# and setting a MWCS coordinate transformation.  If not dispersion corrected
+# the dispersion function is found in the database for the reference
+# spectra and set in the SMW.
+
+procedure dc_gec (spec, im, smw, stp, ap, fd1, fd2)
+
+char	spec[ARB]	#I Spectrum name
+pointer	im		#I IMIO pointer
+pointer	smw		#I SMW pointers
+pointer	stp		#I Symbol table
+pointer	ap		#O Aperture data structure
+int	fd1		#I Logfile descriptor
+int	fd2		#I Logfile descriptor
+
+double	wt1, wt2, dval
+int	i, j, k, l, dc, sfd, naps, ncoeffs, offset, slope
+pointer	sp, str1, str2, coeff, coeffs, ct1, ct2, un1, un2, un3
+pointer	pshift1, pshift2, pshift3, pcoeff1, pcoeff2, pcoeff3
+
+bool	un_compare()
+double	smw_c1trand()
+int	imaccf(), nscan(), stropen()
+pointer	smw_sctran(), un_open()
+errchk	dc_gecdb, imgstr, smw_sctran, un_open
+
+define	done_	90
+
+begin
+	call smark (sp)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (str2, SZ_LINE, TY_CHAR)
+	coeff = NULL
+
+	# Set WCS attributes
+	naps = IM_LEN(im,2)
+	call calloc (ap, LEN_AP(naps), TY_STRUCT)
+	do i = 1, naps {
+	    DC_PL(ap,i) = i
+	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), coeff)
+	    if (i == 1) {
+		iferr (call mw_gwattrs (SMW_MW(smw,0), 1, "units", Memc[str1],
+		    SZ_LINE))
+		    Memc[str1] = EOS
+		DC_UN(ap,i) = un_open (Memc[str1])
+	    }
+	    dc = DC_DT(ap,i)
+	}
+
+	# Check if the spectra have been dispersion corrected
+	# by an earlier version of DISPCOR.  If so then don't allow
+	# another database dispersion correction.  This assumes all
+	# spectra have the same dispersion type.  Check for a
+	# reference spectrum.
+
+	#if ((imaccf (im, "REFSPEC1") == NO) ||
+	#    (dc > -1 && imaccf (im, "DCLOG1") == NO)) {
+	if (imaccf (im, "REFSPEC1") == NO) {
+	    if (fd1 != NULL) {
+		call fprintf (fd1,
+		    "%s: Resampling using current coordinate system\n")
+		    call pargstr (spec)
+	    }
+	    if (fd2 != NULL) {
+		call fprintf (fd2,
+		    "%s: Resampling using current coordinate system\n")
+		    call pargstr (spec)
+	    }
+	    goto done_
+	}
+
+	# Get the reference spectra dispersion function from the database
+	# and determine a reference shift.
+
+	iferr {
+	    call imgstr (im, "REFSPEC1", Memc[str1], SZ_LINE)
+	    call sscan (Memc[str1])
+	    call gargwrd (Memc[str1], SZ_LINE)
+	    call gargd (wt1)
+	    if (nscan() == 1)
+		wt1 = 1.
+	} then {
+	    call strcpy (spec, Memc[str1], SZ_LINE)
+	    wt1 = 1.
+	}
+	call salloc (pshift1, naps, TY_DOUBLE)
+	call salloc (pcoeff1, naps, TY_POINTER)
+	slope = 0
+	iferr (call dc_gecdb (Memc[str1], stp, ap, un1, Memd[pshift1],
+	    Memi[pcoeff1], naps, offset, slope)) {
+	    call sfree (sp)
+	    call erract (EA_ERROR)
+	}
+	if (fd1 != NULL) {
+	    call fprintf (fd1, "%s: REFSPEC1 = '%s %.8g'\n")
+		call pargstr (spec)
+		call pargstr (Memc[str1])
+		call pargd (wt1)
+	}
+	if (fd2 != NULL) {
+	    call fprintf (fd2, "%s: REFSPEC1 = '%s %.8g'\n")
+		call pargstr (spec)
+		call pargstr (Memc[str1])
+		call pargd (wt1)
+	}
+
+        iferr {
+            call imgstr (im, "refshft1", Memc[str1], SZ_LINE)
+	    call salloc (pshift3, naps, TY_DOUBLE)
+	    call salloc (pcoeff3, naps, TY_POINTER)
+	    call dc_gecdb (Memc[str1], stp, ap, un3, Memd[pshift3],
+		Memi[pcoeff3], naps, offset, slope)
+            if (fd1 != NULL) {
+                call fprintf (fd1, "%s: REFSHFT1 = '%s', shift = %.6g\n")
+                    call pargstr (spec)
+                    call pargstr (Memc[str1])
+                    call pargd (Memd[pshift3])
+            }
+            if (fd2 != NULL) {
+                call fprintf (fd2, "%s: REFSHFT1 = '%s', shift = %.6g\n")
+                    call pargstr (spec)
+                    call pargstr (Memc[str1])
+                    call pargd (Memd[pshift3])
+            }
+	    call aaddd (Memd[pshift1], Memd[pshift3], Memd[pshift1], naps)
+        } then
+            ;
+
+	iferr {
+            call imgstr (im, "REFSPEC2", Memc[str1], SZ_LINE)
+            call sscan (Memc[str1])
+            call gargwrd (Memc[str1], SZ_LINE)
+            call gargd (wt2)
+            if (nscan() == 1)
+                wt2 = 1.
+	    call salloc (pshift2, naps, TY_DOUBLE)
+	    call salloc (pcoeff2, naps, TY_POINTER)
+	    call dc_gecdb (Memc[str1], stp, ap, un2, Memd[pshift2],
+		Memi[pcoeff2], naps, offset, slope)
+	    if (fd1 != NULL) {
+		call fprintf (fd1, "%s: REFSPEC2 = '%s %.8g'\n")
+		    call pargstr (spec)
+		    call pargstr (Memc[str1])
+		    call pargd (wt2)
+	    }
+	    if (fd2 != NULL) {
+		call fprintf (fd2, "%s: REFSPEC2 = '%s %.8g'\n")
+		    call pargstr (spec)
+		    call pargstr (Memc[str1])
+		    call pargd (wt2)
+	    }
+
+	    iferr {
+		call imgstr (im, "refshft2", Memc[str1], SZ_LINE)
+		call salloc (pshift3, naps, TY_DOUBLE)
+		call salloc (pcoeff3, naps, TY_POINTER)
+		call dc_gecdb (Memc[str1], stp, ap, un3, Memd[pshift3],
+		    Memi[pcoeff3], naps, offset, slope)
+		if (fd1 != NULL) {
+		    call fprintf (fd1, "%s: REFSHFT2 = '%s', shift = %.6g\n")
+			call pargstr (spec)
+                    call pargstr (Memc[str1])
+                    call pargd (Memd[pshift3])
+		}
+		if (fd2 != NULL) {
+		    call fprintf (fd2, "%s: REFSHFT2 = '%s', shift = %.6g\n")
+			call pargstr (spec)
+                    call pargstr (Memc[str1])
+                    call pargd (Memd[pshift3])
+		}
+		call aaddd (Memd[pshift1], Memd[pshift3], Memd[pshift1], naps)
+	    } then
+		;
+	} then
+	    wt2 = 0.
+
+	# Adjust weights to unit sum.
+	dval = wt1 + wt2
+	wt1 = wt1 / dval
+	wt2 = wt2 / dval
+
+	# Enter dispersion function in the MWCS.
+	do i = 1, naps {
+	    coeffs = Memi[pcoeff1+i-1]
+	    ncoeffs = nint (Memd[coeffs])
+	    l = 20 * (ncoeffs + 2)
+	    if (wt2 > 0.)
+		l = 2 * l
+	    call realloc (coeff, l, TY_CHAR)
+	    call aclrc (Memc[coeff], l)
+	    sfd = stropen (Memc[coeff], l, NEW_FILE)
+	    call fprintf (sfd, "%.8g %g")
+		call pargd (wt1)
+		call pargd (Memd[pshift1+i-1])
+
+	    # The following assumes some knowledge of the data structure in
+	    # order to shortten the the attribute string.
+
+	    call fprintf (sfd, " %d %d %.8g %.8g")
+		call pargi (nint (Memd[coeffs+1]))
+		call pargi (nint (Memd[coeffs+2]))
+		call pargd (Memd[coeffs+3])
+		call pargd (Memd[coeffs+4])
+	    do j = 5, ncoeffs {
+		call fprintf (sfd, " %.15g")
+		    call pargd (Memd[coeffs+j])
+	    }
+
+	    if (wt2 > 0.) {
+		coeffs = Memi[pcoeff2+i-1]
+		ncoeffs = nint (Memd[coeffs])
+		call fprintf (sfd, "%.8g %g")
+		    call pargd (wt2)
+		    call pargd (Memd[pshift2+i-1])
+		call fprintf (sfd, " %d %d %.8g %.8g")
+		    call pargi (nint (Memd[coeffs+1]))
+		    call pargi (nint (Memd[coeffs+2]))
+		    call pargd (Memd[coeffs+3])
+		    call pargd (Memd[coeffs+4])
+		do j = 5, ncoeffs {
+		    call fprintf (sfd, " %.15g")
+			call pargd (Memd[coeffs+j])
+		}
+		if (!un_compare (un1, un2)) {
+		    call sfree (sp)
+		    call error (2,
+			"Can't combine references with different units")
+		}
+	    }
+
+	    if (i == 1) {
+		if (UN_LABEL(un1) != EOS)
+		    call mw_swattrs (SMW_MW(smw,0), 1, "label", UN_LABEL(un1))
+		if (UN_UNITS(un1) != EOS)
+		    call mw_swattrs (SMW_MW(smw,0), 1, "units", UN_UNITS(un1))
+		call un_close (DC_UN(ap,i))
+		DC_UN(ap,i) = un1
+	    }
+	    DC_DT(ap,i) = 2
+	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), Memc[coeff])
+	    call strclose (sfd)
+	}
+
+	# Update the linear part of WCS.
+	ct1 = smw_sctran (smw, "logical", "physical", 2)
+	ct2 = smw_sctran (smw, "physical", "world", 3)
+	do i = 1, naps {
+	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), coeff)
+	    wt1 = nint (smw_c1trand (ct1, double(i)))
+	    call smw_c2trand (ct2, 1D0, wt1, DC_W1(ap,i), wt2)
+	    call smw_c2trand (ct2, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
+	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
+	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
+		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
+		DC_LW(ap,i), DC_UP(ap,i), Memc[coeff])
+	}
+	call smw_ctfree (ct1)
+	call smw_ctfree (ct2)
+
+done_	# Set aperture parameters in terms of logical image.
+	ct1 = smw_sctran (smw, "physical", "logical", 1)
+	j = nint (smw_c1trand (ct1, 1D0))
+	do i = 1, naps {
+	    k = nint (smw_c1trand (ct1, double(DC_NW(ap,i))))
+	    DC_NW(ap,i) = min (IM_LEN(im,1), max (j, k))
+	}
+	call smw_ctfree (ct1)
+
+	ct1 = smw_sctran (smw, "logical", "world", 3)
+	do i = 1, naps {
+	    wt1 = i
+	    call smw_c2trand (ct1, 1D0, wt1, DC_W1(ap,i), wt2)
+	    call smw_c2trand (ct1, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
+	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
+	}
+	call smw_ctfree (ct1)
+
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end
+
+
+# DC_GECDB -- Get a dispersion database entry.
+# The database entry is read only once from the database and stored in a
+# symbol table keyed by the spectrum name.  Subsequent requests for the
+# reference spectrum returns the data from the symbol table.
+
+procedure dc_gecdb (spec, stp, ap, un, shifts, pcoeff, naps, offset, slope)
+
+char	spec[ARB]	# Spectrum image name
+pointer	stp		# Symbol table pointer
+pointer	ap		# Aperture data structure
+pointer	un		# Units
+double	shifts[naps]	# Shifts
+pointer	pcoeff[naps]	# Pointer to coefficients
+int	naps		# Number of apertures
+int     offset          # Aperture to order offset
+int     slope           # Aperture to order slope
+
+double  shift
+real    dtgetr()
+int     i, rec, offst, slpe, n, dtlocate(), dtgeti()
+pointer sp, str, coeffs, sym, db, dt
+pointer	stfind(), stenter(), strefsbuf(), dtmap1(), un_open()
+errchk  dtmap1, dtlocate, dtgeti, dtgad, un_open
+
+begin
+	# Check if dispersion solution is in the symbol table from a previous
+	# call.  If not in the symbol table get it from the database and
+	# store it in the symbol table.
+
+	sym = stfind (stp, spec)
+	if (sym == NULL) {
+	    call smark (sp)
+	    call salloc (str, SZ_LINE, TY_CHAR)
+	    call strcpy ("ec", Memc[str], SZ_LINE)
+	    call imgcluster (spec, Memc[str+2], SZ_LINE-2)
+	    call xt_imroot (Memc[str+2], Memc[str+2], SZ_LINE-2)
+	    db = strefsbuf (stp, Memi[stfind (stp, "database")])
+	    dt = dtmap1 (Memc[db], Memc[str], READ_ONLY)
+
+            call sprintf (Memc[str], SZ_LINE, "ecidentify %s")
+                call pargstr (spec)
+            iferr (rec = dtlocate (dt, Memc[str])) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "DISPCOR: Echelle dispersion function not found (%s/%s)")
+		    call pargstr (DT_DNAME(dt))
+		    call pargstr (DT_FNAME(dt))
+                call fatal (0, Memc[str])
+	    }
+
+	    iferr (call dtgstr (dt, rec, "units", Memc[str], SZ_LINE))
+		call strcpy ("Angstroms", Memc[str], SZ_LINE)
+	    un = un_open (Memc[str])
+            iferr (offst = dtgeti (dt, rec, "offset"))
+                offst = 0
+            iferr (slpe = dtgeti (dt, rec, "slope"))
+                slpe = 1
+            iferr (shift = dtgetr (dt, rec, "shift"))
+                shift = 0.
+            n = dtgeti (dt, rec, "coefficients")
+            call malloc (coeffs, n, TY_DOUBLE)
+            call dtgad (dt, rec, "coefficients", Memd[coeffs], n, n)
+
+            sym = stenter (stp, spec, LEN_DC)
+	    DC_FORMAT(sym) = 2
+	    DC_PUN(sym) = un
+            DC_OFFSET(sym) = offst
+            DC_SLOPE(sym) = slpe
+            DC_SHIFT(sym) = shift
+            DC_COEFFS(sym) = coeffs
+
+	    call dtunmap (dt)
+	    call sfree (sp)
+	} else {
+	    if (DC_FORMAT(sym) != 2)
+		call error (1, "Not an echelle dispersion function")
+	    un = DC_PUN(sym)
+            offst = DC_OFFSET(sym)
+            slpe = DC_SLOPE(sym)
+            coeffs = DC_COEFFS(sym)
+            shift = DC_SHIFT(sym)
+	}
+
+	# Check aperture to order parameters.
+	if (slope == 0) {
+	    offset = offst
+	    slope = slpe
+	} else if (offset != offst || slope != slpe) {
+	    call eprintf (
+		"WARNING: Echelle order offsets/slopes are not the same.\n")
+	}
+
+	# Convert to multispec coefficients
+	do i = 1, naps {
+	    DC_BM(ap,i) = offset + slope * DC_AP(ap,i)
+	    call dc_ecms (DC_BM(ap,i), Memd[coeffs], pcoeff[i])
+	    shifts[i] = shift / DC_BM(ap,i)
+	}
+end
+
+
+# DC_ECMS -- Convert echelle dispersion coefficients to multispec coefficients
+
+procedure dc_ecms (order, eccoeff, mscoeff)
+
+int	order			# Echelle order
+double	eccoeff[ARB]		# Echelle dispersion coefficients
+pointer	mscoeff			# Pointer to multispec coefficients
+
+int	i, j, k, type, xorder, yorder
+double	xmin, xmax, ymin, ymax, ymaxmin, yrange, y, coeff, a, b, c
+
+begin
+	type = nint (eccoeff[1])
+	xorder = nint (eccoeff[2])
+	yorder = nint (eccoeff[3])
+	xmin = eccoeff[5]
+	xmax = eccoeff[6]
+	ymin = eccoeff[7]
+	ymax = eccoeff[8]
+
+	yrange = 2. / (ymax - ymin)
+	ymaxmin = (ymax + ymin) / 2
+	y = (order - ymaxmin) * yrange
+
+	call malloc (mscoeff, 5+xorder, TY_DOUBLE)
+	Memd[mscoeff] = 4+xorder
+	Memd[mscoeff+1] = type
+	Memd[mscoeff+2] = xorder
+	Memd[mscoeff+3] = xmin
+	Memd[mscoeff+4] = xmax
+
+	switch (type) {
+	case 1:
+	    do k = 1, xorder {
+		j = 9 + k - 1
+		coeff = eccoeff[j]
+		if (yorder > 1) {
+		    j = j + xorder
+		    coeff = coeff + eccoeff[j] * y
+		}
+		if (yorder > 2) {
+		    a = 1
+		    b = y
+		    do i = 3, yorder {
+			c = 2 * y * b - a
+			j = j + xorder
+			coeff = coeff + eccoeff[j] * c
+			a = b
+			b = c
+		    }
+		}
+		Memd[mscoeff+4+k] = coeff / order
+	    }
+	case 2:
+	    do k = 1, xorder {
+		j = 9 + k - 1
+		coeff = eccoeff[j]
+		if (yorder > 1) {
+		    j = j + xorder
+		    coeff = coeff + eccoeff[j] * y
+		}
+		if (yorder > 2) {
+		    a = 1
+		    b = y
+		    do i = 3, yorder {
+			c = ((2 * i - 3) * y * b - (i - 2) * a) / (i - 1)
+			j = j + xorder
+			coeff = coeff + eccoeff[j] * c
+			a = b
+			b = c
+		    }
+		}
+		Memd[mscoeff+4+k] = coeff / order
+	    }
+	}
+end
diff --git a/noao/onedspec/dispcor/dctable.h b/noao/onedspec/dispcor/dctable.h
new file mode 100644
index 00000000..4cf3657a
--- /dev/null
+++ b/noao/onedspec/dispcor/dctable.h
@@ -0,0 +1,11 @@
+# Wavelength table structure
+define	TBL_LEN		14
+define	TBL_W1		Memd[P2D($1)]	# Starting wavelength
+define	TBL_W2		Memd[P2D($1+2)]	# Ending wavelength
+define	TBL_DW		Memd[P2D($1+4)]	# Wavelength interval
+define	TBL_WMIN	Memd[P2D($1+6)]	# Minimum wavelength for global
+define	TBL_WMAX	Memd[P2D($1+8)]	# Maximum wavelength for global
+define	TBL_AP		Memi[$1+10]	# Aperture
+define	TBL_NW		Memi[$1+11]	# Number of points
+define	TBL_NWMAX	Memi[$1+12]	# Maximum number of points for global
+define	TBL_CONFIRM	Memi[$1+13]	# Confirm?
diff --git a/noao/onedspec/dispcor/dctable.x b/noao/onedspec/dispcor/dctable.x
new file mode 100644
index 00000000..93f27531
--- /dev/null
+++ b/noao/onedspec/dispcor/dctable.x
@@ -0,0 +1,145 @@
+include	<imhdr.h>
+include	<mach.h>
+include	"dctable.h"
+include	<smw.h>
+
+
+# DC_TABLE -- Set default wavelengths.
+# This may be specified by the task parameters alone, from a reference image,
+# or from a text table.  A reference image or table allows separate
+# wavelength parameters for each aperture.  The text table columns are the
+# aperture number, starting wavelength, ending wavelength, wavelength
+# interval per pixel, and number of pixels.  Any of these values may be
+# INDEF.
+
+procedure dc_table (table, naps)
+
+pointer	table			# Table pointer (returned)
+int	naps			# Number of apertures (returned)
+
+int	i, j, ap, nw, fd, clgeti(), open(), fscan(), nscan(), btoi(), nowhite()
+double	ws, we, dw, clgetd()
+pointer	sp, fname, tbl, mw, sh, immap(), smw_openim()
+bool	clgetb()
+errchk	smw_openim(), shdr_open()
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+	call clgstr ("table", Memc[fname], SZ_FNAME)
+
+	# Set defaults.
+	naps = 0
+	call malloc (table, 10, TY_INT)
+	call malloc (Memi[table], TBL_LEN, TY_STRUCT)
+	tbl= Memi[table]
+	TBL_W1(tbl) = clgetd ("w1")
+	TBL_W2(tbl) = clgetd ("w2")
+	TBL_DW(tbl) = clgetd ("dw")
+	TBL_NW(tbl) = clgeti ("nw")
+	TBL_WMIN(tbl) = MAX_REAL
+	TBL_WMAX(tbl) = -MAX_REAL
+	TBL_NWMAX(tbl) = 0
+	TBL_CONFIRM(tbl) = btoi (clgetb ("confirm"))
+
+	# Read a reference image or table if specified and add entries to
+	# the table array.
+
+	if (nowhite (Memc[fname], Memc[fname], SZ_FNAME) > 0) {
+	    ifnoerr (fd = immap (Memc[fname], READ_ONLY, 0)) {
+		mw = smw_openim (fd)
+		call shdr_open (fd, mw, 1, 1, INDEFI, SHHDR, sh)
+		if (DC(sh) == DCLINEAR || DC(sh) == DCLOG) {
+		    do j = 1, IM_LEN(fd,2) {
+			call shdr_open (fd, mw, j, 1, INDEFI, SHHDR, sh)
+			call dc_getentry (false, AP(sh), table, naps, i)
+			tbl = Memi[table+i]
+			TBL_AP(tbl) = AP(sh)
+			TBL_NW(tbl) = SN(sh)
+			TBL_W1(tbl) = W0(sh)
+			TBL_W2(tbl) = W1(sh)
+			TBL_DW(tbl) = WP(sh)
+		    }
+		}
+		call shdr_close (sh)
+		call smw_close (mw)
+		call imunmap (fd)
+	    } else {
+		ifnoerr (fd = open (Memc[fname], READ_ONLY, TEXT_FILE)) {
+		    while (fscan (fd) != EOF) {
+			call gargi (ap)
+			call gargd (ws)
+			call gargd (we)
+			call gargd (dw)
+			call gargi (nw)
+			if (nscan() < 5)
+			    next
+
+			call dc_getentry (false, ap, table, naps, i)
+			tbl = Memi[table+i]
+			TBL_AP(tbl) = ap
+			TBL_W1(tbl) = ws
+			TBL_W2(tbl) = we
+			TBL_DW(tbl) = dw
+			TBL_NW(tbl) = nw
+		    }
+		    call close (fd)
+		} else
+		    call error (1, "Can't access wavelength table")
+	    }
+	}
+
+	# If ignoreaps=yes then replace INDEFs in the default entry with
+	# the first non-INDEF entry.
+
+	if (clgetb ("ignoreaps") && naps > 0) {
+	    tbl= Memi[table]
+	    if (IS_INDEFD(TBL_W1(tbl)))
+		TBL_W1(tbl) = TBL_W1(Memi[table+1])
+	    if (IS_INDEFD(TBL_W2(tbl)))
+		TBL_W2(tbl) = TBL_W2(Memi[table+1])
+	    if (IS_INDEFD(TBL_DW(tbl)))
+		TBL_DW(tbl) = TBL_DW(Memi[table+1])
+	    if (IS_INDEFI(TBL_NW(tbl)))
+		TBL_NW(tbl) = TBL_NW(Memi[table+1])
+	}
+
+	call sfree (sp)
+end
+
+
+# DC_GETENTRY -- Get entry from wavelength table.  Return the index.  Allocate
+# a new entry if needed.
+
+procedure dc_getentry (apflag, ap, table, naps, index)
+
+bool	apflag		# Ignore aperture numbers?
+int	ap		# Aperture
+pointer	table		# Wavelength table
+int	naps		# Number of apertures
+int	index		# Table index of entry
+
+pointer	tbl
+
+begin
+	for (index=1; index<=naps; index=index+1)
+	    if (apflag || TBL_AP(Memi[table+index]) == ap)
+		return
+
+	naps = naps + 1
+	if (mod (naps, 10) == 0)
+	    call realloc (table, naps+10, TY_INT)
+	call malloc (Memi[table+naps], TBL_LEN, TY_STRUCT)
+
+	index = naps
+	tbl = Memi[table+index]
+	TBL_AP(tbl) = ap
+	TBL_W1(tbl) = TBL_W1(Memi[table])
+	TBL_W2(tbl) = TBL_W2(Memi[table])
+	TBL_DW(tbl) = TBL_DW(Memi[table])
+	TBL_NW(tbl) = TBL_NW(Memi[table])
+	TBL_WMIN(tbl) = TBL_WMIN(Memi[table])
+	TBL_WMAX(tbl) = TBL_WMAX(Memi[table])
+	TBL_NWMAX(tbl) = TBL_NWMAX(Memi[table])
+	TBL_CONFIRM(tbl) = TBL_CONFIRM(Memi[table])
+end
diff --git a/noao/onedspec/dispcor/dispcor.h b/noao/onedspec/dispcor/dispcor.h
new file mode 100644
index 00000000..51167973
--- /dev/null
+++ b/noao/onedspec/dispcor/dispcor.h
@@ -0,0 +1,16 @@
+# Aperture data structure
+
+define	LEN_AP	($1*20)				# Length of DC data structure
+define	DC_PL	Memi[$1+($2-1)*20+1]		# Physical line number
+define	DC_AP	Memi[$1+($2-1)*20+2]		# Aperture number
+define	DC_BM	Memi[$1+($2-1)*20+3]		# Beam number
+define	DC_DT	Memi[$1+($2-1)*20+4]		# Dispersion type
+define	DC_NW	Memi[$1+($2-1)*20+5]		# Number of pixels in spectrum
+define	DC_W1	Memd[P2D($1+($2-1)*20+6)]	# Wavelength of first pixel
+define	DC_W2	Memd[P2D($1+($2-1)*20+8)]	# Wavelength of last pixel
+define	DC_DW	Memd[P2D($1+($2-1)*20+10)]	# Wavelength interval per pixel
+define	DC_Z	Memd[P2D($1+($2-1)*20+12)]	# Redshift
+define	DC_LW	Memr[P2R($1+($2-1)*20+14)]	# Aperture lower limit (2)
+define	DC_UP	Memr[P2R($1+($2-1)*20+16)]	# Aperture upper limit (2)
+define	DC_CF	Memi[$1+($2-1)*20+18]		# Pointer to coefficients
+define	DC_UN	Memi[$1+($2-1)*20+19]		# Units
diff --git a/noao/onedspec/dispcor/dispcor.x b/noao/onedspec/dispcor/dispcor.x
new file mode 100644
index 00000000..7f2c32a8
--- /dev/null
+++ b/noao/onedspec/dispcor/dispcor.x
@@ -0,0 +1,233 @@
+include	<math/iminterp.h>
+
+# DISPCOR -- Dispersion correct input spectrum to output spectrum.
+# This procedure uses the MWCS forward and inverse transformations
+# and interpolate the input data, conserving flux if desired.  Image
+# interpolation uses the image interpolation package and flux conservation
+# integrates the interpolation function across the output pixel.  This
+# procedure does some CLIO to get the interpolation function and to
+# query whether to conserve flux.
+
+procedure dispcor (cti, linei, cto, lineo, in, npts, out, nw, flux)
+
+pointer	cti			#I MWCS input inverse transformation
+int	linei			#I Spectrum line
+pointer	cto			#I MWCS output forward transformation
+int	lineo			#I Spectrum line
+real	in[npts]		#I Input spectrum
+int	npts			#I Number of input pixels
+real	out[nw]			#O Output spectrum
+int	nw			#I Number of output pixels
+bool	flux			#I Conserve flux
+
+char	interp[10]
+bool	ofb_a, ofb_b
+int	i, j, ia, ib, clgwrd()
+real	a, b, sum, asieval(), asigrl()
+double	x, xmin, xmax, w, y1, y2, smw_c1trand()
+pointer	asi, temp
+
+begin
+	# Get the image buffers fit the interpolation function to the
+	# input spectrum.  Extend the interpolation by one pixel at each end.
+
+	call malloc (temp, npts+2, TY_REAL)
+	call amovr (in, Memr[temp+1], npts)
+	Memr[temp] = in[1]
+	Memr[temp+npts+1] = in[npts]
+
+	call asiinit (asi, clgwrd ("interp", interp, 10, II_FUNCTIONS))
+	call asifit (asi, Memr[temp], npts+2)
+
+	call mfree (temp, TY_REAL)
+
+	# Determine edges of output pixels in input spectrum and integrate
+	# using ASIGRL.  If not flux conserving take the average.
+
+	xmin = 0.5
+	xmax = npts + 0.5
+
+	x = 0.5
+	if (IS_INDEFI(lineo))
+	    w = smw_c1trand (cto, x)
+	else {
+	    y1 = lineo
+	    call smw_c2trand (cto, x, y1, w, y2)
+	}
+	if (IS_INDEFI(linei))
+	    x = smw_c1trand (cti, w)
+	else {
+	    #y2 = linei
+	    call smw_c2trand (cti, w, y2, x, y1)
+	}
+	ofb_b = (x < xmin || x > xmax)
+	b = max (xmin, min (xmax, x)) + 1
+	do i = 1, nw {
+	    ofb_a = ofb_b
+	    a = b
+	    x = i + 0.5
+	    if (IS_INDEFI(lineo))
+		w = smw_c1trand (cto, x)
+	    else {
+		y1 = lineo
+		call smw_c2trand (cto, x, y1, w, y2)
+	    }
+	    if (IS_INDEFI(linei))
+		x = smw_c1trand (cti, w)
+	    else {
+		#y2 = linei
+		call smw_c2trand (cti, w, y2, x, y1)
+	    }
+	    ofb_b = (x < xmin || x > xmax)
+	    b = max (xmin, min (xmax, x)) + 1
+	    if (ofb_a && ofb_b)
+		out[i] = 0.
+	    else if (a <= b) {
+		ia = nint (a + 0.5)
+		ib = nint (b - 0.5)
+		if (abs (a+0.5-ia) < 0.00001 && abs (b-0.5-ib) < 0.00001) {
+		    sum = 0.
+		    do j = ia, ib
+			sum = sum + asieval (asi, real(j))
+		    out[i] = sum
+		} else
+		    out[i] = asigrl (asi, a, b)
+		if (!flux)
+		    out[i] = out[i] / max (b - a, 1e-4)
+	    } else {
+		ib = nint (b + 0.5)
+		ia = nint (a - 0.5)
+		if (abs (a-0.5-ia) < 0.00001 && abs (b+0.5-ib) < 0.00001) {
+		    sum = 0.
+		    do j = ib, ia
+			sum = sum + asieval (asi, real(j))
+		    out[i] = sum
+		} else
+		    out[i] = asigrl (asi, b, a)
+		if (!flux)
+		    out[i] = out[i] / max (a - b, 1e-4)
+	    }
+	}
+
+	call asifree (asi)
+end
+
+
+# DISPCORA -- Dispersion correct input spectrum to output spectrum.
+# This procedure uses the MWCS forward and inverse transformations
+# and interpolate the input data, conserving flux if desired.  Image
+# interpolation uses the image interpolation package and flux conservation
+# integrates the interpolation function across the output pixel.  This
+# procedure does some CLIO to get the interpolation function and to
+# query whether to conserve flux.
+# 
+# This differs from DISPCOR by the "blank" argument.
+
+procedure dispcora (cti, linei, cto, lineo, in, npts, out, nw, flux, blank)
+
+pointer	cti			#I MWCS input inverse transformation
+int	linei			#I Spectrum line
+pointer	cto			#I MWCS output forward transformation
+int	lineo			#I Spectrum line
+real	in[npts]		#I Input spectrum
+int	npts			#I Number of input pixels
+real	out[nw]			#O Output spectrum
+int	nw			#I Number of output pixels
+bool	flux			#I Conserve flux
+real	blank			#I Out of bounds value (INDEF to leave unchanged
+
+char	interp[10]
+bool	ofb_a, ofb_b
+int	i, j, ia, ib, clgwrd()
+real	a, b, sum, asieval(), asigrl()
+double	x, xmin, xmax, w, y1, y2, smw_c1trand()
+pointer	asi, temp
+
+begin
+	# Get the image buffers fit the interpolation function to the
+	# input spectrum.  Extend the interpolation by one pixel at each end.
+
+	call malloc (temp, npts+2, TY_REAL)
+	call amovr (in, Memr[temp+1], npts)
+	Memr[temp] = in[1]
+	Memr[temp+npts+1] = in[npts]
+
+	call asiinit (asi, clgwrd ("interp", interp, 10, II_FUNCTIONS))
+	call asifit (asi, Memr[temp], npts+2)
+
+	call mfree (temp, TY_REAL)
+
+	# Determine edges of output pixels in input spectrum and integrate
+	# using ASIGRL.  If not flux conserving take the average.
+
+	xmin = 0.5
+	xmax = npts + 0.5
+
+	x = 0.5
+	if (IS_INDEFI(lineo))
+	    w = smw_c1trand (cto, x)
+	else {
+	    y1 = lineo
+	    call smw_c2trand (cto, x, y1, w, y2)
+	}
+	if (IS_INDEFI(linei))
+	    x = smw_c1trand (cti, w)
+	else {
+	    #y2 = linei
+	    call smw_c2trand (cti, w, y2, x, y1)
+	}
+	ofb_b = (x < xmin || x > xmax)
+	b = max (xmin, min (xmax, x)) + 1
+	do i = 1, nw {
+	    ofb_a = ofb_b
+	    a = b
+	    x = i + 0.5
+	    if (IS_INDEFI(lineo))
+		w = smw_c1trand (cto, x)
+	    else {
+		y1 = lineo
+		call smw_c2trand (cto, x, y1, w, y2)
+	    }
+	    if (IS_INDEFI(linei))
+		x = smw_c1trand (cti, w)
+	    else {
+		#y2 = linei
+		call smw_c2trand (cti, w, y2, x, y1)
+	    }
+	    ofb_b = (x < xmin || x > xmax)
+	    b = max (xmin, min (xmax, x)) + 1
+	    if (ofb_a && ofb_b) {
+		if (!IS_INDEFR(blank))
+		    out[i] = blank
+	    } else if (a == b) {
+		if (!IS_INDEFR(blank))
+		    out[i] = blank
+	    } else if (a < b) {
+		ia = nint (a + 0.5)
+		ib = nint (b - 0.5)
+		if (abs (a+0.5-ia) < 0.00001 && abs (b-0.5-ib) < 0.00001) {
+		    sum = 0.
+		    do j = ia, ib
+			sum = sum + asieval (asi, real(j))
+		    out[i] = sum
+		} else
+		    out[i] = asigrl (asi, a, b)
+		if (!flux)
+		    out[i] = out[i] / max (b - a, 1e-4)
+	    } else {
+		ib = nint (b + 0.5)
+		ia = nint (a - 0.5)
+		if (abs (a-0.5-ia) < 0.00001 && abs (b+0.5-ib) < 0.00001) {
+		    sum = 0.
+		    do j = ib, ia
+			sum = sum + asieval (asi, real(j))
+		    out[i] = sum
+		} else
+		    out[i] = asigrl (asi, b, a)
+		if (!flux)
+		    out[i] = out[i] / max (a - b, 1e-4)
+	    }
+	}
+
+	call asifree (asi)
+end
diff --git a/noao/onedspec/dispcor/mkpkg b/noao/onedspec/dispcor/mkpkg
new file mode 100644
index 00000000..e609106e
--- /dev/null
+++ b/noao/onedspec/dispcor/mkpkg
@@ -0,0 +1,28 @@
+# DISPCOR Task
+
+$checkout libpkg.a ..
+$update   libpkg.a
+$checkin  libpkg.a ..
+$exit
+
+libpkg.a:
+	dcio.x		dispcor.h <error.h> <imhdr.h> <imset.h> <pkg/dttext.h>\
+			<smw.h> <units.h>
+	dctable.x	dctable.h <imhdr.h> <mach.h> <smw.h>
+	dispcor.x	<math/iminterp.h>
+	ranges.x	<ctype.h> <mach.h>
+	refaverage.x	refspectra.h
+	reffollow.x	refspectra.h <mach.h>
+	refgspec.x	refspectra.com refspectra.h <error.h>
+	refinterp.x	refspectra.h <mach.h>
+	refmatch.x	refspectra.h
+	refmsgs.x	refspectra.com refspectra.h
+	refnearest.x	refspectra.h <mach.h>
+	refnoextn.x	
+	refprecede.x	refspectra.h <mach.h>
+	refspectra.x	refspectra.com refspectra.h
+	reftable.x	refspectra.h <error.h>
+	t_dispcor.x	dctable.h dispcor.h <error.h> <imhdr.h> <imio.h>\
+			<mach.h> <mwset.h> <smw.h> <units.h>
+	t_disptrans.x	<error.h> <imhdr.h> <math/curfit.h> <smw.h> <units.h>
+	;
diff --git a/noao/onedspec/dispcor/ranges.x b/noao/onedspec/dispcor/ranges.x
new file mode 100644
index 00000000..403b81f7
--- /dev/null
+++ b/noao/onedspec/dispcor/ranges.x
@@ -0,0 +1,239 @@
+include	<mach.h>
+include	<ctype.h>
+
+define	FIRST	0		# Default starting range
+define	LAST	MAX_INT		# Default ending range
+define	STEP	1		# Default step
+define	EOLIST	-1		# End of list
+
+# DECODE_RANGES -- Parse a string containing a list of integer numbers or
+# ranges, delimited by either spaces or commas.  Return as output a list
+# of ranges defining a list of numbers, and the count of list numbers.
+# Range limits must be positive nonnegative integers.  ERR is returned as
+# the function value if a conversion error occurs.  The list of ranges is
+# delimited by EOLIST.
+
+int procedure decode_ranges (range_string, ranges, max_ranges, nvalues)
+
+char	range_string[ARB]	# Range string to be decoded
+int	ranges[3, max_ranges]	# Range array
+int	max_ranges		# Maximum number of ranges
+int	nvalues			# The number of values in the ranges
+
+int	ip, nrange, first, last, step, ctoi()
+
+begin
+	ip = 1
+	nvalues = 0
+
+	do nrange = 1, max_ranges - 1 {
+	    # Defaults to all nonnegative integers
+	    first = FIRST
+	    last = LAST
+	    step = STEP
+
+	    # Skip delimiters
+	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		ip = ip + 1
+
+	    # Get first limit.
+	    # Must be a number, '-', 'x', or EOS.  If not return ERR.
+	    if (range_string[ip] == EOS) {			# end of list
+		if (nrange == 1) {
+		    # Null string defaults
+		    ranges[1, 1] = first
+		    ranges[2, 1] = last
+		    ranges[3, 1] = step
+		    ranges[1, 2] = EOLIST
+	    	    nvalues = MAX_INT
+		    return (OK)
+		} else {
+		    ranges[1, nrange] = EOLIST
+		    return (OK)
+		}
+	    } else if (range_string[ip] == '-')
+		;
+	    else if (range_string[ip] == 'x')
+		;
+	    else if (IS_DIGIT(range_string[ip])) {		# ,n..
+		if (ctoi (range_string, ip, first) == 0)
+		    return (ERR)
+	    } else
+		return (ERR)
+
+	    # Skip delimiters
+	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		ip = ip + 1
+
+	    # Get last limit
+	    # Must be '-', or 'x' otherwise last = first.
+	    if (range_string[ip] == 'x')
+		;
+	    else if (range_string[ip] == '-') {
+		ip = ip + 1
+	        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		    ip = ip + 1
+		if (range_string[ip] == EOS)
+		    ;
+		else if (IS_DIGIT(range_string[ip])) {
+		    if (ctoi (range_string, ip, last) == 0)
+		        return (ERR)
+		} else if (range_string[ip] == 'x')
+		    ;
+		else
+		    return (ERR)
+	    } else
+		last = first
+
+	    # Skip delimiters
+	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		ip = ip + 1
+
+	    # Get step.
+	    # Must be 'x' or assume default step.
+	    if (range_string[ip] == 'x') {
+		ip = ip + 1
+	        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		    ip = ip + 1
+		if (range_string[ip] == EOS)
+		    ;
+		else if (IS_DIGIT(range_string[ip])) {
+		    if (ctoi (range_string, ip, step) == 0)
+		        ;
+		} else if (range_string[ip] == '-')
+		    ;
+		else
+		    return (ERR)
+	    }
+
+	    # Output the range triple.
+	    ranges[1, nrange] = first
+	    ranges[2, nrange] = last
+	    ranges[3, nrange] = step
+	    nvalues = nvalues + abs (last-first) / step + 1
+	}
+
+	return (ERR)					# ran out of space
+end
+
+
+# GET_NEXT_NUMBER -- Given a list of ranges and the current file number,
+# find and return the next file number.  Selection is done in such a way
+# that list numbers are always returned in monotonically increasing order,
+# regardless of the order in which the ranges are given.  Duplicate entries
+# are ignored.  EOF is returned at the end of the list.
+
+int procedure get_next_number (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Both input and output parameter
+
+int	ip, first, last, step, next_number, remainder
+
+begin
+	# If number+1 is anywhere in the list, that is the next number,
+	# otherwise the next number is the smallest number in the list which
+	# is greater than number+1.
+
+	number = number + 1
+	next_number = MAX_INT
+
+	for (ip=1;  ranges[ip] != EOLIST;  ip=ip+3) {
+	    first = min (ranges[ip], ranges[ip+1])
+	    last = max (ranges[ip], ranges[ip+1])
+	    step = ranges[ip+2]
+	    if (number >= first && number <= last) {
+		remainder = mod (number - first, step)
+		if (remainder == 0)
+		    return (number)
+		if (number - remainder + step <= last)
+		    next_number = number - remainder + step
+	    } else if (first > number)
+		next_number = min (next_number, first)
+	}
+
+	if (next_number == MAX_INT)
+	    return (EOF)
+	else {
+	    number = next_number
+	    return (number)
+	}
+end
+
+
+# GET_PREVIOUS_NUMBER -- Given a list of ranges and the current file number,
+# find and return the previous file number.  Selection is done in such a way
+# that list numbers are always returned in monotonically decreasing order,
+# regardless of the order in which the ranges are given.  Duplicate entries
+# are ignored.  EOF is returned at the end of the list.
+
+int procedure get_previous_number (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Both input and output parameter
+
+int	ip, first, last, step, next_number, remainder
+
+begin
+	# If number-1 is anywhere in the list, that is the previous number,
+	# otherwise the previous number is the largest number in the list which
+	# is less than number-1.
+
+	number = number - 1
+	next_number = 0
+
+	for (ip=1;  ranges[ip] != EOLIST;  ip=ip+3) {
+	    first = min (ranges[ip], ranges[ip+1])
+	    last = max (ranges[ip], ranges[ip+1])
+	    step = ranges[ip+2]
+	    if (number >= first && number <= last) {
+		remainder = mod (number - first, step)
+		if (remainder == 0)
+		    return (number)
+		if (number - remainder >= first)
+		    next_number = number - remainder
+	    } else if (last < number) {
+		remainder = mod (last - first, step)
+		if (remainder == 0)
+		    next_number = max (next_number, last)
+		else if (last - remainder >= first)
+		    next_number = max (next_number, last - remainder)
+	    }
+	}
+
+	if (next_number == 0)
+	    return (EOF)
+	else {
+	    number = next_number
+	    return (number)
+	}
+end
+
+
+# IS_IN_RANGE -- Test number to see if it is in range.
+# If the number is INDEFI then it is mapped to the maximum integer.
+
+bool procedure is_in_range (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Number to be tested against ranges
+
+int	ip, first, last, step, num
+
+begin
+	if (IS_INDEFI (number))
+	    num = MAX_INT
+	else
+	    num = number
+
+	for (ip=1;  ranges[ip] != EOLIST;  ip=ip+3) {
+	    first = min (ranges[ip], ranges[ip+1])
+	    last = max (ranges[ip], ranges[ip+1])
+	    step = ranges[ip+2]
+	    if (num >= first && num <= last)
+		if (mod (num - first, step) == 0)
+		    return (true)
+	}
+
+	return (false)
+end
diff --git a/noao/onedspec/dispcor/refaverage.x b/noao/onedspec/dispcor/refaverage.x
new file mode 100644
index 00000000..e25866c4
--- /dev/null
+++ b/noao/onedspec/dispcor/refaverage.x
@@ -0,0 +1,84 @@
+include	"refspectra.h"
+
+# REFAVERAGE -- Assign reference spectrum by averageing reference list.
+# In earlier version the reference apertures were always set to all
+
+procedure refaverage (input, refs)
+
+pointer	input			# List of input spectra
+pointer	refs			# List of reference spectra
+
+int	ap
+double	sortval
+real	wt1, wt2
+pointer	sp, image, ref1, ref2, gval
+
+bool	refgref(), refginput()
+int	imtgetim(), imtlen()
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (ref1, SZ_FNAME, TY_CHAR)
+	call salloc (ref2, SZ_FNAME, TY_CHAR)
+
+	# Get reference spectra to average.
+	switch (imtlen (refs)) {
+	case 0:
+	    call error (0, "No reference spectra specified")
+	case 1:
+	    ap = imtgetim (refs, Memc[ref1], SZ_FNAME)
+	    call refnoextn (Memc[ref1])
+	    if (!refgref (Memc[ref1], ap, sortval, gval)) {
+		call sfree (sp)
+		return
+	    }
+	    wt1 = 1.
+	    wt2 = 0.
+	case 2:
+	    ap = imtgetim (refs, Memc[ref1], SZ_FNAME)
+	    ap = imtgetim (refs, Memc[ref2], SZ_FNAME)
+	    call refnoextn (Memc[ref1])
+	    call refnoextn (Memc[ref2])
+	    if (!refgref (Memc[ref1], ap, sortval, gval)) {
+		call sfree (sp)
+		return
+	    }
+	    if (!refgref (Memc[ref2], ap, sortval, gval)) {
+		call sfree (sp)
+		return
+	    }
+	    wt1 = 0.5
+	    wt2 = 0.5
+	default:
+	    ap = imtgetim (refs, Memc[ref1], SZ_FNAME)
+	    ap = imtgetim (refs, Memc[ref2], SZ_FNAME)
+	    call refnoextn (Memc[ref1])
+	    call refnoextn (Memc[ref2])
+	    if (!refgref (Memc[ref1], ap, sortval, gval)) {
+		call sfree (sp)
+		return
+	    }
+	    if (!refgref (Memc[ref2], ap, sortval, gval)) {
+		call sfree (sp)
+		return
+	    }
+	    wt1 = 0.5
+	    wt2 = 0.5
+	    call eprintf ("WARNING: Averaging only first two reference spectra")
+	}
+
+	# Assign reference spectra to each input spectrum.
+	# Skip spectra which are not of the appropriate aperture
+	# or have been assigned previously (unless overriding).
+
+	while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refginput (Memc[image], ap, sortval, gval))
+		next
+
+	    call refspectra (Memc[image], Memc[ref1], wt1, Memc[ref2], wt2)
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/reffollow.x b/noao/onedspec/dispcor/reffollow.x
new file mode 100644
index 00000000..a320c504
--- /dev/null
+++ b/noao/onedspec/dispcor/reffollow.x
@@ -0,0 +1,114 @@
+include	<mach.h>
+include	"refspectra.h"
+
+
+# REFFOLLOW -- Assign following reference spectrum based on sort key.
+# If there is no following spectrum assign the nearest preceding spectrum.
+
+procedure reffollow (input, refs)
+
+pointer	input			# List of input spectra
+pointer	refs			# List of reference spectra
+
+bool	ignoreaps		# Ignore apertures?
+
+int	i, i1, i2, nrefs, ap
+double	sortval, d, d1, d2
+pointer	sp, image, gval, refimages, refaps, refvals, refgvals
+
+bool	clgetb(), streq(), refginput(), refgref()
+int	imtgetim(), imtlen()
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+
+	# Task parameters
+	ignoreaps = clgetb ("ignoreaps")
+
+	# Tabulate reference spectra.  This expands the reference list,
+	# checks the spectrum is a reference spectrum of the appropriate
+	# aperture.
+
+	call salloc (refimages, imtlen (refs), TY_INT)
+	call salloc (refaps, imtlen (refs), TY_INT)
+	call salloc (refvals, imtlen (refs), TY_DOUBLE)
+	call salloc (refgvals, imtlen (refs), TY_INT)
+	nrefs = 0
+	while (imtgetim (refs, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refgref (Memc[image], ap, sortval, gval))
+		next
+
+	    for (i=0; i<nrefs; i=i+1)
+		if (streq (Memc[image], Memc[Memi[refimages+i]]))
+		    break
+	    if (i == nrefs) {
+		call salloc (Memi[refimages+nrefs], SZ_FNAME, TY_CHAR)
+		call salloc (Memi[refgvals+nrefs], SZ_FNAME, TY_CHAR)
+		call strcpy (Memc[image], Memc[Memi[refimages+i]], SZ_FNAME)
+		Memi[refaps+i] = ap
+		Memd[refvals+i] = sortval
+		call strcpy (Memc[gval], Memc[Memi[refgvals+i]], SZ_FNAME)
+		nrefs = i + 1
+	    }
+	}
+	if (nrefs < 1)
+	    call error (0, "No reference images specified")
+
+	# Assign following reference spectra to each input spectrum.
+	# Skip input spectra which are not of the appropriate aperture
+	# or have been assigned previously (unless overriding).
+
+	while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refginput (Memc[image], ap, sortval, gval))
+		next
+
+	    i1 = 0
+	    i2 = 0
+	    d1 = MAX_REAL
+	    d2 = -MAX_REAL
+	    do i = 1, nrefs {
+		if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]]))
+		    next
+	        if (!ignoreaps  && ap != Memi[refaps+i-1])
+		    next
+	        d = sortval - Memd[refvals+i-1]
+	        if ((d > 0.) && (d < d1)) {
+		    i1 = i
+		    d1 = d
+	        }
+	        if ((d <= 0.) && (d > d2)) {
+		    i2 = i
+		    d2 = d
+	        }
+	    }
+
+	    if (i2 > 0)		# Nearest following spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i2-1]], 1.,
+		    Memc[Memi[refimages+i2-1]], 0.)
+	    else if (i1 > 0)	# Nearest preceding spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], 1.,
+		    Memc[Memi[refimages+i1-1]], 0.)
+	    else {		# No reference spectrum found
+		call refprint (STDERR, NO_REFSPEC, Memc[image], "", "", "",
+		    ap, 0, "")
+		do i = 1, nrefs {
+		    if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]])) {
+			call refprint (STDERR, REF_GROUP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		    }
+		    if (!ignoreaps  && ap != Memi[refaps+i-1])
+			call refprint (STDERR, REF_AP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/refgspec.x b/noao/onedspec/dispcor/refgspec.x
new file mode 100644
index 00000000..bb851307
--- /dev/null
+++ b/noao/onedspec/dispcor/refgspec.x
@@ -0,0 +1,268 @@
+include	<error.h>
+include	"refspectra.h"
+
+# REFOPEN   -- Set verbose and log file descriptors and open symbol table.
+# REFCLOSE  -- Close file descriptors and symbol table
+# REFGSPEC  -- Get a spectrum from the symbol table.  Map it only once.
+# REFGINPUT -- Get input spectrum.  Apply various checks.
+# REFGREF   -- Get reference spectrum.  Apply various checks.
+
+define	REF_LEN		6		# Length of reference structure
+define	REF_SORTVAL	Memd[P2D($1)]	# Sort value
+define	REF_AP		Memi[$1+2]	# Aperture number
+define	REF_GVAL	Memi[$1+3]	# Sort value
+define	REF_SPEC1	Memi[$1+4]	# Offset for reference spectrum 1
+define	REF_SPEC2	Memi[$1+5]	# Offset for reference spectrum 2
+
+
+# REFOPEN  -- Set verbose and log file descriptors and open symbol table.
+# The file descriptors and symbol table pointer are in common.  A null
+# file descriptor indicates no output.
+
+procedure refopen ()
+
+bool	clgetb()
+real	clgetr()
+pointer	rng_open(), stopen()
+int	fd, btoi(), clpopnu(), clgfil(), open(), nowhite()
+errchk	open()
+
+include	"refspectra.com"
+
+begin
+	call malloc (sort, SZ_FNAME, TY_CHAR)
+	call malloc (group, SZ_FNAME, TY_CHAR)
+
+	# Check log files
+	logfiles = clpopnu ("logfiles")
+	while (clgfil (logfiles, Memc[sort], SZ_FNAME) != EOF) {
+	    fd = open (Memc[sort], APPEND, TEXT_FILE)
+	    call close (fd)
+	}
+	call clprew (logfiles)
+
+	# Get other parameters
+	call clgstr ("apertures", Memc[sort], SZ_FNAME)
+	iferr (aps = rng_open (Memc[sort], INDEF, INDEF, INDEF))
+	    call error (0, "Bad aperture list")
+	call clgstr ("refaps", Memc[sort], SZ_FNAME)
+	iferr (raps = rng_open (Memc[sort], INDEF, INDEF, INDEF))
+	    call error (0, "Bad reference aperture list")
+	call clgstr ("sort", Memc[sort], SZ_FNAME)
+	call clgstr ("group", Memc[group], SZ_FNAME)
+	time = btoi (clgetb ("time"))
+	timewrap = clgetr ("timewrap")
+	verbose = btoi (clgetb ("verbose"))
+
+	fd = nowhite (Memc[sort], Memc[sort], SZ_FNAME)
+	fd = nowhite (Memc[group], Memc[group], SZ_FNAME)
+
+	# Open symbol table.
+	stp = stopen ("refspectra", 10, 20, 10*SZ_FNAME)
+end
+
+
+# REFCLOSE  -- Finish up
+
+procedure refclose ()
+
+include	"refspectra.com"
+
+begin
+	call mfree (sort, TY_CHAR)
+	call mfree (group, TY_CHAR)
+	call clpcls (logfiles)
+	call stclose (stp)
+	call rng_close (raps)
+	call rng_close (aps)
+end
+
+
+# REFGSPEC  -- Get a spectrum from the symbol table.  Map it only once.
+# All access to spectra is through this routine.  It returns header parameters.
+# Because the spectra may be accessed in very random order and many times
+# the information is stored in a symbol table keyed on the spectrum name.
+# The spectrum need be mapped only once!  Any error from IMMAP is returned.
+
+procedure refgspec (spec, ap, sortval, gval, ref1, ref2)
+
+char	spec[ARB]	# Spectrum image name
+int	ap		# Spectrum aperture number
+double	sortval		# Spectrum sort value
+pointer	gval		# Group string
+pointer	ref1		# Reference spectrum 1
+pointer	ref2		# Reference spectrum 2
+
+pointer	sym, stfind(), stenter(), stpstr(), strefsbuf()
+pointer	im, str, immap()
+bool	streq()
+int	imgeti(), strlen()
+double	imgetd()
+errchk	immap, imgetd, imgstr
+
+include	"refspectra.com"
+
+begin
+	# Check if spectrum is in the symbol table from a previous call.
+	# If not in the symbol table map the image, get the header parameters,
+	# and store them in the symbol table.
+
+	sym = stfind (stp, spec)
+	if (sym == NULL) {
+	    im = immap (spec, READ_ONLY, 0)
+	    iferr (ap = imgeti (im, "BEAM-NUM"))
+		ap = 1
+
+	    # Failure to find a specified keyword is a fatal error.
+	    iferr {
+		if (Memc[sort] == EOS || streq (Memc[sort], "none") ||
+		    select == MATCH || select == AVERAGE)
+		    sortval = INDEFD
+		else {
+		    sortval = imgetd (im, Memc[sort])
+		    if (time == YES)
+			sortval = mod (sortval + 24. - timewrap, 24.0D0)
+		}
+
+		call malloc (str, SZ_FNAME, TY_CHAR)
+		if (Memc[group] == EOS || streq (Memc[group], "none") ||
+		    select == MATCH || select == AVERAGE)
+		    Memc[str] = EOS
+		else
+		    call imgstr (im, Memc[group], Memc[str], SZ_FNAME)
+		gval = stpstr (stp, Memc[str], strlen (Memc[str])+1)
+	    } then
+		call erract (EA_FATAL)
+
+	    iferr (call imgstr (im, "refspec1", Memc[str], SZ_FNAME))
+		Memc[str] = EOS
+	    ref1 = stpstr (stp, Memc[str], strlen (Memc[str])+1)
+	    iferr (call imgstr (im, "refspec2", Memc[str], SZ_FNAME))
+		Memc[str] = EOS
+	    ref2 = stpstr (stp, Memc[str], strlen (Memc[str])+1)
+	    call mfree (str, TY_CHAR)
+
+	    call imunmap (im)
+
+	    sym = stenter (stp, spec, REF_LEN)
+	    REF_AP(sym) = ap
+	    REF_SORTVAL(sym) = sortval
+	    REF_GVAL(sym) = gval
+	    REF_SPEC1(sym) = ref1
+	    REF_SPEC2(sym) = ref2
+	}
+	ap = REF_AP(sym)
+	sortval = REF_SORTVAL(sym)
+	gval = strefsbuf (stp, REF_GVAL(sym))
+	ref1 = strefsbuf (stp, REF_SPEC1(sym))
+	ref2 = strefsbuf (stp, REF_SPEC2(sym))
+end
+
+
+# REFGINPUT -- Get input spectrum.  Apply various checks.
+# This calls REFGSPEC and then checks:
+#	1. The spectrum is found.
+#	2. The spectrum has not been assigned reference spectra previously.
+#	   If it has then determine whether to override the assignment.
+#	3. Check if the aperture is correct.
+# Return true if the spectrum is acceptable and false if not.
+
+bool procedure refginput (spec, ap, val, gval)
+
+char	spec[ARB]	# Spectrum image name
+int	ap		# Spectrum aperture number (returned)
+double	val		# Spectrum sort value (returned)
+pointer	gval		# Spectrum group value (returned)
+
+bool	clgetb(), rng_elementi()
+pointer	ref1, ref2
+errchk	refgspec
+
+include	"refspectra.com"
+
+define	err_	99
+
+begin
+	# Get the spectrum from the symbol table.
+	iferr (call refgspec (spec, ap, val, gval, ref1, ref2)) {
+	    call refmsgs (NO_SPEC, spec, "", "", "", ap, 0, "")
+	    goto err_
+	}
+
+	# Check if it has a previous reference spectrum.  Override if desired.
+	if (Memc[ref1] != EOS) {
+	    if (!clgetb ("override")) {
+		call refmsgs (DEF_REFSPEC, spec, Memc[ref1], "", "", ap, 0,
+		    Memc[ref2])
+	       goto err_
+	    } else {
+		call refmsgs (OVR_REFSPEC, spec, Memc[ref1], "", "", ap, 0,
+		    Memc[ref2])
+	    }
+	}
+	
+	# Check aperture numbers.
+	if (aps != NULL) {
+	    if (!rng_elementi (aps, ap)) {
+	        call refmsgs (BAD_AP, spec, "", "", "", ap, 0, "")
+	        goto err_
+	    }
+	}
+
+	return (true)
+
+err_
+	return (false)
+end
+
+
+# REFGREF   -- Get reference spectrum.  Apply various checks.
+# This calls REFGSPEC and then checks:
+#	1. The spectrum is found.
+#	2. The spectrum is a reference spectrum, i.e. has an IDENTIFY
+#	   record.  This is signaled by having a reference equivalent to
+#	   itself.
+#	3. Check if the aperture is correct.
+# Return true if the spectrum is acceptable and false if not.
+
+bool procedure refgref (spec, ap, val, gval)
+
+char	spec[ARB]	# Spectrum image name
+int	ap		# Spectrum aperture number (returned)
+double	val		# Spectrum sort value (returned)
+pointer	gval		# Spectrum group value (returned)
+
+bool	strne(), rng_elementi()
+pointer	ref1, ref2
+errchk	refgspec
+
+include	"refspectra.com"
+
+define	err_	99
+
+begin
+	# Get spectrum from symbol table.
+	iferr (call refgspec (spec, ap, val, gval, ref1, ref2)) {
+	    call refmsgs (NO_REF, spec, "", "", "", ap, 0, "")
+	    goto err_
+	}
+
+	# Check if spectrum is a reference spectrum.
+	if (strne (spec, Memc[ref1])) {
+	    call refmsgs (NOT_REFSPEC, spec, "", "", "", ap, 0, "")
+	    goto err_
+	}
+	
+	# Check aperture numbers.
+	if (raps != NULL) {
+	    if (!rng_elementi (raps, ap)) {
+	        call refmsgs (BAD_REFAP, spec, "", "", "", ap, 0, "")
+	        goto err_
+	    }
+	}
+
+	return (true)
+
+err_
+	return (false)
+end
diff --git a/noao/onedspec/dispcor/refinterp.x b/noao/onedspec/dispcor/refinterp.x
new file mode 100644
index 00000000..b074c053
--- /dev/null
+++ b/noao/onedspec/dispcor/refinterp.x
@@ -0,0 +1,127 @@
+include	<mach.h>
+include	"refspectra.h"
+
+
+# REFINTERP -- Assign reference spectra to interpolate between based on sort
+# key.  The nearest preceding and following spectra are assigned weights based
+# on their distance.  If there is no preceding and following spectrum then
+# the nearest spectrum is assigned.
+
+procedure refinterp (input, refs)
+
+pointer	input			# List of input spectra
+pointer	refs			# List of reference spectra
+
+bool	ignoreaps		# Ignore apertures?
+
+int	i, i1, i2, nrefs, ap
+double	sortval, d, d1, d2
+real	wt1, wt2
+pointer	sp, image, gval, refimages, refaps, refvals, refgvals
+
+bool	clgetb(), streq(), refginput(), refgref()
+int	imtgetim(), imtlen()
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+
+	# Task parameters
+	ignoreaps = clgetb ("ignoreaps")
+
+	# Tabulate reference spectra.  This expands the reference list,
+	# checks the spectrum is a reference spectrum of the appropriate
+	# aperture.
+
+	call salloc (refimages, imtlen (refs), TY_INT)
+	call salloc (refaps, imtlen (refs), TY_INT)
+	call salloc (refvals, imtlen (refs), TY_DOUBLE)
+	call salloc (refgvals, imtlen (refs), TY_INT)
+	nrefs = 0
+	while (imtgetim (refs, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refgref (Memc[image], ap, sortval, gval))
+		next
+
+	    for (i=0; i<nrefs; i=i+1)
+		if (streq (Memc[Memi[refimages+i]], Memc[image]))
+		    break
+	    if (i == nrefs) {
+		call salloc (Memi[refimages+nrefs], SZ_FNAME, TY_CHAR)
+		call salloc (Memi[refgvals+nrefs], SZ_FNAME, TY_CHAR)
+		call strcpy (Memc[image], Memc[Memi[refimages+i]], SZ_FNAME)
+		Memi[refaps+i] = ap
+		Memd[refvals+i] = sortval
+		call strcpy (Memc[gval], Memc[Memi[refgvals+i]], SZ_FNAME)
+		nrefs = i + 1
+	    }
+	}
+	if (nrefs < 1)
+	    call error (0, "No reference images specified")
+
+
+	# Assign following reference spectra to each input spectrum.
+	# Skip input spectra which are not of the appropriate aperture
+	# or have been assigned previously (unless overriding).
+
+	while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refginput (Memc[image], ap, sortval, gval))
+		next
+
+	    i1 = 0
+	    i2 = 0
+	    d1 = MAX_REAL
+	    d2 = -MAX_REAL
+	    do i = 1, nrefs {
+		if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]]))
+		    next
+	        if (!ignoreaps  && ap != Memi[refaps+i-1])
+		    next
+	        d = sortval - Memd[refvals+i-1]
+	        if ((d >= 0.) && (d < d1)) {
+		    i1 = i
+		    d1 = d
+	        } else if ((d <= 0.) && (d > d2)) {
+		    i2 = i
+		    d2 = d
+	        }
+	    }
+
+	    if (i1 > 0 && i2 > 0) {	# Weight spectra
+		if (d1 - d2 == 0.) {
+		    wt1 = 0.5
+		    wt2 = 0.5
+		} else {
+		    wt1 = -d2 / (d1 - d2)
+		    wt2 =  d1 / (d1 - d2)
+		}
+		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], wt1,
+		    Memc[Memi[refimages+i2-1]], wt2)
+	    } else if (i1 > 0)	# Nearest preceding spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], 1.,
+		    Memc[Memi[refimages+i1-1]], 0.)
+	    else if (i2 > 0)	# Nearest following spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i2-1]], 1.,
+		    Memc[Memi[refimages+i2-1]], 0.)
+	    else {		# No reference spectrum found
+		call refprint (STDERR, NO_REFSPEC, Memc[image], "", "", "",
+		    ap, 0, "")
+		do i = 1, nrefs {
+		    if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]])) {
+			call refprint (STDERR, REF_GROUP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		    }
+		    if (!ignoreaps  && ap != Memi[refaps+i-1])
+			call refprint (STDERR, REF_AP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/refmatch.x b/noao/onedspec/dispcor/refmatch.x
new file mode 100644
index 00000000..c94a7113
--- /dev/null
+++ b/noao/onedspec/dispcor/refmatch.x
@@ -0,0 +1,43 @@
+include	"refspectra.h"
+
+# REFMATCH -- Assign reference spectrum by match against reference list.
+
+procedure refmatch (input, refs)
+
+pointer	input			# List of input spectra
+pointer	refs			# List of reference spectra
+
+int	ap
+double	sortval
+pointer	sp, image, refimage, gval
+
+bool	refgref(), refginput()
+int	imtgetim(), imtlen()
+
+begin
+	if (imtlen (input) != imtlen (refs))
+	    call error (0, "Input and reference list have different lengths")
+
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (refimage, SZ_FNAME, TY_CHAR)
+
+	# Assign reference spectra to each input spectrum.
+	# Skip spectra which are not of the appropriate aperture
+	# or have been assigned previously (unless overriding).
+
+	while ((imtgetim (input, Memc[image], SZ_FNAME) != EOF) &&
+	       (imtgetim (refs, Memc[refimage], SZ_FNAME) != EOF)) {
+	    call refnoextn (Memc[image])
+	    call refnoextn (Memc[refimage])
+	    if (!refginput (Memc[image], ap, sortval, gval))
+		next
+	    if (!refgref (Memc[refimage], ap, sortval, gval))
+		next
+
+	    call refspectra (Memc[image], Memc[refimage], 1.,
+		    Memc[refimage], 0.)
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/refmsgs.x b/noao/onedspec/dispcor/refmsgs.x
new file mode 100644
index 00000000..a374088e
--- /dev/null
+++ b/noao/onedspec/dispcor/refmsgs.x
@@ -0,0 +1,108 @@
+include	"refspectra.h"
+
+
+# REFMSGS -- Print any verbose messages to log files.  All messages
+# except the assignments go through this procedure.  It calls REFPRINT with
+# each output stream.
+
+procedure refmsgs (msg, spec, ref, gval, gvalref, ap, apref, ref2)
+
+int	msg		# Message code
+char	spec[ARB]	# Spectrum
+char	ref[ARB]	# Reference spectrum
+char	gval[ARB]	# Group value
+char	gvalref[ARB]	# Group value in reference
+int	ap		# Aperture
+int	apref		# Aperture in reference
+char	ref2[ARB]	# Reference spectrum 2
+
+int	fd, clgfil(), open()
+pointer	sp, logfile
+include	"refspectra.com"
+
+begin
+	if (verbose == NO)
+	    return
+
+	call smark (sp)
+	call salloc (logfile, SZ_FNAME, TY_CHAR)
+	while (clgfil (logfiles, Memc[logfile], SZ_FNAME) != EOF) {
+	    fd = open (Memc[logfile], APPEND, TEXT_FILE)
+	    call refprint (fd, msg, spec, ref, gval, gvalref, ap, apref, ref2)
+	    call close (fd)
+	}
+	call clprew (logfiles)
+
+	call sfree (sp)
+end
+
+
+# REFPRINT -- Print requested message with appropriate parameters if non-null
+# stream is specified.
+
+procedure refprint (fd, msg, spec, ref, gval, gvalref, ap, apref, ref2)
+
+int	fd		# File descriptor
+int	msg		# Message code
+char	spec[ARB]	# Spectrum
+char	ref[ARB]	# Reference spectrum
+char	gval[ARB]	# Group value
+char	gvalref[ARB]	# Group value in reference
+int	ap		# Aperture
+int	apref		# Aperture in reference
+char	ref2[ARB]	# Reference spectrum 2
+
+include	"refspectra.com"
+
+begin
+	if (fd == NULL)
+	    return
+
+	switch (msg) {
+	case NO_SPEC:
+	    call fprintf (fd, "[%s] Spectrum not found\n")
+		call pargstr (spec)
+	case NO_REF:
+	    call fprintf (fd, "[%s] Reference spectrum not found\n")
+		call pargstr (spec)
+	case NOT_REFSPEC:
+	    call fprintf (fd, "[%s] Not a reference spectrum\n")
+		call pargstr (spec)
+	case NO_REFSPEC:
+	    call fprintf (fd, "[%s] No reference spectrum found\n")
+		call pargstr (spec)
+	case DEF_REFSPEC:
+	    call fprintf (fd, "[%s] Reference spectra already defined: %s %s\n")
+		call pargstr (spec)
+		call pargstr (ref)
+		call pargstr (ref2)
+	case OVR_REFSPEC:
+	    call fprintf (fd, 
+		"[%s] Overriding previous reference spectra: %s %s\n")
+		call pargstr (spec)
+		call pargstr (ref)
+		call pargstr (ref2)
+	case BAD_AP:
+	    call fprintf (fd, "[%s] Wrong aperture: %d\n")
+		call pargstr (spec)
+		call pargi (ap)
+	case BAD_REFAP:
+	    call fprintf (fd, "[%s] Wrong reference aperture: %d\n")
+		call pargstr (spec)
+		call pargi (ap)
+	case REF_GROUP:
+	    call fprintf (fd, "Input [%s] %s = %s : Ref [%s] %s = %s\n")
+		call pargstr (spec)
+		call pargstr (Memc[group])
+		call pargstr (gval)
+		call pargstr (ref)
+		call pargstr (Memc[group])
+		call pargstr (gvalref)
+	case REF_AP:
+	    call fprintf (fd, "Input [%s] ap = %d : Ref [%s] ap = %d\n")
+		call pargstr (spec)
+		call pargi (ap)
+		call pargstr (ref)
+		call pargi (apref)
+	}
+end
diff --git a/noao/onedspec/dispcor/refnearest.x b/noao/onedspec/dispcor/refnearest.x
new file mode 100644
index 00000000..78ebb62b
--- /dev/null
+++ b/noao/onedspec/dispcor/refnearest.x
@@ -0,0 +1,104 @@
+include	<mach.h>
+include	"refspectra.h"
+
+
+# REFNEAREST -- Assign nearest reference spectrum based on sort key.
+
+procedure refnearest (input, refs)
+
+pointer	input			# List of input spectra
+pointer	refs			# List of reference spectra
+
+bool	ignoreaps		# Ignore apertures?
+
+int	i, i1, nrefs, ap
+double	sortval, d, d1
+pointer	sp, image, gval, refimages, refaps, refvals, refgvals
+
+bool	clgetb(), streq(), refginput(), refgref()
+int	imtgetim(), imtlen()
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+
+	# Task parameters
+	ignoreaps = clgetb ("ignoreaps")
+
+	# Tabulate reference spectra.  This expands the reference list,
+	# checks the spectrum is a reference spectrum of the appropriate
+	# aperture.
+
+	call salloc (refimages, imtlen (refs), TY_POINTER)
+	call salloc (refaps, imtlen (refs), TY_INT)
+	call salloc (refvals, imtlen (refs), TY_DOUBLE)
+	call salloc (refgvals, imtlen (refs), TY_POINTER)
+	nrefs = 0
+	while (imtgetim (refs, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refgref (Memc[image], ap, sortval, gval))
+		next
+
+	    for (i=0; i<nrefs; i=i+1)
+		if (streq (Memc[image], Memc[Memi[refimages+i]]))
+		    break
+	    if (i == nrefs) {
+		call salloc (Memi[refimages+nrefs], SZ_FNAME, TY_CHAR)
+		call salloc (Memi[refgvals+nrefs], SZ_FNAME, TY_CHAR)
+		call strcpy (Memc[image], Memc[Memi[refimages+i]], SZ_FNAME)
+		Memi[refaps+i] = ap
+		Memd[refvals+i] = sortval
+		call strcpy (Memc[gval], Memc[Memi[refgvals+i]], SZ_FNAME)
+		nrefs = i + 1
+	    }
+	}
+	if (nrefs < 1)
+	    call error (0, "No reference images specified")
+
+
+	# Assign nearest reference spectra to each input spectrum.
+	# Skip input spectra which are not of the appropriate aperture
+
+	while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refginput (Memc[image], ap, sortval, gval))
+		next
+
+	    i1 = 0
+	    d1 = MAX_REAL
+	    do i = 1, nrefs {
+		if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]]))
+		    next
+	        if (!ignoreaps && ap != Memi[refaps+i-1])
+		    next
+	        d = abs (sortval - Memd[refvals+i-1])
+	        if (d < d1) {
+		    i1 = i
+		    d1 = d
+	        }
+	    }
+
+	    if (i1 > 0)		# Assign nearest reference spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], 1.,
+		    Memc[Memi[refimages+i1-1]], 0.)
+	    else {		# No reference spectrum found
+		call refprint (STDERR, NO_REFSPEC, Memc[image], "", "", "",
+		    ap, 0, "")
+		do i = 1, nrefs {
+		    if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]])) {
+			call refprint (STDERR, REF_GROUP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		    }
+		    if (!ignoreaps  && ap != Memi[refaps+i-1])
+			call refprint (STDERR, REF_AP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/refnoextn.x b/noao/onedspec/dispcor/refnoextn.x
new file mode 100644
index 00000000..4c48b194
--- /dev/null
+++ b/noao/onedspec/dispcor/refnoextn.x
@@ -0,0 +1,29 @@
+# REFNOEXTN -- Strip any image extensions
+
+procedure refnoextn (spec)
+
+char	spec[ARB]		# Image name
+
+int	i, strlen()
+bool	streq()
+
+begin
+	i = strlen (spec)
+	call imgimage (spec, spec, i)
+
+	i = strlen (spec)
+	switch (spec[i]) {
+	case 'h':
+	    if (i > 3 && spec[i-3] == '.')
+		spec[i-3] = EOS
+	case 'l':
+	    if (i > 2 && streq (spec[i-2], ".pl"))
+		spec[i-2] = EOS
+	case 's':
+	    if (i > 4 && streq (spec[i-4], ".fits"))
+		spec[i-4] = EOS
+	case 't':
+	    if (i > 3 && streq (spec[i-3], ".fit"))
+		spec[i-3] = EOS
+	}
+end
diff --git a/noao/onedspec/dispcor/refprecede.x b/noao/onedspec/dispcor/refprecede.x
new file mode 100644
index 00000000..c2ba0467
--- /dev/null
+++ b/noao/onedspec/dispcor/refprecede.x
@@ -0,0 +1,114 @@
+include	<mach.h>
+include	"refspectra.h"
+
+
+# REFPRECEDE -- Assign preceding reference spectrum based on sort key.
+# If there is no preceding spectrum assign the nearest following spectrum.
+
+procedure refprecede (input, refs)
+
+pointer	input			# List of input spectra
+pointer	refs			# List of reference spectra
+
+bool	ignoreaps		# Ignore aperture numbers?
+
+int	i, i1, i2, nrefs, ap
+double	sortval, d, d1, d2
+pointer	sp, image, gval, refimages, refaps, refvals, refgvals
+
+bool	clgetb(), streq(), refginput(), refgref()
+int	imtgetim(), imtlen()
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+
+	# Task parameters
+	ignoreaps = clgetb ("ignoreaps")
+
+	# Tabulate reference spectra.  This expands the reference list,
+	# checks the spectrum is a reference spectrum of the appropriate
+	# aperture.
+
+	call salloc (refimages, imtlen (refs), TY_INT)
+	call salloc (refaps, imtlen (refs), TY_INT)
+	call salloc (refvals, imtlen (refs), TY_DOUBLE)
+	call salloc (refgvals, imtlen (refs), TY_INT)
+	nrefs = 0
+	while (imtgetim (refs, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refgref (Memc[image], ap, sortval, gval))
+		next
+
+	    for (i=0; i<nrefs; i=i+1)
+		if (streq (Memc[image], Memc[Memi[refimages+i]]))
+		    break
+	    if (i == nrefs) {
+		call salloc (Memi[refimages+nrefs], SZ_FNAME, TY_CHAR)
+		call salloc (Memi[refgvals+nrefs], SZ_FNAME, TY_CHAR)
+		call strcpy (Memc[image], Memc[Memi[refimages+i]], SZ_FNAME)
+		Memi[refaps+i] = ap
+		Memd[refvals+i] = sortval
+		call strcpy (Memc[gval], Memc[Memi[refgvals+i]], SZ_FNAME)
+		nrefs = i + 1
+	    }
+	}
+	if (nrefs < 1)
+	    call error (0, "No reference images specified")
+
+	# Assign preceding reference spectra to each input spectrum.
+	# Skip input spectra which are not of the appropriate aperture
+	# or have been assigned previously (unless overriding).
+
+	while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    if (!refginput (Memc[image], ap, sortval, gval))
+		next
+
+	    i1 = 0
+	    i2 = 0
+	    d1 = MAX_REAL
+	    d2 = -MAX_REAL
+	    do i = 1, nrefs {
+		if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]]))
+		    next
+	        if (!ignoreaps && ap != Memi[refaps+i-1])
+		    next
+	        d = sortval - Memd[refvals+i-1]
+	        if ((d >= 0.) && (d < d1)) {
+		    i1 = i
+		    d1 = d
+	        }
+	        if ((d < 0.) && (d < d2)) {
+		    i2 = i
+		    d2 = d
+	        }
+	    }
+
+	    if (i1 > 0)		# Nearest preceding spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], 1.,
+		    Memc[Memi[refimages+i1-1]], 0.)
+	    else if (i2 > 0)	# Nearest following spectrum
+		call refspectra (Memc[image], Memc[Memi[refimages+i2-1]], 1.,
+		    Memc[Memi[refimages+i2-1]], 0.)
+	    else {		# No reference spectrum found
+		call refprint (STDERR, NO_REFSPEC, Memc[image], "", "", "",
+		    ap, 0, "")
+		do i = 1, nrefs {
+		    if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]])) {
+			call refprint (STDERR, REF_GROUP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		    }
+		    if (!ignoreaps  && ap != Memi[refaps+i-1])
+			call refprint (STDERR, REF_AP, Memc[image],
+			    Memc[Memi[refimages+i-1]], Memc[gval],
+			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
+			next
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/refspectra.com b/noao/onedspec/dispcor/refspectra.com
new file mode 100644
index 00000000..7c5dc622
--- /dev/null
+++ b/noao/onedspec/dispcor/refspectra.com
@@ -0,0 +1,15 @@
+# Common parameters for logging and the spectrum symbol table.
+
+pointer	aps		# Pointer to aperture list
+pointer	raps		# Pointer to reference aperture list
+pointer	sort		# Pointer to sort keyword
+pointer	group		# Pointer to group keyword
+int	select		# Selection type
+int	time		# Is sort keyword a time?
+real	timewrap	# Timewrap parameter
+int	verbose		# Verbose output?
+int	logfiles	# List of log files
+pointer	stp		# Symbol table for previously mapped spectra
+
+common	/refcom/ aps, raps, sort, group, select, time, timewrap, verbose,
+		logfiles, stp
diff --git a/noao/onedspec/dispcor/refspectra.h b/noao/onedspec/dispcor/refspectra.h
new file mode 100644
index 00000000..c7ba397a
--- /dev/null
+++ b/noao/onedspec/dispcor/refspectra.h
@@ -0,0 +1,30 @@
+# Selection method keywords and codes.
+
+define	SELECT	"|match|nearest|preceding|following|interp|average|"
+define	MATCH		1	# Match input and reference lists
+define	NEAREST		2	# Nearest reference
+define	PRECEDING	3	# Preceding reference
+define	FOLLOWING	4	# Following reference
+define	INTERP		5	# Interpolate between nearest references
+define	AVERAGE		6	# Average first two reference spectra
+
+# Reference list types.
+
+define	LIST		1	# References are an image list
+define	TABLE		2	# Referenece are a table
+
+# Maximum number of aperture ranges.
+define	NRANGES		100
+
+# Message codes (see procedure refprint)
+
+define	NO_SPEC		1	# Spectrum not found (immap failed)
+define	NO_REF		2	# Reference spectrum not found (immap failed)
+define	NOT_REFSPEC	3	# Not a reference spectrum
+define	NO_REFSPEC	4	# No reference spectrum found
+define	DEF_REFSPEC	5	# Reference spectra already defined
+define	OVR_REFSPEC	6	# Override reference spectra
+define	BAD_AP		7	# Bad aperture
+define	BAD_REFAP	8	# Bad reference aperture
+define	REF_GROUP	9	# Group
+define	REF_AP		10	# Aperture
diff --git a/noao/onedspec/dispcor/refspectra.x b/noao/onedspec/dispcor/refspectra.x
new file mode 100644
index 00000000..57a18e43
--- /dev/null
+++ b/noao/onedspec/dispcor/refspectra.x
@@ -0,0 +1,186 @@
+include	"refspectra.h"
+
+
+# T_REFSPECTRA -- Assign reference spectra.
+# Reference spectra are assigned to input spectra from a specified list of
+# reference spectra with various criteria.  This procedure only gets some
+# of the task parameters and switches to separate procedures for each
+# implemented assignment method.  The reference spectra may be specified by
+# and image list or a lookup table.  The difference is determined by attempting
+# to map the first reference element in the list as an image.
+
+procedure t_refspectra ()
+
+pointer	input			# List of input images
+pointer	refs			# List of reference images
+#int	select			# Selection method for reference spectra
+int	type			# Type of reference specification
+
+int	clgwrd(), imtgetim()
+pointer	sp, ref, im, imtopenp(), immap()
+errchk	immap
+
+include	"refspectra.com"
+
+begin
+	call smark (sp)
+	call salloc (ref, SZ_LINE, TY_CHAR)
+
+	# Get input and reference spectra lists.  Determine selection method.
+	input = imtopenp ("input")
+	call clgstr ("records", Memc[ref], SZ_LINE)
+	call odr_openp (input, Memc[ref])
+	refs = imtopenp ("references")
+	select = clgwrd ("select", Memc[ref], SZ_FNAME, SELECT)
+
+	# Determine if reference list is a table.
+	if (imtgetim (refs, Memc[ref], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[ref])
+	    iferr {
+		im = immap (Memc[ref], READ_ONLY, 0)
+	        call imunmap (im)
+		type = LIST
+	    } then
+		type = TABLE
+	} else
+	    call error (0, "No reference spectra specified")
+	call imtrew (refs)
+
+	# Initialize confirm flag, symbol table and logging streams.
+	call refconfirm1 ()
+	call refopen ()
+
+	# Switch of reference list type and selection method.
+	if (type == LIST) {
+	    switch (select) {
+	    case MATCH:
+	        call refmatch(input, refs)
+	    case NEAREST:
+	        call refnearest (input, refs)
+	    case PRECEDING:
+	        call refprecede (input, refs)
+	    case FOLLOWING:
+	        call reffollow (input, refs)
+	    case INTERP:
+	        call refinterp (input, refs)
+	    case AVERAGE:
+	        call refaverage (input, refs)
+	    }
+	} else
+	    call reftable (input, Memc[ref], select)
+
+	call refclose ()
+	call imtclose (input)
+	call imtclose (refs)
+	call sfree (sp)
+end
+
+
+# REFSPECTRA -- Confirm and set reference spectra in header.
+# 	1.  Confirm assignments if desired.
+#	2.  Log output to logfiles if desired.
+#	3.  Update assignment if desired.
+# Note that if wt1 > 0.995 then only the first reference spectrum is
+# set with no weight specified.  No weight implies no interpolation.
+
+procedure refspectra (image, ref1, wt1, ref2, wt2)
+
+char	image[ARB]		# Spectrum image name
+char	ref1[ARB]		# Reference spectrum image name
+real	wt1			# Weight
+char	ref2[ARB]		# Reference spectrum image name
+real	wt2			# Weight
+bool	confirm			# Confirm assignments?
+
+int	fd, clgfil(), open(), clgwrd()
+bool	clgetb(), streq()
+pointer	im, sp, str, immap()
+errchk	immap
+
+include	"refspectra.com"
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Confirm assignments.
+	if (confirm) {
+	    if (wt1 < 0.995) {
+	        call printf ("[%s] refspec1='%s %.8g'\n")
+	            call pargstr (image)
+	            call pargstr (ref1)
+		    call pargr (wt1)
+	        call printf ("[%s] refspec2='%s %.8g'   ")
+	            call pargstr (image)
+	            call pargstr (ref2)
+		    call pargr (wt2)
+	    } else {
+	        call printf ("[%s] refspec1='%s'   ")
+	            call pargstr (image)
+	            call pargstr (ref1)
+	    }
+	    call flush (STDOUT)
+	    fd = clgwrd ("answer", Memc[str], SZ_LINE, "|no|yes|YES|")
+	    switch (fd) {
+	    case 1:
+		call sfree (sp)
+		return
+	    case 3:
+		confirm = false
+	    }
+	}
+
+	# Log output.
+	while (clgfil (logfiles, Memc[str], SZ_LINE) != EOF) {
+	    if (streq (Memc[str], "STDOUT") && confirm)
+		next
+	    fd = open (Memc[str], APPEND, TEXT_FILE)
+	    if (wt1 < 0.995) {
+	        call fprintf (fd, "[%s] refspec1='%s %.8g'\n")
+	            call pargstr (image)
+	            call pargstr (ref1)
+		    call pargr (wt1)
+	        call fprintf (fd, "[%s] refspec2='%s %.8g'\n")
+	            call pargstr (image)
+	            call pargstr (ref2)
+		    call pargr (wt2)
+	    } else {
+	        call fprintf (fd, "[%s] refspec1='%s'\n")
+	            call pargstr (image)
+	            call pargstr (ref1)
+	    }
+	    call close (fd)
+	}
+	call clprew (logfiles)
+
+	# If updating the assigments map the spectrum READ_WRITE and set
+	# the keywords REFSPEC1 and REFSPEC2.  REFSPEC2 is not set if not
+	# interpolating.
+
+	if (clgetb ("assign")) {
+	    im = immap (image, READ_WRITE, 0)
+	    if (wt1 < 0.9999995D0) {
+	        call sprintf (Memc[str], SZ_LINE, "%s %.8g")
+	            call pargstr (ref1)
+		    call pargr (wt1)
+	        call imastr (im, "refspec1", Memc[str])
+	        call sprintf (Memc[str], SZ_LINE, "%s %.8g")
+	            call pargstr (ref2)
+		    call pargr (wt2)
+	        call imastr (im, "refspec2", Memc[str])
+	    } else {
+	        call imastr (im, "refspec1", ref1)
+	        iferr (call imdelf (im, "refspec2"))
+		    ;
+	    }
+	    call imunmap (im)
+	}
+
+	call sfree (sp)
+	return
+
+entry refconfirm1 ()
+
+	confirm = clgetb ("confirm")
+
+end
diff --git a/noao/onedspec/dispcor/reftable.x b/noao/onedspec/dispcor/reftable.x
new file mode 100644
index 00000000..abdf1f4a
--- /dev/null
+++ b/noao/onedspec/dispcor/reftable.x
@@ -0,0 +1,109 @@
+include	<error.h>
+include	"refspectra.h"
+
+
+# REFTABLE -- For each input image select reference spectrum list from a table.
+# The table is read from the file and stored in a simple symbol table.
+#
+# The table consists of pairs of words.  The first word is a list of spectra
+# and the second word is the reference spectrum list to be used for each
+# spectrum in the first list.  Note that the first list is  not an input
+# list.  As a convenience if a reference list is missing the preceding list
+# is implied.  Some examples follow.
+#
+#	spec1		spec2,spec3,spec4
+#	spec5
+#	spec6,spec7	spect8,spec9
+#	spec10		spec11
+#	spec12		spec13
+#	spec14		spec15
+
+procedure reftable (list, table, select)
+
+pointer	list			# List of input spectra
+char	table[ARB]		# Reference table
+int	select			# Selection method
+
+int	i, fd, input, refs
+pointer	stp, sym
+pointer	sp, image, ref1, ref2
+
+pointer	stopen(), strefsbuf(), stenter(), stpstr(), stfind(), imtopen()
+int	imtgetim(), open(), fscan(), nscan()
+errchk	open
+
+begin
+	# Read the table.  Return an error if the file can't be opened.
+	# Read each table entry of spectrum list and reference list.
+	# Expand the input list to make a symbol table keyed on the
+	# spectrum with the reference list string as it's value.
+	# As a convenience if a reference list is missing the preceding
+	# list is implied.
+
+	fd = open (table, READ_ONLY, TEXT_FILE)
+
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (ref1, SZ_FNAME, TY_CHAR)
+	call salloc (ref2, SZ_FNAME, TY_CHAR)
+
+	stp = stopen ("table", 10, 10, 20*SZ_FNAME)
+	while (fscan (fd) != EOF) {
+	    call gargwrd (Memc[image], SZ_FNAME)
+	    call gargwrd (Memc[ref1], SZ_FNAME)
+	    if (nscan() < 1)
+		next
+	    if (nscan() < 2)
+		call strcpy (Memc[ref2], Memc[ref1], SZ_FNAME)
+	    else
+		call strcpy (Memc[ref1], Memc[ref2], SZ_FNAME)
+
+	    i = stpstr (stp, Memc[ref1], SZ_FNAME)
+
+	    input = imtopen (Memc[image])
+	    while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
+		call refnoextn (Memc[image])
+	        sym = stenter (stp, Memc[image], 1)
+	        Memi[sym] = i
+	    }
+	    call imtclose (input)
+	}
+	call close (fd)
+
+	# For each input spectrum find the appropriate reference spectrum list.
+	# If no list is found print a message and continue.  Switch on the
+	# selection method.
+
+	while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+	    call refnoextn (Memc[image])
+	    sym = stfind (stp, Memc[image])
+	    if (sym == NULL) {
+		call refmsgs (NO_REFSPEC, Memc[image], "", "", "", 0, 0, "")
+		next
+	    }
+
+	    input = imtopen (Memc[image])
+	    refs = imtopen (Memc[strefsbuf (stp, Memi[sym])])
+
+	    switch (select) {
+	    case MATCH:
+	        call refmatch(input, refs)
+	    case NEAREST:
+	        call refnearest (input, refs)
+	    case PRECEDING:
+	        call refprecede (input, refs)
+	    case FOLLOWING:
+	        call reffollow (input, refs)
+	    case INTERP:
+	        call refinterp (input, refs)
+	    case AVERAGE:
+	        call refaverage (input, refs)
+	    }
+
+	    call imtclose (input)
+	    call imtclose (refs)
+	}
+
+	call stclose (stp)
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/t_dispcor.x b/noao/onedspec/dispcor/t_dispcor.x
new file mode 100644
index 00000000..94aeba44
--- /dev/null
+++ b/noao/onedspec/dispcor/t_dispcor.x
@@ -0,0 +1,1336 @@
+include	<error.h>
+include	<imhdr.h>
+include	<imio.h>
+include	<mach.h>
+include	<mwset.h>
+include	"dispcor.h"
+include	"dctable.h"
+include	<smw.h>
+include	<units.h>
+ 
+# Dispersion types.
+define	MULTISPEC	1
+define	ECHELLE		2
+
+ 
+# T_DISPCOR -- Dispersion correct spectra.
+ 
+procedure t_dispcor ()
+ 
+int	in			# List of input spectra
+int	out			# List of output spectra
+bool	linearize		# Linearize spectra?
+bool	log			# Log scale?
+bool	flux			# Conserve flux?
+real	blank			# Blank value
+int	ignoreaps		# Ignore aperture numbers?
+int	fd1			# Log file descriptor
+int	fd2			# Log file descriptor
+ 
+int	i, format, naps
+int	open(), nowhite(), imtopenp(), imtgetim(), errcode(), btoi()
+pointer	sp, input, output, str, err, stp, table
+pointer	im, im1, smw, smw1, ap, immap(), smw_openim()
+bool	clgetb()
+real	clgetr()
+errchk	open, immap, smw_openim, dc_gms, dc_gec, dc_multispec, dc_echelle
+ 
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (err, SZ_LINE, TY_CHAR)
+ 
+	# Task parameters
+	in = imtopenp ("input")
+	out = imtopenp ("output")
+	call clgstr ("records", Memc[str], SZ_LINE)
+	call odr_openp (in, Memc[str])
+	call odr_openp (out, Memc[str])
+	call clgstr ("database", Memc[str], SZ_FNAME)
+	call clgstr ("logfile", Memc[err], SZ_LINE)
+	linearize = clgetb ("linearize")
+	ignoreaps = btoi (clgetb ("ignoreaps"))
+ 
+	# Initialize the database cacheing and wavelength table.
+	call dc_open (stp, Memc[str])
+	if (linearize) {
+	    log = clgetb ("log")
+	    flux = clgetb ("flux")
+	    blank = clgetr ("blank")
+
+	    call dc_table (table, naps)
+	    if (clgetb ("global")) {
+		if (clgetb ("samedisp"))
+		    call dc_global1 (in, stp, log, table, naps)
+		else
+		    call dc_global (in, stp, log, table, naps)
+	    }
+	}
+
+	# Open logfile if specified.
+	if (clgetb ("verbose"))
+	    fd1 = STDOUT
+	if (nowhite (Memc[err], Memc[err], SZ_LINE) != 0)
+	    fd2 = open (Memc[err], APPEND, TEXT_FILE)
+	else
+	    fd2 = NULL
+
+	# Loop through each input image.  Do the dispersion correction
+	# in place if no output spectrum list is given or if the input
+	# and output spectra names are the same.
+ 
+	while (imtgetim (in, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (out, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+		
+	    iferr {
+		im = NULL; im1 = NULL
+		smw = NULL; smw1 = NULL
+		ap = NULL
+
+		i = immap (Memc[input], READ_ONLY, 0); im = i
+		i = smw_openim (im); smw = i
+
+		switch (SMW_FORMAT(smw)) {
+		case SMW_ND:
+		    # Use first line for reference.
+		    switch (SMW_LDIM(smw)) {
+		    case 1:
+			call strcpy (Memc[input], Memc[str], SZ_LINE)
+		    case 2:
+			switch (SMW_LAXIS(smw,1)) {
+			case 1:
+			    call sprintf (Memc[str], SZ_LINE, "%s[*,1]")
+				call pargstr (Memc[input])
+			case 2:
+			    call sprintf (Memc[str], SZ_LINE, "%s[1,*]")
+				call pargstr (Memc[input])
+			}
+		    case 3:
+			switch (SMW_LAXIS(smw,1)) {
+			case 1:
+			    call sprintf (Memc[str], SZ_LINE, "%s[*,1,1]")
+				call pargstr (Memc[input])
+			case 2:
+			    call sprintf (Memc[str], SZ_LINE, "%s[1,*,1]")
+				call pargstr (Memc[input])
+			case 3:
+			    call sprintf (Memc[str], SZ_LINE, "%s[*,1,1]")
+				call pargstr (Memc[input])
+			}
+		    }
+		    im1 = immap (Memc[str], READ_ONLY, 0)
+		    smw1 = smw_openim (im1)
+		    call smw_ndes (im1, smw1)
+		    if (SMW_PDIM(smw1) == 1)
+			call smw_esms (smw1)
+
+		    call dc_gms (Memc[input], im1, smw1, stp, YES, ap, fd1, fd2)
+		    call dc_ndspec (im, smw, smw1, ap, Memc[input],
+			Memc[output], linearize, log, flux, blank, table, naps,
+			fd1, fd2)
+		default:
+		    # Get dispersion functions.  Determine type of dispersion
+		    # by the error return.
+
+		    format = MULTISPEC
+		    iferr (call dc_gms (Memc[input], im, smw, stp, ignoreaps,
+			ap, fd1, fd2)) {
+			if (errcode() > 1 && errcode() < 100)
+			    call erract (EA_ERROR)
+			format = ECHELLE
+			iferr (call dc_gec (Memc[input], im, smw, stp, ap,
+			    fd1, fd2)) {
+			    if (errcode() > 1 && errcode() < 100)
+				call erract (EA_ERROR)
+			    call erract (EA_WARN)
+			    iferr (call dc_gms (Memc[input], im, smw, stp,
+				ignoreaps, ap, fd1, fd2))
+				call erract (EA_WARN)
+			    call sprintf (Memc[err], SZ_LINE,
+				"%s: Dispersion data not found")
+				call pargstr (Memc[input])
+			    call error (1, Memc[err])
+			}
+		    }
+
+		    switch (format) {
+		    case MULTISPEC:
+			call dc_multispec (im, smw, ap, Memc[input],
+			    Memc[output], linearize, log, flux, blank, table,
+			    naps, fd1, fd2)
+		    case ECHELLE:
+			call dc_echelle (im, smw, ap, Memc[input],
+			    Memc[output], linearize, log, flux, blank, table,
+			    naps, fd1, fd2)
+		    }
+		}
+	    } then
+		call erract (EA_WARN)
+
+	    if (ap != NULL)
+		call mfree (ap, TY_STRUCT)
+	    if (smw1 != NULL)
+		call smw_close (smw1)
+	    if (im1 != NULL)
+		call imunmap (im1)
+	    if (smw != NULL)
+		call smw_close (smw)
+	    if (im != NULL)
+		call imunmap (im)
+	}
+
+	# Finish up.
+	if (linearize)
+	    do i = 0, naps
+		call mfree (Memi[table+i], TY_STRUCT)
+	call mfree (table, TY_INT)
+	call dc_close (stp)
+	call imtclose (in)
+	call imtclose (out)
+	if (fd1 != NULL)
+	    call close (fd1)
+	if (fd2 != NULL)
+	    call close (fd2)
+	call sfree (sp)
+end
+ 
+ 
+# DC_NDSPEC -- Dispersion correct N-dimensional spectrum.
+ 
+procedure dc_ndspec (in, smw, smw1, ap, input, output, linearize, log, flux,
+	blank, table, naps, fd1, fd2)
+ 
+pointer	in			# Input IMIO pointer
+pointer	smw			# SMW pointer
+pointer	smw1			# SMW pointer
+pointer	ap			# Aperture pointer
+char	input[ARB]		# Input multispec spectrum
+char	output[ARB]		# Output root name
+bool	linearize		# Linearize?
+bool	log			# Log wavelength parameters?
+bool	flux			# Conserve flux?
+real	blank			# Blank value
+pointer	table			# Wavelength table
+int	naps			# Number of apertures
+int	fd1			# Log file descriptor
+int	fd2			# Log file descriptor
+
+int	i, j, nin, ndim, dispaxis, n1, n2, n3
+pointer	sp, temp, str, out, mwout, cti, cto, indata, outdata
+pointer	immap(), imgs3r(), imps3r(), mw_open(), smw_sctran()
+bool	clgetb(), streq()
+errchk	immap, mw_open, smw_open, dispcor, imgs3r, imps3r
+
+begin
+	# Determine the wavelength parameters.
+	call dc_wavelengths (in, ap, output, log, table, naps, 1,
+	    DC_AP(ap,1), DC_W1(ap,1), DC_W2(ap,1), DC_DW(ap,1), DC_NW(ap,1))
+	DC_Z(ap,1) = 0.
+	if (log)
+	    DC_DT(ap,1) = 1
+	else
+	    DC_DT(ap,1) = 0
+	
+	call dc_log (fd1, output, ap, 1, log)
+	call dc_log (fd2, output, ap, 1, log)
+
+	if (clgetb ("listonly"))
+	    return
+
+	call smark (sp)
+	call salloc (temp, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+ 
+	# Open output image.  Use temp. image if output is the same as input.
+	if (streq (input, output)) {
+	    call mktemp ("temp", Memc[temp], SZ_LINE)
+	    out = immap (Memc[temp], NEW_COPY, in)
+	    if (IM_PIXTYPE(out) != TY_DOUBLE)
+		IM_PIXTYPE(out) = TY_REAL
+	} else {
+	    out = immap (output, NEW_COPY, in)
+	    if (IM_PIXTYPE(out) != TY_DOUBLE)
+		IM_PIXTYPE(out) = TY_REAL
+	}
+
+	# Set dimensions.
+	ndim = SMW_LDIM(smw)
+	dispaxis = SMW_LAXIS(smw,1)
+	n1 = DC_NW(ap,1)
+	n2 = SMW_LLEN(smw,2)
+	n3 = SMW_LLEN(smw,3)
+	nin = IM_LEN(in,dispaxis)
+	IM_LEN(out,dispaxis) = n1
+
+	# Set WCS header.
+	mwout = mw_open (NULL, ndim)
+	call mw_newsystem (mwout, "world", ndim)
+	do i = 1, ndim
+	    call mw_swtype (mwout, i, 1, "linear", "")
+	if (UN_LABEL(DC_UN(ap,1)) != EOS)
+	    call mw_swattrs (mwout, dispaxis, "label", UN_LABEL(DC_UN(ap,1)))
+	if (UN_UNITS(DC_UN(ap,1)) != EOS)
+	    call mw_swattrs (mwout, dispaxis, "units", UN_UNITS(DC_UN(ap,1)))
+	call smw_open (mwout, NULL, out)
+	call smw_swattrs (mwout, INDEFI, INDEFI, INDEFI, INDEFI, DC_DT(ap,1),
+	    DC_W1(ap,1), DC_DW(ap,1), DC_NW(ap,1), DC_Z(ap,1), INDEFR, INDEFR,
+	    "")
+
+	# Set WCS transformations.
+	cti = smw_sctran (smw1, "world", "logical", 3)
+	switch (dispaxis) {
+	case 1:
+	    cto = smw_sctran (mwout, "logical", "world", 1)
+	case 2:
+	    cto = smw_sctran (mwout, "logical", "world", 2)
+	case 3:
+	    cto = smw_sctran (mwout, "logical", "world", 4)
+	}
+
+	# Dispersion correct.
+	do j = 1, n3 {
+	    do i = 1, n2 {
+		switch (dispaxis) {
+		case 1:
+		    indata = imgs3r (in, 1, nin, i, i, j, j)
+		    outdata = imps3r (out, 1, n1, i, i, j, j)
+		case 2:
+		    indata = imgs3r (in, i, i, 1, nin, j, j)
+		    outdata = imps3r (out, i, i, 1, n1, j, j)
+		case 3:
+		    indata = imgs3r (in, i, i, j, j, 1, nin)
+		    outdata = imps3r (out, i, i, j, j, 1, n1)
+		}
+
+		call aclrr (Memr[outdata], n1)
+		call dispcora (cti, 1, cto, INDEFI, Memr[indata], nin,
+		    Memr[outdata], n1, flux, blank)
+	    }
+	}
+
+	# Save REFSPEC keywords if present.
+	call dc_refspec (out)
+ 
+	# Finish up.  Replace input by output if needed.
+	call smw_ctfree (cti)
+	call smw_ctfree (cto)
+	call smw_saveim (mwout, out)
+	call smw_close (mwout)
+	call imunmap (out)
+	call imunmap (in)
+	if (streq (input, output)) {
+	    call imdelete (input)
+	    call imrename (Memc[temp], output)
+	}
+ 
+	call sfree (sp)
+end
+ 
+ 
+# DC_MULTISPEC -- Linearize multispec apertures into an MULTISPEC format
+# spectrum.  The number of pixels in each image line is the maximum
+# required to contain the longest spectrum.
+ 
+procedure dc_multispec (in, smw, ap, input, output, linearize, log, flux,
+	blank, table, naps, fd1, fd2)
+ 
+pointer	in			# Input IMIO pointer
+pointer	smw			# SMW pointer
+pointer	ap			# Aperture pointer
+char	input[ARB]		# Input multispec spectrum
+char	output[ARB]		# Output root name
+bool	linearize		# Linearize?
+bool	log			# Log wavelength parameters?
+bool	flux			# Conserve flux?
+real	blank			# Blank value
+pointer	table			# Wavelength table
+int	naps			# Number of apertures
+int	fd1			# Log file descriptor
+int	fd2			# Log file descriptor
+ 
+int	i, j, nc, nl, nb, axis[2]
+pointer	sp, temp, str, out, mwout, cti, cto, indata, outdata
+pointer	immap(), imgl3r(), impl3r()
+pointer	mw_open(), smw_sctran()
+bool	clgetb(), streq()
+errchk	immap, mw_open, smw_open, dispcor, imgl3r, impl3r
+
+data	axis/1,2/
+ 
+begin
+	# Determine the wavelength parameters for each aperture.
+	# The options are to have all apertures have the same dispersion
+	# or have each aperture have independent dispersion.  The global
+	# parameters have already been calculated if needed.
+
+	nc = IM_LEN(in,1)
+	nl = IM_LEN(in,2)
+	nb = IM_LEN(in,3)
+ 
+	if (linearize) {
+	    if (log)
+		DC_DT(ap,1) = 1
+	    else
+		DC_DT(ap,1) = 0
+	    if (clgetb ("samedisp")) {
+		call dc_wavelengths1 (in, smw, ap, output, log, table, naps,
+		    DC_W1(ap,1), DC_W2(ap,1), DC_DW(ap,1), DC_NW(ap,1))
+		if ((DC_DW(ap,1)*(DC_W2(ap,1)-DC_W1(ap,1)) <= 0.) ||
+		    (DC_NW(ap,1) < 1))
+		    call error (1, "Error in wavelength scale")
+		do i = 2, nl {
+		    DC_W1(ap,i) = DC_W1(ap,1)
+		    DC_W2(ap,i) = DC_W2(ap,1)
+		    DC_DW(ap,i) = DC_DW(ap,1)
+		    DC_NW(ap,i) = DC_NW(ap,1)
+		    DC_Z(ap,i) = 0.
+		    DC_DT(ap,i) = DC_DT(ap,1)
+		}
+	    } else {
+		do i = 1, nl {
+		    call dc_wavelengths (in, ap, output, log, table, naps, i,
+			DC_AP(ap,i), DC_W1(ap,i), DC_W2(ap,i), DC_DW(ap,i),
+			DC_NW(ap,i))
+		    DC_Z(ap,i) = 0.
+		    DC_DT(ap,i) = DC_DT(ap,1)
+		}
+	    }
+	}
+	call dc_log (fd1, output, ap, nl, log)
+	call dc_log (fd2, output, ap, nl, log)
+
+	if (clgetb ("listonly"))
+	    return
+
+	call smark (sp)
+	call salloc (temp, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+ 
+	# Use a temporary image if the output has the same name as the input.
+	if (streq (input, output)) {
+	    if (linearize) {
+		call mktemp ("temp", Memc[temp], SZ_LINE)
+		out = immap (Memc[temp], NEW_COPY, in)
+		if (IM_PIXTYPE(out) != TY_DOUBLE)
+		    IM_PIXTYPE(out) = TY_REAL
+	    } else {
+		call imunmap (in)
+		i = immap (input, READ_WRITE, 0)
+		in = i
+		out = i
+	    }
+	} else {
+	    out = immap (output, NEW_COPY, in)
+	    if (IM_PIXTYPE(out) != TY_DOUBLE)
+		IM_PIXTYPE(out) = TY_REAL
+	}
+
+	# Set MWCS or linearize
+	if (!linearize) {
+	    if (out != in)
+		do j = 1, nb
+		    do i = 1, nl
+			call amovr (Memr[imgl3r(in,i,j)], Memr[impl3r(out,i,j)],
+			    IM_LEN(in,1))
+	    call smw_saveim (smw, out)
+	} else {
+	    if (nb > 1)
+		i = 3
+	    else
+		i = 2
+	    mwout = mw_open (NULL, i)
+	    call mw_newsystem (mwout, "multispec", i)
+	    call mw_swtype (mwout, axis, 2, "multispec", "")
+	    if (UN_LABEL(DC_UN(ap,1)) != EOS)
+		call mw_swattrs (mwout, 1, "label", UN_LABEL(DC_UN(ap,1)))
+	    if (UN_UNITS(DC_UN(ap,1)) != EOS)
+		call mw_swattrs (mwout, 1, "units", UN_UNITS(DC_UN(ap,1)))
+	    if (i == 3)
+		call mw_swtype (mwout, 3, 1, "linear", "")
+	    call smw_open (mwout, NULL, out)
+	    do i = 1, nl {
+		call smw_swattrs (mwout, i, 1, DC_AP(ap,i), DC_BM(ap,i),
+		    DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i),
+		    DC_Z(ap,i), DC_LW(ap,i), DC_UP(ap,i), "")
+		call smw_gapid (smw, i, 1, Memc[str], SZ_LINE)
+		call smw_sapid (mwout, i, 1, Memc[str])
+	    }
+
+	    IM_LEN(out,1) = DC_NW(ap,1)
+	    do i = 2, nl
+		IM_LEN(out,1) = max (DC_NW(ap,i), IM_LEN(out,1))
+	    cti = smw_sctran (smw, "world", "logical", 3)
+	    cto = smw_sctran (mwout, "logical", "world", 3)
+	    do j = 1, nb {
+		do i = 1, nl {
+		    indata = imgl3r (in, i, j)
+		    outdata = impl3r (out, i, j)
+		    call aclrr (Memr[outdata], IM_LEN(out,1))
+		    call dispcora (cti, i, cto, i, Memr[indata], nc,
+			Memr[outdata], DC_NW(ap,i), flux, blank)
+		    if (DC_NW(ap,i) < IM_LEN(out,1))
+			call amovkr (Memr[outdata+DC_NW(ap,i)-1],
+			    Memr[outdata+DC_NW(ap,i)],IM_LEN(out,1)-DC_NW(ap,i))
+		}
+	    }
+	    call smw_ctfree (cti)
+	    call smw_ctfree (cto)
+	    call smw_saveim (mwout, out)
+	    call smw_close (mwout)
+	}
+
+	# Save REFSPEC keywords if present.
+	call dc_refspec (out)
+ 
+	# Finish up.  Replace input by output if needed.
+	if (out == in) {
+	    call imunmap (in)
+	} else {
+	    call imunmap (in)
+	    call imunmap (out)
+	    if (streq (input, output)) {
+		call imdelete (input)
+		call imrename (Memc[temp], output)
+	    }
+	}
+ 
+	call sfree (sp)
+end
+ 
+ 
+# DC_ECHELLE -- Linearize echelle orders into an ECHELLE format
+# spectrum.  The number of pixels in each image line is the maximum
+# required to contain the longest spectrum.
+ 
+procedure dc_echelle (in, smw, ap, input, output, linearize, log, flux,
+	blank, table, naps, fd1, fd2)
+ 
+pointer	in			# IMIO pointer
+pointer	smw			# SMW pointers
+pointer	ap			# Aperture pointer
+char	input[ARB]		# Input multispec spectrum
+char	output[ARB]		# Output root name
+bool	linearize		# Linearize?
+bool	log			# Log wavelength parameters?
+bool	flux			# Conserve flux?
+real	blank			# Blank value
+pointer	table			# Wavelength table
+int	naps			# Number of apertures
+int	fd1			# Log file descriptor
+int	fd2			# Log file descriptor
+ 
+int	i, j, nc, nl, nb, axis[2]
+pointer	sp, temp, str, out, mwout, cti, cto, indata, outdata
+pointer	immap(), imgl3r(), impl3r()
+pointer	mw_open(), smw_sctran()
+bool	clgetb(), streq()
+errchk	immap, mw_open, smw_open, dispcor, imgl3r, impl3r
+
+data	axis/1,2/
+ 
+begin
+	# Determine the wavelength parameters for each aperture.
+ 
+	nc = IM_LEN(in,1)
+	nl = IM_LEN(in,2)
+	nb = IM_LEN(in,3)
+ 
+	if (linearize) {
+	    if (log)
+		DC_DT(ap,1) = 1
+	    else
+		DC_DT(ap,1) = 0
+	    do i = 1, nl {
+		call dc_wavelengths (in, ap, output, log, table, naps,
+		    i, DC_AP(ap,i), DC_W1(ap,i), DC_W2(ap,i), DC_DW(ap,i),
+		    DC_NW(ap,i))
+		DC_Z(ap,i) = 0.
+		DC_DT(ap,i) = DC_DT(ap,1)
+	    }
+	}
+	call dc_log (fd1, output, ap, nl, log)
+	call dc_log (fd2, output, ap, nl, log)
+
+	if (clgetb ("listonly"))
+	    return
+
+	call smark (sp)
+	call salloc (temp, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Use a temporary image if the output has the same name as the input.
+	if (streq (input, output)) {
+	    if (linearize) {
+		call mktemp ("temp", Memc[temp], SZ_LINE)
+		out = immap (Memc[temp], NEW_COPY, in)
+		if (IM_PIXTYPE(out) != TY_DOUBLE)
+		    IM_PIXTYPE(out) = TY_REAL
+	    } else {
+		call imunmap (in)
+		i = immap (input, READ_WRITE, 0)
+		in = i
+		out = i
+	    }
+	} else {
+	    out = immap (output, NEW_COPY, in)
+	    if (IM_PIXTYPE(out) != TY_DOUBLE)
+		IM_PIXTYPE(out) = TY_REAL
+	}
+
+	# Set MWCS or linearize
+	if (!linearize) {
+	    if (out != in)
+		do j = 1, nb
+		    do i = 1, nl
+			call amovr (Memr[imgl3r(in,i,j)], Memr[impl3r(out,i,j)],
+			    IM_LEN(in,1))
+	    call smw_saveim (smw, out)
+	} else {
+	    if (nb > 1)
+		i = 3
+	    else
+		i = 2
+	    mwout = mw_open (NULL, i)
+	    call mw_newsystem (mwout, "multispec", i)
+	    call mw_swtype (mwout, axis, 2, "multispec", "")
+	    if (UN_LABEL(DC_UN(ap,1)) != EOS)
+		call mw_swattrs (mwout, 1, "label", UN_LABEL(DC_UN(ap,1)))
+	    if (UN_UNITS(DC_UN(ap,1)) != EOS)
+		call mw_swattrs (mwout, 1, "units", UN_UNITS(DC_UN(ap,1)))
+	    if (i == 3)
+		call mw_swtype (mwout, 3, 1, "linear", "")
+	    call smw_open (mwout, NULL, out)
+	    do i = 1, nl {
+		call smw_swattrs (mwout, i, 1, DC_AP(ap,i), DC_BM(ap,i),
+		    DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i),
+		    DC_Z(ap,i), DC_LW(ap,i), DC_UP(ap,i), "")
+		call smw_gapid (smw, i, 1, Memc[str], SZ_LINE)
+		call smw_sapid (mwout, i, 1, Memc[str])
+	    }
+
+	    IM_LEN(out,1) = DC_NW(ap,1)
+	    do i = 2, nl
+		IM_LEN(out,1) = max (DC_NW(ap,i), IM_LEN(out,1))
+	    cti = smw_sctran (smw, "world", "logical", 3)
+	    cto = smw_sctran (mwout, "logical", "world", 3)
+	    do j = 1, nb {
+		do i = 1, nl {
+		    indata = imgl3r (in, i, j)
+		    outdata = impl3r (out, i, j)
+		    call aclrr (Memr[outdata], IM_LEN(out,1))
+		    call dispcora (cti, i, cto, i, Memr[indata], nc,
+			Memr[outdata], DC_NW(ap,i), flux, blank)
+		    if (DC_NW(ap,i) < IM_LEN(out,1))
+			call amovkr (Memr[outdata+DC_NW(ap,i)-1],
+			    Memr[outdata+DC_NW(ap,i)],IM_LEN(out,1)-DC_NW(ap,i))
+		}
+	    }
+	    call smw_ctfree (cti)
+	    call smw_ctfree (cto)
+	    call smw_saveim (mwout, out)
+	    call smw_close (mwout)
+	}
+
+	# Save REFSPEC keywords if present.
+	call dc_refspec (out)
+ 
+	# Finish up.  Replace input by output if needed.
+	if (out == in) {
+	    call imunmap (in)
+	} else {
+	    call imunmap (in)
+	    call imunmap (out)
+	    if (streq (input, output)) {
+		call imdelete (input)
+		call imrename (Memc[temp], output)
+	    }
+	}
+ 
+	call sfree (sp)
+end
+ 
+ 
+# DC_GLOBAL1 -- Set global wavelength parameters using the minimum and
+# maximum wavelengths and and the minimum dispersion over all apertures.
+ 
+procedure dc_global1 (in, stp, log, table, naps)
+ 
+pointer	in			# Input list
+pointer	stp			# Symbol table
+bool	log			# Logarithmic scale?
+pointer	table			# Wavelength table
+int	naps			# Number of apertures
+ 
+int	i, nwmax, imtgetim()
+double	w1, w2, dw, wmin, wmax, dwmin
+pointer	sp, input, str, im, mw, ap, tbl, immap(), smw_openim()
+errchk	dc_gms, dc_gec, smw_openim
+ 
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+ 
+	# Go through all the reference spectra and determine the
+	# minimum and maximum wavelengths and maximum number of pixels.
+	# If there is no entry in the wavelength table add it.
+
+	wmin = MAX_REAL
+	wmax = -MAX_REAL
+	dwmin = MAX_REAL
+
+	while (imtgetim (in, Memc[input], SZ_FNAME) != EOF) {
+	    iferr (im = immap (Memc[input], READ_ONLY, 0))
+		next
+	    mw = smw_openim (im)
+	    switch (SMW_FORMAT(mw)) {
+	    case SMW_ND:
+		nwmax = SMW_NW(mw)
+		dw = SMW_DW(mw)
+		w1 = SMW_W1(mw)
+		w2 = w1 + dw * (nwmax - 1)
+		wmin = min (wmin, w1, w2)
+		wmax = max (wmax, w1, w2)
+		dwmin = min (dwmin, abs (dw))
+	    default:
+		iferr {
+		    iferr (call dc_gms (Memc[input], im, mw, stp, NO, ap,
+			NULL, NULL)) {
+			iferr (call dc_gec (Memc[input], im, mw, stp, ap,
+			    NULL, NULL)) {
+			    call sprintf (Memc[str], SZ_LINE,
+				"%s: Dispersion data not found")
+				call pargstr (Memc[input])
+			    call error (1, Memc[str])
+			}
+		    }
+     
+		    do i = 1, IM_LEN(im,2) {
+			w1 = DC_W1(ap,i)
+			w2 = DC_W2(ap,i)
+			dw = DC_DW(ap,i)
+			wmin = min (wmin, w1, w2)
+			wmax = max (wmax, w1, w2)
+			dwmin = min (dwmin, abs (dw))
+		    }
+		} then
+		    ;
+	    }
+
+	    call mfree (ap, TY_STRUCT)
+	    call smw_close (mw)
+	    call imunmap (im)
+	}
+	call imtrew (in)
+
+	nwmax = (wmax - wmin) / dwmin + 1.5
+ 
+	# Enter the global entry in the first table entry.
+	tbl = Memi[table]
+	call dc_defaults (wmin, wmax, nwmax,
+	    TBL_W1(tbl), TBL_W2(tbl), TBL_DW(tbl), TBL_NW(tbl))
+ 
+	call sfree (sp)
+end
+ 
+ 
+# DC_GLOBAL -- Set global wavelength parameters.  This is done for each
+# aperture separately.  The wavelength table may be used to specify separate
+# fixed parameters for each aperture.
+ 
+procedure dc_global (in, stp, log, table, naps)
+ 
+pointer	in			# Input list
+pointer	stp			# Symbol table
+bool	log			# Logarithmic scale?
+pointer	table			# Wavelength table
+int	naps			# Number of apertures
+ 
+int	i, j, nw, imtgetim()
+double	w1, w2, dw
+pointer	sp, input, str, im, mw, ap, tbl, immap(), smw_openim()
+errchk	dc_gms, dc_gec, smw_openim
+ 
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+ 
+	# Go through all the reference spectra and determine the
+	# minimum and maximum wavelengths and maximum number of pixels.
+	# Do this by aperture.  If there is no entry in the wavelength
+	# table add it.
+ 
+	while (imtgetim (in, Memc[input], SZ_FNAME) != EOF) {
+	    iferr (im = immap (Memc[input], READ_ONLY, 0))
+		next
+	    mw = smw_openim (im)
+	    switch (SMW_FORMAT(mw)) {
+	    case SMW_ND:
+		tbl = Memi[table]
+		nw = SMW_NW(mw)
+		dw = SMW_DW(mw)
+		w1 = SMW_W1(mw)
+		w2 = w1 + dw * (nw - 1)
+		TBL_WMIN(tbl) = min (TBL_WMIN(tbl), w1, w2)
+		TBL_WMAX(tbl) = max (TBL_WMAX(tbl), w1, w2)
+		TBL_NWMAX(tbl) = max (TBL_NWMAX(tbl), nw)
+	    default:
+		iferr {
+		    iferr (call dc_gms (Memc[input], im, mw, stp, NO, ap,
+			NULL, NULL)) {
+			iferr (call dc_gec (Memc[input], im, mw, stp, ap,
+			    NULL, NULL)) {
+			    call sprintf (Memc[str], SZ_LINE,
+				"%s: Dispersion data not found")
+				call pargstr (Memc[input])
+			    call error (1, Memc[str])
+			}
+		    }
+     
+		    do i = 1, IM_LEN(im,2) {
+			call dc_getentry (false, DC_AP(ap,i), table, naps, j)
+			tbl = Memi[table+j]
+	 
+			nw = DC_NW(ap,i)
+			w1 = DC_W1(ap,i)
+			w2 = DC_W2(ap,i)
+			TBL_WMIN(tbl) = min (TBL_WMIN(tbl), w1, w2)
+			TBL_WMAX(tbl) = max (TBL_WMAX(tbl), w1, w2)
+			TBL_NWMAX(tbl) = max (TBL_NWMAX(tbl), nw)
+		    }
+		} then
+		    ;
+	    }
+
+	    call mfree (ap, TY_STRUCT)
+	    call smw_close (mw)
+	    call imunmap (im)
+	}
+	call imtrew (in)
+ 
+	do i = 0, naps {
+	    tbl = Memi[table+i]
+	    call dc_defaults (TBL_WMIN(tbl), TBL_WMAX(tbl), TBL_NWMAX(tbl),
+		TBL_W1(tbl), TBL_W2(tbl), TBL_DW(tbl), TBL_NW(tbl))
+	}
+
+	call sfree (sp)
+end
+ 
+
+# DC_WAVELENGTHS1 -- Set output wavelength parameters for a spectrum.
+# Fill in any INDEF values using the limits of the dispersion function
+# over all apertures and the minimum dispersion over all apertures.  The
+# user may then confirm and change the wavelength parameters if desired.
+ 
+procedure dc_wavelengths1 (im, smw, ap, output, log, table, naps, w1, w2, dw,nw)
+ 
+pointer	im		# IMIO pointer
+pointer	smw		# SMW pointer
+pointer	ap		# Aperture structure
+char	output[ARB]	# Output image name
+bool	log		# Logarithm wavelength parameters?
+pointer	table		# Wavelength table
+int	naps		# Number of apertures
+double	w1, w2, dw	# Image wavelength parameters
+int	nw		# Image wavelength parameter
+ 
+int	i, n, nwt, clgeti(), clgwrd()
+double	a, b, c, w1t, w2t, dwt, y1, y2, dy, clgetd()
+pointer	sp, key, str, tbl
+bool	clgetb()
+ 
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+ 
+	# Get aperture parameters.
+	tbl = Memi[table]
+	w1t = TBL_W1(tbl)
+	w2t = TBL_W2(tbl)
+	dwt = TBL_DW(tbl)
+	nwt = TBL_NW(tbl)
+ 
+	# If there are undefined wavelength scale parameters get
+	# defaults based on the reference spectrum.
+ 
+	if (IS_INDEFD(w1t)||IS_INDEFD(w2t)||IS_INDEFD(dwt)||IS_INDEFD(nwt)) {
+	    a = MAX_REAL
+	    b = -MAX_REAL
+	    c = MAX_REAL
+ 
+	    do i = 1, IM_LEN(im,2) {
+		n = DC_NW(ap,i)
+		y1 = DC_W1(ap,i)
+		y2 = DC_W2(ap,i)
+		dy = DC_DW(ap,i)
+		a = min (a, y1, y2)
+		b = max (b, y1, y2)
+		c = min (c, dy)
+	    }
+	    n = (b - a) / c + 1.5
+	}
+ 
+	call dc_defaults (a, b, n, w1t, w2t, dwt, nwt)
+	w1 = w1t
+	w2 = w2t
+	dw = dwt
+	nw = nwt
+ 
+	# Print the wavelength scale and allow the user to confirm and
+	# change the wavelength scale.  A test is done to check which
+	# parameters the user changes and give them priority in filling
+	# in the remaining parameters.
+ 
+	if (TBL_CONFIRM(tbl) == YES) {
+	    repeat {
+		call printf ("%s: w1 = %g, w2 = %g, dw = %g, nw = %d\n")
+		    call pargstr (output)
+		    call pargd (w1)
+		    call pargd (w2)
+		    call pargd (dw)
+		    call pargi (nw)
+ 
+		i = clgwrd ("dispcor1.change", Memc[str],SZ_LINE, "|yes|no|NO|")
+		switch (i) {
+		case 2:
+		    break
+		case 3:
+		    TBL_CONFIRM(tbl) = NO
+		    break
+		}
+		call clputd ("dispcor1.w1", w1)
+		call clputd ("dispcor1.w2", w2)
+		call clputd ("dispcor1.dw", dw)
+		call clputi ("dispcor1.nw", nw)
+		a = w1
+		b = w2
+		c = dw
+		n = nw
+		w1 = clgetd ("dispcor1.w1")
+		w2 = clgetd ("dispcor1.w2")
+		dw = clgetd ("dispcor1.dw")
+		nw = clgeti ("dispcor1.nw")
+ 
+		# If no INDEF's set unchanged parameters to INDEF.
+		i = 0
+		if (IS_INDEFD(w1))
+		    i = i + 1
+		if (IS_INDEFD(w2))
+		    i = i + 1
+		if (IS_INDEFD(dw))
+		    i = i + 1
+		if (IS_INDEFI(nw))
+		    i = i + 1
+		if (i == 0) {
+		    if (w1 == a)
+		        w1 = INDEFD
+		    if (w2 == b)
+		        w2 = INDEFD
+		    if (dw == c)
+		        dw = INDEFD
+		    if (nw == n)
+		        nw = INDEFI
+		}
+ 
+	        call dc_defaults (a, b, n, w1, w2, dw, nw)
+ 
+		if (clgetb ("global")) {
+		    TBL_W1(tbl) = w1
+		    TBL_W2(tbl) = w2
+		    TBL_DW(tbl) = dw
+		    TBL_NW(tbl) = nw
+		}
+	    }
+	}
+	call sfree (sp)
+end
+ 
+ 
+# DC_WAVELENGTHS -- Set output wavelength parameters for a spectrum for
+# each aperture.  The fixed parameters are given in the wavelength table.
+# If there is no entry in the table for an aperture use the global
+# default (entry 0).  Fill in INDEF values using the limits and number
+# of pixels for the aperture.  The user may then confirm and change the
+# wavelength parameters if desired.
+ 
+procedure dc_wavelengths (im, ap, output, log, table, naps, line, apnum,
+	w1, w2, dw, nw)
+ 
+pointer	im		# IMIO pointer
+pointer	ap		# Aperture structure
+char	output[ARB]	# Output image name
+bool	log		# Logarithm wavelength parameters?
+pointer	table		# Wavelength table
+int	naps		# Number of apertures
+int	line		# Line
+int	apnum		# Aperture number
+double	w1, w2, dw	# Image wavelength parameters
+int	nw		# Image wavelength parameter
+ 
+int	i, n, nwt, clgeti(), clgwrd()
+double	a, b, c, w1t, w2t, dwt, clgetd()
+pointer	sp, str, tbl
+bool	clgetb()
+ 
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+ 
+	# Get aperture parameters.
+	call dc_getentry (false, apnum, table, naps, i)
+	tbl = Memi[table+i]
+ 
+	w1t = TBL_W1(tbl)
+	w2t = TBL_W2(tbl)
+	dwt = TBL_DW(tbl)
+	nwt = TBL_NW(tbl)
+ 
+	# If there are undefined wavelength scale parameters get
+	# defaults based on the reference spectrum.
+ 
+	if (IS_INDEFD(w1t)||IS_INDEFD(w2t)||IS_INDEFD(dwt)||IS_INDEFI(nwt)) {
+	    a = DC_W1(ap,line)
+	    b = DC_W2(ap,line)
+	    n = DC_NW(ap,line)
+	}
+ 
+	call dc_defaults (a, b, n, w1t, w2t, dwt, nwt)
+	w1 = w1t
+	w2 = w2t
+	dw = dwt
+	nw = nwt
+ 
+	# Print the wavelength scale and allow the user to confirm and
+	# change the wavelength scale.  A test is done to check which
+	# parameters the user changes and give them priority in filling
+	# in the remaining parameters.
+ 
+	if (TBL_CONFIRM(tbl) == YES) {
+	    repeat {
+		call printf (
+		    "%s: ap = %d, w1 = %g, w2 = %g, dw = %g, nw = %d\n")
+		    call pargstr (output)
+		    call pargi (apnum)
+		    call pargd (w1)
+		    call pargd (w2)
+		    call pargd (dw)
+		    call pargi (nw)
+		i = clgwrd ("dispcor1.change", Memc[str],SZ_LINE, "|yes|no|NO|")
+		switch (i) {
+		case 2:
+		    break
+		case 3:
+		    TBL_CONFIRM(tbl) = NO
+		    break
+		}
+		call clputd ("dispcor1.w1", w1)
+		call clputd ("dispcor1.w2", w2)
+		call clputd ("dispcor1.dw", dw)
+		call clputi ("dispcor1.nw", nw)
+		a = w1
+		b = w2
+		c = dw
+		n = nw
+		w1 = clgetd ("dispcor1.w1")
+		w2 = clgetd ("dispcor1.w2")
+		dw = clgetd ("dispcor1.dw")
+		nw = clgeti ("dispcor1.nw")
+ 
+		# If no INDEF's set unchanged parameters to INDEF.
+		i = 0
+		if (IS_INDEFD(w1))
+		    i = i + 1
+		if (IS_INDEFD(w2))
+		    i = i + 1
+		if (IS_INDEFD(dw))
+		    i = i + 1
+		if (IS_INDEFI(nw))
+		    i = i + 1
+		if (i == 0) {
+		    if (w1 == a)
+		        w1 = INDEFD
+		    if (w2 == b)
+		        w2 = INDEFD
+		    if (dw == c)
+		        dw = INDEFD
+		    if (nw == n)
+		        nw = INDEFI
+		}
+ 
+	        call dc_defaults (a, b, n, w1, w2, dw, nw)
+ 
+		if (clgetb ("global")) {
+		    TBL_W1(tbl) = w1
+		    TBL_W2(tbl) = w2
+		    TBL_DW(tbl) = dw
+		    TBL_NW(tbl) = nw
+		}
+	    }
+	}
+	call sfree (sp)
+end
+
+
+# DC_DEFAULTS -- Given some set of wavelength scale with others undefined
+# (INDEF) plus some defaults fill in the undefined parameters and make
+# the wavelength scale consistent.  The logic of this task is complex
+# and is meant to provide an "intelligent" result based on what users
+# want.
+
+procedure dc_defaults (a, b, n, w1, w2, dw, nw)
+
+double	a		# Default wavelength endpoint
+double	b		# Default wavelength endpoint
+int	n		# Default number of pixels
+double	w1		# Starting wavelength
+double	w2		# Ending wavelength
+double	dw		# Wavelength interval
+int	nw		# Number of pixels
+
+int	nindef
+
+begin
+	# Determine how many input parameters are specfied.
+	nindef = 0
+	if (IS_INDEFD(w1))
+	    nindef = nindef + 1
+	if (IS_INDEFD(w2))
+	    nindef = nindef + 1
+	if (IS_INDEFD(dw))
+	    nindef = nindef + 1
+	if (IS_INDEFI(nw))
+	    nindef = nindef + 1
+
+	# Depending on how many parameters are specified fill in the
+	# INDEF parameters.
+
+	switch (nindef) {
+	case 0:
+	    # All parameters specified.  First round NW to be consistent with
+	    # w1, w2, and dw.  Then adjust w2 to nearest pixel.  It is possible
+	    # that nw will be negative.  Checks for this should be made by the
+	    # call in program.
+
+	    nw = (w2 - w1) / dw + 1.5
+	    w2 = w1 + dw * (nw - 1)
+	case 1:
+	    # Find the unspecified parameter and compute it from the other
+	    # three specified parameters.  For nw need to adjust w2 to
+	    # agree with a pixel.
+
+	    if (IS_INDEFD(w1))
+		w1 = w2 - dw * (nw - 1)
+	    if (IS_INDEFD(w2))
+		w2 = w1 + dw * (nw - 1)
+	    if (IS_INDEFD(dw))
+		dw = (w2 - w1) / (nw - 1)
+	    if (IS_INDEFI(nw)) {
+	        nw = (w2 - w1) / dw + 1.5
+		w2 = w1 + dw * (nw - 1)
+	    }
+	case 2:
+	    # Fill in two unspecified parameters using the defaults.
+	    # This is tricky.
+
+	    if (IS_INDEFD(dw)) {
+		if (IS_INDEFD(w1)) {
+		    if (abs (w2 - a) > abs (w2 - b))
+			w1 = a
+		    else
+			w1 = b
+		    dw = (w2 - w1) / (nw - 1)
+		} else if (IS_INDEFD(w2)) {
+		    if (abs (w1 - a) > abs (w1 - b))
+			w2 = a
+		    else
+			w2 = b
+		    dw = (w2 - w1) / (nw - 1)
+		} else {
+		    dw = (b - a) / n
+		    nw = abs ((w2 - w1) / dw) + 1.5
+		    dw = (w2 - w1) / (nw - 1)
+		}
+	    } else if (IS_INDEFI(nw)) {
+		if (IS_INDEFD(w1)) {
+		    if (dw > 0.)
+			w1 = min (a, b)
+		    else
+			w1 = max (a, b)
+		    nw = (w2 - w1) / dw + 1.5
+		    w1 = w2 - dw * (nw - 1)
+		} else {
+		    if (dw > 0.)
+			w2 = max (a, b)
+		    else
+			w2 = min (a, b)
+		    nw = (w2 - w1) / dw + 1.5
+		    w2 = w1 + dw * (nw - 1)
+		}
+	    } else {
+		if (dw > 0.)
+		    w1 = min (a, b)
+		else
+		    w1 = max (a, b)
+		w2 = w1 + dw * (nw - 1)
+	    }
+	case 3:
+	    # Find the one specfied parameter and compute the others using
+	    # the supplied defaults.
+
+	    if (!IS_INDEFD(w1)) {
+		if (abs (w1 - a) > abs (w1 - b))
+		    w2 = a
+		else
+		    w2 = b
+		dw = (b - a) / n
+		nw = abs ((w2 - w1) / dw) + 1.5
+		dw = (w2 - w1) / (nw - 1)
+	    } else if (!IS_INDEFD(w2)) {
+		if (abs (w2 - a) > abs (w2 - b))
+		    w1 = a
+		else
+		    w1 = b
+		dw = (b - a) / n
+		nw = abs ((w2 - w1) / dw) + 1.5
+		dw = (w2 - w1) / (nw - 1)
+	    } else if (!IS_INDEFI(nw)) {
+		w1 = min (a, b)
+		w2 = max (a, b)
+	        dw = (w2 - w1) / (nw - 1)
+	    } else if (dw < 0.) {
+		w1 = max (a, b)
+		w2 = min (a, b)
+		nw = (w2 - w1) / dw + 1.5
+		w2 = w1 + dw * (nw - 1)
+	    } else {
+		w1 = min (a, b)
+		w2 = max (a, b)
+		nw = (w2 - w1) / dw + 1.5
+		w2 = w1 + dw * (nw - 1)
+	    }
+	case 4:
+	    # Given only defaults compute a wavelength scale.  The dispersion
+	    # is kept close to the default.
+	    w1 = min (a, b)
+	    w2 = max (a, b)
+	    dw = (b - a) / (n - 1)
+	    nw = abs ((w2 - w1) / dw) + 1.5
+	    dw = (w2 - w1) / (nw - 1)
+	}
+end
+
+
+# DC_LOG -- Print log of wavlength paramters
+
+procedure dc_log (fd, output, ap, naps, log)
+
+int	fd		# Output file descriptor
+char	output[ARB]	# Output image name
+pointer	ap		# Aperture structure
+int	naps		# Number of apertures
+bool	log		# Log dispersion?
+
+int	i
+
+begin
+	if (fd == NULL)
+	    return
+
+	for (i=2; i<=naps; i=i+1) {
+	    if (DC_W1(ap,i) != DC_W1(ap,1))
+		break
+	    if (DC_W2(ap,i) != DC_W2(ap,1))
+		break
+	    if (DC_DW(ap,i) != DC_DW(ap,1))
+		break
+	    if (DC_NW(ap,i) != DC_NW(ap,1))
+		break
+	}
+
+	if (naps == 1 || i <= naps) {
+	    do i = 1, naps {
+		call fprintf (fd,
+		    "%s: ap = %d, w1 = %8g, w2 = %8g, dw = %8g, nw = %d")
+		    call pargstr (output)
+		    call pargi (DC_AP(ap,i))
+		    call pargd (DC_W1(ap,i))
+		    call pargd (DC_W2(ap,i))
+		    call pargd (DC_DW(ap,i))
+		    call pargi (DC_NW(ap,i))
+		if (log) {
+		    call fprintf (fd, ", log = %b")
+			call pargb (log)
+		}
+		call fprintf (fd, "\n")
+	    }
+	} else {
+	    call fprintf (fd,
+		"%s: w1 = %8g, w2 = %8g, dw = %8g, nw = %d")
+		call pargstr (output)
+		call pargd (DC_W1(ap,1))
+		call pargd (DC_W2(ap,1))
+		call pargd (DC_DW(ap,1))
+		call pargi (DC_NW(ap,1))
+	    if (log) {
+		call fprintf (fd, ", log = %b")
+		    call pargb (log)
+	    }
+	    call fprintf (fd, "\n")
+	}
+	call flush (fd)
+end
+
+
+# DC_REFSPEC -- Save REFSPEC keywords in DCLOG keywords.
+
+procedure dc_refspec (im)
+
+pointer	im			#U IMIO pointer
+
+int	i, j, imaccf()
+pointer	sp, dckey, dcstr, refkey, refstr
+
+begin
+	call smark (sp)
+	call salloc (dckey, SZ_FNAME, TY_CHAR)
+	call salloc (dcstr, SZ_LINE, TY_CHAR)
+	call salloc (refkey, SZ_FNAME, TY_CHAR)
+	call salloc (refstr, SZ_LINE, TY_CHAR)
+
+	for (i=1;; i=i+1) {
+	    call sprintf (Memc[dckey], SZ_FNAME, "DCLOG%d")
+		call pargi (i)
+	    if (imaccf (im, Memc[dckey]) == NO)
+		break
+	}
+
+	do j = 1, 4 {
+	    if (j == 1)
+		call strcpy ("REFSPEC1", Memc[refkey], SZ_FNAME)
+	    else if (j == 2)
+		call strcpy ("REFSPEC2", Memc[refkey], SZ_FNAME)
+	    else if (j == 3)
+		call strcpy ("REFSHFT1", Memc[refkey], SZ_FNAME)
+	    else if (j == 4)
+		call strcpy ("REFSHFT2", Memc[refkey], SZ_FNAME)
+
+	    ifnoerr (call imgstr (im, Memc[refkey], Memc[refstr], SZ_LINE)) {
+		call sprintf (Memc[dckey], SZ_FNAME, "DCLOG%d")
+		    call pargi (i)
+		call sprintf (Memc[dcstr], SZ_LINE, "%s = %s")
+		    call pargstr (Memc[refkey])
+		    call pargstr (Memc[refstr])
+		call imastr (im, Memc[dckey], Memc[dcstr])
+		call imdelf (im, Memc[refkey])
+		i = i + 1
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/dispcor/t_disptrans.x b/noao/onedspec/dispcor/t_disptrans.x
new file mode 100644
index 00000000..ee108472
--- /dev/null
+++ b/noao/onedspec/dispcor/t_disptrans.x
@@ -0,0 +1,413 @@
+include	<error.h>
+include	<imhdr.h>
+include	<math/curfit.h>
+include	<smw.h>
+include	<units.h>
+
+define	AIRVAC	"|none|air2vac|vac2air|"
+define	NONE	1		# No correction
+define	AIR2VAC	2		# Correct air to vacuum
+define	VAC2AIR	3		# Correct vacuum to air
+
+
+# T_DISPTRANS -- Tranform dispersion systems and apply air-vac conversion.
+# This task uses the UNITS package to convert the input dispersion
+# coordinates to the desired output coordinates.  An air to vacuum or
+# vacuum to air correction is made.  Since the input and output units
+# may not be linearly related and the MWCS supports only polynomial
+# representations a cubic splines are fit to the desired output coordinates
+# until an error tolerance is reached.  The user may then select to
+# store the new WCS as either the spline approximation or to linearize
+# the coordinates by resampling the data.  Note that if the input and
+# output units ARE linearly related and there is no air/vacuum conversion
+# then linearization or storing of a nonlinear dispersion function is
+# skipped.  The operations are done in double precision.
+
+procedure t_disptrans ()
+
+int	inlist			# List of input spectra
+int	outlist			# List of output spectra
+pointer	units			# Output dispersion units
+double	maxerr			# Maximum error (in pixels)
+bool	linearize		# Linearize ouptut dispersion?
+bool	verbose			# Verbose?
+
+int	air			# Air-vacuum conversion?
+double	t			# Temperture in degrees C
+double	p			# Pressure in mmHg
+double	f			# Water vapour pressure in mmHg
+
+int	i, j, n, nw, format, dtype, dtype1, axis[2]
+double	err, w1, dw
+pointer	ptr, in, out, mwin, mwout, ctin, ctout, sh, cv, inbuf, outbuf
+pointer	sp, input, output, title, coeff, x, y, w, nx
+
+bool	clgetb(), streq()
+int	clgwrd(), imtopenp(), imtgetim()
+double	clgetd(), shdr_lw(), dcveval()
+pointer	immap(), smw_openim(), smw_sctran(), mw_open(), imgl3r(), impl3r()
+errchk	immap, impl3r
+errchk	smw_openim, smw_gwattrs, shdr_open, mw_open
+errchk	dt_airvac, dt_cvfit, dt_setwcs, dispcor
+
+data	axis/1,2/
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (units, SZ_FNAME, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	coeff = NULL
+
+	# Parameters
+	inlist = imtopenp ("input")
+	outlist = imtopenp ("output")
+	call clgstr ("units", Memc[units], SZ_FNAME)
+	maxerr = clgetd ("error")
+	linearize = clgetb ("linearize")
+	verbose = clgetb ("verbose")
+	air = clgwrd ("air", Memc[input], SZ_FNAME, AIRVAC)
+	t = clgetd ("t")
+	p = clgetd ("p")
+	f = clgetd ("f")
+
+	# Loop over input images.
+	while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+
+	    iferr {
+		in = NULL
+		out = NULL
+		mwin = NULL
+		mwout = NULL
+		ctin = NULL
+		ctout = NULL
+		sh = NULL
+		cv = NULL
+
+		# Open input and output images and wcs.
+		if (streq (Memc[input], Memc[output]))
+		    ptr = immap (Memc[input], READ_WRITE, 0)
+		else
+		    ptr = immap (Memc[input], READ_ONLY, 0)
+		in = ptr
+
+		if (streq (Memc[input], Memc[output]))
+		    ptr = in
+		else
+		    ptr = immap (Memc[output], NEW_COPY, in)
+		out = ptr
+
+		ptr = smw_openim (in); mwin = ptr
+		format = SMW_FORMAT(mwin)
+		switch (format) {
+		case SMW_ND:
+		    call error (1,
+			"DISPTRANS does not apply to 2D and 3D images")
+		case SMW_ES:
+		    call smw_esms (mwin)
+		}
+
+		if (IM_NDIM(out) == 3 && IM_LEN(out,3) == 1)
+		    IM_NDIM(out) = 2
+		i = max (2, IM_NDIM(out))
+		ptr = mw_open (NULL, i); mwout = ptr
+		call mw_newsystem (mwout, "multispec", i)
+		call mw_swtype (mwout, axis, 2, "multispec", "")
+		if (i == 3)
+		    call mw_swtype (mwout, 3, 1, "linear", "")
+		call smw_open (mwout, NULL, out)
+		    
+		# Allocate and set arrays.
+		call malloc (x, SMW_LLEN(mwin,1), TY_DOUBLE)
+		call malloc (y, SMW_LLEN(mwin,1), TY_DOUBLE)
+		call malloc (w, SMW_LLEN(mwin,1), TY_DOUBLE)
+		call malloc (nx, SMW_NSPEC(mwin), TY_INT)
+		do i = 1, SMW_LLEN(mwin,1)
+		    Memd[x+i-1] = i
+
+		# Set the output MWCS dispersion function.
+		# Only compute new coordinates once if possible.
+
+		dtype = DCLINEAR
+		do i = 1, SMW_NSPEC(mwin) {
+		    if (format == SMW_MS || i == 1) {
+			call shdr_open (in, mwin, i, 1, INDEFI, SHDATA, sh)
+			call shdr_units (sh, Memc[units])
+			n = SN(sh)
+			do j = 1, n
+			    Memd[y+j-1] = shdr_lw (sh, Memd[x+j-1])
+			call dt_airvac (sh, Memd[y], n, air, t, p, f)
+
+			# Fit dispersion function.
+			dtype1 = DCLINEAR
+			call dt_cvfit (cv, CHEBYSHEV, 2, Memd[x], Memd[y],
+			    Memd[w], n, err)
+			if (err > maxerr) {
+			    dtype1 = DCFUNC
+			    do j = 1, n-4 {
+				call dt_cvfit (cv, SPLINE3, j, Memd[x],
+				    Memd[y], Memd[w], n, err)
+				if (err <= maxerr)
+				    break
+			    }
+			}
+
+			w1 = dcveval (cv, 1D0)
+			dw = (dcveval (cv, double(n)) - w1) / (n - 1)
+		    }
+		    if (linearize) {
+			call dt_setwcs (cv, mwin, mwin, i, dtype1, w1, dw)
+			call dt_setwcs (cv, mwin, mwout, i, DCLINEAR, w1, dw)
+			if (dtype1 != DCLINEAR)
+			    dtype = dtype1
+		    } else
+			call dt_setwcs (cv, mwin, mwout, i, dtype1, w1, dw)
+		    Memi[nx+i-1] = n
+		}
+		call dcvfree (cv)
+
+		# Set label and units.  The check on unit class is done
+		# so that if not a velocity the dictionary expansion
+		# unit is used.  However for velocity the units do not
+		# include the reference coordinate so the user string
+		# is used.
+
+		call mw_swattrs (SMW_MW(mwout,0), 1, "label", LABEL(sh))
+		if (UN_CLASS(UN(sh)) != UN_VEL) {
+		    call mw_swattrs (SMW_MW(mwout,0), 1, "units", UNITS(sh))
+		    call mw_swattrs (SMW_MW(mwout,0), 1, "units_display",
+			UNITS(sh))
+		} else {
+		    call mw_swattrs (SMW_MW(mwout,0), 1, "units", Memc[units])
+		    call mw_swattrs (SMW_MW(mwout,0), 1, "units_display",
+			Memc[units])
+		}
+
+		# Linearize or copy the pixels as requested.
+		if (linearize && dtype != DCLINEAR) {
+		    ptr = smw_sctran (mwin, "world", "logical", 3); ctin = ptr
+		    ptr = smw_sctran (mwout, "logical", "world", 3); ctout = ptr
+		    n = IM_LEN(in,1)
+		    do j = 1, IM_LEN(out,3) {
+			do i = 1, IM_LEN(out,2) {
+			    nw = Memi[nx+i-1]
+			    inbuf = imgl3r (in, i, j)
+			    outbuf = impl3r (out, i, j)
+			    call dispcor (ctin, i, ctout, i, Memr[inbuf], n,
+				Memr[outbuf], nw, NO)
+			    if (nw < n)
+				call amovkr (Memr[outbuf+nw-1], Memr[outbuf+nw],
+				    n-nw)
+			}
+		    }
+		    call smw_ctfree (ctin)
+		    call smw_ctfree (ctout)
+		} else if (in != out) {
+		    n = IM_LEN(in,1)
+		    do j = 1, IM_LEN(out,3) {
+			do i = 1, IM_LEN(out,2) {
+			    inbuf = imgl3r (in, i, j)
+			    outbuf = impl3r (out, i, j)
+			    call amovr (Memr[inbuf], Memr[outbuf], n)
+			}
+		    }
+		}
+
+		# Verbose output
+		if (verbose) {
+		    call printf ("%s: Dispersion transformed to %s")
+			call pargstr (Memc[output])
+			call pargstr (UNITS(sh))
+		    switch (air) {
+		    case 1:
+			call printf (".\n")
+		    case 2:
+			call printf (" in vacuum with\n")
+			call printf (
+			    "  t = %.4g C, p = %.6g mmHg, f = %.4g mmHg.\n")
+			    call pargd (t)
+			    call pargd (p)
+			    call pargd (f)
+		    case 3:
+			call printf (" in air with\n")
+			call printf (
+			    "  t = %.4g C, p = %.6g mmHg, f = %.4g mmHg.\n")
+			    call pargd (t)
+			    call pargd (p)
+			    call pargd (f)
+		    }
+		    call flush (STDOUT)
+		}
+
+	    } then {
+		if (out != NULL && out != in) {
+		    call imunmap (out)
+		    call imdelete (Memc[output])
+		}
+		call erract (EA_WARN)
+	    }
+
+	    # Finish up.
+	    call mfree (x, TY_DOUBLE)
+	    call mfree (y, TY_DOUBLE)
+	    call mfree (w, TY_DOUBLE)
+	    call mfree (nx, TY_INT)
+	    if (mwout != NULL && out != NULL)
+		call smw_saveim (mwout, out)
+	    if (sh != NULL)
+		call shdr_close (sh)
+	    if (ctin != NULL)
+		call smw_ctfree (ctin)
+	    if (ctout != NULL)
+		call smw_ctfree (ctout)
+	    if (mwin != NULL)
+		call smw_close (mwin)
+	    if (mwout != NULL)
+		call smw_close (mwout)
+	    if (out != NULL && out != in)
+		call imunmap (out)
+	    if (in != NULL)
+		call imunmap (in)
+	}
+
+	call imtclose (inlist)
+	call imtclose (outlist)
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end
+
+
+# DT_AIRVAC -- Convert dispersion coordinates to air or vacuum values.
+# The coordinates are first transformed to microns since that is what
+# the formulas expect.  After correction they are transformed back to the
+# original units.  The index of refraction formulas used are from
+# Allen's Astrophysical Quantities (1973).
+
+procedure dt_airvac (sh, x, n, air, t, p, f)
+
+pointer	sh		#I Spectrum pointer
+double	x[n]		#U Dispersion vector
+int	n		#I Number of pixels
+int	air		#I Correction type
+double	t		#I Temperture in deg C
+double	p		#I Total pressure in mmHg
+double	f		#I Water vapour pressure in mmHg
+
+int	i
+double	x2, a
+pointer	un, un_open()
+errchk	un_open, un_ctrand
+
+begin
+	if (air == NONE)
+	    return
+
+	un = un_open ("microns")
+	call un_ctrand (UN(sh), un, x, x, n)
+	do i = 1, n {
+	    x2 = 1 / x[i] **2
+	    a = 64.328 + 29498.1 / (146 - x2) + 255.4 / (41 - x2)
+	    a = a * p * (1 + (1.049 - 0.0157 * t) * 1e-6 * p) /
+		(720.883 * (1 + 0.003661 * t))
+	    a = a - (0.0624 - 0.000680 * x2) / (1 + 0.003661 * t) * f
+	    a = 1 + a / 1e6
+	    switch (air) {
+	    case AIR2VAC:
+		x[i] = a * x[i]
+	    case VAC2AIR:
+		x[i] = x[i] / a
+	    }
+	}
+	call un_ctrand (un, UN(sh), x, x, n)
+	call un_close (un)
+end
+
+
+# DT_CVFIT -- Fit a dispersion function and return the curfit pointer and
+# maximum error in pixels.
+
+procedure dt_cvfit (cv, func, order, x, y, w, n, maxerr)
+
+pointer	cv		#O Fitted dispersion function
+int	func		#I Dispersion function type
+int	order		#I Dispersion function order
+double	x[n]		#I Pixel coordinates
+double	y[n]		#I Desired world coordinates
+double	w[n]		#O Errors in pixels
+int	n		#I Number of pixels
+double	maxerr		#O Maximum error
+
+int	i
+double	minerr, dcveval()
+
+begin
+	if (cv != NULL)
+	    call dcvfree (cv)
+	call dcvinit (cv, func, order, x[1], x[n])
+	call dcvfit (cv, x, y, w, n, WTS_UNIFORM, i)
+	do i = 2, n-1
+	    w[i] = abs ((y[i] - dcveval (cv, x[i])) / ((y[i+1] - y[i-1]) / 2))
+	w[1] = abs ((y[1] - dcveval (cv, x[1])) / (y[2] - y[1]))
+	w[n] = abs ((y[n] - dcveval (cv, x[n])) / (y[n] - y[n-1]))
+	call alimd (w, n, minerr, maxerr)
+end
+
+
+# DT_SETWCS -- Set the multispec WCS.  If the type is nonlinear then
+# the fitted function is stored.
+
+procedure dt_setwcs (cv, mwin, mwout, l, dtype, w1, dw)
+
+pointer	cv		#I Dispersion function
+pointer mwin		#I Input SMW pointer
+pointer mwout		#I Output, SMW pointer
+int	l		#I Image line
+int	dtype		#I Dispersion function type
+double	w1		#I Coordinate of first pixel
+double	dw		#I Coordinate interval at first physical pixel
+
+int	i, ap, bm, dt, nw, n, fd, dcvstati(), stropen()
+double	a, b, z, lw, up
+pointer	sp, title, coeff, coeffs
+
+begin
+	call smark (sp)
+	call salloc (title, SZ_LINE, TY_CHAR)
+
+	coeff = NULL
+	call smw_gwattrs (mwin, l, 1, ap, bm, dt, a, b, nw, z, lw, up, coeff)
+	call smw_gapid (mwin, l, 1, Memc[title], SZ_LINE)
+
+	switch (dtype) {
+	case DCFUNC:
+	    n = dcvstati (cv, CVNSAVE)
+	    call malloc (coeffs, n, TY_DOUBLE)
+	    call dcvsave (cv, Memd[coeffs])
+	    call realloc (coeff, 20*(n+2), TY_CHAR)
+	    fd = stropen (Memc[coeff], 20*(n+2), NEW_FILE)
+	    call fprintf (fd, "1 0 %d %d")
+	    call pargi (nint (Memd[coeffs]))
+	    call pargi (nint (Memd[coeffs+1]))
+	    do i = 2, n-1 {
+		call fprintf (fd, " %g")
+		    call pargd (Memd[coeffs+i])
+	    }
+	    call close (fd)
+	    call mfree (coeffs, TY_DOUBLE)
+	default:
+	    Memc[coeff] = EOS
+	}
+	dt = dtype
+	a = w1
+	b = dw
+	z = 0.
+	call smw_swattrs (mwout, l, 1, ap, bm, dt, a, b, nw, z, lw, up,
+	    Memc[coeff])
+	call smw_sapid (mwout, l, 1, Memc[title])
+
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end