Initial commit

31 files changed, 4809 insertions, 0 deletions

diff --git a/noao/onedspec/smw/README b/noao/onedspec/smw/README
new file mode 100644
index 00000000..2f2c27b0
--- /dev/null
+++ b/noao/onedspec/smw/README
@@ -0,0 +1,6 @@
+This directory contains interface routines for the spectral world
+coordinate systems.  The interface has two functions.  The first is to
+convert various input formats, including old formats, to one of the WCS
+formats used by the ONEDSPEC package.  These are MULTISPEC, EQUISPEC, and
+NDSPEC.  The second is to split large numbers of spectra which exceed the
+limits of the MWCS for a single WCS into a number of smaller WCS.
diff --git a/noao/onedspec/smw/funits.x b/noao/onedspec/smw/funits.x
new file mode 100644
index 00000000..e3c8ddf5
--- /dev/null
+++ b/noao/onedspec/smw/funits.x
@@ -0,0 +1,445 @@
+include	<ctype.h>
+include	<error.h>
+include	<funits.h>
+
+
+# FUN_OPEN -- Open funits package
+# It is allowed to open an unknown funit type
+
+pointer procedure fun_open (funits)
+
+char	funits[ARB]		# Units string
+pointer	fun			# Units pointer returned
+
+begin
+	call calloc (fun, FUN_LEN, TY_STRUCT)
+	iferr (call fun_decode (fun, funits)) {
+	    call fun_close (fun)
+	    call erract (EA_ERROR)
+	}
+	return (fun)
+end
+
+
+# FUN_CLOSE -- Close funits package
+
+procedure fun_close (fun)
+
+pointer	fun			# Units pointer
+
+begin
+	call mfree (fun, TY_STRUCT)
+end
+
+
+# FUN_COPY -- Copy funits pointer
+
+procedure fun_copy (fun1, fun2)
+
+pointer	fun1, fun2		# Units pointers
+
+begin
+	if (fun2 == NULL)
+	    call malloc (fun2, FUN_LEN, TY_STRUCT)
+	call amovi (Memi[fun1], Memi[fun2], FUN_LEN)
+end
+
+
+# FUN_DECODE -- Decode funits string and set up funits structure.
+# The main work is done in FUN_DECODE1 so that the funits string may
+# be recursive; i.e. the funits string may contain other funits strings.
+
+procedure fun_decode (fun, funits)
+
+pointer	fun			# Units pointer
+char	funits[ARB]		# Units string
+
+bool	streq()
+pointer	sp, funits1, temp
+errchk	fun_decode1, fun_ctranr
+
+begin
+	if (streq (funits, FUN_USER(fun)))
+	    return
+
+	call smark (sp)
+	call salloc (funits1, SZ_LINE, TY_CHAR)
+	call salloc (temp, FUN_LEN, TY_STRUCT)
+
+	# Save a copy to restore in case of an error.
+	call fun_copy (fun, temp)
+
+	iferr (call fun_decode1 (fun, funits, Memc[funits1], SZ_LINE)) {
+	    call fun_copy (temp, fun)
+	    call sfree (sp)
+	    call erract (EA_ERROR)
+	}
+
+	call sfree (sp)
+end
+
+
+# FUN_DECODE1 -- Decode funits string and set up funits structure.
+# Unknown funit strings are allowed.
+
+procedure fun_decode1 (fun, funits, funits1, sz_funits1)
+
+pointer	fun			# Units pointer
+char	funits[ARB]		# Units string
+char	funits1[sz_funits1]	# Secondary funits string to return
+int	sz_funits1		# Size of secondary funits string
+
+int	funmod, funtype
+int	i, j, k, nscan(), strdic(), strlen()
+real	funscale
+pointer	sp, str
+
+int	class[FUN_NUNITS]
+real	scale[FUN_NUNITS]
+data	class /FUN_FREQ,FUN_FREQ,FUN_FREQ,FUN_WAVE/
+data	scale /FUN_J,FUN_FU,FUN_CGSH,FUN_CGSA/
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	call strcpy (funits, Memc[str], SZ_FNAME)
+	call strlwr (Memc[str])
+	call sscan (Memc[str])
+	funtype = 0
+	funmod = 0
+	do i = 1, 2 {
+	    call gargwrd (Memc[str], SZ_FNAME)
+	    if (nscan() != i)
+		break
+
+	    j = strdic (Memc[str], Memc[str], SZ_FNAME, FUN_DIC)
+	    for (k=strlen(Memc[str]); k>0 &&
+		(IS_WHITE(Memc[str+k-1]) || Memc[str+k-1]=='\n'); k=k-1)
+		Memc[str+k-1] = EOS
+
+	    if (j > FUN_NUNITS) {
+		if (funmod != 0)
+		    break
+		funmod = j - FUN_NUNITS
+	    } else {
+		funtype = j
+		break
+	    }
+	}
+	i = nscan()
+	call gargr (funscale)
+	if (nscan() != i+1)
+	    funscale = 1
+
+	if (funtype == 0) {
+	    FUN_TYPE(fun) = 0
+	    FUN_CLASS(fun) = FUN_UNKNOWN
+	    FUN_LABEL(fun) = EOS
+	    call strcpy (funits, FUN_UNITS(fun), SZ_UNITS)
+	} else {
+	    FUN_TYPE(fun) = funtype
+	    FUN_CLASS(fun) = class[funtype]
+	    FUN_MOD(fun) = funmod
+	    FUN_SCALE(fun) = scale[funtype] * funscale
+	    FUN_LABEL(fun) = EOS
+	    FUN_UNITS(fun) = EOS
+	    call strcpy (funits, FUN_USER(fun), SZ_UNITS)
+	    switch (funmod) {
+	    case FUN_LOG:
+		call strcat ("Log ", FUN_LABEL(fun), SZ_UNITS)
+	    case FUN_MAG:
+		call strcat ("Mag ", FUN_LABEL(fun), SZ_UNITS)
+	    }
+	    call strcat ("Flux", FUN_LABEL(fun), SZ_UNITS)
+	    if (funscale != 1) {
+		call sprintf (FUN_UNITS(fun), SZ_UNITS, "%sx%.1g")
+		    call pargstr (Memc[str])
+		    call pargr (funscale)
+	    } else {
+		call sprintf (FUN_UNITS(fun), SZ_UNITS, "%s")
+		    call pargstr (Memc[str])
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# FUN_COMPARE -- Compare two funits
+
+bool procedure fun_compare (fun1, fun2)
+
+pointer	fun1, fun2		# Units pointers to compare
+bool	strne()
+
+begin
+	if (strne (FUN_UNITS(fun1), FUN_UNITS(fun2)))
+	    return (false)
+	if (strne (FUN_LABEL(fun1), FUN_LABEL(fun2)))
+	    return (false)
+	return (true)
+end
+
+
+# FUN_CTRANR -- Transform funits
+# Error is returned if the transform cannot be made
+
+procedure fun_ctranr (fun1, fun2, dun, dval, fval1, fval2, nvals)
+
+pointer	fun1			# Input funits pointer
+pointer	fun2			# Output funits pointer
+pointer	dun			# Input units pointer
+real	dval[nvals]		# Input dispersion values
+real	fval1[nvals]		# Input flux values
+real	fval2[nvals]		# Output flux values
+int	nvals			# Number of values
+
+int	i
+real	s, lambda
+pointer	ang, un_open()
+bool	fun_compare()
+errchk	un_open, un_ctranr
+
+begin
+	if (fun_compare (fun1, fun2)) {
+	    call amovr (fval1, fval2, nvals)
+	    return
+	}
+
+	if (FUN_CLASS(fun1) == FUN_UNKNOWN || FUN_CLASS(fun2) == FUN_UNKNOWN)
+	    call error (1, "Cannot convert between selected funits")
+
+	call amovr (fval1, fval2, nvals)
+
+	s = FUN_SCALE(fun1)
+	switch (FUN_MOD(fun1)) {
+	case FUN_LOG:
+	    do i = 1, nvals
+		fval2[i] = 10. ** fval2[i]
+	case FUN_MAG:
+	    do i = 1, nvals
+		fval2[i] = 10. ** (-0.4 * fval2[i])
+	}
+	switch (FUN_CLASS(fun1)) {
+	case FUN_FREQ:
+	    do i = 1, nvals
+		fval2[i] = fval2[i] / s
+	case FUN_WAVE:
+	    if (FUN_CLASS(fun2) != FUN_WAVE) {
+		s = s * FUN_VLIGHT
+		ang = un_open ("angstroms")
+		do i = 1, nvals {
+		    call un_ctranr (dun, ang, dval[i], lambda, 1)
+		    fval2[i] = fval2[i] / s * lambda**2
+		}
+		call un_close (ang)
+	    } else {
+		do i = 1, nvals
+		    fval2[i] = fval2[i] / s
+	    }
+	}
+
+	s = FUN_SCALE(fun2)
+	switch (FUN_CLASS(fun2)) {
+	case FUN_FREQ:
+	    do i = 1, nvals
+		fval2[i] = fval2[i] * s
+	case FUN_WAVE:
+	    if (FUN_CLASS(fun1) != FUN_WAVE) {
+		s = s * FUN_VLIGHT
+		ang = un_open ("angstroms")
+		do i = 1, nvals {
+		    call un_ctranr (dun, ang, dval[i], lambda, 1)
+		    fval2[i] = fval2[i] * s / lambda**2
+		}
+		call un_close (ang)
+	    } else {
+		do i = 1, nvals
+		    fval2[i] = fval2[i] * s
+	    }
+	}
+	switch (FUN_MOD(fun2)) {
+	case FUN_LOG:
+	    do i = 1, nvals
+		fval2[i] = log10 (fval2[i])
+	case FUN_MAG:
+	    do i = 1, nvals
+		fval2[i] = -2.5 * log10 (fval2[i])
+	}
+end
+
+
+# FUN_CHANGER -- Change funits
+# Error is returned if the conversion cannot be made
+
+procedure fun_changer (fun, funits, dun, dvals, fvals, nvals, update)
+
+pointer	fun			# Units pointer (may be changed)
+char	funits[ARB]		# Desired funits
+pointer	dun			# Dispersion units pointer
+real	dvals[nvals]		# Dispersion values
+real	fvals[nvals]		# Flux Values
+int	nvals			# Number of values
+int	update			# Update funits pointer?
+
+bool	streq(), fun_compare()
+pointer	fun1, fun_open()
+errchk	fun_open, fun_ctranr
+
+begin
+
+	# Check for same funit string
+	if (streq (funits, FUN_USER(fun)))
+	    return
+
+	# Check for error in funits string, or the same funits.
+	fun1 = fun_open (funits)
+	if (fun_compare (fun1, fun)) {
+	    call strcpy (funits, FUN_USER(fun), SZ_UNITS)
+	    call fun_close (fun1)
+	    return
+	}
+
+	iferr {
+	    call fun_ctranr (fun, fun1, dun, dvals, fvals, fvals, nvals)
+	    if (update == YES)
+		call fun_copy (fun1, fun)
+	    call fun_close(fun1)
+	} then {
+	    call fun_close(fun1)
+	    call erract (EA_ERROR)
+	}
+end
+
+
+# FUN_CTRAND -- Transform funits
+# Error is returned if the transform cannot be made
+
+procedure fun_ctrand (fun1, fun2, dun, dval, fval1, fval2, nvals)
+
+pointer	fun1			# Input funits pointer
+pointer	fun2			# Output funits pointer
+pointer	dun			# Input dispersion units pointer
+double	dval[nvals]		# Input dispersion values
+double	fval1[nvals]		# Input flux values
+double	fval2[nvals]		# Output flux values
+int	nvals			# Number of values
+
+int	i
+double	s, lambda
+pointer	ang, un_open()
+bool	fun_compare()
+errchk	un_open, un_ctrand
+
+begin
+	if (fun_compare (fun1, fun2)) {
+	    call amovd (fval1, fval2, nvals)
+	    return
+	}
+
+	if (FUN_CLASS(fun1) == FUN_UNKNOWN || FUN_CLASS(fun2) == FUN_UNKNOWN)
+	    call error (1, "Cannot convert between selected funits")
+
+	call amovd (fval1, fval2, nvals)
+
+	s = FUN_SCALE(fun1)
+	switch (FUN_MOD(fun1)) {
+	case FUN_LOG:
+	    do i = 1, nvals
+		fval2[i] = 10. ** fval2[i]
+	case FUN_MAG:
+	    do i = 1, nvals
+		fval2[i] = 10. ** (-0.4 * fval2[i])
+	}
+	switch (FUN_CLASS(fun1)) {
+	case FUN_FREQ:
+	    do i = 1, nvals
+		fval2[i] = fval2[i] / s
+	case FUN_WAVE:
+	    if (FUN_CLASS(fun2) != FUN_WAVE) {
+		s = s * FUN_VLIGHT
+		ang = un_open ("angstroms")
+		do i = 1, nvals {
+		    call un_ctrand (dun, ang, dval[i], lambda, 1)
+		    fval2[i] = fval2[i] / s * lambda**2
+		}
+		call un_close (ang)
+	    } else {
+		do i = 1, nvals
+		    fval2[i] = fval2[i] / s
+	    }
+	}
+
+	s = FUN_SCALE(fun2)
+	switch (FUN_CLASS(fun2)) {
+	case FUN_FREQ:
+	    do i = 1, nvals
+		fval2[i] = fval2[i] * s
+	case FUN_WAVE:
+	    if (FUN_CLASS(fun1) != FUN_WAVE) {
+		s = s * FUN_VLIGHT
+		ang = un_open ("angstroms")
+		do i = 1, nvals {
+		    call un_ctrand (dun, ang, dval[i], lambda, 1)
+		    fval2[i] = fval2[i] * s / lambda**2
+		}
+		call un_close (ang)
+	    } else {
+		do i = 1, nvals
+		    fval2[i] = fval2[i] * s
+	    }
+	}
+	switch (FUN_MOD(fun2)) {
+	case FUN_LOG:
+	    do i = 1, nvals
+		fval2[i] = log10 (fval2[i])
+	case FUN_MAG:
+	    do i = 1, nvals
+		fval2[i] = -2.5 * log10 (fval2[i])
+	}
+
+end
+
+
+# FUN_CHANGED -- Change funits
+# Error is returned if the conversion cannot be made
+
+procedure fun_changed (fun, funits, dun, dvals, fvals, nvals, update)
+
+pointer	fun			# Units pointer (may be changed)
+char	funits[ARB]		# Desired funits
+pointer	dun			# Input dispersion pointer
+double	dvals[nvals]		# Input dispersion values
+double	fvals[nvals]		# Flux values
+int	nvals			# Number of values
+int	update			# Update funits pointer?
+
+bool	streq(), fun_compare()
+pointer	fun1, fun_open()
+errchk	fun_open, fun_ctrand
+
+begin
+
+	# Check for same funit string
+	if (streq (funits, FUN_USER(fun)))
+	    return
+
+	# Check for error in funits string, or the same funits.
+	fun1 = fun_open (funits)
+	if (fun_compare (fun1, fun)) {
+	    call strcpy (funits, FUN_USER(fun), SZ_UNITS)
+	    call fun_close (fun1)
+	    return
+	}
+
+	iferr {
+	    call fun_ctrand (fun, fun1, dun, dvals, fvals, fvals, nvals)
+	    if (update == YES)
+		call fun_copy (fun1, fun)
+	    call fun_close(fun1)
+	} then {
+	    call fun_close(fun1)
+	    call erract (EA_ERROR)
+	}
+end
diff --git a/noao/onedspec/smw/mkpkg b/noao/onedspec/smw/mkpkg
new file mode 100644
index 00000000..64326969
--- /dev/null
+++ b/noao/onedspec/smw/mkpkg
@@ -0,0 +1,48 @@
+# SMW/SHDR Interface
+
+update:
+	$checkout libsmw.a noaolib$
+	$update   libsmw.a
+	$checkin  libsmw.a noaolib$
+	;
+
+generic:
+	$set	GEN = "$$generic -k"
+
+	$ifolder (smwctran.x, smwctran.gx)
+	    $(GEN) smwctran.gx -o smwctran.x $endif
+	;
+
+libsmw.a:
+	$ifeq (USE_GENERIC, yes) $call generic $endif
+
+	funits.x	<ctype.h> <error.h> <funits.h>
+	shdr.x		<error.h> <funits.h> <imset.h> <math/iminterp.h>\
+			<smw.h> <units.h> <imhdr.h>
+	smwclose.x	<smw.h>
+	smwct.x	<smw.h>
+	smwctfree.x	<smw.h>
+	smwctran.x	<smw.h>
+	smwdaxis.x	<smw.h>
+	smwequispec.x	<mwset.h> <smw.h> <imhdr.h>
+	smwesms.x	<mwset.h> <smw.h>
+	smwgapid.x	<smw.h>
+	smwgwattrs.x	<error.h> <smw.h>
+	smwmerge.x	<mwset.h> <smw.h>
+	smwmultispec.x	<smw.h>
+	smwmw.x		<smw.h>
+	smwnd.x		<imhdr.h> <smw.h>
+	smwndes.x	<imhdr.h> <smw.h>
+	smwnewcopy.x	<smw.h>
+	smwoldms.x	<mwset.h> <smw.h>
+	smwonedspec.x	<smw.h> <imhdr.h>
+	smwopen.x	<smw.h>
+	smwopenim.x	<imio.h> <mwset.h> <imhdr.h>
+	smwsapid.x	<smw.h>
+	smwsaveim.x	<imio.h> <smw.h> <imhdr.h>
+	smwsaxes.x	<imhdr.h> <mwset.h> <smw.h>
+	smwsctran.x	<smw.h>
+	smwsmw.x	<smw.h>
+	smwswattrs.x	<error.h> <smw.h>
+	units.x		<ctype.h> <error.h> <units.h>
+	;
diff --git a/noao/onedspec/smw/shdr.x b/noao/onedspec/smw/shdr.x
new file mode 100644
index 00000000..bdcc4b95
--- /dev/null
+++ b/noao/onedspec/smw/shdr.x
@@ -0,0 +1,1269 @@
+include	<error.h>
+include <imhdr.h>
+include	<imset.h>
+include	<smw.h>
+include	<units.h>
+include	<funits.h>
+include	<math/iminterp.h>
+
+
+# SHDR_OPEN    -- Open the SHDR spectrum header structure.
+# SHDR_TYPE    -- Determine spectrum type.
+# SHDR_GTYPE   -- Get the selected spectrum type.
+# SHDR_CLOSE   -- Close and free the SHDR structure.
+# SHDR_COPY    -- Make a copy of an SHDR structure.
+# SHDR_SYSTEM  -- Set or change the WCS system.
+# SHDR_UNITS   -- Set or change user units.
+# SHDR_LW      -- Logical to world coordinate transformation.
+# SHDR_WL      -- World to logical coordinate transformation.
+# SHDR_REBIN   -- Rebin spectrum to dispersion of reference spectrum.
+# SHDR_LINEAR  -- Rebin spectrum to linear dispersion.
+# SHDR_EXTRACT -- Extract a specific wavelength region.
+# SHDR_GI      -- Load an integer value from the header.
+# SHDR_GR      -- Load a real value from the header.
+
+
+# SHDR_OPEN -- Open SHDR spectrum header structure.
+#
+# This routine sets header information, WCS transformations, and extracts the
+# spectrum from EQUISPEC, MULTISPEC, and NDSPEC format images.  The spectrum
+# from a 2D/3D format is specified by a logical line and band number.
+# Optionally an EQUISPEC or MULTISPEC spectrum may be selected by it's
+# aperture number.  The access modes are header only or header and data.
+# Special checks are made to avoid repeated setting of the header and WCS
+# information common to all spectra in an image provided the previously set
+# structure is input.  Note that the logical to world and world to logical
+# transformations require that the MWCS pointer not be closed.
+
+procedure shdr_open (im, smw, index1, index2, ap, mode, sh)
+
+pointer	im			# IMIO pointer
+pointer	smw			# SMW pointer
+int	index1			# Image index desired
+int	index2			# Image index desired
+int	ap			# Aperture number desired
+int	mode			# Access mode
+pointer	sh			# SHDR pointer
+
+int	i, j, k, l, n, np, np1, np2, aplow[2], aphigh[2], strncmp()
+real	smw_c1tranr(), asumr()
+double	dval, shdr_lw()
+bool	newim, streq()
+pointer	sp, key, str, coeff, mw, ct, buf
+pointer	smw_sctran(), imgs3r(), un_open(), fun_open()
+errchk	smw_sctran, imgstr, imgeti, imgetr, smw_gwattrs
+errchk	un_open, fun_open, fun_ctranr, imgs3r, shdr_gtype
+
+define	data_	90
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Allocate basic structure or check if the same spectrum is requested.
+	if (sh == NULL) {
+	    call calloc (sh, LEN_SHDR, TY_STRUCT)
+	    call calloc (SID(sh,1), LEN_SHDRS, TY_CHAR)
+	    newim = true
+	} else {
+	    call imstats (im, IM_IMAGENAME, Memc[str], SZ_LINE)
+	    newim = !streq (Memc[str], IMNAME(sh))
+	    if (!newim) {
+		if (LINDEX(sh,1)==index1 && LINDEX(sh,2)==index2) {
+		    if (IS_INDEFI(ap) || AP(sh)==ap) {
+			np1 = NP1(sh)
+			np2 = NP2(sh)
+			np = np2 - np1 + 1
+			if (CTLW(sh) == NULL || CTWL(sh) == NULL)
+			    goto data_
+			if (mode == SHHDR) {
+			    do i = 1, SH_NTYPES
+				call mfree (SPEC(sh,i), TY_REAL)
+			} else {
+			    switch (SMW_FORMAT(smw)) {
+			    case SMW_ND:
+				if (mode == SHDATA && SPEC(sh,mode) == NULL)
+				    goto data_
+			    case SMW_ES, SMW_MS:
+				if (SPEC(sh,mode) == NULL)
+				    goto data_
+			    }
+			}
+			call sfree (sp)
+			return
+		    }
+		}
+	    }
+	}
+
+	# Set parameters common to an entire image.
+	if (newim) {
+	    call imstats (im, IM_IMAGENAME, IMNAME(sh), LEN_SHDRS)
+	    IM(sh) = im
+	    MW(sh) = smw
+
+	    # Get standard parameters.
+	    call shdr_gi (im, "OFLAG", OBJECT, OFLAG(sh))
+	    call shdr_gr (im, "EXPOSURE", INDEF, IT(sh))
+	    call shdr_gr (im, "ITIME", IT(sh), IT(sh))
+	    call shdr_gr (im, "EXPTIME", IT(sh), IT(sh))
+	    call shdr_gr (im, "RA", INDEF, RA(sh))
+	    call shdr_gr (im, "DEC", INDEF, DEC(sh))
+	    call shdr_gr (im, "UT", INDEF, UT(sh))
+	    call shdr_gr (im, "ST", INDEF, ST(sh))
+	    call shdr_gr (im, "HA", INDEF, HA(sh))
+	    call shdr_gr (im, "AIRMASS", INDEF, AM(sh))
+	    call shdr_gi (im, "DC-FLAG", DCNO, DC(sh))
+	    call shdr_gi (im, "EX-FLAG", ECNO, EC(sh))
+	    call shdr_gi (im, "CA-FLAG", FCNO, FC(sh))
+	    iferr (call imgstr (im, "DEREDDEN", RC(sh), LEN_SHDRS))
+		RC(sh) = EOS
+
+	    # Flag bad airmass value; i.e. 0.
+	    if (!IS_INDEF (AM(sh)) && AM(sh) < 1.)
+		AM(sh) = INDEF
+
+	    # Set the SMW information.
+	    if (SMW_FORMAT(smw) == SMW_MS)
+		i = 3B
+	    else
+		i = 2 ** (SMW_PAXIS(smw,1) - 1)
+	    CTLW1(sh) = smw_sctran (smw, "logical", "world", i)
+	    CTWL1(sh) = smw_sctran (smw, "world", "logical", i)
+
+	    # Set the units.
+	    mw = SMW_MW(smw,0)
+	    i = SMW_PAXIS(smw,1)
+	    iferr (call mw_gwattrs (mw, i, "label", LABEL(sh),LEN_SHDRS))
+		call strcpy ("", LABEL(sh), LEN_SHDRS)
+	    if (streq (LABEL(sh), "equispec") || streq (LABEL(sh), "multispe"))
+		call strcpy ("", LABEL(sh), LEN_SHDRS)
+	    iferr (call mw_gwattrs (mw, i, "units", UNITS(sh),LEN_SHDRS)) {
+		call sprintf (Memc[key], SZ_FNAME, "cunit%d")
+		    call pargi (i)
+		iferr (call imgstr (im, Memc[key], UNITS(sh), LEN_SHDRS)) {
+		    call strlwr (LABEL(sh))
+		    if (LABEL(sh) == EOS)
+			call strcpy ("", UNITS(sh), LEN_SHDRS)
+		    else if (streq (LABEL(sh), "lambda"))
+			call strcpy ("angstroms", UNITS(sh), LEN_SHDRS)
+		    else if (streq (LABEL(sh), "freq"))
+			call strcpy ("hertz", UNITS(sh), LEN_SHDRS)
+		    else if (strncmp (LABEL(sh), "velo", 4) == 0)
+			call strcpy ("m/s", UNITS(sh), LEN_SHDRS)
+		    else if (streq (LABEL(sh), "waveleng"))
+			call strcpy ("angstroms", UNITS(sh), LEN_SHDRS)
+		    else
+			call strcpy ("", UNITS(sh), LEN_SHDRS)
+		}
+		if (strncmp (LABEL(sh), "velo", 4) == 0)
+		    call strcat (" 21 centimeters", UNITS(sh), LEN_SHDRS)
+	    }
+	    if (UNITS(sh) == EOS && DC(sh) != DCNO)
+		call strcpy ("Angstroms", UNITS(sh), LEN_SHDRS)
+	    MWUN(sh) = un_open (UNITS(sh))
+	    call un_copy (MWUN(sh), UN(sh))
+
+	    iferr (call imgstr (im, "bunit", Memc[str], SZ_LINE))
+		call strcpy ("", Memc[str], SZ_LINE)
+	    FUNIM(sh) = fun_open (Memc[str])
+	    if (FUN_CLASS(FUNIM(sh)) != FUN_UNKNOWN)
+		FC(sh) = FCYES
+
+	    call fun_copy (FUNIM(sh), FUN(sh))
+	    call strcpy (FUN_LABEL(FUN(sh)), FLABEL(sh), LEN_SHDRS)
+	    call strcpy (FUN_UNITS(FUN(sh)), FUNITS(sh), LEN_SHDRS)
+	}
+
+	# Set WCS parameters for spectrum type.
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    # Set physical and logical indices.
+	    if (!IS_INDEFI (ap)) {
+		i = max (1, min (SMW_NSPEC(smw), ap))
+		j = 1
+	    } else {
+		i = max (1, index1)
+		j = max (1, index2)
+	    }
+	    call smw_mw (smw, i, j, mw, k, l)
+
+	    LINDEX(sh,1) = max (1, min (SMW_LLEN(smw,2), k))
+	    LINDEX(sh,2) = max (1, min (SMW_LLEN(smw,3), l))
+	    PINDEX(sh,1) = LINDEX(sh,1)
+	    PINDEX(sh,2) = LINDEX(sh,2)
+	    APINDEX(sh) = LINDEX(sh,1)
+
+	    # Set aperture information.  Note the use of the logical index.
+	    np1 = 1
+	    call smw_gwattrs (smw, i, j, AP(sh), BEAM(sh), DC(sh),
+		dval, dval, np2, dval, APLOW(sh,1), APHIGH(sh,1), coeff) 
+
+	    call smw_gapid (smw, i, j, TITLE(sh), LEN_SHDRS)
+	    Memc[SID(sh,1)] = EOS
+
+	    switch (SMW_LDIM(smw)) {
+	    case 1:
+		IMSEC(sh) = EOS
+	    case 2:
+		if (APLOW(sh,1) == APHIGH(sh,1)) {
+		    if (SMW_LAXIS(smw,1) == 1)
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d]")
+		    else
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,*]")
+		    call pargi (nint (APLOW(sh,1)))
+		} else {
+		    if (SMW_LAXIS(smw,1) == 1)
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d:%d]")
+		    else
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,*]")
+		    call pargi (nint (APLOW(sh,1)))
+		    call pargi (nint (APHIGH(sh,1)))
+		}
+	    case 3:
+		if (APLOW(sh,1)==APHIGH(sh,1) && APLOW(sh,2)==APHIGH(sh,2)) {
+		    switch (SMW_LAXIS(smw,1)) {
+		    case 1:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d,%d]")
+		    case 2:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,*,%d]")
+		    case 3:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,%d,*]")
+		    }
+		    call pargi (nint (APLOW(sh,1)))
+		    call pargi (nint (APLOW(sh,2)))
+		} else if (APLOW(sh,1) == APHIGH(sh,1)) {
+		    switch (SMW_LAXIS(smw,1)) {
+		    case 1:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d,%d:%d]")
+		    case 2:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,*,%d:%d]")
+		    case 3:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,%d:%d,*]")
+		    }
+		    call pargi (nint (APLOW(sh,1)))
+		    call pargi (nint (APLOW(sh,2)))
+		    call pargi (nint (APHIGH(sh,2)))
+		} else if (APLOW(sh,2) == APHIGH(sh,2)) {
+		    switch (SMW_LAXIS(smw,1)) {
+		    case 1:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d:%d,%d]")
+		    case 2:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,*,%d]")
+		    case 3:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,%d,*]")
+		    }
+		    call pargi (nint (APLOW(sh,1)))
+		    call pargi (nint (APHIGH(sh,1)))
+		    call pargi (nint (APLOW(sh,2)))
+		} else {
+		    switch (SMW_LAXIS(smw,1)) {
+		    case 1:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d:%d,%d:%d]")
+		    case 2:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,*,%d:%d]")
+		    case 3:
+			call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,%d:%d,*]")
+		    }
+		    call pargi (nint (APLOW(sh,1)))
+		    call pargi (nint (APHIGH(sh,1)))
+		    call pargi (nint (APLOW(sh,2)))
+		    call pargi (nint (APHIGH(sh,2)))
+		}
+	    }
+
+	case SMW_ES, SMW_MS:
+	    # Set the image and aperture indices.
+	    if (SMW_PAXIS(smw,2) != 3) {
+		PINDEX(sh,1) = max (1, min (SMW_LLEN(smw,2), index1))
+		PINDEX(sh,2) = max (1, min (SMW_LLEN(smw,3), index2))
+		LINDEX(sh,1) = PINDEX(sh,1)
+		LINDEX(sh,2) = PINDEX(sh,2)
+		APINDEX(sh) = LINDEX(sh,1)
+	    } else {
+		PINDEX(sh,1) = 1
+		PINDEX(sh,2) = max (1, min (SMW_LLEN(smw,2), index2))
+		LINDEX(sh,1) = PINDEX(sh,2)
+		LINDEX(sh,2) = 1
+		APINDEX(sh) = 1
+	    }
+
+	    # If an aperture is specified first try and find it.
+	    # If it is not specified or found then use the physical index.
+
+	    coeff = NULL
+	    AP(sh) = 0
+	    if (!IS_INDEFI(ap)) {
+		do i = 1, SMW_NSPEC(smw) {
+		    call smw_gwattrs (smw, i, 1, AP(sh), BEAM(sh), DC(sh),
+			dval, dval, np2, dval, APLOW(sh,1), APHIGH(sh,1), coeff)
+		    if (AP(sh) == ap && SMW_PAXIS(smw,2) != 3) {
+			PINDEX(sh,1) = i
+			LINDEX(sh,1) = i
+			APINDEX(sh) = i
+			break
+		    }
+		}
+	    }
+	    if (AP(sh) != ap)
+		call smw_gwattrs (smw, APINDEX(sh), 1, AP(sh), BEAM(sh), DC(sh),
+		    dval, dval, np2, dval, APLOW(sh,1), APHIGH(sh,1), coeff) 
+	    call mfree (coeff, TY_CHAR)
+
+	    np1 = 1
+	    if (SMW_PDIM(smw) > 1) {
+		ct = smw_sctran (smw, "logical", "physical", 2)
+		PINDEX(sh,1) = nint (smw_c1tranr (ct, real(PINDEX(sh,1))))
+		call smw_ctfree (ct)
+	    }
+	    if (SMW_PDIM(smw) > 2) {
+		ct = smw_sctran (smw, "logical", "physical", 4)
+		PINDEX(sh,2) = nint (smw_c1tranr (ct, real(PINDEX(sh,2))))
+		call smw_ctfree (ct)
+	    }
+
+	    call smw_gapid (smw, APINDEX(sh), 1, TITLE(sh), LEN_SHDRS)
+	    call shdr_type (sh, 1, PINDEX(sh,2))
+
+	    switch (SMW_LDIM(smw)) {
+	    case 1:
+		IMSEC(sh) = EOS
+	    case 2:
+		call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d]")
+		call pargi (APINDEX(sh))
+	    case 3:
+		call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d,%d]")
+		call pargi (APINDEX(sh))
+		call pargi (LINDEX(sh,2))
+	    }
+	}
+	
+	# Set NP1 and NP2 in logical coordinates.
+	i = 2 ** (SMW_PAXIS(smw,1) - 1)
+	ct = smw_sctran (smw, "physical", "logical", i)
+	i = max (1, min (int (smw_c1tranr (ct, real (np1))), SMW_LLEN(smw,1)))
+	j = max (1, min (int (smw_c1tranr (ct, real (np2))), SMW_LLEN(smw,1)))
+	call smw_ctfree (ct)
+	np1 = min (i, j)
+	np2 = max (i, j)
+	np = np2 - np1 + 1
+
+	NP1(sh) = np1
+	NP2(sh) = np2
+	SN(sh) = np
+
+
+data_	# Set the coordinate and data arrays if desired otherwise free them.
+	CTLW(sh) = CTLW1(sh)
+	CTWL(sh) = CTWL1(sh)
+
+	# Set linear dispersion terms.
+	W0(sh) = shdr_lw (sh, double(1))
+	W1(sh) = shdr_lw (sh, double(np))
+	WP(sh) = (W1(sh) - W0(sh)) / (np - 1)
+	SN(sh) = np
+
+	if (mode == SHHDR) {
+	    do i = 1, SH_NTYPES
+		call mfree (SPEC(sh,i), TY_REAL)
+	    call sfree (sp)
+	    return
+	}
+
+	# Set WCS array
+	if (SX(sh) == NULL)
+	    call malloc (SX(sh), np, TY_REAL)
+	else
+	    call realloc (SX(sh), np, TY_REAL)
+	do i = 1, np
+	    Memr[SX(sh)+i-1] = shdr_lw (sh, double(i))
+
+	# Set spectrum array in most efficient way.
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    if (mode == SHDATA || SY(sh) == NULL) {
+		if (SY(sh) == NULL)
+		    call malloc (SY(sh), np, TY_REAL)
+		else
+		    call realloc (SY(sh), np, TY_REAL)
+		call aclrr (Memr[SY(sh)], np)
+		if (IS_INDEF(APLOW(sh,1)))
+		    aplow[1] = 1
+		else
+		    aplow[1] = nint (APLOW(sh,1))
+		if (IS_INDEF(APHIGH(sh,1)))
+		    aphigh[1] = 1
+		else
+		    aphigh[1] = nint (APHIGH(sh,1))
+		if (IS_INDEF(APLOW(sh,2)))
+		    aplow[2] = 1
+		else
+		    aplow[2] = nint (APLOW(sh,2))
+		if (IS_INDEF(APHIGH(sh,2)))
+		    aphigh[2] = 1
+		else
+		    aphigh[2] = nint (APHIGH(sh,2))
+		k = aplow[1]
+		l = aphigh[1]
+		n = aphigh[1] - aplow[1] + 1
+		if (SMW_LAXIS(smw,1) == 1) {
+		    do j = aplow[2], aphigh[2] {
+			do i = aplow[1], aphigh[1] {
+			    buf = imgs3r (im, np1, np2, i, i, j, j)
+			    call aaddr (Memr[buf], Memr[SY(sh)],
+				 Memr[SY(sh)], np)
+			}
+		    }
+		} else if (SMW_LAXIS(smw,1) == 2) {
+		    do j = aplow[2], aphigh[2] {
+			do i = np1, np2 {
+			    buf = imgs3r (im, k, l, i, i, j, j)
+			    Memr[SY(sh)+i-np1] = Memr[SY(sh)+i-np1] +
+				asumr (Memr[buf], n)
+			}
+		    }
+		} else {
+		    do i = np1, np2 {
+			do j = aplow[2], aphigh[2] {
+			    buf = imgs3r (im, k, l, j, j, i, i)
+			    Memr[SY(sh)+i-np1] = Memr[SY(sh)+i-np1] +
+				asumr (Memr[buf], n)
+			}
+		    }
+		}
+	    }
+	case SMW_ES, SMW_MS:
+	    if (mode == SHDATA || SY(sh) == NULL) {
+		if (SY(sh) == NULL)
+		    call malloc (SY(sh), np, TY_REAL)
+		else
+		    call realloc (SY(sh), np, TY_REAL)
+		i = LINDEX(sh,1)
+		j = LINDEX(sh,2)
+		buf = imgs3r (im, np1, np2, i, i, j, j)
+		call amovr (Memr[buf], Memr[SY(sh)], np)
+	    }
+
+	    if (mode > SHDATA)
+		call shdr_gtype (sh, mode)
+	}
+
+	# Guess flux scale if necessary.
+	if (FC(sh) == FCYES && FUN_CLASS(FUNIM(sh)) == FUN_UNKNOWN) {
+	    if (Memr[SY(sh)+np/2] < 1e-18)
+		call strcpy ("erg/cm2/s/Hz", Memc[str], SZ_LINE)
+	    else if (Memr[SY(sh)+np/2] < 1e-5)
+		call strcpy ("erg/cm2/s/A", Memc[str], SZ_LINE)
+	    call fun_close (FUNIM(sh))
+	    FUNIM(sh) = fun_open (Memc[str])
+	    if (FUN_CLASS(FUN(sh)) == FUN_UNKNOWN) {
+		call fun_copy (FUNIM(sh), FUN(sh))
+		call strcpy (FUN_LABEL(FUN(sh)), FLABEL(sh), LEN_SHDRS)
+		call strcpy (FUN_UNITS(FUN(sh)), FUNITS(sh), LEN_SHDRS)
+	    }
+	}
+	if (SPEC(sh,mode) != 0)
+	    iferr (call fun_ctranr (FUNIM(sh), FUN(sh), UN(sh), Memr[SX(sh)],
+		Memr[SPEC(sh,mode)], Memr[SPEC(sh,mode)], SN(sh)))
+		;
+
+	call sfree (sp)
+end
+
+
+# SHDR_GTYPE -- Get the selected spectrum type.
+# Currently this only works for multispec data.
+
+procedure shdr_gtype (sh, type)
+
+pointer	sh			#I SHDR pointer
+int	type			#I Spectrum type
+
+int	i, j, ctowrd(), strdic()
+pointer	sp, key, str, im, smw, ct, buf, smw_sctran(), imgs3r()
+real	smw_c1tranr()
+
+begin
+	im = IM(sh)
+	smw = MW(sh)
+
+	if (SMW_FORMAT(smw) == SMW_ND)
+	    return
+	if (SMW_PDIM(smw) < 3) {
+	    if (type != SHDATA && type != SHRAW) {
+		if (SID(sh,type) != NULL)
+		    call mfree (SID(sh,type), TY_CHAR)
+		if (SPEC(sh,type) != NULL)
+		    call mfree (SPEC(sh,type), TY_REAL)
+	    }
+	    return
+	}
+
+	# Find the band.
+	call smark (sp)
+	call salloc (key, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	do i = 1, 5 {
+	    call sprintf (Memc[key], SZ_LINE, "BANDID%d")
+		call pargi (i)
+	    ifnoerr (call imgstr (im, Memc[key], Memc[str], SZ_LINE)) {
+		j = 1
+		if (ctowrd (Memc[str], j, Memc[key], SZ_LINE) == 0)
+		    next
+		if (strdic (Memc[key], Memc[key], SZ_LINE, STYPES) != type)
+		    next
+		if (SID(sh,type) == NULL)
+		    call malloc (SID(sh,type), LEN_SHDRS, TY_CHAR)
+		call strcpy (Memc[str], Memc[SID(sh,type)], LEN_SHDRS)
+		STYPE(sh,type) = type
+		break
+	    }
+	}
+	call sfree (sp)
+	if (i == 6) {
+	    if (SID(sh,type) != NULL)
+		call mfree (SID(sh,type), TY_CHAR)
+	    if (SPEC(sh,type) != NULL)
+		call mfree (SPEC(sh,type), TY_REAL)
+	    return
+	}
+
+	# Map the physical band to logical vector.
+	ct = smw_sctran (smw, "physical", "logical", 4)
+	i = nint (smw_c1tranr (ct, real(i)))
+	call smw_ctfree (ct)
+	if (SMW_PAXIS(smw,2) != 3) {
+	    if (i > SMW_LLEN(smw,3))
+		return
+	    j = i
+	    i = LINDEX(sh,1)
+	} else {
+	    if (i > SMW_LLEN(smw,2))
+		return
+	    j = 1
+	}
+
+	# Get the spectrum.
+	if (SPEC(sh,type) == NULL)
+	    call malloc (SPEC(sh,type), SN(sh), TY_REAL)
+	else
+	    call realloc (SPEC(sh,type), SN(sh), TY_REAL)
+	buf = imgs3r (im, NP1(sh), NP2(sh), i, i, j, j)
+	call amovr (Memr[buf], Memr[SPEC(sh,type)], SN(sh))
+end
+
+
+# SHDR_TYPE -- Determine the spectrum type.
+# Currently this only works for multispec data.
+
+procedure shdr_type (sh, index, band)
+
+pointer	sh			#I SHDR pointer
+int	index			#I Index
+int	band			#I Physical band
+
+int	i, ctowrd(), strdic()
+pointer	sp, key
+
+begin
+	if (SMW_FORMAT(MW(sh)) == SMW_ND)
+	    return
+
+	call smark (sp)
+	call salloc (key, SZ_LINE, TY_CHAR)
+
+	if (SID(sh,index) == NULL)
+	    call malloc (SID(sh,index), LEN_SHDRS, TY_CHAR)
+
+	call sprintf (Memc[key], SZ_FNAME, "BANDID%d")
+	    call pargi (band)
+	 iferr (call imgstr (IM(sh), Memc[key], Memc[SID(sh,index)], LEN_SHDRS))
+	    Memc[SID(sh,index)] = EOS
+
+	i = 1
+	if (ctowrd (Memc[SID(sh,index)], i, Memc[key], SZ_LINE) > 0)
+	    STYPE(sh,index) = strdic (Memc[key], Memc[key], SZ_LINE, STYPES)
+	else
+	    STYPE(sh,index) = 0
+
+	call sfree (sp)
+end
+
+
+# SHDR_CLOSE -- Close and free the SHDR structure.
+
+procedure shdr_close (sh)
+
+pointer	sh			# SHDR structure
+int	i
+
+begin
+	if (sh == NULL)
+	    return
+	do i = 1, SH_NTYPES {
+	    call mfree (SPEC(sh,i), TY_REAL)
+	    call mfree (SID(sh,i), TY_CHAR)
+	}
+	call un_close (UN(sh))
+	call un_close (MWUN(sh))
+	call fun_close (FUN(sh))
+	call fun_close (FUNIM(sh))
+	if (MW(sh) != NULL) {
+	    call smw_ctfree (CTLW1(sh))
+	    call smw_ctfree (CTWL1(sh))
+	}
+	call mfree (sh, TY_STRUCT)
+end
+
+
+# SHDR_COPY -- Make a copy of an SHDR structure.
+# The image pointer is not copied and the MWCS pointer and transform pointers
+# may or may not be copied .  The uncopied pointers mean that they will be
+# shared by multiple spectrum structures but it also means that when they are
+# closed the structures will have invalid pointers.  The advantage of not
+# copying is that many spectra may come from the same image and the overhead
+# of having copies of the IMIO and MWCS pointers can be avoided.
+
+procedure shdr_copy (sh1, sh2, wcs)
+
+pointer	sh1		# SHDR structure to copy
+pointer	sh2		# SHDR structure copy
+int	wcs		# Make copy of wcs?
+
+int	i
+pointer	un, mwun, fun, funim, spec[SH_NTYPES], sid[SH_NTYPES], smw_newcopy()
+errchk	shdr_system
+
+begin
+	if (sh2 == NULL) {
+	    call calloc (sh2, LEN_SHDR, TY_STRUCT)
+	    call calloc (SID(sh2,1), LEN_SHDRS, TY_CHAR)
+	}
+
+	un = UN(sh2)
+	mwun = MWUN(sh2)
+	fun = FUN(sh2)
+	funim = FUNIM(sh2)
+	call amovi (SPEC(sh2,1), spec, SH_NTYPES)
+	call amovi (SID(sh2,1), sid, SH_NTYPES)
+	call amovi (Memi[sh1], Memi[sh2], LEN_SHDR)
+	call amovi (spec, SPEC(sh2,1), SH_NTYPES)
+	call amovi (sid, SID(sh2,1), SH_NTYPES)
+	UN(sh2) = un
+	MWUN(sh2) = mwun
+	FUN(sh2) = fun
+	FUNIM(sh2) = funim
+	call un_copy (UN(sh1), UN(sh2))
+	call un_copy (MWUN(sh1), MWUN(sh2))
+	call fun_copy (FUN(sh1), FUN(sh2))
+	call fun_copy (FUNIM(sh1), FUNIM(sh2))
+	do i = 1, SH_NTYPES {
+	    if (SPEC(sh1,i) != NULL) {
+		if (SPEC(sh2,i) == NULL)
+		    call malloc (SPEC(sh2,i), SN(sh1), TY_REAL)
+		else
+		    call realloc (SPEC(sh2,i), SN(sh1), TY_REAL)
+		call amovr (Memr[SPEC(sh1,i)], Memr[SPEC(sh2,i)], SN(sh1))
+	    }
+	}
+
+	if (wcs == YES && MW(sh1) != NULL) {
+	    MW(sh2) = smw_newcopy (MW(sh1))
+	    CTLW1(sh2) = NULL
+	    CTWL1(sh2) = NULL
+	    call shdr_system (sh2, "world")
+	}
+end
+
+
+# SHDR_SYSTEM -- Set or change the WCS system.
+
+procedure shdr_system (sh, system)
+
+pointer	sh			# SHDR pointer
+char	system[ARB]		# System
+
+int	i, sn
+bool	streq()
+double	shdr_lw()
+pointer	smw, mw, smw_sctran(), un_open()
+errchk	smw_sctran, un_open
+
+begin
+	smw = MW(sh)
+	if (smw == NULL)
+	    call error (1, "shdr_system: MWCS not defined")
+
+	call smw_ctfree (CTLW1(sh))
+	call smw_ctfree (CTWL1(sh))
+
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND, SMW_ES:
+	    i = 2 ** (SMW_PAXIS(smw,1) - 1)
+	    CTLW1(sh) = smw_sctran (smw, "logical", system, i)
+	    CTWL1(sh) = smw_sctran (smw, system, "logical", i)
+	case SMW_MS:
+	    CTLW1(sh) = smw_sctran (smw, "logical", system, 3B)
+	    CTWL1(sh) = smw_sctran (smw, system, "logical", 3B)
+	}
+	CTLW(sh) = CTLW1(sh)
+	CTWL(sh) = CTWL1(sh)
+
+	# Set labels and units
+	call un_close (MWUN(sh))
+	if (streq (system, "physical")) {
+	    call strcpy ("Pixel", LABEL(sh), LEN_SHDRS)
+	    call strcpy ("", UNITS(sh), LEN_SHDRS)
+	    MWUN(sh) = un_open (UNITS(sh))
+	} else {
+	    call smw_mw (smw, 1, 1, mw, i, i)
+	    iferr (call mw_gwattrs (mw, SMW_PAXIS(smw,1), "label", LABEL(sh),
+		LEN_SHDRS))
+		call strcpy ("", LABEL(sh), LEN_SHDRS)
+	    if (streq (LABEL(sh), "equispec") || streq (LABEL(sh), "multispe"))
+		call strcpy ("", LABEL(sh), LEN_SHDRS)
+	    iferr (call mw_gwattrs (mw, SMW_PAXIS(smw,1), "units", UNITS(sh),
+		LEN_SHDRS))
+		call strcpy ("", UNITS(sh), LEN_SHDRS)
+	    MWUN(sh) = un_open (UNITS(sh))
+	    call strcpy (UN_LABEL(UN(sh)), LABEL(sh), LEN_SHDRS)
+	    call strcpy (UN_UNITS(UN(sh)), UNITS(sh), LEN_SHDRS)
+	}
+
+	sn = SN(sh)
+	W0(sh) = shdr_lw (sh, double(1))
+	W1(sh) = shdr_lw (sh, double(sn))
+	WP(sh) = (W1(sh) - W0(sh)) / (sn - 1)
+	if (SX(sh) != NULL)
+	    do i = 1, sn
+		Memr[SX(sh)+i-1] = shdr_lw (sh, double(i))
+end
+
+
+# SHDR_UNITS -- Set or change the WCS system.
+# This changes W0, W1, WP, and SX.
+
+procedure shdr_units (sh, units)
+
+pointer	sh			# SHDR pointer
+char	units[ARB]		# Units
+
+int	i, sn
+bool	streq()
+double	shdr_lw()
+pointer	str, un, un_open()
+errchk	un_open
+
+begin
+	# Check for unknown units.
+	if (streq (units, "display")) {
+	    call malloc (str, SZ_LINE, TY_CHAR)
+	    iferr (call mw_gwattrs (SMW_MW(MW(sh),0), SMW_PAXIS(MW(sh),1),
+		"units_display", Memc[str], SZ_LINE)) {
+		un = NULL
+		call un_copy (MWUN(sh), un)
+	    } else
+		un = un_open (Memc[str])
+	    call mfree (str, TY_CHAR)
+	} else if (streq (units, "default")) {
+	    un = NULL
+	    call un_copy (MWUN(sh), un)
+	} else
+	    un = un_open (units)
+	if (UN_CLASS(un) == UN_UNKNOWN || UN_CLASS(MWUN(sh)) == UN_UNKNOWN) {
+	    call un_close (un)
+	    call error (1, "Cannot convert to specified units")
+	}
+
+	# Update the coordinates.
+	call un_close (UN(sh))
+	UN(sh) = un
+
+	call strcpy (UN_LABEL(UN(sh)), LABEL(sh), LEN_SHDRS)
+	call strcpy (UN_UNITS(UN(sh)), UNITS(sh), LEN_SHDRS)
+
+	sn = SN(sh)
+	W0(sh) = shdr_lw (sh, double(1))
+	W1(sh) = shdr_lw (sh, double(sn))
+	WP(sh) = (W1(sh) - W0(sh)) / (sn - 1)
+	if (SX(sh) != NULL)
+	    do i = 1, sn
+		Memr[SX(sh)+i-1] = shdr_lw (sh, double(i))
+end
+
+
+# SHDR_LW -- Logical to world coordinate transformation.
+# The transformation pointer is generally NULL only after SHDR_LINEAR
+
+double procedure shdr_lw (sh, l)
+
+pointer	sh			# SHDR pointer
+double	l			# Logical coordinate
+double	w			# World coordinate
+
+double	l0, l1, l2, w1, smw_c1trand()
+
+begin
+	l0 = l + NP1(sh) - 1
+	if (CTLW(sh) != NULL) {
+	    switch (SMW_FORMAT(MW(sh))) {
+	    case SMW_ND, SMW_ES:
+		w = smw_c1trand (CTLW(sh), l0)
+	    case SMW_MS:
+		call smw_c2trand (CTLW(sh), l0, double (APINDEX(sh)), w, w1)
+	    }
+	} else {
+	    switch (DC(sh)) {
+	    case DCLOG:
+		w = W0(sh) * 10. ** (log10(W1(sh)/W0(sh)) * (l0-1) / (SN(sh)-1))
+	    case DCFUNC:
+		w = W0(sh)
+		call smw_c2trand (CTWL1(sh), w, double (AP(sh)), l1, w1)
+		w = W1(sh)
+		call smw_c2trand (CTWL1(sh), w, double (AP(sh)), l2, w1)
+		if (SN(sh) > 1)
+		    l1 = (l2 - l1) / (SN(sh) - 1) * (l0 - 1) + l1
+		else
+		    l1 = l0 - 1 + l1
+		call smw_c2trand (CTLW1(sh), l1, double (APINDEX(sh)), w, w1)
+	    default:
+		w = W0(sh) + (l0 - 1) * WP(sh)
+	    }
+	}
+
+	iferr (call un_ctrand (MWUN(sh), UN(sh), w, w, 1))
+	    ;
+	return (w)
+end
+
+
+# SHDR_WL -- World to logical coordinate transformation.
+# The transformation pointer is generally NULL only after SHDR_LINEAR
+
+double procedure shdr_wl (sh, w)
+
+pointer	sh			# SHDR pointer
+double	w			# World coordinate
+double	l			# Logical coordinate
+
+double	w1, l1, l2, smw_c1trand()
+
+begin
+	iferr (call un_ctrand (UN(sh), MWUN(sh), w, w1, 1))
+	    w1 = w
+
+	if (CTWL(sh) != NULL) {
+	    switch (SMW_FORMAT(MW(sh))) {
+	    case SMW_ND, SMW_ES:
+		l = smw_c1trand (CTWL(sh), w1)
+	    case SMW_MS:
+		call smw_c2trand (CTWL(sh), w1, double (AP(sh)),l,l1)
+	    }
+	} else {
+	    switch (DC(sh)) {
+	    case DCLOG:
+		l = log10(w1/W0(sh)) / log10(W1(sh)/W0(sh)) * (SN(sh)-1) + 1
+	    case DCFUNC:
+		call smw_c2trand (CTWL1(sh), w1, double (AP(sh)), l, l1)
+
+		w1 = W0(sh)
+		call smw_c2trand (CTWL1(sh), w1, double (AP(sh)), l1, w1)
+		w1 = W1(sh)
+		call smw_c2trand (CTWL1(sh), w1, double (AP(sh)), l2, w1)
+		if (l1 != l2)
+		    l = (SN(sh) - 1) / (l2 - l1) * (l - l1) + 1
+		else
+		    l = l - l1 + 1
+	    default:
+		l = (w1 - W0(sh)) / WP(sh) + 1
+	    }
+	}
+
+	return (l-NP1(sh)+1)
+end
+
+
+# SHDR_REBIN -- Rebin spectrum to dispersion of reference spectrum.
+# The interpolation function is set by ONEDINTERP.
+
+procedure shdr_rebin (sh, shref)
+
+pointer	sh		# Spectrum to be rebinned
+pointer	shref		# Reference spectrum
+
+char	interp[10]
+int	i, j, type, ia, ib, n, clgwrd()
+real	a, b, sum, asieval(), asigrl()
+double	x, w, xmin, xmax, shdr_lw(), shdr_wl()
+pointer	unref, unsave, asi, spec
+bool	fp_equalr()
+
+begin
+	# Check if rebinning is needed
+	if (DC(sh) == DC(shref) && DC(sh) != DCFUNC &&
+	    fp_equalr (W0(sh), W0(shref)) && fp_equalr(WP(sh), WP(shref)) &&
+	    SN(sh) == SN(shref))
+	    return
+
+	# Do everything in units of reference WCS.
+	unref = UN(shref)
+	unsave = UN(sh)
+	UN(SH) = unref
+
+	call asiinit (asi, clgwrd ("interp", interp, 10, II_FUNCTIONS))
+	do type = 1, SH_NTYPES {
+	    if (SPEC(sh,type) == NULL)
+		next
+
+	    # Fit the interpolation function to the spectrum.
+	    # Extend the interpolation by one pixel at each end.
+
+	    n = SN(sh)
+	    call malloc (spec, n+2, TY_REAL)
+	    call amovr (Memr[SPEC(sh,type)], Memr[spec+1], n)
+	    Memr[spec] = Memr[SPEC(sh,type)]
+	    Memr[spec+n+1] = Memr[SPEC(sh,type)+n-1]
+	    call asifit (asi, Memr[spec], n+2)
+	    call mfree (spec, TY_REAL)
+
+	    xmin = 0.5
+	    xmax = n + 0.5
+
+	    # Reallocate data array
+	    if (n != SN(shref)) {
+		n = SN(shref)
+		call realloc (SPEC(sh,type), n, TY_REAL)
+		call aclrr (Memr[SPEC(sh,type)], n)
+	    }
+	    spec = SPEC(sh,type)
+
+	    # Compute the average flux in each output pixel.
+
+	    x = 0.5
+	    w = shdr_lw (shref, x)
+	    x = shdr_wl (sh, w)
+	    b = max (xmin, min (xmax, x)) + 1
+	    do i = 1, n {
+		x = i + 0.5
+		w = shdr_lw (shref, x)
+		x = shdr_wl (sh, w)
+		a = b
+		b = max (xmin, min (xmax, x)) + 1
+		if (a <= b) {
+		    ia = nint (a + 0.5)
+		    ib = nint (b - 0.5)
+		    if (abs (a+0.5-ia) < .00001 && abs (b-0.5-ib) < .00001) {
+			sum = 0.
+			do j = ia, ib
+			    sum = sum + asieval (asi, real(j))
+			if (ib - ia > 0)
+			    sum = sum / (ib - ia)
+		    } else {
+			sum = asigrl (asi, a, b)
+			if (b - a > 0.)
+			    sum = sum / (b - a)
+		    }
+		} else {
+		    ib = nint (b + 0.5)
+		    ia = nint (a - 0.5)
+		    if (abs (a-0.5-ia) < .00001 && abs (b+0.5-ib) < .00001) {
+			sum = 0.
+			do j = ib, ia
+			    sum = sum + asieval (asi, real(j))
+			if (ia - ib > 0)
+			    sum = sum / (ia - ib)
+		    } else {
+			sum = asigrl (asi, b, a)
+			if (a - b > 0.)
+			    sum = sum / (a - b)
+		    }
+		}
+
+		Memr[spec] = sum
+		spec = spec + 1
+	    }
+	}
+	call asifree (asi)
+
+	# Set the rest of the header.  The coordinate transformations are
+	# canceled to indicate they are not valid for the data.  They
+	# are not freed because the same pointer may be used in other
+	# spectra from the same image.
+
+	if (SN(sh) != n)
+	    call realloc (SX(sh), n, TY_REAL)
+	call amovr (Memr[SX(shref)], Memr[SX(sh)], n)
+	DC(sh) = DC(shref)
+	W0(sh) = W0(shref)
+	W1(sh) = W1(shref)
+	WP(sh) = WP(shref)
+	SN(sh) = SN(shref)
+
+	CTLW(sh) = NULL
+	CTWL(sh) = NULL
+
+	# Restore original units
+	UN(sh) = unsave
+	iferr (call un_ctranr (unref, UN(sh), Memr[SX(sh)], Memr[SX(sh)],
+	    SN(sh)))
+	    ;
+end
+
+
+# SHDR_LINEAR -- Rebin spectrum to linear dispersion.
+# The interpolation function is set by ONEDINTERP
+
+procedure shdr_linear (sh, w0, w1, sn, dc)
+
+pointer	sh		# Spectrum to be rebinned
+real	w0		# Wavelength of first logical pixel
+real	w1		# Wavelength of last logical pixel
+int	sn		# Number of pixels
+int	dc		# Dispersion type (DCLINEAR | DCLOG)
+
+char	interp[10]
+int	i, j, type, ia, ib, n, clgwrd()
+real	w0mw, w1mw, a, b, sum, asieval(), asigrl()
+double	x, w, w0l, wp, xmin, xmax, shdr_wl()
+pointer	unsave, asi, spec
+bool	fp_equalr()
+
+begin
+	# Check if rebinning is needed
+	if (DC(sh) == dc && fp_equalr (W0(sh), w0) &&
+	    fp_equalr (W1(sh), w1) && SN(sh) == sn)
+	    return
+
+	# Do everything in units of MWCS.
+	call un_ctranr (UN(sh), MWUN(sh), w0, w0mw, 1)
+	call un_ctranr (UN(sh), MWUN(sh), w1, w1mw, 1)
+	unsave = UN(sh)
+	UN(SH) = MWUN(sh)
+
+	call asiinit (asi, clgwrd ("interp", interp, 10, II_FUNCTIONS))
+	do type = 1, SH_NTYPES {
+	    if (SPEC(sh,type) == NULL)
+		next
+
+	    # Fit the interpolation function to the spectrum.
+	    # Extend the interpolation by one pixel at each end.
+
+	    n = SN(sh)
+	    call malloc (spec, n+2, TY_REAL)
+	    call amovr (Memr[SPEC(sh,type)], Memr[spec+1], n)
+	    Memr[spec] = Memr[SPEC(sh,type)]
+	    Memr[spec+n+1] = Memr[SPEC(sh,type)+n-1]
+	    call asifit (asi, Memr[spec], n+2)
+	    call mfree (spec, TY_REAL)
+
+	    xmin = 0.5
+	    xmax = n + 0.5
+
+	    # Reallocate spectrum data array
+	    if (n != sn) {
+		n = sn
+		call realloc (SPEC(sh,type), n, TY_REAL)
+	    }
+	    spec = SPEC(sh,type)
+
+	    # Integrate across pixels using ASIGRL.
+
+	    x = 0.5
+	    if (dc == DCLOG) {
+		w0l = log10 (w0mw)
+		wp = (log10 (w1mw) - log10(w0mw)) / (n - 1)
+		w = 10. ** (w0l+(x-1)*wp)
+	    } else {
+		wp = (w1mw - w0mw) / (n - 1)
+		w = w0mw + (x - 1) * wp
+	    }
+	    x = shdr_wl (sh, w)
+	    b = max (xmin, min (xmax, x)) + 1
+	    do i = 1, n {
+		x = i + 0.5
+		if (dc == DCLOG)
+		    w = 10. ** (w0l + (x - 1) * wp)
+		else
+		    w = w0mw + (x - 1) * wp
+		x = shdr_wl (sh, w)
+		a = b
+		b = max (xmin, min (xmax, x)) + 1
+		if (a <= b) {
+		    ia = nint (a + 0.5)
+		    ib = nint (b - 0.5)
+		    if (abs (a+0.5-ia) < .00001 && abs (b-0.5-ib) < .00001) {
+			sum = 0.
+			do j = ia, ib
+			    sum = sum + asieval (asi, real(j))
+			if (ib - ia > 0)
+			    sum = sum / (ib - ia)
+		    } else {
+			sum = asigrl (asi, a, b)
+			if (b - a > 0.)
+			    sum = sum / (b - a)
+		    }
+		} else {
+		    ib = nint (b + 0.5)
+		    ia = nint (a - 0.5)
+		    if (abs (a-0.5-ia) < .00001 && abs (b+0.5-ib) < .00001) {
+			sum = 0.
+			do j = ib, ia
+			    sum = sum + asieval (asi, real(j))
+			if (ia - ib > 0)
+			    sum = sum / (ia - ib)
+		    } else {
+			sum = asigrl (asi, b, a)
+			if (a - b > 0.)
+			    sum = sum / (a - b)
+		    }
+		}
+		Memr[spec] = sum
+		spec = spec + 1
+	    }
+	}
+	call asifree (asi)
+
+	# Set the rest of the header.  The coordinate transformations are
+	# canceled to indicate they are not valid for the data.  They
+	# are not freed because the same pointer may be used in other
+	# spectra from the same image.
+
+	if (SN(sh) != n)
+	    call realloc (SX(sh), n, TY_REAL)
+	do i = 0, n-1 {
+	    if (dc == DCLOG)
+		w = 10. ** (w0l + i * wp)
+	    else
+		w = w0mw + i * wp
+	    Memr[SX(sh)+i] = w
+	}
+	W0(sh) = w0
+	W1(sh) = w1
+	WP(sh) = (w1 - w0) / (sn - 1)
+	SN(sh) = sn
+	NP1(sh) = 1
+	NP2(sh) = sn
+	DC(sh) = dc
+
+	CTLW(sh) = NULL
+	CTWL(sh) = NULL
+
+	# Restore original units
+	UN(sh) = unsave
+	iferr (call un_ctranr (MWUN(sh), UN(sh), Memr[SX(sh)], Memr[SX(sh)],
+	    sn))
+	    ;
+end
+
+
+# SHDR_EXTRACT -- Extract a specific wavelength region.
+
+procedure shdr_extract (sh, w1, w2, rebin)
+
+pointer	sh			# SHDR structure
+real	w1			# Starting wavelength
+real	w2			# Ending wavelength
+bool	rebin			# Rebin wavelength region?
+
+int	i, j, i1, i2, n
+double	l1, l2
+pointer	buf
+bool	fp_equald()
+double	shdr_wl(), shdr_lw()
+errchk	shdr_linear, shdr_lw, shdr_wl
+
+begin
+	l1 = shdr_wl (sh, double (w1))
+	l2 = shdr_wl (sh, double (w2))
+	if (fp_equald(l1,l2) || max(l1,l2) < 1 || min (l1,l2) > SN(sh))
+	    call error (1, "No pixels to extract")
+	l1 = max (1D0, min (double (SN(sh)), l1))
+	l2 = max (1D0, min (double (SN(sh)), l2))
+	i1 = nint (l1)
+	i2 = nint (l2)
+	n = abs (i2 - i1) + 1
+
+	if (rebin) {
+	    l1 = shdr_lw (sh, l1)
+	    l2 = shdr_lw (sh, l2)
+	    if (DC(sh) == DCFUNC)
+		call shdr_linear (sh, real (l1), real (l2), n, DCLINEAR)
+	    else
+		call shdr_linear (sh, real (l1), real (l2), n, DC(sh))
+	} else {
+	    if (i1 == 1 && i2 == SN(sh))
+		return
+
+	    if (i1 <= i2) {
+		do j = 1, SH_NTYPES
+		    if (SPEC(sh,j) != NULL)
+			call amovr (Memr[SPEC(sh,j)+i1-1], Memr[SPEC(sh,j)], n)
+	    } else {
+		call malloc (buf, n, TY_REAL)
+		do j = 1, SH_NTYPES {
+		    if (SPEC(sh,j) != NULL) {
+			do i = i1, i2, -1
+			    Memr[buf+i1-i] = Memr[SPEC(sh,j)+i-1]
+			call amovr (Memr[buf], Memr[SPEC(sh,j)], n)
+		    }
+		}
+		call mfree (buf, TY_REAL)
+	    }
+	    W0(sh) = Memr[SX(sh)]
+	    W1(sh) = Memr[SX(sh)+n-1]
+	    SN(sh) = n
+	    NP1(sh) = 1
+	    NP2(sh) = n
+	    if (n > 1)
+		WP(sh) = (W1(sh) - W0(sh)) / (SN(sh) - 1)
+	    CTLW(sh) = NULL
+	    CTWL(sh) = NULL
+	}
+end
+
+
+# SHDR_GI -- Load an integer value from the header.
+
+procedure shdr_gi (im, field, default, ival)
+
+pointer	im
+char	field[ARB]
+int	default
+int	ival
+
+int	dummy, imaccf(), imgeti()
+
+begin
+	ival = default
+	if (imaccf (im, field) == YES) {
+	    iferr (dummy = imgeti (im, field))
+		call erract (EA_WARN)
+	    else
+		ival = dummy
+	}
+end
+
+
+# SHDR_GR -- Load a real value from the header.
+
+procedure shdr_gr (im, field, default, rval)
+
+pointer	im
+char	field[ARB]
+real	default
+real	rval
+
+int	imaccf()
+real	dummy, imgetr()
+
+begin
+	rval = default
+	if (imaccf (im, field) == YES) {
+	    iferr (dummy = imgetr (im, field))
+		call erract (EA_WARN)
+	    else
+		rval = dummy
+	}
+end
diff --git a/noao/onedspec/smw/smwclose.x b/noao/onedspec/smw/smwclose.x
new file mode 100644
index 00000000..339ebd98
--- /dev/null
+++ b/noao/onedspec/smw/smwclose.x
@@ -0,0 +1,46 @@
+include	<smw.h>
+
+
+# SMW_CLOSE -- Close the SMW data structure.
+# This includes closing the MWCS pointers.
+
+procedure smw_close (smw)
+
+pointer	smw		# SMW pointer
+
+int	i
+pointer	apids
+
+begin
+	if (smw == NULL)
+	    return
+
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    call mfree (SMW_APID(smw), TY_CHAR)
+	    call mw_close (SMW_MW(smw,0))
+	case SMW_ES:
+	    call mfree (SMW_APS(smw), TY_INT)
+	    call mfree (SMW_BEAMS(smw), TY_INT)
+	    call mfree (SMW_APLOW(smw), TY_REAL)
+	    call mfree (SMW_APHIGH(smw), TY_REAL)
+	    call mfree (SMW_APID(smw), TY_CHAR)
+	    apids = SMW_APIDS(smw) - 1
+	    do i = 1, SMW_NSPEC(smw)
+		call mfree (Memi[apids+i], TY_CHAR)
+	    call mfree (SMW_APIDS(smw), TY_POINTER)
+	    call mw_close (SMW_MW(smw,0))
+	case SMW_MS:
+	    call mfree (SMW_APS(smw), TY_INT)
+	    call mfree (SMW_BEAMS(smw), TY_INT)
+	    call mfree (SMW_APLOW(smw), TY_REAL)
+	    call mfree (SMW_APHIGH(smw), TY_REAL)
+	    call mfree (SMW_APID(smw), TY_CHAR)
+	    apids = SMW_APIDS(smw) - 1
+	    do i = 1, SMW_NSPEC(smw)
+		call mfree (Memi[apids+i], TY_CHAR)
+	    do i = 0, SMW_NMW(smw)-1
+		call mw_close (SMW_MW(smw,i))
+	}
+	call mfree (smw, TY_STRUCT)
+end
diff --git a/noao/onedspec/smw/smwct.x b/noao/onedspec/smw/smwct.x
new file mode 100644
index 00000000..b568f759
--- /dev/null
+++ b/noao/onedspec/smw/smwct.x
@@ -0,0 +1,19 @@
+include	<smw.h>
+
+
+# SMW_CT -- Get MCWS CT pointer for the specified physical line.
+
+pointer procedure smw_ct (sct, line)
+
+pointer	sct		#I SMW pointer
+int	line		#I Physical line
+
+begin
+	if (SMW_NCT(sct) == 1)
+	    return (SMW_CT(sct,0))
+
+	if (line < 1 || line > SMW_NSPEC(SMW_SMW(sct)))
+	    call error (1, "smw_ct: aperture not found")
+
+	return (SMW_CT(sct,(line-1)/SMW_NSPLIT))
+end
diff --git a/noao/onedspec/smw/smwctfree.x b/noao/onedspec/smw/smwctfree.x
new file mode 100644
index 00000000..90a506d7
--- /dev/null
+++ b/noao/onedspec/smw/smwctfree.x
@@ -0,0 +1,19 @@
+include	<smw.h>
+
+
+# SMW_CTFREE -- Free a spectral SMW coordinate transform pointer.
+
+procedure smw_ctfree (ct)
+
+pointer	ct		# SMW CT pointer
+int	i
+
+begin
+	if (ct == NULL)
+	    return
+
+	do i = 0, SMW_NCT(ct)-1
+	    call mw_ctfree (SMW_CT(ct,i))
+	call mw_ctfree (SMW_CTL(ct))
+	call mfree (ct, TY_STRUCT)
+end
diff --git a/noao/onedspec/smw/smwctran.gx b/noao/onedspec/smw/smwctran.gx
new file mode 100644
index 00000000..4029aaab
--- /dev/null
+++ b/noao/onedspec/smw/smwctran.gx
@@ -0,0 +1,166 @@
+include	<smw.h>
+
+
+# Evaluate SMW coordinate transform.  These procedures call the
+# MWCS procedures unless the WCS is a split MULTISPEC WCS.  In that
+# case the appropriate piece needs to be determined and the physical
+# line numbers manipulated.  For log sampled spectra conversions
+# must be made for EQUISPEC/NDSPEC.  The convention is that coordinates
+# are always input and output and linear.  Note that the MULTISPEC
+# function driver already takes care of this.
+#
+# SMW_CTRANR -- N dimensional real coordinate transformation.
+# SMW_CTRAND -- N dimensional double coordinate transformation.
+# SMW_C1TRANR -- One dimensional real coordinate transformation.
+# SMW_C1TRAND -- One dimensional double coordinate transformation.
+# SMW_C2TRANR -- Two dimensional real coordinate transformation.
+# SMW_C2TRAND -- Two dimensional double coordinate transformation.
+
+
+$for (rd)
+# SMW_CTRAN -- N dimensional coordinate transformation.
+
+procedure smw_ctran$t (ct, p1, p2, ndim)
+
+pointer	ct		#I SMW CT pointer
+PIXEL	p1[ndim]	#I Input coordinate
+PIXEL	p2[ndim]	#O Output coordinate
+int	ndim		#I Dimensionality
+
+int	i, j, format, daxis, aaxis, dtype, naps
+pointer	smw, aps
+errchk	mw_ctran$t
+
+begin
+	if (SMW_NCT(ct) != 1)
+	    call error (1, "SMW_CTRAN: Wrong WCS type")
+
+	call amov$t (p1, p2, ndim)
+
+	smw = SMW_SMW(ct)
+	format = SMW_FORMAT(smw)
+	daxis = SMW_DAXIS(ct)
+	aaxis = SMW_AAXIS(ct)
+	dtype = SMW_DTYPE(smw)
+	naps = SMW_NSPEC(smw)
+	aps = SMW_APS(smw)
+	switch (format) {
+	case SMW_ND, SMW_ES:
+	    switch (SMW_CTTYPE(ct)) {
+	    case SMW_LW, SMW_PW:
+		call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim)
+		if (daxis > 0 && dtype == DCLOG)
+		    p2[daxis] = 10. ** max (-20$F, min (20$F, p2[daxis]))
+		if (aaxis > 0 && format == SMW_ES) {
+		    i = max (1, min (naps, nint (p2[aaxis])))
+		    p2[aaxis] = Memi[aps+i-1]
+		}
+	    case SMW_WL, SMW_WP:
+		if (daxis > 0 && dtype == DCLOG)
+		    p2[daxis] = log10 (p2[daxis])
+		if (aaxis > 0 && format == SMW_ES) {
+		    j = nint (p2[aaxis])
+		    p2[aaxis] = 1
+		    do i = 1, naps {
+			if (j == Memi[aps+i-1]) {
+			    p2[aaxis] = i
+			    break
+			}
+		    }
+		}
+		call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim)
+	    default:
+		call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim)
+	    }
+	case SMW_MS:
+	    call mw_ctran$t (SMW_CT(ct,0), p1, p2, ndim)
+	}
+end
+
+
+# SMW_C1TRAN -- One dimensional coordinate transformation.
+
+PIXEL procedure smw_c1tran$t (ct, x1)
+
+pointer	ct		#I SMW CT pointer
+PIXEL	x1		#I Input coordinate
+PIXEL	x2		#O Output coordinate
+
+errchk	mw_ctran$t
+
+begin
+	call smw_ctran$t (ct, x1, x2, 1)
+	return (x2)
+end
+
+
+# SMW_C2TRAN -- Two dimensional coordinate transformation.
+
+procedure smw_c2tran$t (ct, x1, y1, x2, y2)
+
+pointer	ct		#I SMW CT pointer
+PIXEL	x1, y1		#I Input coordinates
+PIXEL	x2, y2		#O Output coordinates
+
+PIXEL	p1[2], p2[2]
+int	i, j, naps
+PIXEL	xp, yp
+pointer	aps, smw_ct()
+errchk	smw_ct, mw_c2tran$t
+
+begin
+	# Unsplit WCS.
+	if (SMW_NCT(ct) == 1) {
+	    p1[1] = x1
+	    p1[2] = y1
+	    call smw_ctran$t (ct, p1, p2, 2)
+	    x2 = p2[1]
+	    y2 = p2[2]
+	    return
+	}
+
+	# If we get here then we are dealing with a split MULTISPEC WCS.
+	# Depending on the systems being transformed there may need to
+	# transformation made on the physical coordinate system.
+
+	switch (SMW_CTTYPE(ct)) {
+	case SMW_LW:
+	    call mw_c2tran$t (SMW_CTL(ct), x1, y1, xp, yp)
+	    i = nint (yp)
+	    yp = mod (i-1, SMW_NSPLIT) + 1
+	    call mw_c2tran$t (smw_ct(ct,i), xp, yp, x2, y2)
+	case SMW_PW:
+	    i = nint (y1)
+	    yp = mod (i-1, SMW_NSPLIT) + 1
+	    call mw_c2tran$t (smw_ct(ct,i), x1, yp, x2, y2)
+	case SMW_WL:
+	    aps = SMW_APS(SMW_SMW(ct))
+	    naps = SMW_NSPEC(SMW_SMW(ct))
+	    j = nint (y1)
+	    do i = 1, naps {
+		if (Memi[aps+i-1] == j) {
+		    call mw_c2tran$t (smw_ct(ct,i), x1, y1, xp, yp)
+		    yp = i
+		    call mw_c2tran$t (SMW_CTL(ct), xp, yp, x2, y2)
+		    return
+		}
+	    }
+	    call error (1, "Aperture not found")
+	case SMW_WP:
+	    aps = SMW_APS(SMW_SMW(ct))
+	    naps = SMW_NSPEC(SMW_SMW(ct))
+	    j = nint (y1)
+	    do i = 1, naps {
+		if (Memi[aps+i-1] == j) {
+		    call mw_c2tran$t (smw_ct(ct,i), x1, y1, x2, y2)
+		    y2 = i
+		    return
+		}
+	    }
+	    call error (1, "Aperture not found")
+	default:
+	    x2 = x1
+	    y2 = y1
+	}
+end
+$endfor
diff --git a/noao/onedspec/smw/smwctran.x b/noao/onedspec/smw/smwctran.x
new file mode 100644
index 00000000..38967be2
--- /dev/null
+++ b/noao/onedspec/smw/smwctran.x
@@ -0,0 +1,312 @@
+include	<smw.h>
+
+
+# Evaluate SMW coordinate transform.  These procedures call the
+# MWCS procedures unless the WCS is a split MULTISPEC WCS.  In that
+# case the appropriate piece needs to be determined and the physical
+# line numbers manipulated.  For log sampled spectra conversions
+# must be made for EQUISPEC/NDSPEC.  The convention is that coordinates
+# are always input and output and linear.  Note that the MULTISPEC
+# function driver already takes care of this.
+#
+# SMW_CTRANR -- N dimensional real coordinate transformation.
+# SMW_CTRAND -- N dimensional double coordinate transformation.
+# SMW_C1TRANR -- One dimensional real coordinate transformation.
+# SMW_C1TRAND -- One dimensional double coordinate transformation.
+# SMW_C2TRANR -- Two dimensional real coordinate transformation.
+# SMW_C2TRAND -- Two dimensional double coordinate transformation.
+
+
+
+# SMW_CTRAN -- N dimensional coordinate transformation.
+
+procedure smw_ctranr (ct, p1, p2, ndim)
+
+pointer	ct		#I SMW CT pointer
+real	p1[ndim]	#I Input coordinate
+real	p2[ndim]	#O Output coordinate
+int	ndim		#I Dimensionality
+
+int	i, j, format, daxis, aaxis, dtype, naps
+pointer	smw, aps
+errchk	mw_ctranr
+
+begin
+	if (SMW_NCT(ct) != 1)
+	    call error (1, "SMW_CTRAN: Wrong WCS type")
+
+	call amovr (p1, p2, ndim)
+
+	smw = SMW_SMW(ct)
+	format = SMW_FORMAT(smw)
+	daxis = SMW_DAXIS(ct)
+	aaxis = SMW_AAXIS(ct)
+	dtype = SMW_DTYPE(smw)
+	naps = SMW_NSPEC(smw)
+	aps = SMW_APS(smw)
+	switch (format) {
+	case SMW_ND, SMW_ES:
+	    switch (SMW_CTTYPE(ct)) {
+	    case SMW_LW, SMW_PW:
+		call mw_ctranr (SMW_CT(ct,0), p2, p2, ndim)
+		if (daxis > 0 && dtype == DCLOG)
+		    p2[daxis] = 10. ** max (-20.0, min (20.0, p2[daxis]))
+		if (aaxis > 0 && format == SMW_ES) {
+		    i = max (1, min (naps, nint (p2[aaxis])))
+		    p2[aaxis] = Memi[aps+i-1]
+		}
+	    case SMW_WL, SMW_WP:
+		if (daxis > 0 && dtype == DCLOG)
+		    p2[daxis] = log10 (p2[daxis])
+		if (aaxis > 0 && format == SMW_ES) {
+		    j = nint (p2[aaxis])
+		    p2[aaxis] = 1
+		    do i = 1, naps {
+			if (j == Memi[aps+i-1]) {
+			    p2[aaxis] = i
+			    break
+			}
+		    }
+		}
+		call mw_ctranr (SMW_CT(ct,0), p2, p2, ndim)
+	    default:
+		call mw_ctranr (SMW_CT(ct,0), p2, p2, ndim)
+	    }
+	case SMW_MS:
+	    call mw_ctranr (SMW_CT(ct,0), p1, p2, ndim)
+	}
+end
+
+
+# SMW_C1TRAN -- One dimensional coordinate transformation.
+
+real procedure smw_c1tranr (ct, x1)
+
+pointer	ct		#I SMW CT pointer
+real	x1		#I Input coordinate
+real	x2		#O Output coordinate
+
+errchk	mw_ctranr
+
+begin
+	call smw_ctranr (ct, x1, x2, 1)
+	return (x2)
+end
+
+
+# SMW_C2TRAN -- Two dimensional coordinate transformation.
+
+procedure smw_c2tranr (ct, x1, y1, x2, y2)
+
+pointer	ct		#I SMW CT pointer
+real	x1, y1		#I Input coordinates
+real	x2, y2		#O Output coordinates
+
+real	p1[2], p2[2]
+int	i, j, naps
+real	xp, yp
+pointer	aps, smw_ct()
+errchk	smw_ct, mw_c2tranr
+
+begin
+	# Unsplit WCS.
+	if (SMW_NCT(ct) == 1) {
+	    p1[1] = x1
+	    p1[2] = y1
+	    call smw_ctranr (ct, p1, p2, 2)
+	    x2 = p2[1]
+	    y2 = p2[2]
+	    return
+	}
+
+	# If we get here then we are dealing with a split MULTISPEC WCS.
+	# Depending on the systems being transformed there may need to
+	# transformation made on the physical coordinate system.
+
+	switch (SMW_CTTYPE(ct)) {
+	case SMW_LW:
+	    call mw_c2tranr (SMW_CTL(ct), x1, y1, xp, yp)
+	    i = nint (yp)
+	    yp = mod (i-1, SMW_NSPLIT) + 1
+	    call mw_c2tranr (smw_ct(ct,i), xp, yp, x2, y2)
+	case SMW_PW:
+	    i = nint (y1)
+	    yp = mod (i-1, SMW_NSPLIT) + 1
+	    call mw_c2tranr (smw_ct(ct,i), x1, yp, x2, y2)
+	case SMW_WL:
+	    aps = SMW_APS(SMW_SMW(ct))
+	    naps = SMW_NSPEC(SMW_SMW(ct))
+	    j = nint (y1)
+	    do i = 1, naps {
+		if (Memi[aps+i-1] == j) {
+		    call mw_c2tranr (smw_ct(ct,i), x1, y1, xp, yp)
+		    yp = i
+		    call mw_c2tranr (SMW_CTL(ct), xp, yp, x2, y2)
+		    return
+		}
+	    }
+	    call error (1, "Aperture not found")
+	case SMW_WP:
+	    aps = SMW_APS(SMW_SMW(ct))
+	    naps = SMW_NSPEC(SMW_SMW(ct))
+	    j = nint (y1)
+	    do i = 1, naps {
+		if (Memi[aps+i-1] == j) {
+		    call mw_c2tranr (smw_ct(ct,i), x1, y1, x2, y2)
+		    y2 = i
+		    return
+		}
+	    }
+	    call error (1, "Aperture not found")
+	default:
+	    x2 = x1
+	    y2 = y1
+	}
+end
+
+# SMW_CTRAN -- N dimensional coordinate transformation.
+
+procedure smw_ctrand (ct, p1, p2, ndim)
+
+pointer	ct		#I SMW CT pointer
+double	p1[ndim]	#I Input coordinate
+double	p2[ndim]	#O Output coordinate
+int	ndim		#I Dimensionality
+
+int	i, j, format, daxis, aaxis, dtype, naps
+pointer	smw, aps
+errchk	mw_ctrand
+
+begin
+	if (SMW_NCT(ct) != 1)
+	    call error (1, "SMW_CTRAN: Wrong WCS type")
+
+	call amovd (p1, p2, ndim)
+
+	smw = SMW_SMW(ct)
+	format = SMW_FORMAT(smw)
+	daxis = SMW_DAXIS(ct)
+	aaxis = SMW_AAXIS(ct)
+	dtype = SMW_DTYPE(smw)
+	naps = SMW_NSPEC(smw)
+	aps = SMW_APS(smw)
+	switch (format) {
+	case SMW_ND, SMW_ES:
+	    switch (SMW_CTTYPE(ct)) {
+	    case SMW_LW, SMW_PW:
+		call mw_ctrand (SMW_CT(ct,0), p2, p2, ndim)
+		if (daxis > 0 && dtype == DCLOG)
+		    p2[daxis] = 10. ** max (-20.0D0, min (20.0D0, p2[daxis]))
+		if (aaxis > 0 && format == SMW_ES) {
+		    i = max (1, min (naps, nint (p2[aaxis])))
+		    p2[aaxis] = Memi[aps+i-1]
+		}
+	    case SMW_WL, SMW_WP:
+		if (daxis > 0 && dtype == DCLOG)
+		    p2[daxis] = log10 (p2[daxis])
+		if (aaxis > 0 && format == SMW_ES) {
+		    j = nint (p2[aaxis])
+		    p2[aaxis] = 1
+		    do i = 1, naps {
+			if (j == Memi[aps+i-1]) {
+			    p2[aaxis] = i
+			    break
+			}
+		    }
+		}
+		call mw_ctrand (SMW_CT(ct,0), p2, p2, ndim)
+	    default:
+		call mw_ctrand (SMW_CT(ct,0), p2, p2, ndim)
+	    }
+	case SMW_MS:
+	    call mw_ctrand (SMW_CT(ct,0), p1, p2, ndim)
+	}
+end
+
+
+# SMW_C1TRAN -- One dimensional coordinate transformation.
+
+double procedure smw_c1trand (ct, x1)
+
+pointer	ct		#I SMW CT pointer
+double	x1		#I Input coordinate
+double	x2		#O Output coordinate
+
+errchk	mw_ctrand
+
+begin
+	call smw_ctrand (ct, x1, x2, 1)
+	return (x2)
+end
+
+
+# SMW_C2TRAN -- Two dimensional coordinate transformation.
+
+procedure smw_c2trand (ct, x1, y1, x2, y2)
+
+pointer	ct		#I SMW CT pointer
+double	x1, y1		#I Input coordinates
+double	x2, y2		#O Output coordinates
+
+double	p1[2], p2[2]
+int	i, j, naps
+double	xp, yp
+pointer	aps, smw_ct()
+errchk	smw_ct, mw_c2trand
+
+begin
+	# Unsplit WCS.
+	if (SMW_NCT(ct) == 1) {
+	    p1[1] = x1
+	    p1[2] = y1
+	    call smw_ctrand (ct, p1, p2, 2)
+	    x2 = p2[1]
+	    y2 = p2[2]
+	    return
+	}
+
+	# If we get here then we are dealing with a split MULTISPEC WCS.
+	# Depending on the systems being transformed there may need to
+	# transformation made on the physical coordinate system.
+
+	switch (SMW_CTTYPE(ct)) {
+	case SMW_LW:
+	    call mw_c2trand (SMW_CTL(ct), x1, y1, xp, yp)
+	    i = nint (yp)
+	    yp = mod (i-1, SMW_NSPLIT) + 1
+	    call mw_c2trand (smw_ct(ct,i), xp, yp, x2, y2)
+	case SMW_PW:
+	    i = nint (y1)
+	    yp = mod (i-1, SMW_NSPLIT) + 1
+	    call mw_c2trand (smw_ct(ct,i), x1, yp, x2, y2)
+	case SMW_WL:
+	    aps = SMW_APS(SMW_SMW(ct))
+	    naps = SMW_NSPEC(SMW_SMW(ct))
+	    j = nint (y1)
+	    do i = 1, naps {
+		if (Memi[aps+i-1] == j) {
+		    call mw_c2trand (smw_ct(ct,i), x1, y1, xp, yp)
+		    yp = i
+		    call mw_c2trand (SMW_CTL(ct), xp, yp, x2, y2)
+		    return
+		}
+	    }
+	    call error (1, "Aperture not found")
+	case SMW_WP:
+	    aps = SMW_APS(SMW_SMW(ct))
+	    naps = SMW_NSPEC(SMW_SMW(ct))
+	    j = nint (y1)
+	    do i = 1, naps {
+		if (Memi[aps+i-1] == j) {
+		    call mw_c2trand (smw_ct(ct,i), x1, y1, x2, y2)
+		    y2 = i
+		    return
+		}
+	    }
+	    call error (1, "Aperture not found")
+	default:
+	    x2 = x1
+	    y2 = y1
+	}
+end
+
diff --git a/noao/onedspec/smw/smwdaxis.x b/noao/onedspec/smw/smwdaxis.x
new file mode 100644
index 00000000..0bea9375
--- /dev/null
+++ b/noao/onedspec/smw/smwdaxis.x
@@ -0,0 +1,109 @@
+include	<smw.h>
+
+define	CTYPES	"|LAMBDA|FREQ|WAVELENGTH|VELO|VELO-LSR|VELO-HEL|VELO-OBS|"
+
+# SMW_DAXIS -- Set physical dispersion axis and summing factors.
+# A default value of zero for the dispersion axis will cause the dispersion
+# axis to be sought in the image header and, if not found, from the CL
+# "dispaxis" parameter.  A default value of zero for the summing factors will
+# cause them to be queried from the CL "nsum" parameter.  A default value of
+# INDEFI in either parameter will leave the current default unchanged.
+#
+# When this procedure is called with an SMW and IMIO pointer the SMW
+# pointer is updated to desired default dispersion axis and summing
+# parameters.
+
+procedure smw_daxis (smw, im, daxisp, nsum1, nsum2)
+
+pointer	smw			#I SMW pointer
+pointer	im			#I IMIO pointer
+int	daxisp			#I Default dispersion axis
+int	nsum1, nsum2		#I Default summing factors
+
+int	i, da, ns[2]
+int	imgeti(), clgeti(), clscan(), nscan(), nowhite(), strdic()
+pointer	sp, key, val
+data	da/0/, ns/0,0/
+errchk	clgeti
+
+begin
+	# Set defaults.
+	# A value of 0 will use the image DISPAXIS or query the CL and
+	# a value of INDEFI will leave the current default unchanged.
+
+	if (!IS_INDEFI (daxisp))
+	    da = daxisp
+	if (!IS_INDEFI (nsum1))
+	    ns[1] = nsum1
+	if (!IS_INDEFI (nsum2))
+	    ns[2] = nsum2
+
+	if (smw == NULL)
+	    return
+
+	# This procedure is specific to the NDSPEC format.
+	if (SMW_FORMAT(smw) != SMW_ND)
+	    return
+
+	# Set dispersion axis.
+	if (da == 0) {
+	    if (im == NULL)
+		SMW_PAXIS(smw,1) = clgeti ("dispaxis")
+	    else {
+		iferr (SMW_PAXIS(smw,1) = imgeti (im, "DISPAXIS")) {
+		    SMW_PAXIS(smw,1) = clgeti ("dispaxis")
+		    call smark (sp)
+		    call salloc (key, 8, TY_CHAR)
+		    call salloc (val, SZ_FNAME, TY_CHAR)
+		    do i = 1, 7 {
+			call sprintf (Memc[key], 8, "CTYPE%d")
+			    call pargi (i)
+			iferr (call imgstr (im, Memc[key], Memc[val], SZ_FNAME))
+			    break
+			if (nowhite (Memc[val], Memc[val], SZ_FNAME) > 0) {
+			    call strupr (Memc[val])
+			    if (strdic(Memc[val],Memc[val],SZ_FNAME,CTYPES)>0) {
+				SMW_PAXIS(smw,1) = i
+				break
+			    }
+			}
+		    }
+		    call sfree (sp)
+		}
+		if (SMW_PAXIS(smw,1) < 1 || SMW_PAXIS(smw,1) > 7) {
+		    i = SMW_PAXIS(smw,1)
+		    SMW_PAXIS(smw,1) = clgeti ("dispaxis")
+		    call eprintf (
+	"WARNING: Image header dispersion axis %d invalid. Using axis %d.\n")
+			call pargi (i)
+			call pargi (SMW_PAXIS(smw,1))
+		}
+	    }
+	} else
+	    SMW_PAXIS(smw,1) = da
+
+	# Set summing parameters.
+	if (ns[1] == 0 || ns[2] == 0) {
+	    if (clscan("nsum") == EOF)
+		call error (1, "smw_daxis: Error in 'nsum' parameter") 
+	    call gargi (i)
+	    if (ns[1] == 0) {
+		if (nscan() == 1)
+		    SMW_NSUM(smw,1) = max (1, i)
+		else
+		    call error (1, "smw_daxis: Error in 'nsum' parameter")
+	    } else
+		SMW_NSUM(smw,1) = ns[1]
+	    call gargi (i)
+	    if (ns[2] == 0) {
+		if (nscan() == 2)
+		    SMW_NSUM(smw,2) = max (1, i)
+		else
+		    SMW_NSUM(smw,2) = SMW_NSUM(smw,1)
+	    } else
+		SMW_NSUM(smw,2) = ns[2]
+	} else {
+	    SMW_NSUM(smw,1) = ns[1]
+	    SMW_NSUM(smw,2) = ns[2]
+	}
+end
diff --git a/noao/onedspec/smw/smwequispec.x b/noao/onedspec/smw/smwequispec.x
new file mode 100644
index 00000000..5cda57c0
--- /dev/null
+++ b/noao/onedspec/smw/smwequispec.x
@@ -0,0 +1,86 @@
+include <imhdr.h>
+include	<mwset.h>
+include	<smw.h>
+
+
+# SMW_EQUISPEC -- Setup the EQUISPEC SMW parameters.
+# The aperture information is in the APNUM and APID keywords and the
+# WCS information is in the linear MWCS.
+
+procedure smw_equispec (im, smw)
+
+pointer	im			#I IMIO pointer
+pointer	smw			#U MWCS pointer input SMW pointer output
+
+int	i, j, k, nchar, ap, beam, dtype, nw
+double	w1, dw, z
+real	aplow[2], aphigh[2], mw_c1tranr()
+pointer	sp, key, val, wterm, mw, ct, mw_sctran()
+int	imgeti(), mw_stati(), ctoi(), ctor()
+errchk	imgstr, mw_gwtermd, mw_sctran
+errchk	smw_open, smw_saxes, smw_swattrs, smw_sapid
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (val, SZ_LINE, TY_CHAR)
+
+	# Determine the dispersion parameters
+	mw = smw
+	i = mw_stati (mw, MW_NDIM)
+	call salloc (wterm, 2*i+i*i, TY_DOUBLE)
+	call mw_gwtermd (mw, Memd[wterm], Memd[wterm+i], Memd[wterm+2*i], i)
+	w1 = Memd[wterm+i]
+	dw = Memd[wterm+2*i]
+
+	# Determine the number of physical pixels.
+	ct = mw_sctran (mw, "logical", "physical", 1)
+	nw = max (mw_c1tranr (ct, 1.), mw_c1tranr (ct, real(IM_LEN(im,1))))
+	call mw_ctfree (ct)
+
+	# Determine the dispersion type.
+	iferr (dtype = imgeti (im, "DC-FLAG"))
+	    dtype = DCNO
+	if (dtype==DCLOG) {
+	    if (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.)
+		dtype = DCLINEAR
+	    else {
+		w1 = 10D0 ** w1
+		dw = w1 * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1)
+	    }
+	}
+
+	# Set the SMW data structure.
+	call smw_open (smw, NULL, im)
+	do i = 1, SMW_NSPEC(smw) {
+	    call smw_mw (smw, i, 1, mw, j, k)
+	    if (j < 1000)
+		call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
+	    else
+		call sprintf (Memc[key], SZ_FNAME, "AP%d")
+		call pargi (j)
+	    call imgstr (im, Memc[key], Memc[val], SZ_LINE)
+	    k = 1
+	    nchar = ctoi (Memc[val], k, ap)
+	    nchar = ctoi (Memc[val], k, beam)
+	    if (ctor (Memc[val], k, aplow[1]) == 0)
+		aplow[1] = INDEF
+	    if (ctor (Memc[val], k, aphigh[1]) == 0)
+		aphigh[1] = INDEF
+	    if (ctor (Memc[val], k, aplow[2]) == 0)
+		aplow[2] = INDEF
+	    if (ctor (Memc[val], k, aphigh[2]) == 0)
+		aphigh[2] = INDEF
+	    z = 0.
+
+	    call smw_swattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, "")
+
+	    call sprintf (Memc[key], SZ_FNAME, "APID%d")
+		call pargi (j)
+	    ifnoerr (call imgstr (im, Memc[key], Memc[val], SZ_LINE))
+		call smw_sapid (smw, i, 1, Memc[val])
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwesms.x b/noao/onedspec/smw/smwesms.x
new file mode 100644
index 00000000..d1b8a368
--- /dev/null
+++ b/noao/onedspec/smw/smwesms.x
@@ -0,0 +1,96 @@
+include	<mwset.h>
+include	<smw.h>
+
+
+# SMW_ESMS -- Convert EQUISPEC WCS into MULTISPEC WCS.
+# This is called by SMW_SWATTRS when the equal linear coordinate system
+# requirement of the EQUISPEC WCS is violated.
+
+procedure smw_esms (smw)
+
+pointer	smw			#U SMW pointer
+
+int	i, j, k, pdim1, pdim2, ap, beam, dtype, nw, axes[2]
+double	w1, dw, z
+real	aplow, aphigh
+pointer	sp, key, str, lterm, mw, mw1, mw2, apid, mw_open()
+data	axes/1,2/
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (lterm, 12, TY_DOUBLE)
+
+	# Set the basic MWCS
+	mw1 = SMW_MW(smw,0)
+	pdim1 = SMW_PDIM(smw)
+	pdim2 = max (2, pdim1)
+	mw2 = mw_open (NULL, pdim2)
+	call mw_newsystem (mw2, "multispec", pdim2)
+	call mw_swtype (mw2, axes, 2, "multispec", "")
+	if (pdim2 > 2)
+	    call mw_swtype (mw2, 3, 1, "linear", "")
+	call mw_gltermd (mw1, Memd[lterm+pdim2], Memd[lterm], pdim1)
+	if (pdim1 == 1) {
+	    Memd[lterm+1] = 0.
+	    Memd[lterm+3] = 0.
+	    Memd[lterm+4] = 0.
+	    Memd[lterm+5] = 1.
+	}
+	call mw_sltermd (mw2, Memd[lterm+pdim2], Memd[lterm], pdim2)
+	ifnoerr (call mw_gwattrs (mw1, 1, "label", Memc[str], SZ_LINE))
+	    call mw_swattrs (mw2, 1, "label", Memc[str])
+	ifnoerr (call mw_gwattrs (mw1, 1, "units", Memc[str], SZ_LINE))
+	    call mw_swattrs (mw2, 1, "units", Memc[str])
+	ifnoerr (call mw_gwattrs (mw1, 1, "units_display", Memc[str], SZ_LINE))
+	    call mw_swattrs (mw2, 1, "units_display", Memc[str])
+
+	 # Write the MULTISPEC WCS
+	dtype = SMW_DTYPE(smw)
+	w1 = SMW_W1(smw)
+	dw = SMW_DW(smw)
+	nw = SMW_NW(smw)
+	z = SMW_Z(smw)
+	if (dtype == DCLOG) {
+	    dw = log10 ((w1+(nw-1)*dw)/w1)/(nw-1)
+	    w1 = log10 (w1)
+	}
+
+         call smw_open (mw2, smw, 0)
+         do i = 1, SMW_NSPEC(smw) {
+	    ap = Memi[SMW_APS(smw)+i-1]
+	    beam = Memi[SMW_BEAMS(smw)+i-1]
+	    aplow = Memr[SMW_APLOW(smw)+2*i-2]
+	    aphigh = Memr[SMW_APHIGH(smw)+2*i-2]
+	    apid = Memi[SMW_APIDS(smw)+i-1]
+	    
+	    call smw_mw (mw2, i, 1, mw, j, k)
+	    call sprintf (Memc[key], SZ_FNAME, "spec%d")
+		call pargi (j)
+	    call sprintf (Memc[str], SZ_LINE,
+		"%d %d %d %.14g %.14g %d %.14g %.2f %.2f")
+		call pargi (ap)
+		call pargi (beam)
+		call pargi (dtype)
+		call pargd (w1)
+		call pargd (dw)
+		call pargi (nw)
+		call pargd (z)
+		call pargr (aplow)
+		call pargr (aphigh)
+	    call mw_swattrs (mw, 2, Memc[key], Memc[str])
+
+	    if (apid != NULL)
+		call smw_sapid (mw2, i, 1, Memc[apid])
+
+	    # This is used if there is a split MULTISPEC WCS.
+	    if (SMW_APS(mw2) != NULL)
+		Memi[SMW_APS(mw2)+i-1] = ap
+	}
+
+	call smw_close (smw)
+	smw = mw2
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwgapid.x b/noao/onedspec/smw/smwgapid.x
new file mode 100644
index 00000000..214be533
--- /dev/null
+++ b/noao/onedspec/smw/smwgapid.x
@@ -0,0 +1,30 @@
+include	<smw.h>
+
+
+# SMW_GAPID -- Get aperture id.
+
+procedure smw_gapid (smw, index1, index2, apid, maxchar)
+
+pointer	smw				#I SMW pointer
+int	index1				#I Spectrum index
+int	index2				#I Spectrum index
+char	apid[maxchar]			#O Aperture id
+int	maxchar				#I Maximum number of characters
+
+pointer	ptr
+
+begin
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    call strcpy (Memc[SMW_APID(smw)], apid, maxchar)
+	case SMW_ES, SMW_MS:
+	    if (index1 < 0 || index1 > SMW_NSPEC(smw))
+		call error (1, "smw_gapid: index out of bounds")
+
+	    ptr = Memi[SMW_APIDS(smw)+index1-1]
+	    if (index1 == 0 || ptr == NULL)
+		call strcpy (Memc[SMW_APID(smw)], apid, maxchar)
+	    else
+		call strcpy (Memc[ptr], apid, maxchar)
+	}
+end
diff --git a/noao/onedspec/smw/smwgwattrs.x b/noao/onedspec/smw/smwgwattrs.x
new file mode 100644
index 00000000..4084fd4c
--- /dev/null
+++ b/noao/onedspec/smw/smwgwattrs.x
@@ -0,0 +1,134 @@
+include	<error.h>
+include	<smw.h>
+
+
+# SMW_GWATTRS -- Get spectrum attribute parameters.
+# BE CAREFUL OF OUTPUT VARIABLES BEING THE SAME MEMORY ADDRESS!
+
+procedure smw_gwattrs (smw, index1, index2, ap, beam, dtype, w1, dw, nw, z,
+	aplow, aphigh, coeff)
+
+pointer	smw				# SMW pointer
+int	index1				# Spectrum index
+int	index2				# Spectrum index
+int	ap				# Aperture number
+int	beam				# Beam number
+int	dtype				# Dispersion type
+double	w1				# Starting coordinate
+double	dw				# Coordinate interval
+int	nw				# Number of valid pixels
+double	z				# Redshift factor
+real	aplow[2], aphigh[2]		# Aperture limits
+pointer	coeff				# Nonlinear coeff string (input/output)
+
+int	i, j, n, ip, sz_coeff, strlen(), ctoi(), ctor(), ctod()
+double	a, b
+pointer	sp, key, mw
+errchk	smw_mw, mw_gwattrs
+
+data	sz_coeff /SZ_LINE/
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+
+	if (coeff == NULL)
+	    call malloc (coeff, sz_coeff, TY_CHAR)
+	else
+	    call realloc (coeff, sz_coeff, TY_CHAR)
+
+	# Determine parameters based on the SMW format.
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    call smw_mw (smw, index1, index2, mw, i, j)
+
+	    dtype = SMW_DTYPE(smw)
+	    nw = SMW_NW(smw)
+	    w1 = SMW_W1(smw)
+	    dw = SMW_DW(smw)
+	    z = SMW_Z(smw)
+
+	    ap = index1
+	    beam = 0
+	    aplow[1] = 1
+	    aphigh[1] = 1
+	    aplow[2] = 1
+	    aphigh[2] = 1
+	    if (SMW_LDIM(smw) > 1) {
+		aplow[1] = i - (SMW_NSUM(smw,1)-1) / 2
+		aphigh[1] = nint (aplow[1]) + SMW_NSUM(smw,1) - 1
+		aplow[1] = max (1, nint (aplow[1]))
+		aphigh[1] = min (SMW_LLEN(smw,2), nint (aphigh[1]))
+	    }
+	    if (SMW_LDIM(smw) > 2) {
+		aplow[2] = j - (SMW_NSUM(smw,2)-1) / 2
+		aphigh[2] = nint (aplow[2]) + SMW_NSUM(smw,2) - 1
+		aplow[2] = max (1, nint (aplow[2]))
+		aphigh[2] = min (SMW_LLEN(smw,3), nint (aphigh[2]))
+	    }
+
+	    Memc[coeff] = EOS
+	case SMW_ES:
+	    call smw_mw (smw, index1, index2, mw, i, j)
+
+	    dtype = SMW_DTYPE(smw)
+	    nw = SMW_NW(smw)
+	    w1 = SMW_W1(smw)
+	    dw = SMW_DW(smw)
+	    z = SMW_Z(smw)
+
+	    ap = Memi[SMW_APS(smw)+index1-1]
+	    beam = Memi[SMW_BEAMS(smw)+index1-1]
+	    aplow[1] = Memr[SMW_APLOW(smw)+2*index1-2]
+	    aphigh[1] = Memr[SMW_APHIGH(smw)+2*index1-2]
+	    aplow[2] = Memr[SMW_APLOW(smw)+2*index1-1]
+	    aphigh[2] = Memr[SMW_APHIGH(smw)+2*index1-1]
+
+	    Memc[coeff] = EOS
+	case SMW_MS:
+	    call smw_mw (smw, index1, index2, mw, i, j)
+
+	    call sprintf (Memc[key], SZ_FNAME, "spec%d")
+		call pargi (i)
+
+	    call mw_gwattrs (mw, 2, Memc[key], Memc[coeff], sz_coeff)
+	    while (strlen (Memc[coeff]) == sz_coeff) {
+		sz_coeff = 2 * sz_coeff
+		call realloc (coeff, sz_coeff, TY_CHAR)
+		call mw_gwattrs (mw, 2, Memc[key], Memc[coeff], sz_coeff)
+	    }
+
+	    ip = 1
+	    i = ctoi (Memc[coeff], ip, ap)
+	    i = ctoi (Memc[coeff], ip, beam)
+	    i = ctoi (Memc[coeff], ip, j)
+	    i = ctod (Memc[coeff], ip, a)
+	    i = ctod (Memc[coeff], ip, b)
+	    i = ctoi (Memc[coeff], ip, n)
+	    i = ctod (Memc[coeff], ip, z)
+	    i = ctor (Memc[coeff], ip, aplow[1])
+	    i = ctor (Memc[coeff], ip, aphigh[1])
+	    aplow[2] = INDEF
+	    aphigh[2] = INDEF
+	    if (Memc[coeff+ip-1] != EOS)
+		call strcpy (Memc[coeff+ip], Memc[coeff], sz_coeff)
+	    else
+		Memc[coeff] = EOS
+
+	    if (j==DCLOG) {
+		if (abs(a)>20. || abs(a+(n-1)*b)>20.)
+		    j = DCLINEAR
+		else {
+		    a = 10D0 ** a
+		    b = a * (10D0 ** ((n-1)*b) - 1) / (n - 1)
+		}
+	    }
+
+	    dtype = j
+	    w1 = a
+	    dw = b
+	    nw = n
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwmerge.x b/noao/onedspec/smw/smwmerge.x
new file mode 100644
index 00000000..d3e09bd1
--- /dev/null
+++ b/noao/onedspec/smw/smwmerge.x
@@ -0,0 +1,102 @@
+include	<mwset.h>
+include	<smw.h>
+
+
+# SMW_MERGE -- Merge split MWCS array to a single MWCS.
+
+procedure smw_merge (smw)
+
+pointer	smw			#U Input split WCS, output single WCS
+
+int	i, pdim, naps, format, beam, dtype, dtype1, nw, nw1
+int	ap, axes[3]
+double	w1, dw, z, w11, dw1, z1
+real	aplow[2], aphigh[2]
+pointer	sp, key, val, term, coeff, mw, mw1, mw_open()
+data	axes/1,2,3/
+
+begin
+	if (SMW_NMW(smw) == 1)
+	    return
+
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (val, SZ_LINE, TY_CHAR)
+	call salloc (term, 15, TY_DOUBLE)
+	coeff = NULL
+
+	pdim = SMW_PDIM(smw)
+	naps = SMW_NSPEC(smw)
+	mw1 = SMW_MW(smw,0)
+
+	# Determine output WCS format.
+	format = SMW_ES
+	do i = 1, naps {
+	    call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, coeff)
+	    if (i == 1) {
+		dtype1 = dtype
+		w11 = w1
+		dw1 = dw
+		z1 = z
+		nw1 = nw
+	    }
+	    if (dtype>1||dtype!=dtype1||w1!=w11||dw!=dw1||nw!=nw1||z!=z1) {
+		format = SMW_MS
+		break
+	    }
+	}
+
+	# Setup WCS.
+	switch (format) {
+	case SMW_ES:
+	    mw = mw_open (NULL, pdim)
+	    call mw_newsystem (mw, "equispec", pdim)
+	    call mw_swtype (mw, axes, pdim, "linear", "")
+
+	case SMW_MS:
+	    mw = mw_open (NULL, pdim)
+	    call mw_newsystem (mw, "multispec", pdim)
+	    call mw_swtype (mw, axes, pdim, "multispec", "")
+	    if (pdim > 2)
+		call mw_swtype (mw, 3, 1, "linear", "")
+	}
+
+	ifnoerr (call mw_gwattrs (mw1, 1, "label", Memc[val], SZ_LINE))
+	    call mw_swattrs (mw, 1, "label", Memc[val])
+	ifnoerr (call mw_gwattrs (mw1, 1, "units", Memc[val], SZ_LINE))
+	    call mw_swattrs (mw, 1, "units", Memc[val])
+	ifnoerr (call mw_gwattrs (mw1, 1, "units_display", Memc[val], SZ_LINE))
+	    call mw_swattrs (mw, 1, "units_display", Memc[val])
+	call mw_gltermd (mw1, Memd[term+pdim], Memd[term], pdim)
+	call mw_sltermd (mw, Memd[term+pdim], Memd[term], pdim)
+	call mw_gwtermd (mw1, Memd[term], Memd[term+pdim],
+	    Memd[term+2*pdim], pdim)
+	Memd[term] = 1.
+	Memd[term+pdim] = w1 / (1 + z)
+	Memd[term+2*pdim] = dw / (1 + z)
+	call mw_swtermd (mw, Memd[term], Memd[term+pdim],
+	    Memd[term+2*pdim], pdim)
+
+	# Set the SMW structure.
+	call smw_open (mw, smw, NULL)
+	if (format == SMW_MS) {
+	    do i = 1, SMW_NMW(mw) - 1
+		call mw_close (SMW_MW(mw,i))
+	    SMW_NMW(mw) = 1
+	}
+	do i = 1, naps {
+	    call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, coeff)
+	    call smw_swattrs (mw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, Memc[coeff])
+	    call smw_gapid (smw, i, 1, Memc[val], SZ_LINE)
+	    call smw_sapid (mw, i, 1, Memc[val])
+	}
+
+	call smw_close (smw)
+	smw = mw
+
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwmultispec.x b/noao/onedspec/smw/smwmultispec.x
new file mode 100644
index 00000000..18e2dbd0
--- /dev/null
+++ b/noao/onedspec/smw/smwmultispec.x
@@ -0,0 +1,30 @@
+include	<smw.h>
+
+
+# SMW_MULTISPEC -- Setup the MULTISPEC SMW parameters.
+
+procedure smw_multispec (im, smw)
+
+pointer	im			#I IMIO pointer
+pointer	smw			#U MWCS pointer input SMW pointer output
+
+int	i, j, k
+pointer	sp, key, val, mw
+errchk	smw_open, smw_saxes, smw_sapid
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (val, SZ_LINE, TY_CHAR)
+
+	call smw_open (smw, NULL, im)
+	do i = 1, SMW_NSPEC(smw) {
+	    call smw_mw (smw, i, 1, mw, j, k)
+	    call sprintf (Memc[key], SZ_FNAME, "APID%d")
+		call pargi (j)
+	    ifnoerr (call imgstr (im, Memc[key], Memc[val], SZ_LINE))
+		call smw_sapid (smw, i, 1, Memc[val])
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwmw.x b/noao/onedspec/smw/smwmw.x
new file mode 100644
index 00000000..a79aaf98
--- /dev/null
+++ b/noao/onedspec/smw/smwmw.x
@@ -0,0 +1,38 @@
+include	<smw.h>
+
+
+# SMW_MW -- Get MWCS pointer and coordinates from spectrum line and band
+
+procedure smw_mw (smw, line, band, mw, x, y)
+
+pointer	smw		#I SMW pointer
+int	line		#I Spectrum line
+int	band		#I Spectrum band
+pointer	mw		#O MWCS pointer
+int	x, y		#O MWCS coordinates
+
+real	mw_c1tranr()
+
+begin
+	if (line < 1 || line > SMW_NSPEC(smw))
+	    call error (1, "smw_mw: spectrum not found")
+
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    mw = SMW_MW(smw,0)
+	    x = mod (line - 1, SMW_LLEN(smw,2)) + 1
+	    y = (line - 1) / SMW_LLEN(smw,2) + band
+	default:
+	    if (SMW_NMW(smw) == 1) {
+		mw = SMW_MW(smw,0)
+		x = line
+		y = band
+		if (SMW_CTLP(smw) != NULL)
+		    x = nint (mw_c1tranr (SMW_CTLP(smw), real(line)))
+	    } else {
+		mw = SMW_MW(smw,(line-1)/SMW_NSPLIT)
+		x = mod (line - 1, SMW_NSPLIT) + 1
+		y = band
+	    }
+	}
+end
diff --git a/noao/onedspec/smw/smwnd.x b/noao/onedspec/smw/smwnd.x
new file mode 100644
index 00000000..10b48079
--- /dev/null
+++ b/noao/onedspec/smw/smwnd.x
@@ -0,0 +1,19 @@
+include	<imhdr.h>
+include	<smw.h>
+
+
+# SMW_ND -- Setup the NDSPEC SMW.
+# If there is only one spectrum convert it to EQUISPEC if possible.
+
+procedure smw_nd (im, smw)
+
+pointer	im			#I IMIO pointer
+pointer	smw			#U MWCS pointer input SMW pointer output
+
+errchk	smw_open, smw_daxis, smw_ndes
+
+begin
+	call smw_open (smw, NULL, im)
+	if (SMW_NSPEC(smw) == 1)
+	    call smw_ndes (im, smw)
+end
diff --git a/noao/onedspec/smw/smwndes.x b/noao/onedspec/smw/smwndes.x
new file mode 100644
index 00000000..9d35ca6d
--- /dev/null
+++ b/noao/onedspec/smw/smwndes.x
@@ -0,0 +1,82 @@
+include	<imhdr.h>
+include	<smw.h>
+
+
+# SMW_NDES -- Convert NDSPEC WCS into EQUISPEC WCS.
+# This requires that the logical dispersion axis be 1.
+
+procedure smw_ndes (im, smw)
+
+pointer	im			#I IMIO pointer
+pointer	smw			#U Input NDSPEC SMW, output EQUISPEC SMW
+
+int	i, pdim1, pdim2, daxis, ap, beam, dtype, nw, axes[2]
+real	aplow[2], aphigh[2]
+double	w1, dw, z
+pointer	sp, key, str, lterm1, lterm2, coeff, mw1, mw2, mw_open()
+errchk	mw_open, mw_gltermd, mw_gwtermd, smw_open, smw_saxes, smw_gwattrs
+data	axes/1,2/, coeff/NULL/
+
+begin
+	# Require the dispersion to be along the first logical axis.
+	if (SMW_LAXIS(smw,1) != 1)
+	    return
+
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (lterm1, 15, TY_DOUBLE)
+	call salloc (lterm2, 15, TY_DOUBLE)
+
+	# Set the MWCS.  Only the logical and world transformations along
+	# the dispersion axis are transfered.
+
+	pdim1 = SMW_PDIM(smw)
+	pdim2 = IM_NDIM(im)
+	daxis = SMW_PAXIS(smw,1)
+	mw1 = SMW_MW(smw,0)
+
+	mw2 = mw_open (NULL, pdim2)
+	call mw_newsystem (mw2, "equispec", pdim2)
+	call mw_swtype (mw2, axes, pdim2, "linear", "")
+	ifnoerr (call mw_gwattrs (mw1, daxis, "label", Memc[str], SZ_LINE))
+	    call mw_swattrs (mw2, 1, "label", Memc[str])
+	ifnoerr (call mw_gwattrs (mw1, daxis, "units", Memc[str], SZ_LINE))
+	    call mw_swattrs (mw2, 1, "units", Memc[str])
+	ifnoerr (call mw_gwattrs (mw1, daxis, "units_display", Memc[str],
+	    SZ_LINE))
+	    call mw_swattrs (mw2, 1, "units_display", Memc[str])
+
+	call mw_gltermd (mw1, Memd[lterm1+pdim1], Memd[lterm1], pdim1)
+	call mw_gltermd (mw2, Memd[lterm2+pdim2], Memd[lterm2], pdim2)
+	Memd[lterm2] = Memd[lterm1+daxis-1]
+	Memd[lterm2+pdim2] = Memd[lterm1+pdim1+(pdim1+1)*(daxis-1)]
+	call mw_sltermd (mw2, Memd[lterm2+pdim2], Memd[lterm2], pdim2)
+
+	call mw_gwtermd (mw1, Memd[lterm1], Memd[lterm1+pdim1],
+	     Memd[lterm1+2*pdim1], pdim1)
+	call mw_gwtermd (mw2, Memd[lterm2], Memd[lterm2+pdim2],
+	     Memd[lterm2+2*pdim2], pdim2)
+	Memd[lterm2] = Memd[lterm1+daxis-1]
+	Memd[lterm2+pdim2] = Memd[lterm1+pdim1+daxis-1]
+	Memd[lterm2+2*pdim2] = Memd[lterm1+2*pdim1+(pdim1+1)*(daxis-1)]
+	call mw_swtermd (mw2, Memd[lterm2], Memd[lterm2+pdim2],
+	     Memd[lterm2+2*pdim2], pdim2)
+
+	# Set the EQUISPEC SMW.
+	IM_LEN(im,2) = SMW_NSPEC(smw)
+	IM_LEN(im,3) = SMW_NBANDS(smw)
+	call smw_open (mw2, NULL, im)
+	do i = 1, SMW_NSPEC(smw) {
+	    call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, coeff)
+	    call smw_swattrs (mw2, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, Memc[coeff])
+	}
+	call mfree (coeff, TY_CHAR)
+
+	call smw_close (smw)
+	smw = mw2
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwnewcopy.x b/noao/onedspec/smw/smwnewcopy.x
new file mode 100644
index 00000000..230ed205
--- /dev/null
+++ b/noao/onedspec/smw/smwnewcopy.x
@@ -0,0 +1,58 @@
+include <smw.h>
+
+
+# SMW_NEWCOPY -- Make a new copy of an SMW structure.
+
+pointer procedure smw_newcopy (smw)
+
+pointer smw             #I SMW pointer to copy
+pointer new 	        #O SMW copy
+
+int     i, nspec
+pointer mw_newcopy(), mw_sctran()
+
+begin
+	call calloc (new, SMW_LEN(SMW_NMW(smw)), TY_STRUCT)
+	call amovi (Memi[smw], Memi[new], SMW_LEN(SMW_NMW(smw)))
+
+	if (SMW_APID(smw) != NULL) {
+	    call malloc (SMW_APID(new), SZ_LINE, TY_CHAR)
+	    call strcpy (Memc[SMW_APID(smw)], Memc[SMW_APID(new)], SZ_LINE)
+	}
+
+	nspec = SMW_NSPEC(smw)
+	if (SMW_APS(smw) != NULL) {
+	    call malloc (SMW_APS(new), nspec, TY_INT)
+	    call amovi (Memi[SMW_APS(smw)], Memi[SMW_APS(new)], nspec)
+	}
+	if (SMW_BEAMS(smw) != NULL) {
+	    call malloc (SMW_BEAMS(new), nspec, TY_INT)
+	    call amovi (Memi[SMW_BEAMS(smw)], Memi[SMW_BEAMS(new)], nspec)
+	}
+	if (SMW_APLOW(smw) != NULL) {
+	    call malloc (SMW_APLOW(new), 2*nspec, TY_REAL)
+	    call amovr (Memr[SMW_APLOW(smw)], Memr[SMW_APLOW(new)], 2*nspec)
+	}
+	if (SMW_APHIGH(smw) != NULL) {
+	    call malloc (SMW_APHIGH(new), 2*nspec, TY_REAL)
+	    call amovr (Memr[SMW_APHIGH(smw)], Memr[SMW_APHIGH(new)], 2*nspec)
+	}
+	if (SMW_APIDS(smw) != NULL) {
+	    call calloc (SMW_APIDS(new), nspec, TY_POINTER)
+	    do i = 0, nspec-1 {
+		if (Memi[SMW_APIDS(smw)+i] != NULL) {
+		    call malloc (Memi[SMW_APIDS(new)+i], SZ_LINE, TY_CHAR)
+		    call strcpy (Memc[Memi[SMW_APIDS(smw)+i]],
+			Memc[Memi[SMW_APIDS(new)+i]], SZ_LINE)
+		}
+	    }
+	}
+
+        do i = 0, SMW_NMW(smw)-1
+            SMW_MW(new,i) = mw_newcopy (SMW_MW(smw,i))
+
+	if (SMW_PDIM(smw) > 1)
+	    SMW_CTLP(new) = mw_sctran (SMW_MW(new,0), "logical", "physical", 2)
+
+	return (new)
+end
diff --git a/noao/onedspec/smw/smwoldms.x b/noao/onedspec/smw/smwoldms.x
new file mode 100644
index 00000000..6640641f
--- /dev/null
+++ b/noao/onedspec/smw/smwoldms.x
@@ -0,0 +1,101 @@
+include	<mwset.h>
+include	<smw.h>
+
+
+# SMW_OLDMS -- Convert old multispec format into MULTISPEC SMW.
+
+procedure smw_oldms (im, smw)
+
+pointer	im			#I IMIO pointer
+pointer	smw			#U Input MWCS pointer, output SMW pointer
+
+int	i, j, k, nchar, ap, beam, dtype, nw, axes[2]
+double	w1, dw, z
+real	aplow[2], aphigh[2]
+pointer	sp, key, val, lterm, mw, mw_open()
+int	imgeti(), mw_stati(), ctoi(), ctor(), ctod(), imofnlu(), imgnfn()
+errchk	imgstr, mw_gltermd, mw_sltermd
+data	axes/1,2/
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (val, SZ_LINE, TY_CHAR)
+	call salloc (lterm, 12, TY_DOUBLE)
+
+	# Set the basic multispec MWCS
+	i = mw_stati (smw, MW_NDIM)
+	j = max (2, i)
+	mw = mw_open (NULL, j)
+	call mw_newsystem (mw, "multispec", j)
+	call mw_swtype (mw, axes, 2, "multispec", "")
+	if (j > 2)
+	    call mw_swtype (mw, 3, 1, "linear", "")
+	call mw_gltermd (smw, Memd[lterm+j], Memd[lterm], i)
+	if (i == 1) {
+	    Memd[lterm+1] = 0.
+	    Memd[lterm+3] = 0.
+	    Memd[lterm+4] = 0.
+	    Memd[lterm+5] = 1.
+	}
+	call mw_sltermd (mw, Memd[lterm+j], Memd[lterm], j)
+
+	iferr (dtype = imgeti (im, "DC-FLAG"))
+	    dtype = -1
+	else {
+	    call mw_swattrs (mw, 1, "label", "Wavelength")
+	    call mw_swattrs (mw, 1, "units", "Angstroms")
+	}
+
+	call mw_close (smw)
+	smw = mw
+
+	# Set the SMW data structure.
+	call smw_open (smw, NULL, im)
+	do i = 1, SMW_NSPEC(smw) {
+	    call smw_mw (smw, i, 1, mw, j, k)
+	    call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
+		call pargi (j)
+	    call imgstr (im, Memc[key], Memc[val], SZ_LINE)
+	    call imdelf (im, Memc[key])
+
+	    k = 1
+	    nchar = ctoi (Memc[val], k, ap)
+	    nchar = ctoi (Memc[val], k, beam)
+	    nchar = ctod (Memc[val], k, w1)
+	    nchar = ctod (Memc[val], k, dw)
+	    nchar = ctoi (Memc[val], k, nw)
+	    if (ctor (Memc[val], k, aplow[1]) == 0)
+		aplow[1] = INDEF
+	    if (ctor (Memc[val], k, aphigh[1]) == 0)
+		aphigh[1] = INDEF
+	    z = 0.
+
+	    k = dtype
+	    if (k==1 && (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.))
+		k = 0
+	    call smw_swattrs (smw, i, 1, ap, beam, k, w1, dw, nw, z,
+		aplow, aphigh, "")
+
+	    call sprintf (Memc[key], SZ_FNAME, "APID%d")
+		call pargi (j)
+	    ifnoerr (call imgstr (im, Memc[key], Memc[val], SZ_LINE)) {
+		call smw_sapid (smw, i, 1, Memc[val])
+		call imdelf (im, Memc[key])
+	    }
+	}
+
+	# Delete old parameters
+	i = imofnlu (im,
+	    "DISPAXIS,APFORMAT,BEAM-NUM,DC-FLAG,W0,WPC,NP1,NP2")
+	while (imgnfn (i, Memc[key], SZ_FNAME) != EOF) {
+	    iferr (call imdelf (im, Memc[key]))
+		;
+	}
+	call imcfnl (i)
+
+	# Update MWCS
+	call smw_saveim (smw, im)
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwonedspec.x b/noao/onedspec/smw/smwonedspec.x
new file mode 100644
index 00000000..b7d2fa6a
--- /dev/null
+++ b/noao/onedspec/smw/smwonedspec.x
@@ -0,0 +1,109 @@
+include <imhdr.h>
+include	<smw.h>
+
+
+# SMW_ONEDSPEC -- Convert old "onedspec" format to EQUISPEC.
+
+procedure smw_onedspec (im, smw)
+
+pointer	im			#I IMIO pointer
+pointer	smw			#U MWCS pointer input SMW pointer output
+
+int	i, dtype, ap, beam, nw, imgeti(), imofnlu(), imgnfn()
+real	aplow[2], aphigh[2], imgetr(), mw_c1tranr()
+double	ltm, ltv, r, w, dw, z, imgetd()
+pointer	sp, key, mw, ct, mw_openim(), mw_sctran()
+bool	fp_equald()
+errchk	smw_open, smw_saxes, mw_gwtermd, mw_sctran
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+
+	# Convert old W0/WPC keywords if needed.
+	mw = smw
+	iferr (w = imgetd (im, "CRVAL1")) {
+	    ifnoerr (w = imgetd (im, "W0")) {
+		dw = imgetd (im, "WPC")
+		iferr (ltm = imgetd (im, "LTM1_1"))
+		    ltm = 1
+		iferr (ltv = imgetd (im, "LTV1"))
+		    ltv = 0
+		r = ltm + ltv
+		dw = dw / ltm
+		call imaddd (im, "CRPIX1", r)
+		call imaddd (im, "CRVAL1", w)
+		call imaddd (im, "CD1_1", dw)
+		call imaddd (im, "CDELT1", dw)
+		call mw_close(mw)
+		mw = mw_openim (im)
+	    }
+	}
+
+	# Get dispersion and determine number of valid pixels.
+	call mw_gwtermd (mw, r, w, dw, 1)
+	w = w - (r - 1) * dw
+	r = 1
+	call mw_swtermd (mw, r, w, dw, 1)
+	ct = mw_sctran (mw, "logical", "physical", 1)
+	nw = max (mw_c1tranr (ct, 1.), mw_c1tranr (ct, real (IM_LEN(im,1))))
+	call mw_ctfree (ct)
+
+	iferr (dtype = imgeti (im, "DC-FLAG")) {
+	    if (fp_equald (1D0, w) || fp_equald (1D0, dw))
+		dtype = DCNO
+	    else
+		dtype = DCLINEAR
+	}
+	if (dtype==DCLOG) {
+	    if (abs(w)>20. || abs(w+(nw-1)*dw)>20.)
+		dtype = DCLINEAR
+	    else {
+		w = 10D0 ** w
+		dw = w * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1)
+	    }
+	}
+
+	# Convert to EQUISPEC system.
+	call mw_swattrs (mw, 0, "system", "equispec")
+	if (dtype != DCNO) {
+	    iferr (call mw_gwattrs (mw, 1, "label", Memc[key], SZ_FNAME)) {
+		iferr (call mw_gwattrs (mw, 1, "units", Memc[key], SZ_FNAME)) {
+		    call mw_swattrs (mw, 1, "units", "angstroms")
+		    call mw_swattrs (mw, 1, "label", "Wavelength")
+		}
+	    }
+	}
+
+	# Set the SMW data structure.
+	call smw_open (smw, NULL, im)
+
+	# Determine the aperture parameters.
+	iferr (beam = imgeti (im, "BEAM-NUM"))
+	    beam = 1
+	iferr (ap = imgeti (im, "APNUM"))
+	    ap = beam
+	iferr (aplow[1] = imgetr (im, "APLOW"))
+	    aplow[1] = INDEF
+	iferr (aphigh[1] = imgetr (im, "APHIGH"))
+	    aphigh[1] = INDEF
+	iferr (z = imgetd (im, "DOPCOR"))
+	    z = 0.
+
+	call smw_swattrs (smw, 1, 1, ap, beam, dtype, w, dw, nw, z,
+	    aplow, aphigh, "")
+
+        # Delete old parameters
+        i = imofnlu (im,
+            "BEAM-NUM,APNUM,APLOW,APHIGH,DOPCOR,DC-FLAG,W0,WPC,NP1,NP2")
+        while (imgnfn (i, Memc[key], SZ_FNAME) != EOF) {
+            iferr (call imdelf (im, Memc[key]))
+		;
+	}
+	call imcfnl (i)
+
+        # Update MWCS
+        call smw_saveim (smw, im)
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwopen.x b/noao/onedspec/smw/smwopen.x
new file mode 100644
index 00000000..782c8749
--- /dev/null
+++ b/noao/onedspec/smw/smwopen.x
@@ -0,0 +1,70 @@
+include	<smw.h>
+
+
+# SMW_OPEN -- Open SMW structure.
+# The basic MWCS pointer and a template SMW pointer or image is input
+# and the SMW pointer is returned in its place.
+
+procedure smw_open (mw, smw1, im)
+
+pointer	mw		#U MWCS pointer input and SMW pointer output
+pointer	smw1		#I Template SMW pointer
+pointer	im		#I Template IMIO pointer
+
+int	i, nspec, nmw, format, strdic()
+pointer	sp, sys, smw, mw_sctran(), mw_newcopy()
+errchk	smw_daxis, smw_saxes, mw_sctran
+
+begin
+	call smark (sp)
+	call salloc (sys, SZ_FNAME, TY_CHAR)
+
+	call mw_gwattrs (mw, 0, "system", Memc[sys], SZ_FNAME)
+	format = strdic (Memc[sys], Memc[sys], SZ_FNAME, SMW_FORMATS)
+
+	call calloc (smw, SMW_LEN(1), TY_STRUCT)
+	call malloc (SMW_APID(smw), SZ_LINE, TY_CHAR)
+	SMW_FORMAT(smw) = format
+	SMW_DTYPE(smw) = INDEFI
+	SMW_NMW(smw) = 1
+	SMW_MW(smw,0) = mw
+
+	switch (format) {
+	case SMW_ND:
+	    call smw_daxis (smw, im, INDEFI, INDEFI, INDEFI)
+	    call smw_saxes (smw, smw1, im)
+
+	case SMW_ES:
+	    call smw_saxes (smw, smw1, im)
+
+	    nspec = SMW_NSPEC(smw)
+	    call calloc (SMW_APS(smw), nspec, TY_INT)
+	    call calloc (SMW_BEAMS(smw), nspec, TY_INT)
+	    call calloc (SMW_APLOW(smw), 2*nspec, TY_REAL)
+	    call calloc (SMW_APHIGH(smw), 2*nspec, TY_REAL)
+	    call calloc (SMW_APIDS(smw), nspec, TY_POINTER)
+	    if (SMW_PDIM(smw) > 1)
+		SMW_CTLP(smw) = mw_sctran (mw, "logical", "physical", 2)
+
+	case SMW_MS:
+	    call smw_saxes (smw, smw1, im)
+
+	    nspec = SMW_NSPEC(smw)
+	    call calloc (SMW_APIDS(smw), nspec, TY_POINTER)
+	    if (SMW_PDIM(smw) > 1)
+		SMW_CTLP(smw) = mw_sctran (mw, "logical", "physical", 2)
+
+	    nmw = 1 + (nspec - 1) / SMW_NSPLIT
+	    if (nmw > 1) {
+		call realloc (smw, SMW_LEN(nmw), TY_STRUCT)
+		call calloc (SMW_APS(smw), nspec, TY_INT)
+	    }
+	    do i = 1, nmw-1
+		SMW_MW(smw,i) = mw_newcopy (mw)
+	    SMW_NMW(smw) = nmw
+	}
+
+	mw = smw
+	    
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwopenim.x b/noao/onedspec/smw/smwopenim.x
new file mode 100644
index 00000000..468f09a7
--- /dev/null
+++ b/noao/onedspec/smw/smwopenim.x
@@ -0,0 +1,69 @@
+include <imhdr.h>
+include	<imio.h>
+include	<mwset.h>
+
+define	SYSTEMS	"|equispec|multispec|physical|image|world|linear|"
+
+
+# SMW_OPENIM -- Open the spectral MWCS for various input formats.
+
+pointer procedure smw_openim (im)
+
+pointer	im		#I Image pointer
+pointer	mw		#O MWCS pointer
+
+pointer	sp, system, mw_openim()
+bool	streq()
+int	i, wcsdim, sys, strdic(), mw_stati()
+errchk	mw_openim, smw_oldms, smw_linear
+
+begin
+	call smark (sp)
+	call salloc (system, SZ_FNAME, TY_CHAR)
+
+	# Workaround for truncation of header during image header copy.
+	IM_HDRLEN(im) = IM_LENHDRMEM(im)
+
+	# Force higher dimensions to length 1.
+	do i = IM_NDIM(im) + 1, 3
+	    IM_LEN(im,i) = 1
+
+	mw = mw_openim (im)
+	call mw_seti (mw, MW_USEAXMAP, NO)
+	wcsdim = mw_stati (mw, MW_NDIM)
+	call mw_gwattrs (mw, 0, "system", Memc[system], SZ_FNAME)
+	sys = strdic (Memc[system], Memc[system], SZ_FNAME, SYSTEMS)
+
+	# Set various input systems.
+	switch (sys) {
+	case 1:
+	    call smw_equispec (im, mw)
+	case 2:
+	    call smw_multispec (im, mw)
+	default:
+	    if (sys == 0) {
+	        call eprintf (
+	    "WARNING: Unknown coordinate system `%s' - assuming `linear'.\n")
+		    call pargstr (Memc[system])
+	    } else if (sys == 3)
+		call mw_newsystem (mw, "image", wcsdim)
+
+	    # Old "multispec" format.
+	    ifnoerr (call imgstr (im, "APFORMAT", Memc[system], SZ_FNAME)) {
+		if (streq (Memc[system], "onedspec"))
+		    call smw_onedspec (im, mw)
+		else
+		    call smw_oldms (im, mw)
+
+	    # Old "onedspec" format or other 1D image.
+	    } else if (wcsdim == 1) {
+		call smw_onedspec (im, mw)
+
+	    # N-dimensional image.
+	    } else
+		call smw_nd (im, mw)
+	}
+
+	call sfree (sp)
+	return (mw)
+end
diff --git a/noao/onedspec/smw/smwsapid.x b/noao/onedspec/smw/smwsapid.x
new file mode 100644
index 00000000..1bdf30a8
--- /dev/null
+++ b/noao/onedspec/smw/smwsapid.x
@@ -0,0 +1,40 @@
+include	<smw.h>
+
+
+# SMW_SAPID -- Set aperture id.
+
+procedure smw_sapid (smw, index1, index2, apid)
+
+pointer	smw				#I SMW pointer
+int	index1				#I Spectrum index
+int	index2				#I Spectrum index
+char	apid[ARB]			#I Aperture id
+
+pointer	ptr
+bool	streq()
+errchk	malloc
+
+begin
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    call strcpy (apid, Memc[SMW_APID(smw)], SZ_LINE)
+	case SMW_ES, SMW_MS:
+	    if (index1 < 0 || index1 > SMW_NSPEC(smw))
+		call error (1, "smw_sapid: index out of bounds")
+
+	    if (index1 == 0)
+		call strcpy (apid, Memc[SMW_APID(smw)], SZ_LINE)
+
+	    else {
+		ptr = Memi[SMW_APIDS(smw)+index1-1]
+		if (streq (apid, Memc[SMW_APID(smw)]))
+		    call mfree (ptr, TY_CHAR) 
+		else {
+		    if (ptr == NULL)
+			call malloc (ptr, SZ_LINE, TY_CHAR)
+		    call strcpy (apid, Memc[ptr], SZ_LINE)
+		}
+		Memi[SMW_APIDS(smw)+index1-1] = ptr
+	    }
+	}
+end
diff --git a/noao/onedspec/smw/smwsaveim.x b/noao/onedspec/smw/smwsaveim.x
new file mode 100644
index 00000000..892a3319
--- /dev/null
+++ b/noao/onedspec/smw/smwsaveim.x
@@ -0,0 +1,251 @@
+include <imhdr.h>
+include	<imio.h>
+include	<smw.h>
+
+
+# SMW_SAVEIM -- Save spectral WCS in image header.
+# The input and output formats are EQUISPEC and MULTISPEC.  A split input
+# MULTISPEC WCS is first merged to a single EQUISPEC or MULTISPEC WCS.
+# An input MULTISPEC WCS is converted to EQUISPEC output if possible.
+
+procedure smw_saveim (smw, im)
+
+pointer	smw			# SMW pointer
+pointer	im			# Image pointer
+
+int	i, j, format, nl, pdim, pdim1, beam, dtype, dtype1, nw, nw1
+int	ap, axes[3]
+real	aplow[2], aphigh[2]
+double	v, m, w1, dw, z, w11, dw1, z1
+pointer	sp, key, str1, str2, axmap, lterm, coeff, mw, mw1
+
+bool	strne(), fp_equald()
+int	imaccf(), imgeti()
+pointer	mw_open()
+errchk	smw_merge, imdelf
+data	axes/1,2,3/
+
+begin
+	call smark (sp)
+	call salloc (key, SZ_FNAME, TY_CHAR)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (str2, SZ_LINE, TY_CHAR)
+	call salloc (axmap, 6, TY_INT)
+	call salloc (lterm, 15, TY_DOUBLE)
+	coeff = NULL
+
+	# Merge split WCS into a single WCS.
+	call smw_merge (smw)
+
+	mw = SMW_MW(smw,0)
+	pdim = SMW_PDIM(smw)
+	format = SMW_FORMAT(smw)
+	if (IM_NDIM(im) == 1)
+	    nl = 1
+	else
+	    nl = IM_LEN(im,2)
+
+	# If writing to an existing image we must follow IM_NPHYSDIM
+	# but in a NEW_COPY header we may really want a lower dimension.
+	# Since IM_NPHYSDIM is outside the interface we only violate
+	# it here and use a temporary keyword to communicate from the
+	# routine setting up the WCS.
+
+	pdim1 = max (IM_NDIM(im), IM_NPHYSDIM(im))
+	ifnoerr (i = imgeti (im, "SMW_NDIM")) {
+	    pdim1 = i
+	    call imdelf (im, "SMW_NDIM")
+	}
+
+	# Check if MULTISPEC WCS can be converted to EQUISPEC.
+	if (format == SMW_MS) {
+	    format = SMW_ES
+	    do i = 1, nl {
+		call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		    aplow, aphigh, coeff)
+		if (i == 1) {
+		    dtype1 = dtype
+		    w11 = w1
+		    dw1 = dw
+		    z1 = z
+		    nw1 = nw
+		}
+		if (dtype>1||dtype!=dtype1||!fp_equald(w1,w11)||
+		    !fp_equald(dw,dw1)||nw!=nw1||!fp_equald(z,z1)) {
+		    format = SMW_MS
+		    break
+		}
+	    }
+	}
+
+	# Save WCS in desired format.
+	switch (format) {
+	case SMW_ND:
+	    if (SMW_DTYPE(smw) != -1)
+		call imaddi (im, "DC-FLAG", SMW_DTYPE(smw))
+	    else if (imaccf (im, "DC-FLAG") == YES)
+		call imdelf (im, "DC-FLAG")
+	    if (imaccf (im, "DISPAXIS") == YES)
+		call imaddi (im, "DISPAXIS", SMW_PAXIS(smw,1))
+
+	    call smw_gapid (smw, 1, 1, IM_TITLE(im), SZ_IMTITLE)
+	    call mw_saveim (mw, im)
+
+	case SMW_ES:
+	    # Save aperture information.
+	    do i = 1, nl {
+		call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		    aplow, aphigh, coeff)
+		if (i < 1000)
+		    call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
+		else
+		    call sprintf (Memc[key], SZ_FNAME, "AP%d")
+		    call pargi (i)
+		call sprintf (Memc[str1], SZ_LINE, "%d %d")
+		    call pargi (ap)
+		    call pargi (beam)
+		if (!IS_INDEF(aplow[1]) || !IS_INDEF(aphigh[1])) {
+		    call sprintf (Memc[str2], SZ_LINE, " %.2f %.2f")
+			call pargr (aplow[1])
+			call pargr (aphigh[1])
+		    call strcat (Memc[str2], Memc[str1], SZ_LINE)
+		    if (!IS_INDEF(aplow[2]) || !IS_INDEF(aphigh[2])) {
+			call sprintf (Memc[str2], SZ_LINE, " %.2f %.2f")
+			    call pargr (aplow[2])
+			    call pargr (aphigh[2])
+			call strcat (Memc[str2], Memc[str1], SZ_LINE)
+		    }
+		}
+		call imastr (im, Memc[key], Memc[str1])
+		if (i == 1) {
+		    iferr (call imdelf (im, "APID1"))
+			;
+		}
+		call smw_gapid (smw, i, 1, Memc[str1], SZ_LINE)
+		if (Memc[str1] != EOS) {
+		    if (strne (Memc[str1], IM_TITLE(im))) {
+			if (nl == 1) {
+			    call imastr (im, "MSTITLE", IM_TITLE(im))
+			    call strcpy (Memc[str1], IM_TITLE(im), SZ_IMTITLE)
+			} else {
+			    call sprintf (Memc[key], SZ_FNAME, "APID%d")
+				call pargi (i)
+			    call imastr (im, Memc[key], Memc[str1])
+			}
+		    }
+		}
+	    }
+
+	    # Delete unnecessary aperture information.
+	    do i = nl+1, ARB {
+		if (i < 1000)
+		    call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
+		else
+		    call sprintf (Memc[key], SZ_FNAME, "AP%d")
+		    call pargi (i)
+		iferr (call imdelf (im, Memc[key]))
+		    break
+		call sprintf (Memc[key], SZ_FNAME, "APID%d")
+		    call pargi (i)
+		iferr (call imdelf (im, Memc[key]))
+		    ;
+	    }
+
+	    # Add dispersion parameters to image.
+	    if (dtype != -1)
+		call imaddi (im, "DC-FLAG", dtype)
+	    else if (imaccf (im, "DC-FLAG") == YES)
+		call imdelf (im, "DC-FLAG")
+	    if (nw < IM_LEN(im,1))
+		call imaddi (im, "NP2", nw)
+	    else if (imaccf (im, "NP2") == YES)
+		call imdelf (im, "NP2")
+
+	    # Setup EQUISPEC WCS.
+
+	    mw1 = mw_open (NULL, pdim1)
+	    call mw_newsystem (mw1, "equispec", pdim1)
+	    call mw_swtype (mw1, axes, pdim1, "linear", "")
+	    ifnoerr (call mw_gwattrs (mw, 1, "label", Memc[str1], SZ_LINE))
+		call mw_swattrs (mw1, 1, "label", Memc[str1])
+	    ifnoerr (call mw_gwattrs (mw, 1, "units", Memc[str1], SZ_LINE))
+		call mw_swattrs (mw1, 1, "units", Memc[str1])
+	    ifnoerr (call mw_gwattrs (mw, 1, "units_display", Memc[str1],
+		SZ_LINE))
+		call mw_swattrs (mw1, 1, "units_display", Memc[str1])
+	    call mw_gltermd (mw, Memd[lterm+pdim], Memd[lterm], pdim)
+	    v = Memd[lterm]
+	    m = Memd[lterm+pdim]
+	    call mw_gltermd (mw1, Memd[lterm+pdim1], Memd[lterm], pdim1)
+	    Memd[lterm] = v
+	    Memd[lterm+pdim1] = m
+	    call mw_sltermd (mw1, Memd[lterm+pdim1], Memd[lterm], pdim1)
+	    call mw_gwtermd (mw1, Memd[lterm], Memd[lterm+pdim1],
+		Memd[lterm+2*pdim1], pdim1)
+	    Memd[lterm] = 1.
+	    w1 = w1 / (1 + z)
+	    dw = dw / (1 + z)
+	    if (dtype == DCLOG) {
+		dw = log10 ((w1 + (nw - 1) * dw) / w1) / (nw - 1)
+		w1 = log10 (w1)
+	    }
+	    Memd[lterm+pdim1] = w1
+	    Memd[lterm+2*pdim1] = dw
+	    call mw_swtermd (mw1, Memd[lterm], Memd[lterm+pdim1],
+		Memd[lterm+2*pdim1], pdim1)
+	    call mw_saveim (mw1, im)
+	    call mw_close (mw1)
+
+	case SMW_MS:
+	    # Delete any APNUM keywords.  If there is only one spectrum
+	    # define the axis mapping.
+
+	    do j = 1, ARB {
+		if (j < 1000)
+		    call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
+		else
+		    call sprintf (Memc[key], SZ_FNAME, "AP%d")
+		    call pargi (j)
+		iferr (call imdelf (im, Memc[key]))
+		    break
+	    }
+	    if (IM_NDIM(im) == 1) {
+		call aclri (Memi[axmap], 2*pdim)
+		Memi[axmap] = 1
+		call mw_saxmap (mw, Memi[axmap], Memi[axmap+pdim], pdim)
+	    }
+
+	    # Set aperture ids.
+	    do i = 1, nl {
+		if (i == 1) {
+		    iferr (call imdelf (im, "APID1"))
+			;
+		}
+		call smw_gapid (smw, i, 1, Memc[str1], SZ_LINE)
+		if (Memc[str1] != EOS) {
+		    if (strne (Memc[str1], IM_TITLE(im))) {
+			if (nl == 1) {
+			    call imastr (im, "MSTITLE", IM_TITLE(im))
+			    call strcpy (Memc[str1], IM_TITLE(im), SZ_IMTITLE)
+			} else {
+			    call sprintf (Memc[key], SZ_FNAME, "APID%d")
+				call pargi (i)
+			    call imastr (im, Memc[key], Memc[str1])
+			}
+		    }
+		}
+	    }
+
+	    do i = nl+1, ARB {
+		call sprintf (Memc[key], SZ_FNAME, "APID%d")
+		    call pargi (i)
+		iferr (call imdelf (im, Memc[key]))
+		    break
+	    }
+
+	    call mw_saveim (mw, im)
+	}
+
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwsaxes.x b/noao/onedspec/smw/smwsaxes.x
new file mode 100644
index 00000000..f5e31b63
--- /dev/null
+++ b/noao/onedspec/smw/smwsaxes.x
@@ -0,0 +1,247 @@
+include	<imhdr.h>
+include	<mwset.h>
+include	<smw.h>
+
+
+# SMW_SAXES -- Set axes parameters based on previously set dispersion axis.
+# If the dispersion axis has been excluded for NDSPEC allow another axis to
+# be chosen with a warning.  For EQUISPEC and MULTISPEC require the dispersion
+# to be 1 and also to be present.
+
+procedure smw_saxes (smw, smw1, im)
+
+pointer	smw			#I SMW pointer
+pointer	smw1			#I Template SMW pointer
+pointer	im			#I Template IMIO pointer
+
+int	i, pdim, ldim, paxis, laxis, nw, dtype, nspec
+real	smw_c1tranr()
+double	w1, dw
+pointer	sp, str, axno, axval, r, w, cd, mw, ct, smw_sctran()
+int	mw_stati(), imgeti()
+bool	streq(), fp_equald()
+errchk	smw_sctran
+
+begin
+	# If a template SMW pointer is specified just copy the axes parameters.
+	if (smw1 != NULL) {
+	    call strcpy (Memc[SMW_APID(smw1)], Memc[SMW_APID(smw)], SZ_LINE)
+	    SMW_NSPEC(smw) = SMW_NSPEC(smw1)
+	    SMW_NBANDS(smw) = SMW_NBANDS(smw1)
+	    SMW_TRANS(smw) = SMW_TRANS(smw1)
+	    call amovi (SMW_PAXIS(smw1,1), SMW_PAXIS(smw,1), 3)
+	    SMW_LDIM(smw) = SMW_LDIM(smw1)
+	    call amovi (SMW_LAXIS(smw1,1), SMW_LAXIS(smw,1), 3)
+	    call amovi (SMW_LLEN(smw1,1), SMW_LLEN(smw,1), 3)
+	    call amovi (SMW_NSUM(smw1,1), SMW_NSUM(smw,1), 2)
+
+	    mw = SMW_MW(smw,0)
+	    SMW_PDIM(smw) = mw_stati (mw, MW_NDIM)
+	    if (SMW_PDIM(smw) > SMW_PDIM(smw1))
+	       do i = SMW_PDIM(smw1)+1, SMW_PDIM(smw)
+		   SMW_PAXIS(smw,i) = i
+
+	    return
+	}
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (axno, 3, TY_INT)
+	call salloc (axval, 3, TY_INT)
+	call aclri (Memi[axno], 3)
+
+	# Determine the dimensions.
+	mw = SMW_MW(smw,0)
+	pdim = mw_stati (mw, MW_NDIM)
+	ldim = IM_NDIM(im)
+	call mw_gaxmap (mw, Memi[axno], Memi[axval], pdim)
+
+	# Set the physical dispersion axis.
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    call salloc (r, pdim, TY_DOUBLE)
+	    call salloc (w, pdim, TY_DOUBLE)
+	    call salloc (cd, pdim*pdim, TY_DOUBLE)
+
+	    # Check for a transposed or rotated 2D image.
+	    SMW_TRANS(smw) = NO
+	    if (pdim == 2) {
+		call mw_gltermd (mw, Memd[cd], Memd[w], pdim)
+		if (Memd[cd] == 0D0 && Memd[cd+3] == 0D0) {
+		    Memd[cd] = Memd[cd+1]
+		    Memd[cd+1] = 0.
+		    Memd[cd+3] = Memd[cd+2]
+		    Memd[cd+2] = 0.
+		    call mw_sltermd (mw, Memd[cd], Memd[w], pdim)
+		    paxis = SMW_PAXIS(smw,1)
+		    if (paxis == 1)
+			SMW_PAXIS(smw,1) = 2
+		    else
+			SMW_PAXIS(smw,1) = 1
+		    SMW_TRANS(smw) = YES
+		} else if (Memd[cd+1] != 0D0 || Memd[cd+2] != 0D0) {
+		    Memd[w] = 0
+		    Memd[w+1] = 0
+		    Memd[cd] = 1
+		    Memd[cd+1] = 0
+		    Memd[cd+2] = 0
+		    Memd[cd+3] = 1
+		    call mw_sltermd (mw, Memd[cd], Memd[w], pdim)
+		}
+	    }
+
+	    # If the dispersion axis is of length 1 or has been excluded find
+	    # the first longer axis and print a warning.
+
+	    paxis = SMW_PAXIS(smw,1)
+	    i = max (1, min (pdim, paxis))
+	    laxis = max (1, Memi[axno+i-1])
+	    if (IM_LEN(im,laxis) == 1)
+		do laxis = 1, ldim
+		    if (IM_LEN(im,laxis) != 1)
+			break
+
+	    # Determine the number of spectra.
+	    nspec = 1
+	    do i = 1, ldim
+		if (i != laxis)
+		    nspec = nspec * IM_LEN(im,i)
+	    SMW_NSPEC(smw) = nspec
+	    SMW_NBANDS(smw) = 1
+
+	    i = paxis
+	    do paxis = 1, pdim
+		if (Memi[axno+paxis-1] == laxis)
+		    break
+
+	    if (i != paxis && nspec > 1) {
+		call eprintf (
+		    "WARNING: Dispersion axis %d not found. Using axis %d.\n")
+		call pargi (i)
+		call pargi (paxis)
+	    }
+
+	    # Set the dispersion system.
+	    call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim)
+	    if (SMW_TRANS(smw) == YES) {
+		Memd[cd] = Memd[cd+1]
+		Memd[cd+1] = 0.
+		Memd[cd+3] = Memd[cd+2]
+		Memd[cd+2] = 0.
+	    }
+	    if (pdim == 2 && (Memd[cd+1] != 0D0 || Memd[cd+2] != 0D0)) {
+		iferr (dtype = imgeti (im, "DC-FLAG"))
+		    dtype = DCNO
+		if (dtype != DCNO) {
+		    call sfree (sp)
+		    call error (1,
+		    "Rotated, dispersion calibrated spectra are not allowed")
+		}
+		Memd[r] = 0
+		Memd[r+1] = 0
+		Memd[w] = 0
+		Memd[w+1] = 0
+		Memd[cd] = 1
+		Memd[cd+1] = 0
+		Memd[cd+2] = 0
+		Memd[cd+3] = 1
+	    }
+	    do i = 0, pdim-1 {
+		dw = Memd[cd+i*(pdim+1)]
+		if (dw == 0D0)
+		    Memd[cd+i*(pdim+1)] = 1D0
+	    }
+	    call mw_swtermd (mw, Memd[r], Memd[w], Memd[cd], pdim)
+
+	    dw = Memd[cd+(paxis-1)*(pdim+1)]
+	    w1 = Memd[w+paxis-1] - (Memd[r+paxis-1] - 1) * dw
+	    nw = IM_LEN(im,laxis)
+
+	    i = 2 ** (paxis - 1)
+	    ct = smw_sctran (smw, "logical", "physical", i)
+	    nw = max (smw_c1tranr (ct, 0.5), smw_c1tranr (ct, nw+0.5))
+	    call smw_ctfree (ct)
+
+	    iferr (dtype = imgeti (im, "DC-FLAG")) {
+		iferr (call mw_gwattrs (mw,paxis,"axtype",Memc[str],SZ_LINE))
+		    Memc[str] = EOS
+		if (streq (Memc[str], "ra") || streq (Memc[str], "dec"))
+		    dtype = DCNO
+		else if (fp_equald (1D0, w1) || fp_equald (1D0, dw))
+		    dtype = DCNO
+		else
+		    dtype = DCLINEAR
+	    }
+	    if (dtype==DCLOG) {
+		if (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.)
+		    dtype = DCLINEAR
+		else {
+		    w1 = 10D0 ** w1
+		    dw = w1 * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1)
+		}
+	    }
+
+	    if (dtype != DCNO) {
+		
+		
+		iferr (call mw_gwattrs (mw,paxis,"label",Memc[str],SZ_LINE)) {
+		iferr (call mw_gwattrs(mw,paxis,"units",Memc[str],SZ_LINE)) {
+			call mw_swattrs (mw, paxis, "units", "angstroms")
+			call mw_swattrs (mw, paxis, "label", "Wavelength")
+		    }
+		}
+	    }
+
+	    SMW_DTYPE(smw) = INDEFI
+	    call smw_swattrs (smw, 1, 1, INDEFI, INDEFI,
+		dtype, w1, dw, nw, 0D0, INDEFR, INDEFR, "")
+	case SMW_ES, SMW_MS:
+	    paxis = 1
+	    i = Memi[axno+1]
+	    if (i == 0)
+		SMW_NSPEC(smw) = 1
+	    else
+		SMW_NSPEC(smw) = IM_LEN(im,i)
+	    i = Memi[axno+2]
+	    if (i == 0)
+		SMW_NBANDS(smw) = 1
+	    else
+		SMW_NBANDS(smw) = IM_LEN(im,i)
+	}
+
+	# Check and set the physical and logical dispersion axes.
+	laxis = Memi[axno+paxis-1]
+	if (laxis == 0) {
+	    if (Memi[axval+paxis-1] == 0)
+		laxis = paxis
+	    else
+		call error (1, "No dispersion axis")
+	}
+
+	SMW_PDIM(smw) = pdim
+	SMW_LDIM(smw) = ldim
+	SMW_PAXIS(smw,1) = paxis
+	SMW_LAXIS(smw,1) = laxis
+	SMW_LLEN(smw,1) = IM_LEN(im,laxis)
+	SMW_LLEN(smw,2) = 1
+	SMW_LLEN(smw,3) = 1
+
+	# Set the spatial axes.
+	i = 2
+	do laxis = 1, ldim {
+	    if (laxis != SMW_LAXIS(smw,1)) {
+		do paxis = 1, pdim
+		    if (Memi[axno+paxis-1] == laxis)
+			break
+		SMW_PAXIS(smw,i) = paxis
+		SMW_LAXIS(smw,i) = laxis
+		SMW_LLEN(smw,i) = IM_LEN(im,laxis)
+		i = i + 1
+	    }
+	}
+
+	# Set the default title.
+	call smw_sapid (smw, 0, 1, IM_TITLE(im))
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/smwsctran.x b/noao/onedspec/smw/smwsctran.x
new file mode 100644
index 00000000..06f240db
--- /dev/null
+++ b/noao/onedspec/smw/smwsctran.x
@@ -0,0 +1,96 @@
+include	<error.h>
+include	<smw.h>
+
+
+# SMW_SCTRAN -- Set the SMW coordinate system transformation.
+# This routine sets up the SMW_CT structure which may consist of multiple
+# pointers if the MWCS is split.  If the MWCS is not split then the MWCS
+# transformation routine is used directly.  However if the MWCS is split then
+# there may need to be an intermediate mapping from the input coordinate to
+# the physical coordinate in which the split MWCS is defined as well as a
+# transformation for each WCS piece.
+
+pointer procedure smw_sctran (smw, system1, system2, axbits)
+
+pointer	smw		#I SMW pointer
+char	system1[ARB]	#I Input coordinate system
+char	system2[ARB]	#I Output coordinate system
+int	axbits		#I Bitmap defining axes to be transformed
+pointer	ct		#O SMW CT pointer
+
+int	i, cttype, nct, axes[3], naxes, strdic()
+pointer	mw_sctran()
+errchk	mw_sctran
+
+begin
+	# Determine the coordinate transformation type and setup the structure.
+	cttype = 10 * (strdic (system1, system1, ARB, SMW_CTTYPES) + 1) +
+	    strdic (system2, system2, ARB, SMW_CTTYPES) + 1
+
+	nct = SMW_NMW(smw)
+	if (cttype == SMW_LP || cttype == SMW_PL)
+	    nct = 1
+
+	call calloc (ct, SMW_CTLEN(nct), TY_STRUCT)
+	SMW_SMW(ct) = smw
+	SMW_CTTYPE(ct) = cttype
+	SMW_NCT(ct) = nct
+
+	# Determine dispersion and aperture axes.
+	call mw_gaxlist (SMW_MW(smw,0), axbits, axes, naxes)
+	do i = 1, naxes {
+	    if (axes[i] == SMW_PAXIS(smw,1))
+		SMW_DAXIS(ct) = i
+	    if (axes[i] == SMW_PAXIS(smw,2))
+		SMW_AAXIS(ct) = i
+	}
+
+	# If the MWCS is not split use the MWCS transformation directly.
+	if (nct == 1) {
+	    iferr (SMW_CT(ct,0) = mw_sctran (SMW_MW(smw,0), system1, system2,
+		axbits)) {
+		switch (cttype) {
+		case SMW_WL, SMW_WP:
+	            SMW_CT(ct,0) = mw_sctran (SMW_MW(smw,0), "physical",
+			system2, axbits)
+		case SMW_LW, SMW_PW:
+	            SMW_CT(ct,0) = mw_sctran (SMW_MW(smw,0), system1,
+			"physical", axbits)
+		default:
+		    call erract (EA_ERROR)
+		}
+	    }
+	    return(ct)
+	}
+
+	# If there is a split MWCS then setup the intermediary transformation.
+	switch (cttype) {
+	case SMW_LW:
+	    SMW_CTL(ct) = mw_sctran (SMW_MW(smw,0), system1, "physical",
+		axbits)
+	    do i = 0, nct-1 {
+		iferr (SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical",
+		    system2, axbits))
+		    SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical",
+			"physical", axbits)
+	    }
+	case SMW_WL:
+	    do i = 0, nct-1 {
+		iferr (SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), system1,
+		    "physical", axbits))
+		    SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical",
+			"physical", axbits)
+	    }
+	    SMW_CTL(ct) = mw_sctran (SMW_MW(smw,0), "physical", system2,
+		axbits)
+	case SMW_PW, SMW_WP:
+	    do i = 0, nct-1 {
+		iferr (SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), system1,
+		    system2, axbits))
+		    SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical",
+			system2, axbits)
+	    }
+	}
+
+	return (ct)
+end
diff --git a/noao/onedspec/smw/smwsmw.x b/noao/onedspec/smw/smwsmw.x
new file mode 100644
index 00000000..c3870c4a
--- /dev/null
+++ b/noao/onedspec/smw/smwsmw.x
@@ -0,0 +1,21 @@
+include	<smw.h>
+
+
+# SMW_SMW -- Set MCWS pointer
+
+procedure smw_smw (smw, line, mw)
+
+pointer	smw		#I SMW pointer
+int	line		#I Physical line
+pointer	mw		#I MWCS pointer
+
+begin
+	if (SMW_NMW(smw) == 1)
+	    SMW_MW(smw,0) = mw
+
+	else {
+	    if (line < 1 || line > SMW_NSPEC(smw))
+		call error (1, "smw_smw: aperture not found")
+	    SMW_MW(smw,(line-1)/SMW_NSPLIT) = mw
+	}
+end
diff --git a/noao/onedspec/smw/smwswattrs.x b/noao/onedspec/smw/smwswattrs.x
new file mode 100644
index 00000000..ff859cfc
--- /dev/null
+++ b/noao/onedspec/smw/smwswattrs.x
@@ -0,0 +1,162 @@
+include	<error.h>
+include	<smw.h>
+
+
+# SMW_SWATTRS -- Set spectrum attribute parameters
+# This routine has the feature that if the coordinate system of a single
+# spectrum in an EQUISPEC WCS is changed then the image WCS is changed
+# to a MULTISPEC WCS.
+
+procedure smw_swattrs (smw, index1, index2, ap, beam, dtype, w1, dw, nw, z,
+	aplow, aphigh, coeff)
+
+pointer	smw				# SMW pointer
+int	index1				# Spectrum index
+int	index2				# Spectrum index
+int	ap				# Aperture number
+int	beam				# Beam number
+int	dtype				# Dispersion type
+double	w1				# Starting coordinate
+double	dw				# Coordinate interval
+int	nw				# Number of valid pixels
+double	z				# Redshift factor
+real	aplow[2], aphigh[2]		# Aperture limits
+char	coeff[ARB]			# Nonlinear coeff string
+
+bool	fp_equald()
+int	i, j, sz_val, strlen()
+double	a, b
+pointer	sp, str, val, mw
+errchk	smw_mw
+
+define	start_	10
+
+begin
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+start_
+	switch (SMW_FORMAT(smw)) {
+	case SMW_ND:
+	    if (!IS_INDEFI(SMW_DTYPE(smw)) && (!fp_equald(w1,SMW_W1(smw)) ||
+		!fp_equald(dw,SMW_DW(smw)) || !fp_equald(z,SMW_Z(smw)))) {
+		call malloc (val, 15, TY_DOUBLE)
+		mw = SMW_MW(smw,0)
+		i = SMW_PDIM(smw)
+		j = SMW_PAXIS(smw,1)
+		call mw_gwtermd (mw, Memd[val], Memd[val+i], Memd[val+2*i], i)
+		Memd[val+j-1] = 1.
+		switch (dtype) {
+		case DCNO, DCLINEAR:
+		    a = w1 / (1 + z)
+		    b = dw / (1 + z)
+		case DCLOG:
+		    a = log10 (w1 / (1 + z))
+		    b = log10 ((w1 + (nw - 1) * dw) / w1) / (nw - 1)
+		case DCFUNC:
+		    call error (1,
+		       "Nonlinear functions not allowed for NSPEC format")
+		}
+		Memd[val+i+j-1] = a
+		Memd[val+2*i+(i+1)*(j-1)] = b
+		call mw_swtermd (mw, Memd[val], Memd[val+i], Memd[val+2*i], i)
+		call mfree (val, TY_DOUBLE)
+	    }
+	    SMW_DTYPE(smw) = dtype
+	    SMW_NW(smw) = nw
+	    SMW_W1(smw) = w1
+	    SMW_DW(smw) = dw
+	    SMW_Z(smw) = z
+
+	case SMW_ES:
+	    # Check for any changes to the dispersion system.
+	    if (dtype == DCFUNC) {
+		call smw_esms(smw)
+		goto start_
+	    }
+	    if (!IS_INDEFI(SMW_DTYPE(smw)) && (dtype != SMW_DTYPE(smw) ||
+		nw != SMW_NW(smw) || !fp_equald(w1,SMW_W1(smw)) ||
+		!fp_equald(dw,SMW_DW(smw)) || !fp_equald(z,SMW_Z(smw)))) {
+		if (SMW_NSPEC(smw) > 1 && index1 > 0) {
+		    call smw_esms(smw)
+		    goto start_
+		}
+		if (!fp_equald(w1,SMW_W1(smw)) || !fp_equald(dw,SMW_DW(smw)) ||
+		    !fp_equald(z,SMW_Z(smw))) {
+		    call malloc (val, 15, TY_DOUBLE)
+		    mw = SMW_MW(smw,0)
+		    i = SMW_PDIM(smw)
+		    j = SMW_PAXIS(smw,1)
+		    call mw_gwtermd (mw, Memd[val], Memd[val+i],
+			Memd[val+2*i], i)
+		    Memd[val+j-1] = 1.
+		    switch (dtype) {
+		    case DCNO, DCLINEAR:
+			a = w1 / (1 + z)
+			b = dw / (1 + z)
+		    case DCLOG:
+			a = log10 (w1 / (1 + z))
+			b = log10 ((w1 + (nw - 1) * dw) / w1) / (nw - 1)
+		    }
+		    Memd[val+i+j-1] = a
+		    Memd[val+2*i+(i+1)*(j-1)] = b
+		    call mw_swtermd (mw, Memd[val], Memd[val+i],
+			Memd[val+2*i], i)
+		    call mfree (val, TY_DOUBLE)
+		}
+	    }
+
+	    SMW_DTYPE(smw) = dtype
+	    SMW_NW(smw) = nw
+	    SMW_W1(smw) = w1
+	    SMW_DW(smw) = dw
+	    SMW_Z(smw) = z
+
+	    if (index1 > 0) {
+		Memi[SMW_APS(smw)+index1-1] = ap
+		Memi[SMW_BEAMS(smw)+index1-1] = beam
+		Memr[SMW_APLOW(smw)+2*index1-2] = aplow[1]
+		Memr[SMW_APHIGH(smw)+2*index1-2] = aphigh[1]
+		Memr[SMW_APLOW(smw)+2*index1-1] = aplow[2]
+		Memr[SMW_APHIGH(smw)+2*index1-1] = aphigh[2]
+	    }
+
+	case SMW_MS:
+	    # We can't use SPRINTF for the whole string because it can only
+	    # handle a limited length and trucates long coefficient strings.
+	    # Use STRCAT instead.
+
+	    call smw_mw (smw, index1, index2, mw, i, j)
+	    sz_val = strlen (coeff) + SZ_LINE
+	    call salloc (val, sz_val, TY_CHAR)
+	    call sprintf (Memc[str], SZ_LINE, "spec%d")
+		call pargi (i)
+	    call sprintf (Memc[val], sz_val,
+		"%d %d %d %.14g %.14g %d %.14g %.2f %.2f")
+		call pargi (ap)
+		call pargi (beam)
+		call pargi (dtype)
+		if (dtype == DCLOG) {
+		    call pargd (log10 (w1))
+		    call pargd (log10 ((w1+(nw-1)*dw)/w1)/(nw-1))
+		} else {
+		    call pargd (w1)
+		    call pargd (dw)
+		}
+		call pargi (nw)
+		call pargd (z)
+		call pargr (aplow[1])
+		call pargr (aphigh[1])
+	    if (coeff[1] != EOS) {
+		call strcat (" ", Memc[val], sz_val)
+		call strcat (coeff, Memc[val], sz_val)
+	    }
+	    call mw_swattrs (mw, 2, Memc[str], Memc[val])
+
+	    if (SMW_APS(smw) != NULL)
+		Memi[SMW_APS(smw)+index1-1] = ap
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/smw/units.x b/noao/onedspec/smw/units.x
new file mode 100644
index 00000000..f44abb57
--- /dev/null
+++ b/noao/onedspec/smw/units.x
@@ -0,0 +1,529 @@
+include	<ctype.h>
+include	<error.h>
+include	<units.h>
+
+
+# UN_OPEN -- Open units package
+# It is allowed to open an unknown unit type
+
+pointer procedure un_open (units)
+
+char	units[ARB]		# Units string
+pointer	un			# Units pointer returned
+
+begin
+	call calloc (un, UN_LEN, TY_STRUCT)
+	iferr (call un_decode (un, units)) {
+	    call un_close (un)
+	    call erract (EA_ERROR)
+	}
+	return (un)
+end
+
+
+# UN_CLOSE -- Close units package
+
+procedure un_close (un)
+
+pointer	un			# Units pointer
+
+begin
+	call mfree (un, TY_STRUCT)
+end
+
+
+# UN_COPY -- Copy units pointer
+
+procedure un_copy (un1, un2)
+
+pointer	un1, un2		# Units pointers
+
+begin
+	if (un2 == NULL)
+	    call malloc (un2, UN_LEN, TY_STRUCT)
+	call amovi (Memi[un1], Memi[un2], UN_LEN)
+end
+
+
+# UN_DECODE -- Decode units string and set up units structure.
+# The main work is done in UN_DECODE1 so that the units string may
+# be recursive; i.e. the units string may contain other units strings.
+# In particular, this is required for the velocity units to specify
+# a reference wavelength.
+
+procedure un_decode (un, units)
+
+pointer	un			# Units pointer
+char	units[ARB]		# Units string
+
+bool	streq()
+pointer	sp, units1, temp, un1, un2
+errchk	un_decode1, un_ctranr
+
+begin
+	if (streq (units, UN_USER(un)))
+	    return
+
+	call smark (sp)
+	call salloc (units1, SZ_LINE, TY_CHAR)
+	call salloc (temp, UN_LEN, TY_STRUCT)
+
+	# Save a copy to restore in case of an error.
+	call un_copy (un, temp)
+
+	iferr {
+	    # Decode the primary units
+	    call un_decode1 (un, units, Memc[units1], SZ_LINE)
+
+	    # Decode velocity reference wavelength if necessary.
+	    if (UN_CLASS(un) == UN_VEL || UN_CLASS(un) == UN_DOP) {
+		call salloc (un1, UN_LEN, TY_STRUCT)
+		call un_decode1 (un1, Memc[units1], Memc[units1], SZ_LINE)
+		if (UN_CLASS(un1) == UN_VEL || UN_CLASS(un1) == UN_DOP)
+		    call error (1,
+			"Velocity reference units may not be velocity")
+		call salloc (un2, UN_LEN, TY_STRUCT)
+		call un_decode1 (un2, "angstroms", Memc[units1], SZ_LINE)
+		call un_ctranr (un1, un2, UN_VREF(un), UN_VREF(un), 1)
+	    }
+	} then {
+	    call un_copy (temp, un)
+	    call sfree (sp)
+	    call erract (EA_ERROR)
+	}
+
+	call sfree (sp)
+end
+
+
+# UN_DECODE1 -- Decode units string and set up units structure.
+# Return any secondary units string.  Unknown unit strings are allowed.
+
+procedure un_decode1 (un, units, units1, sz_units1)
+
+pointer	un			# Units pointer
+char	units[ARB]		# Units string
+char	units1[sz_units1]	# Secondary units string to return
+int	sz_units1		# Size of secondary units string
+
+int	unlog, uninv, untype
+int	i, j, k, nscan(), strdic(), strlen()
+pointer	sp, str
+pointer	stp, sym, stfind(), strefsbuf()
+
+int	class[UN_NUNITS]
+real	scale[UN_NUNITS]
+data	stp/NULL/
+data	class /UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,
+		UN_FREQ,UN_FREQ,UN_FREQ,UN_FREQ,UN_VEL,UN_VEL,
+		UN_ENERGY,UN_ENERGY,UN_ENERGY,UN_DOP/
+data	scale /UN_ANG,UN_NM,UN_MMIC,UN_MIC,UN_MM,UN_CM,UN_M,UN_HZ,UN_KHZ,
+		UN_MHZ,UN_GHZ,UN_MPS,UN_KPS,UN_EV,UN_KEV,UN_MEV,UN_Z/
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	iferr (call un_abbr (stp))
+	    ;
+
+	call strcpy (units, Memc[str], SZ_FNAME)
+	if (stp != NULL) {
+	    sym = stfind (stp, Memc[str])
+	    if (sym != NULL)
+		call strcpy (Memc[strefsbuf(stp,Memi[sym])],
+		    Memc[str], SZ_FNAME)
+	}
+	call strlwr (Memc[str])
+	call sscan (Memc[str])
+	untype = 0
+	unlog = NO
+	uninv = NO
+	do i = 1, 3 {
+	    call gargwrd (Memc[str], SZ_FNAME)
+	    if (nscan() != i)
+		break
+
+	    j = strdic (Memc[str], Memc[str], SZ_FNAME, UN_DIC)
+	    for (k=strlen(Memc[str]); k>0 &&
+		(IS_WHITE(Memc[str+k-1]) || Memc[str+k-1]=='\n'); k=k-1)
+		Memc[str+k-1] = EOS
+
+	    if (j > UN_NUNITS) {
+	        j = j - UN_NUNITS
+	        if (j == 1) {
+		    if (unlog == YES)
+			break
+		    unlog = YES
+		} else if (j == 2) {
+		    if (uninv == YES)
+			break
+		    uninv = YES
+		}
+	    } else {
+		if (class[j] == UN_VEL || class[j] == UN_DOP) {
+		    call gargr (UN_VREF(un))
+		    call gargstr (units1, sz_units1)
+		    if (nscan() != i+2)
+			call error (1, "Error in velocity reference wavelength")
+		} else
+		    UN_VREF(un) = 0.
+		untype = j
+		break
+	    }
+	}
+
+	if (untype == 0) {
+	    UN_TYPE(un) = 0
+	    UN_CLASS(un) = UN_UNKNOWN
+	    UN_LABEL(un) = EOS
+	    call strcpy (units, UN_UNITS(un), SZ_UNITS)
+	} else {
+	    UN_TYPE(un) = untype
+	    UN_CLASS(un) = class[untype]
+	    UN_LOG(un) = unlog
+	    UN_INV(un) = uninv
+	    UN_SCALE(un) = scale[untype]
+	    UN_LABEL(un) = EOS
+	    UN_UNITS(un) = EOS
+	    call strcpy (units, UN_USER(un), SZ_UNITS)
+
+	    if (unlog == YES)
+		call strcat ("Log ", UN_LABEL(un), SZ_UNITS)
+	    if (uninv == YES)
+		call strcat ("inverse ", UN_UNITS(un), SZ_UNITS)
+	    call strcat (Memc[str], UN_UNITS(un), SZ_UNITS)
+	    switch (class[j]) {
+	    case UN_WAVE:
+		if (uninv == NO)
+		    call strcat ("Wavelength", UN_LABEL(un), SZ_UNITS)
+		else
+		    call strcat ("Wavenumber", UN_LABEL(un), SZ_UNITS)
+	    case UN_FREQ:
+		call strcat ("Frequency", UN_LABEL(un), SZ_UNITS)
+	    case UN_VEL:
+		call strcat ("Velocity", UN_LABEL(un), SZ_UNITS)
+	    case UN_ENERGY:
+		call strcat ("Energy", UN_LABEL(un), SZ_UNITS)
+	    case UN_DOP:
+		call strcat ("Redshift", UN_LABEL(un), SZ_UNITS)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# UN_COMPARE -- Compare two units
+
+bool procedure un_compare (un1, un2)
+
+pointer	un1, un2		# Units pointers to compare
+bool	strne()
+
+begin
+	if (strne (UN_UNITS(un1), UN_UNITS(un2)))
+	    return (false)
+	if (strne (UN_LABEL(un1), UN_LABEL(un2)))
+	    return (false)
+	if (UN_VREF(un1) != UN_VREF(un2))
+	    return (false)
+	return (true)
+end
+
+
+# UN_CTRANR -- Transform units
+# Error is returned if the transform cannot be made
+
+procedure un_ctranr (un1, un2, val1, val2, nvals)
+
+pointer	un1			# Input units pointer
+pointer	un2			# Output units pointer
+real	val1[nvals]		# Input values
+real	val2[nvals]		# Output values
+int	nvals			# Number of values
+
+int	i
+real	s, v, z
+bool	un_compare()
+
+begin
+	if (un_compare (un1, un2)) {
+	    call amovr (val1, val2, nvals)
+	    return
+	}
+
+	if (UN_CLASS(un1) == UN_UNKNOWN || UN_CLASS(un2) == UN_UNKNOWN)
+	    call error (1, "Cannot convert between selected units")
+
+	call amovr (val1, val2, nvals)
+
+	s = UN_SCALE(un1)
+	if (UN_LOG(un1) == YES)
+	    do i = 1, nvals
+		val2[i] = 10. ** val2[i]
+	if (UN_INV(un1) == YES)
+	    do i = 1, nvals
+		val2[i] = 1. / val2[i]
+	switch (UN_CLASS(un1)) {
+	case UN_WAVE:
+	    do i = 1, nvals
+		val2[i] = val2[i] / s
+	case UN_FREQ:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_VEL:
+	    v = UN_VREF(un1)
+	    do i = 1, nvals {
+		z = val2[i] / s
+		val2[i] = sqrt ((1 + z) / (1 - z)) * v
+	    }
+	case UN_ENERGY:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_DOP:
+	    v = UN_VREF(un1)
+	    do i = 1, nvals
+		val2[i] = (val2[i] / s + 1) * v
+	}
+
+	s = UN_SCALE(un2)
+	switch (UN_CLASS(un2)) {
+	case UN_WAVE:
+	    do i = 1, nvals
+		val2[i] = val2[i] * s
+	case UN_FREQ:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_VEL:
+	    v = UN_VREF(un2)
+	    do i = 1, nvals {
+		z = (val2[i] / v) ** 2
+		val2[i] = (z - 1) / (z + 1) * s
+	    }
+	case UN_ENERGY:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_DOP:
+	    v = UN_VREF(un2)
+	    do i = 1, nvals
+		val2[i] = (val2[i] / v - 1) * s
+	}
+	if (UN_INV(un2) == YES)
+	    do i = 1, nvals
+		val2[i] = 1. / val2[i]
+	if (UN_LOG(un2) == YES)
+	    do i = 1, nvals
+		val2[i] = log10 (val2[i])
+end
+
+
+# UN_CHANGER -- Change units
+# Error is returned if the conversion cannot be made
+
+procedure un_changer (un, units, vals, nvals, update)
+
+pointer	un			# Units pointer (may be changed)
+char	units[ARB]		# Desired units
+real	vals[nvals]		# Values
+int	nvals			# Number of values
+int	update			# Update units pointer?
+
+bool	streq(), un_compare()
+pointer	un1, un_open()
+errchk	un_open, un_ctranr
+
+begin
+
+	# Check for same unit string
+	if (streq (units, UN_USER(un)))
+	    return
+
+	# Check for error in units string, or the same units.
+	un1 = un_open (units)
+	if (un_compare (un1, un)) {
+	    call strcpy (units, UN_USER(un), SZ_UNITS)
+	    call un_close (un1)
+	    return
+	}
+
+	iferr {
+	    call un_ctranr (un, un1, vals, vals, nvals)
+	    if (update == YES)
+		call un_copy (un1, un)
+	    call un_close(un1)
+	} then {
+	    call un_close(un1)
+	    call erract (EA_ERROR)
+	}
+end
+
+
+# UN_CTRAND -- Transform units
+# Error is returned if the transform cannot be made
+
+procedure un_ctrand (un1, un2, val1, val2, nvals)
+
+pointer	un1			# Input units pointer
+pointer	un2			# Output units pointer
+double	val1[nvals]		# Input values
+double	val2[nvals]		# Output values
+int	nvals			# Number of values
+
+int	i
+double	s, v, z
+bool	un_compare()
+
+begin
+	if (un_compare (un1, un2)) {
+	    call amovd (val1, val2, nvals)
+	    return
+	}
+
+	if (UN_CLASS(un1) == UN_UNKNOWN || UN_CLASS(un2) == UN_UNKNOWN)
+	    call error (1, "Cannot convert between selected units")
+
+	call amovd (val1, val2, nvals)
+
+	s = UN_SCALE(un1)
+	if (UN_LOG(un1) == YES)
+	    do i = 1, nvals
+		val2[i] = 10. ** val2[i]
+	if (UN_INV(un1) == YES)
+	    do i = 1, nvals
+		val2[i] = 1. / val2[i]
+	switch (UN_CLASS(un1)) {
+	case UN_WAVE:
+	    do i = 1, nvals
+		val2[i] = val2[i] / s
+	case UN_FREQ:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_VEL:
+	    v = UN_VREF(un1)
+	    do i = 1, nvals {
+		z = val2[i] / s
+		val2[i] = sqrt ((1 + z) / (1 - z)) * v
+	    }
+	case UN_ENERGY:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_DOP:
+	    v = UN_VREF(un1)
+	    do i = 1, nvals
+		val2[i] = (val2[i] / s + 1) * v
+	}
+
+	s = UN_SCALE(un2)
+	switch (UN_CLASS(un2)) {
+	case UN_WAVE:
+	    do i = 1, nvals
+		val2[i] = val2[i] * s
+	case UN_FREQ:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_VEL:
+	    v = UN_VREF(un2)
+	    do i = 1, nvals {
+		z = (val2[i] / v) ** 2
+		val2[i] = (z - 1) / (z + 1) * s
+	    }
+	case UN_ENERGY:
+	    do i = 1, nvals
+		val2[i] = s / val2[i]
+	case UN_DOP:
+	    v = UN_VREF(un2)
+	    do i = 1, nvals
+		val2[i] = (val2[i] / v - 1) * s
+	}
+	if (UN_INV(un2) == YES)
+	    do i = 1, nvals
+		val2[i] = 1. / val2[i]
+	if (UN_LOG(un2) == YES)
+	    do i = 1, nvals
+		val2[i] = log10 (val2[i])
+end
+
+
+# UN_CHANGED -- Change units
+# Error is returned if the conversion cannot be made
+
+procedure un_changed (un, units, vals, nvals, update)
+
+pointer	un			# Units pointer (may be changed)
+char	units[ARB]		# Desired units
+double	vals[nvals]		# Values
+int	nvals			# Number of values
+int	update			# Update units pointer?
+
+bool	streq(), un_compare()
+pointer	un1, un_open()
+errchk	un_open, un_ctrand
+
+begin
+
+	# Check for same unit string
+	if (streq (units, UN_USER(un)))
+	    return
+
+	# Check for error in units string, or the same units.
+	un1 = un_open (units)
+	if (un_compare (un1, un)) {
+	    call strcpy (units, UN_USER(un), SZ_UNITS)
+	    call un_close (un1)
+	    return
+	}
+
+	iferr {
+	    call un_ctrand (un, un1, vals, vals, nvals)
+	    if (update == YES)
+		call un_copy (un1, un)
+	    call un_close(un1)
+	} then {
+	    call un_close(un1)
+	    call erract (EA_ERROR)
+	}
+end
+
+
+# UN_ABBR -- Load abbreviations into a symbol table.
+
+procedure un_abbr (stp)
+
+pointer	stp		#U Symbol table
+
+int	fd, open(), fscan(), nscan(), stpstr()
+pointer	sp, key, val
+pointer	sym, stopen(), stfind(), stenter(), strefsbuf()
+errchk	open
+
+begin
+	if (stp != NULL)
+	    return
+
+	fd = open (ABBREVIATIONS, READ_ONLY, TEXT_FILE)
+	stp = stopen ("unabbr", 20, 20, 40*SZ_LINE)
+
+	call smark (sp)
+	call salloc (key, SZ_LINE, TY_CHAR)
+	call salloc (val, SZ_LINE, TY_CHAR)
+
+	while (fscan (fd) != EOF) {
+	    call gargwrd (Memc[key], SZ_LINE)
+	    call gargwrd (Memc[val], SZ_LINE)
+	    if (nscan() != 2)
+		next
+	    if (Memc[key] == '#')
+		next
+
+	    sym = stfind (stp, Memc[key])
+	    if (sym == NULL) {
+		sym = stenter (stp, Memc[key], 1)
+		Memi[sym] = stpstr (stp, Memc[val], SZ_LINE)
+	    } else
+		call strcpy (Memc[val], Memc[strefsbuf(stp,Memi[sym])], SZ_LINE)
+	}
+
+	call close (fd)
+	call sfree (sp)
+end