Repatch (from linux) of OSX IRAF

1 files changed, 1834 insertions, 0 deletions

diff --git a/noao/twodspec/apextract/apextract.x b/noao/twodspec/apextract/apextract.x
new file mode 100644
index 00000000..176ab251
--- /dev/null
+++ b/noao/twodspec/apextract/apextract.x
@@ -0,0 +1,1834 @@
+include	<error.h>
+include	<imhdr.h>
+include	<mach.h>
+include	<math/iminterp.h>
+include	<pkg/gtools.h>
+include	"apertures.h"
+
+# Background fitting types
+define	BACKGROUND	"|none|average|median|minimum|fit|"
+define	B_NONE		1
+define	B_AVERAGE	2
+define	B_MEDIAN	3
+define	B_MINIMUM	4
+define	B_FIT		5
+
+# Weight types
+define	WEIGHTS		"|none|variance|"
+define	W_NONE		1
+define	W_VARIANCE	2
+
+# Profile fitting algorithms
+define	P_FIT		"|fit1d|fit2d|"
+define	P_FIT1D		1
+define	P_FIT2D		2
+
+# Output formats
+define	FORMATS		"|onedspec|multispec|echelle|strip|normalize|flatten\
+			|ratio|difference|fit|noise|"
+define	ONEDSPEC	1	# Individual 1D spectra
+define	MULTISPEC	2	# Multiple spectra
+define	ECHELLE		3	# Echelle spectra
+define	STRIP		4	# Strip spectra
+define	NORM		5	# Normalized spectra
+define	FLAT		6	# Flat spectra
+define	RATIO		7	# Ratio of data to model
+define	DIFF		8	# Difference of data and model
+define	FIT		9	# Model
+define	NOISE		10	# Noise calculation
+
+
+# AP_EXTRACT -- Extract spectra by a weighted sum across the apertures.
+# 
+# This routine does clobber checks on the output images, manages the I/O
+# from the input image in as big of pieces as possible, and loops through
+# each aperture calling routines to determine the sky, do any fitting and
+# extraction, and output the spectra.
+# The extraction may be either a simple, unweighted extraction
+# which is very fast or a weighted extraction using CCD noise
+# parameters.  The weights require dividing out the basic spectrum and
+# smoothing the 2D spectral profile.  The general approach of variance
+# weighting is described by K. Horne (PASP V98, P609, 1986).  The
+# smoothing has two algorithms, fitting columns or lines parallel to the
+# dispersion axis for nearly aligned spectra or fitting a 2D function
+# using a method given by T. Marsh (PASP V101, P1032, 1989).  The profile
+# may also be used to reject cosmic rays by iteration.
+#
+# The extractions require enough memory to get at least one aperture plus
+# background (if needed) into memory.  If possible the region containing
+# all the apertures is read into memory.  The target maximum amount of
+# memory is set by the maxmimum size returned by BEGMEM and the
+# appropriate working set size is requested.  The optimal size can be
+# tuned through BEGMEM, which references a machine dependent include
+# file, if needed.  The algorithm should work well (minimize I/O as well
+# as paging) in all cases but very large image formats with highly tilted
+# spectra (where aperture extraction along the image axes is not really
+# appropriate).  These memory requirements were chosen to minimize image
+# I/O and because the variance weighted algorithms need to make multiple
+# passes through the image.  In principle simple, unweighted extractions
+# with no sky smoothing can be done sequentially but this was not done in
+# order to use nearly the same code for both weighted and unweighted
+# cases.
+#
+# If using variance weighting and a profile image is given then it is used
+# to determine the profile which is then applied to the target image
+# during the final extraction.  If the same profile image is used multiple
+# times it would be more efficient to store the profile but then issues
+# of consistency arise.  For now this possible feature is not implemented.
+
+procedure ap_extract (input, output, format, profiles, aps, naps)
+
+char	input[SZ_FNAME]		# Input image
+char	output[SZ_FNAME]	# Output image (optional root name)
+char	format[SZ_LINE]		# Output format
+char	profiles[SZ_FNAME]	# Profile filename (optional)
+pointer	aps[ARB]		# Apertures
+int	naps			# Number of apertures
+
+# CL parameters
+int	fmt			# Output format
+int	bkg			# Background type
+int	weights			# Extraction weights
+int	pfit			# Profile fitting algorithm
+bool	clean			# Reject cosmic rays?
+real	gain			# Photon/DN gain
+real	rdnoise			# Read out noise
+int	nsubaps			# Number of subapertures
+int	interptype		# Edge interpolation type
+
+int	i, j, k, napsex, aaxis, baxis, namax, na, nb, na1, interpbuf
+int	amin, amax, bmin, bmax
+int	new_size, old_size, max_size, best_size
+real	cmin, cmax, xmin, xmax, shift
+pointer	sp, str, bkgstr, wtstr, cleanstr, apsex
+pointer	a, b, c, astart, spec, specsky, specsig, raw, profile
+pointer	a1, a2, b1, b2, c1, c2, im, pim, ap, cv, ic, dbuf, pbuf, sbuf, svar, ptr
+pointer	asi
+
+bool	clgetb(), apgetb(), strne()
+int	apgeti(), apgwrd(), begmem(), ap_check()
+real	apgimr(), ap_cveval(), ic_getr()
+pointer	ap_immap(), imgs2r(), imgl2r()
+errchk	salloc, malloc, ap_immap, imgs2r, imgl2r, asiinit
+errchk	ap_check, ap_skyeval, ap_profile, ap_variance, ap_output, apgimr
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	napsex = 0
+	do i = 1, naps
+	    if (AP_SELECT(aps[i]) == YES)
+		napsex = napsex + 1
+	if (napsex == 0) {
+	    call sprintf (Memc[str], SZ_LINE,
+	        "EXTRACT - No apertures defined for %s")
+		call pargstr (input)
+	    call ap_log (Memc[str], YES, NO, YES)
+	    call sfree (sp)
+	    return
+	}
+
+	call salloc (bkgstr, SZ_FNAME, TY_CHAR)
+	call salloc (wtstr, SZ_FNAME, TY_CHAR)
+	call salloc (cleanstr, SZ_FNAME, TY_CHAR)
+	call salloc (apsex, napsex, TY_POINTER)
+
+	# Select apertures to extract and fix possible limit error.
+	napsex = 0
+	do i = 1, naps {
+	    if (AP_LOW(aps[i],1) > AP_HIGH(aps[i],1)) {
+		xmax = AP_LOW(aps[i],1)
+		AP_LOW(aps[i],1) = AP_HIGH(aps[i],1)
+		AP_HIGH(aps[i],1) = xmax
+	    }
+	    if (AP_LOW(aps[i],2) > AP_HIGH(aps[i],2)) {
+		xmax = AP_LOW(aps[i],2)
+		AP_LOW(aps[i],2) = AP_HIGH(aps[i],2)
+		AP_HIGH(aps[i],2) = xmax
+	    }
+	    if (AP_SELECT(aps[i]) == NO)
+		next
+	    Memi[apsex+napsex] = aps[i]
+	    napsex = napsex + 1
+	}
+
+	# Get CL parameters
+	bkg = apgwrd ("background", Memc[bkgstr], SZ_FNAME, BACKGROUND)
+	pfit = apgwrd ("pfit", Memc[str], SZ_LINE, P_FIT)
+	clean = apgetb ("clean")
+	if (clean)
+	    call strcpy ("yes", Memc[cleanstr], SZ_FNAME)
+	else
+	    call strcpy ("no", Memc[cleanstr], SZ_FNAME)
+	nsubaps = apgeti ("nsubaps")
+	interptype = II_LINEAR
+
+	# Do clobber checking.  Return if output exists and not clobbering.
+	call apgstr ("ansclobber", Memc[str], SZ_LINE)
+	call appstr ("ansclobber1", Memc[str])
+	fmt = ap_check (input, output, format, Memi[apsex], napsex, nsubaps)
+	if (fmt == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Force weights depending on format or cleaning.
+	switch (fmt) {
+	case FLAT, RATIO, DIFF, FIT, NOISE:
+	    weights = W_VARIANCE
+	default:
+	    if (clean) {
+		call strcpy ("variance", Memc[wtstr], SZ_FNAME)
+		weights = W_VARIANCE
+	    } else
+	        weights = apgwrd ("weights", Memc[wtstr], SZ_FNAME, WEIGHTS)
+	}
+
+	if (clgetb ("verbose")) {
+	    call printf ("Extracting apertures ...\n")
+	    call flush (STDOUT)
+	}
+	    
+	# Open input image and profile image if given.  Set axis parameters
+	# where 'a' is the aperture axis across the dispersion and 'b' is
+	# along the dispersion.
+
+	im = ap_immap (input, aaxis, baxis)
+	namax = IM_LEN(im, aaxis)
+	nb = IM_LEN(im, baxis)
+
+	pim = NULL
+	if (strne(profiles,input) && weights==W_VARIANCE && profiles[1]!=EOS) {
+	    pim = ap_immap (profiles, i, j)
+	    if (i!=aaxis||j!=baxis||IM_LEN(pim,i)!=namax||IM_LEN(pim,j)!=nb) {
+		call imunmap (pim)
+		call imunmap (im)
+		call sfree (sp)
+		call error (1,
+		    "Input image and profile image are not compatible")
+	    }
+	    call sprintf (Memc[str], SZ_LINE,
+	        "EXTRACT - Using profile image %s for %s")
+		call pargstr (profiles)
+		call pargstr (input)
+	    call ap_log (Memc[str], YES, YES, NO)
+	}
+
+	# Determine limits of apertures for use in defining memory requirements
+	# and I/O.
+
+	call salloc (a, 2 * napsex, TY_INT)
+	call salloc (b, 2 * napsex, TY_INT)
+	call salloc (c, 2 * napsex, TY_REAL)
+	a1 = a - 1
+	a2 = a1 + napsex
+	b1 = b - 1
+	b2 = b1 + napsex
+	c1 = c - 1
+	c2 = c1 + napsex
+
+	# Initialize image interpolator for edge pixel weighting.
+	switch (interptype) {
+	case II_LINEAR:
+	    interpbuf = 2
+	case II_POLY3:
+	    interpbuf = 3
+	case II_SINC:
+	    interpbuf = 16
+	default:
+	    interpbuf = 0
+	}
+	if (interptype > 0)
+	    call asiinit (asi, interptype)
+	else
+	    asi = NULL
+
+	na1 = 0
+	do i = 1, napsex {
+	    ap = Memi[apsex+i-1]
+	    cv = AP_CV(ap)
+	    ic = AP_IC(ap)
+
+	    # Dispersion axis limits
+	    bmin = min (nb, max (1, nint (AP_CEN(ap,baxis)+AP_LOW(ap,baxis))))
+	    bmax = max (1, min (nb, nint (AP_CEN(ap,baxis)+AP_HIGH(ap,baxis))))
+
+	    # Aperture axis shifts
+	    if (cv != NULL) {
+		cmin = MAX_REAL
+	        cmax = -MAX_REAL
+	        do j = bmin, bmax {
+	           shift = ap_cveval (cv, real (j))
+	           cmin = min (cmin, shift)
+	           cmax = max (cmax, shift)
+	        }
+	    } else {
+		cmin = 0.
+		cmax = 0.
+	    }
+
+	    # Background region limits.
+	    xmin = AP_LOW(ap,aaxis)
+	    xmax = AP_HIGH(ap,aaxis)
+	    if (weights == W_VARIANCE) {
+		xmin = xmin - 2
+		xmax = xmax + 2
+	    }
+	    xmin = xmin - interpbuf
+	    xmax = xmax + interpbuf
+	    if (bkg != B_NONE && AP_IC(ap) != NULL) {
+		xmin = min (xmin, ic_getr (ic, "xmin"))
+		xmax = max (xmax, ic_getr (ic, "xmax"))
+	    }
+
+	    Memi[a1+i] = min (namax, max (1, nint (AP_CEN(ap,aaxis)+xmin+cmin)))
+	    Memi[a2+i] = max (1, min (namax, nint (AP_CEN(ap,aaxis)+xmax+cmax)))
+	    Memi[b1+i] = bmin
+	    Memi[b2+i] = bmax
+	    Memr[c1+i] = cmin
+	    Memr[c2+i] = cmax
+	}
+	call alimi (Memi[a], 2*napsex, amin, amax)
+	call alimi (Memi[b], 2*napsex, bmin, bmax)
+
+	# The maximum size of the image in memory is 80% of the maximum
+	# working set size returned by begmem or 40% if a profile image
+	# is used.  Later I/O may exceed this since at least one
+	# aperture + background is needed in memory.
+
+	new_size = begmem (0, old_size, max_size)
+	namax = (amax - amin + 1)
+	nb = (bmax - bmin + 1)
+	if (pim == NULL)
+	    namax = min (namax, int (0.8 * max_size / SZ_REAL / nb))
+	else
+	    namax = min (namax, int (0.8 * max_size / SZ_REAL / nb / 2))
+	best_size = 1.2 * namax * nb * SZ_REAL
+	new_size = begmem (best_size, old_size, max_size)
+
+	# Allocate auxilary memory.  Some memory is only dependent on the
+	# number of dispersion points and subapertures and is the same for
+	# all apertures.  Other memory, such as the sky and profile depend on
+	# the aperture widths and tilts which may vary.  The input data is
+	# expected to have the aperture axis along the first dimension.  If
+	# the image is in this orientation then the IMIO buffer is used.
+	# Otherwise sequential I/O is used and transposed into the allocated
+	# memory.
+
+	iferr {
+	    call salloc (astart, nb, TY_INT)
+	    call salloc (spec, nsubaps * nb, TY_REAL)
+	    if (weights == W_VARIANCE) {
+		call salloc (raw, nsubaps * nb, TY_REAL)
+		call salloc (specsig, nsubaps * nb, TY_REAL)
+	    } else {
+		raw = NULL
+		specsig = NULL
+	    }
+	    profile = NULL
+	    if (aaxis == 2) {
+		call calloc (dbuf, namax * nb, TY_REAL)
+		if (pim != NULL)
+		    call calloc (pbuf, namax * nb, TY_REAL)
+	    }
+
+	    # For variance weighting the computations are done in photon units.
+	    if (weights == W_VARIANCE) {
+		gain = apgimr ("gain", im)
+		rdnoise = apgimr ("readnoise", im)
+	    } else {
+		gain = 1
+		rdnoise = 0
+	    }
+
+	    # Loop through each aperture doing the extractions.
+	    amax = 0
+	    do i = 1, napsex {
+		ap = Memi[apsex+i-1]
+
+		# Check if a new input data buffer is needed.  As many apertures
+		# as possible are read at once within the given memory limits
+		# though at least one aperture must be read.  Do a transpose if
+		# needed.
+
+		if (Memi[a1+i] < amin || Memi[a2+i] > amax) {
+		    amin = Memi[a1+i]
+		    amax = Memi[a2+i]
+		    do j = i,  napsex {
+			amin = min (amin, Memi[a1+j]) 
+			amax = max (amax, Memi[a2+j]) 
+			na = amax - amin + 1
+			if (na > namax)
+			    break
+		    }
+
+		    if (aaxis == 1) {
+			if (fmt == DIFF) {
+			    call mfree (dbuf, TY_REAL)
+			    call malloc (dbuf, na*nb, TY_REAL)
+			    call amovr (Memr[imgs2r(im,amin,amax,bmin,bmax)],
+				Memr[dbuf], na*nb)
+			} else
+			    dbuf = imgs2r (im, amin, amax, bmin, bmax)
+		    } else {
+			if (na > namax) {
+			    call mfree (dbuf, TY_REAL)
+			    namax = na
+			    call calloc (dbuf, namax * nb, TY_REAL)
+			}
+			do j = amin, amax {
+			    sbuf = imgl2r (im, j)
+			    sbuf = sbuf + bmin - 1
+			    ptr = dbuf + j - amin
+			    do k = bmin, bmax {
+				Memr[ptr] = Memr[sbuf]
+				sbuf = sbuf + 1
+				ptr = ptr + na
+			    }
+			}
+		    }
+		    if (pim != NULL) {
+			if (aaxis == 1)
+			    pbuf = imgs2r (pim, amin, amax, bmin, bmax)
+			else {
+			    if (na > namax) {
+				call mfree (pbuf, TY_REAL)
+				namax = na
+				call calloc (pbuf, namax * nb, TY_REAL)
+			    }
+			    do j = amin, amax {
+				sbuf = imgl2r (pim, j)
+				sbuf = sbuf + bmin - 1
+				ptr = pbuf + j - amin
+				do k = bmin, bmax {
+				    Memr[ptr] = Memr[sbuf]
+				    sbuf = sbuf + 1
+				    ptr = ptr + na
+				}
+			    }
+			}
+		    }
+		    if (weights == W_VARIANCE && gain != 1.) {
+			j = na * nb
+			call amulkr (Memr[dbuf], gain, Memr[dbuf], j)
+			if (pim != NULL)
+			    call amulkr (Memr[pbuf], gain, Memr[pbuf], j)
+		    }
+		}
+
+		# To minimize memory a variable integer offset is used to
+		# accomodate the aperture tilts.  The offsets are stored in
+		# the astart array and the width of any one line determined.
+		# If a stored profile is used it is read and it is ASSUMED to
+		# be valid for the input aperture with the same ID.  If no
+		# stored profile is found the profile fitting algorithm
+		# parameter determines whether to fit 1D function along the
+		# image axes (in which case all the profile offsets are the
+		# same) or if the Marsh algorithm for tilted spectra is
+		# used.  In the latter the offsets can be adjusted to mimize
+		# memory and a buffer of two pixels around the aperture is
+		# required by the algorithm.
+
+		if (weights == W_NONE) {
+		    xmin = AP_CEN(ap,aaxis) + AP_LOW(ap,aaxis)
+		    xmax = AP_CEN(ap,aaxis) + AP_HIGH(ap,aaxis)
+		    xmin = xmin - interpbuf
+		    xmax = xmax + interpbuf
+		    na1 = nint (xmax) - nint (xmin) + 1
+		    cv = AP_CV(ap)
+		    do j = bmin, bmax {
+			shift = ap_cveval (cv, real (j))
+			Memi[astart+j-bmin] = nint (xmin + shift)
+		    }
+		} else {
+		    if (pfit == P_FIT1D) {
+			xmin = AP_CEN(ap,aaxis) + AP_LOW(ap,aaxis) + Memr[c1+i]
+			xmax = AP_CEN(ap,aaxis) + AP_HIGH(ap,aaxis) + Memr[c2+i]
+			xmin = xmin - interpbuf
+			xmax = xmax + interpbuf
+			na1 = nint (xmax) - nint (xmin) + 1
+			call amovki (nint (xmin), Memi[astart], nb)
+		    } else if (pfit == P_FIT2D) {
+			xmin = AP_CEN(ap,aaxis) + AP_LOW(ap,aaxis) - 2
+			xmax = AP_CEN(ap,aaxis) + AP_HIGH(ap,aaxis) + 2
+			xmin = xmin - interpbuf
+			xmax = xmax + interpbuf
+			na1 = nint (xmax) - nint (xmin) + 1
+			cv = AP_CV(ap)
+			do j = bmin, bmax {
+			    shift = ap_cveval (cv, real (j))
+			    Memi[astart+j-bmin] = nint (xmin + shift)
+			}
+		    }
+		}
+
+		# Do the sky or background determination if needed.  An array
+		# of the same size as the 2D aperture is returned as well as
+		# a single estimate of the variance in the sky value at each
+		# line based on the fit.  If a profile image is used then the
+		# sky is for the profile image and the object sky is
+		# determined later in order to reuse the sky buffers.
+
+		if (bkg != B_NONE && AP_IC(ap) != NULL) {
+		    call malloc (sbuf, na1 * nb, TY_REAL)
+		    call malloc (svar, nb, TY_REAL)
+		    call malloc (specsky, nsubaps * nb, TY_REAL)
+		    if (pim == NULL)
+			call ap_skyeval (im, ap, dbuf, na, nb, amin, 1,
+			    Memr[sbuf], Memr[svar], Memr[specsky], na1, nb,
+			    Memi[astart], 1, nsubaps, rdnoise)
+		    else
+			call ap_skyeval (pim, ap, pbuf, na, nb, amin, 1,
+			    Memr[sbuf], Memr[svar], Memr[specsky], na1, nb,
+			    Memi[astart], 1, nsubaps, rdnoise)
+		} else {
+		    sbuf = NULL
+		    svar = NULL
+		    specsky = NULL 
+		}
+
+		# Use a quick sum for unweighted extraction.  For weighed
+		# extractions we use either a previously determined profile
+		# or call the profile routine.  If desired the profile is
+		# stored for later use.  Then the variance weighted
+		# extraction routine is called.
+
+		if (weights == W_NONE)
+		    call ap_sum (ap, dbuf, na, nb, amin, 1, sbuf, na1, nb,
+			Memi[astart], 1, Memr[spec], nsubaps, asi)
+		else {
+		    call malloc (profile, na1 * nb, TY_REAL)
+		    if (pim == NULL)
+			call ap_profile (im, ap, dbuf, na, nb, amin, 1, sbuf,
+			    svar, Memr[profile], na1, nb, Memi[astart], 1,
+			    asi)
+		    else {
+			call ap_profile (pim, ap, pbuf, na, nb, amin, 1, sbuf,
+			    svar, Memr[profile], na1, nb, Memi[astart], 1,
+			    asi)
+			if (sbuf != NULL)
+			    call ap_skyeval (im, ap, dbuf, na, nb, amin, 1,
+				Memr[sbuf], Memr[svar], Memr[specsky], na1, nb,
+				Memi[astart], 1, nsubaps, rdnoise)
+		    }
+
+		    call ap_variance (im, ap, dbuf, na, nb, amin, 1, sbuf, svar,
+			Memr[profile], na1, nb, Memi[astart], 1, Memr[spec],
+			Memr[raw], Memr[specsig], nsubaps, asi)
+		}
+
+		# Output the extracted spectrum.  The extras of sky, sigma,
+		# and unweighted spectrum may also be stored.  If the extra
+		# information is not available the pointers will be NULL.
+
+		if (weights == W_VARIANCE && gain != 1.) {
+		    call adivkr (Memr[spec], gain, Memr[spec], nb)
+		    if (raw != NULL)
+			call adivkr (Memr[raw], gain, Memr[raw], nb)
+		    if (specsky != NULL)
+			call adivkr (Memr[specsky], gain, Memr[specsky], nb)
+		    if (specsig != NULL)
+			call adivkr (Memr[specsig], gain, Memr[specsig], nb)
+		    call amulkr (Memr[profile], gain, Memr[profile], nb*na1)
+		}
+
+		call ap_output (input, output, format, Memc[bkgstr],
+		    Memc[wtstr], Memc[cleanstr], gain, im, Memi[apsex], napsex,
+		    i, nsubaps, spec, raw, specsky, specsig, dbuf, na, nb, amin,
+		    1, sbuf, profile, na1, nb, Memi[astart], 1)
+
+		call mfree (profile, TY_REAL)
+		call mfree (sbuf, TY_REAL)
+		call mfree (svar, TY_REAL)
+		call mfree (specsky, TY_REAL)
+	    }
+
+	    # Finish up and restore the working set size.
+	    if (asi != NULL)
+		call asifree (asi)
+	    if (pim != NULL) {
+		if (aaxis == 2)
+		    call mfree (pbuf, TY_REAL)
+		call imunmap (pim)
+	    }
+	    if (aaxis == 2)
+		call mfree (dbuf, TY_REAL)
+	    call imunmap (im)
+	    call fixmem (old_size)
+	    call sfree (sp)
+
+	} then {
+	    call mfree (profile, TY_REAL)
+	    call mfree (sbuf, TY_REAL)
+	    call mfree (svar, TY_REAL)
+	    call mfree (specsky, TY_REAL)
+
+	    if (asi != NULL)
+		call asifree (asi)
+	    if (pim != NULL) {
+		if (aaxis == 2)
+		    call mfree (pbuf, TY_REAL)
+		call imunmap (pim)
+	    }
+	    if (aaxis == 2)
+		call mfree (dbuf, TY_REAL)
+	    call imunmap (im)
+	    call fixmem (old_size)
+	    call sfree (sp)
+
+	    call erract (EA_ERROR)
+	}
+end
+
+
+# AP_CHECK -- Check if output spectra exist.  If the user allows clobbering,
+# delete the spectra.  Return the format.
+
+int procedure ap_check (input, output, format, aps, naps, nsubaps)
+
+char	input[ARB]			# Input image name
+char	output[ARB]			# Output root name
+char	format[ARB]			# Output format
+pointer	aps[naps]			# Apertures
+int	naps				# Number of apertures
+int	nsubaps				# Number of subapertures
+
+int	i, j, fmt
+pointer	sp, name, name1, input1, ksection, ans
+
+int	strdic(), imaccess(), stridxs()
+bool	streq(), ap_answer()
+
+begin
+	call smark (sp)
+	call salloc (name, SZ_LINE, TY_CHAR)
+	call salloc (name1, SZ_LINE, TY_CHAR)
+	call salloc (input1, SZ_LINE, TY_CHAR)
+	call salloc (ksection, SZ_LINE, TY_CHAR)
+	call salloc (ans, SZ_LINE, TY_CHAR)
+
+	fmt = strdic (format, format, SZ_LINE, FORMATS)
+	call imgimage (input, Memc[input1], SZ_LINE)
+
+	switch (fmt) {
+	case MULTISPEC, NORM, FLAT, RATIO, DIFF, FIT:
+	    i = stridxs ("[", Memc[input1])
+	    if (i > 0) {
+		call strcpy (Memc[input1+i-1], Memc[ksection], SZ_LINE)
+		Memc[input1+i-1] = EOS
+	    } else
+		Memc[ksection] = EOS
+	    if (output[1] == EOS)
+	        call strcpy (Memc[input1], Memc[name], SZ_LINE)
+	    else
+	        call strcpy (output, Memc[name], SZ_LINE)
+
+	    switch (fmt) {
+	    case MULTISPEC:
+		if (streq (Memc[input1], Memc[name])) {
+	            call strcat (".ms", Memc[name], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name], SZ_LINE)
+		}
+	    case NORM:
+		if (streq (Memc[input1], Memc[name])) {
+	            call strcat (".norm", Memc[name], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name], SZ_LINE)
+		}
+	    case FLAT:
+		if (streq (Memc[input1], Memc[name])) {
+	            call strcat (".flat", Memc[name], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name], SZ_LINE)
+		}
+	    case RATIO:
+		if (streq (Memc[input1], Memc[name])) {
+	            call strcat (".ratio", Memc[name], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name], SZ_LINE)
+		}
+	    case DIFF:
+		if (streq (Memc[input1], Memc[name])) {
+	            call strcat (".diff", Memc[name], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name], SZ_LINE)
+		}
+	    case FIT:
+		if (streq (Memc[input1], Memc[name])) {
+	            call strcat (".fit", Memc[name], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name], SZ_LINE)
+		}
+	    }
+	    if (imaccess (Memc[name], 0) == YES) {
+		call sprintf (Memc[ans], SZ_LINE,
+		    "Clobber existing output image %s?")
+		    call pargstr (Memc[name])
+		if (ap_answer ("ansclobber1", Memc[ans]))
+		    call imdelete (Memc[name])
+		else {
+		    call sprintf (Memc[ans], SZ_LINE,
+			"EXTRACT - Output spectrum %s already exists")
+			call pargstr (Memc[name])
+		    call ap_log (Memc[ans], YES, NO, YES)
+		    fmt = 0
+		}
+	    }
+	case ECHELLE:
+	    if (output[1] == EOS)
+	        call strcpy (Memc[input1], Memc[name], SZ_LINE)
+	    else
+	        call strcpy (output, Memc[name], SZ_LINE)
+
+	    do i = 1, nsubaps {
+		if (nsubaps == 1)
+		    call strcpy (Memc[name], Memc[name1], SZ_LINE)
+		else {
+		    call sprintf (Memc[name1], SZ_LINE, "%s%0*d")
+			call pargstr (Memc[name])
+			call pargi (int(log10(real(nsubaps)))+1)
+			call pargi (i)
+		}
+		if (streq (Memc[input1], Memc[name])) {
+		    call strcat (".ec", Memc[name1], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name1], SZ_LINE)
+		}
+
+		if (imaccess (Memc[name1], 0) == YES) {
+		    call sprintf (Memc[ans], SZ_LINE,
+			"Clobber existing output image %s?")
+			call pargstr (Memc[name1])
+		    if (ap_answer ("ansclobber1", Memc[ans]))
+			call imdelete (Memc[name1])
+		    else {
+			call sprintf (Memc[ans], SZ_LINE,
+			    "EXTRACT - Output spectrum %s already exists")
+			    call pargstr (Memc[name1])
+			call ap_log (Memc[ans], YES, NO, YES)
+			fmt = 0
+		    }
+		}
+	    }
+	case ONEDSPEC, STRIP:
+	    do i = 1, naps {
+		do j = 1, nsubaps {
+		    call sprintf (Memc[name], SZ_LINE, "%s.%0*d")
+			if (output[1] == EOS)
+			    call pargstr (Memc[input1])
+			else
+			    call pargstr (output)
+			call pargi (int(log10(real(nsubaps)))+4)
+			call pargi (AP_ID(aps[i])+(j-1)*1000)
+		    if (imaccess (Memc[name], 0) == YES) {
+			call sprintf (Memc[ans], SZ_LINE,
+			    "Clobber existing output image %s?")
+			    call pargstr (Memc[name])
+			if (ap_answer ("ansclobber1", Memc[ans]))
+			    call imdelete (Memc[name])
+			else {
+			    call sprintf (Memc[ans], SZ_LINE,
+				"EXTRACT - Output spectrum %s already exists")
+				call pargstr (Memc[name])
+			    call ap_log (Memc[ans], YES, NO, YES)
+			    fmt = 0
+			}
+		    }
+		}
+	    }
+	case NOISE:
+	    ;
+	default:
+	    call sfree (sp)
+	    call error (1, "EXTRACT - Unknown output format")
+	}
+
+	call sfree (sp)
+	return (fmt)
+end
+
+
+# AP_OUTPUT -- Review the extracted spectra and write them to an image.
+# This routine determines the output format and whether to also output sky
+# unweighted, and sigma spectra.  The appropriate header keywords have
+# to be added.
+
+procedure ap_output (image, output, format, bkg, wt, clean, gain, in, aps,
+	naps, iap, nsubaps, spec, raw sky, sig, dbuf, nc, nl, c1, l1, sbuf,
+	profile, nx, ny, xs, ys)
+
+char	image[ARB]		# Input image name
+char	output[ARB]		# Output root name
+char	format[ARB]		# Output format
+char	bkg[ARB]		# Background type
+char	wt[ARB]			# Weight type
+char	clean[ARB]		# Clean?
+real	gain			# Gain
+pointer	in			# Input IMIO pointer
+pointer	aps[naps]		# Apertures
+int	naps			# Number of apertures
+int	iap			# Aperture
+int	nsubaps			# Number of subapertures
+pointer	spec			# Output spectrum
+pointer	raw			# Output raw spectrum
+pointer	sky			# Output sky
+pointer	sig			# Output sigma
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+pointer	sbuf			# Sky values (NULL if none)
+pointer	profile			# Profile (NULL if none)
+int	nx, ny			# Size of sky and profile array
+int	xs[ny], ys		# Origin of sky and profile array
+
+int	fmt			# Output format
+bool	extras			# Include raw spectrum, sky, and sigma
+
+real	low, high, step
+int	i, k, l, m, apid, apaxis, dispaxis
+pointer	sp, str, str1, name, name1, input, ksection
+pointer	ap, out, outsave, gt, apmw, buf
+pointer	sum2, sum4, nsum
+
+real	clgetr()
+int	scan(), strdic(), imaccf(), stridxs()
+bool	streq(), ap_answer(), apgetb()
+pointer	immap(), imgl2r(), impl2r(), impl3r()
+pointer	gt_init(), apmw_open()
+errchk	immap, impl2r, impl3r, imps2r, ap_strip, ap_pstrip, apmw_open
+errchk	ap_fitspec, ap_lnorm, ap_cnorm, ap_lflat, ap_cflat
+
+begin
+	# Allocate string and file name arrays.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (name, SZ_LINE, TY_CHAR)
+	call salloc (name1, SZ_LINE, TY_CHAR)
+	call salloc (input, SZ_LINE, TY_CHAR)
+	call salloc (ksection, SZ_LINE, TY_CHAR)
+
+	fmt = strdic (format, format, SZ_LINE, FORMATS)
+	extras = apgetb ("extras")
+
+	ap = aps[iap]
+	apaxis = AP_AXIS(ap)
+	dispaxis = mod (apaxis, 2) + 1
+
+	# Set output name.
+	call imgimage (image, Memc[input], SZ_LINE)
+	i = stridxs ("[", Memc[input])
+	if (i > 0) {
+	    call strcpy (Memc[input+i-1], Memc[ksection], SZ_LINE)
+	    Memc[input+i-1] = EOS
+	    i = stridxs ("]", Memc[ksection])
+	    call strcpy (",append]", Memc[ksection+i-1], SZ_LINE)
+	} else
+	    Memc[ksection] = EOS
+	if (output[1] == EOS)
+	    call strcpy (Memc[input], Memc[name], SZ_LINE)
+	else
+	    call strcpy (output, Memc[name], SZ_LINE)
+
+	switch (fmt) {
+	case ECHELLE:
+	    ;
+	case MULTISPEC:
+	    if (streq (Memc[input], Memc[name])) {
+		call strcat (".ms", Memc[name], SZ_LINE)
+		call strcat (Memc[ksection], Memc[name], SZ_LINE)
+	    }
+	case NORM:
+	    if (streq (Memc[input], Memc[name])) {
+	        call strcat (".norm", Memc[name], SZ_LINE)
+		call strcat (Memc[ksection], Memc[name], SZ_LINE)
+	    }
+	case FLAT:
+	    if (streq (Memc[input], Memc[name])) {
+	        call strcat (".flat", Memc[name], SZ_LINE)
+		call strcat (Memc[ksection], Memc[name], SZ_LINE)
+	    }
+	case RATIO:
+	    if (streq (Memc[input], Memc[name])) {
+	        call strcat (".ratio", Memc[name], SZ_LINE)
+		call strcat (Memc[ksection], Memc[name], SZ_LINE)
+	    }
+	case DIFF:
+	    if (streq (Memc[input], Memc[name])) {
+	        call strcat (".diff", Memc[name], SZ_LINE)
+		call strcat (Memc[ksection], Memc[name], SZ_LINE)
+	    }
+	case FIT:
+	    if (streq (Memc[input], Memc[name])) {
+	        call strcat (".fit", Memc[name], SZ_LINE)
+		call strcat (Memc[ksection], Memc[name], SZ_LINE)
+	    }
+	case NOISE:
+	    Memc[name] = EOS
+	}
+
+
+	# Set the review graph title.
+	call sprintf (Memc[str], SZ_LINE, "%s: %s - Aperture %s")
+	    call pargstr (image)
+	    call pargstr (IM_TITLE(in))
+	    call pargi (AP_ID(ap))
+
+	gt = gt_init ()
+	call gt_sets (gt, GTTITLE, Memc[str])
+
+	# Query the user whether to review the extraction.
+	call sprintf (Memc[str], SZ_LINE,
+	    "Review extracted spectrum for aperture %d from %s?")
+	    call pargi (AP_ID(ap))
+	    call pargstr (image)
+
+	# If reviewing graph the spectrum, do a cursor loop, and allow
+	# the user to skip the output or define a new output image.
+	if (ap_answer ("ansreview1", Memc[str])) {
+	    call ap_graph1 (gt, Memr[spec], ny, nsubaps)
+	    
+	    if (fmt == ONEDSPEC && nsubaps == 1) {
+	        call printf (
+		    "Output image name [use # to skip output] (%s): ")
+		    call pargstr (Memc[name])
+		call flush (STDOUT)
+	        if (scan() != EOF) {
+	            call gargwrd (Memc[str], SZ_LINE)
+	            if (Memc[str] == '#') {
+			call gt_free (gt)
+			call sfree (sp)
+		        return
+		    }
+	            if (Memc[str] != EOS)
+		        call strcpy (Memc[str], Memc[name], SZ_LINE)
+		}
+	    }
+	}
+
+	# Output the image.
+	switch (fmt) {
+	case MULTISPEC:
+	    if (iap == 1) {
+		out = immap (Memc[name], NEW_COPY, in)
+
+		IM_PIXTYPE(out) = TY_REAL
+		IM_NDIM(out) = 1
+		IM_LEN(out, 1) = ny
+		IM_LEN(out, 2) = nsubaps * naps
+		IM_LEN(out, 3) = 1
+		if (extras) {
+		    if (sky != NULL)
+		        IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		    if (raw != NULL)
+		        IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		    if (sig != NULL)
+		        IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		}
+		if (IM_LEN(out, 2) > 1)
+		    IM_NDIM(out) = 2
+		if (IM_LEN(out, 3) > 1)
+		    IM_NDIM(out) = 3
+
+		apmw = apmw_open (in, out, dispaxis, nsubaps*naps, ny)
+
+		# Write BAND IDs.
+		k = 1
+		call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+		    call pargi (k)
+		call sprintf (Memc[str], SZ_LINE,
+		    "spectrum - background %s, weights %s, clean %s")
+		    call pargstr (bkg)
+		    call pargstr (wt)
+		    call pargstr (clean)
+		call imastr (out, Memc[str1], Memc[str])
+		k = k + 1
+		if (extras) {
+		    if (raw != NULL) {
+			call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			    call pargi (k)
+			call sprintf (Memc[str], SZ_LINE,
+			    "raw - background %s, weights none, clean no")
+			    call pargstr (bkg)
+			call imastr (out, Memc[str1], Memc[str])
+			k = k + 1
+		    }
+		    if (sky != NULL) {
+			call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			    call pargi (k)
+			call sprintf (Memc[str], SZ_LINE,
+			    "background - background %s")
+			    call pargstr (bkg)
+			call imastr (out, Memc[str1], Memc[str])
+			k = k + 1
+		    }
+		    if (sig != NULL) {
+			call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			    call pargi (k)
+			call sprintf (Memc[str], SZ_LINE,
+			    "sigma - background %s, weights %s, clean %s")
+			    call pargstr (bkg)
+			    call pargstr (wt)
+			    call pargstr (clean)
+			call imastr (out, Memc[str1], Memc[str])
+		    }
+		}
+
+		do k = 1, naps {
+		    low = AP_CEN(aps[k],apaxis) + AP_LOW(aps[k],apaxis)
+		    high = AP_CEN(aps[k],apaxis) + AP_HIGH(aps[k],apaxis)
+		    step = (high - low) / nsubaps
+		    low = low - step
+		    do l = 1, nsubaps {
+			apid = AP_ID(aps[k]) + (l - 1) * 1000
+			low = low + step
+			high = low + step
+			call apmw_setap (apmw, (k-1)*nsubaps+l,
+			    apid, AP_BEAM(aps[k]), low, high)
+		    }
+		}
+		do k = 1, naps {
+		    if (AP_TITLE(aps[k]) != NULL) {
+			do l = 1, nsubaps {
+			    call sprintf (Memc[str], SZ_LINE, "APID%d")
+				call pargi ((k-1)*nsubaps+l)
+		    	    call imastr (out, Memc[str],
+				Memc[AP_TITLE(aps[k])])
+			}
+		    }
+		}
+	    }
+	
+	    do l = 1, nsubaps {
+		k = (iap - 1) * nsubaps + l
+		buf = impl2r (out, k)
+	        call amovr (Memr[spec+(l-1)*ny], Memr[buf], ny)
+		if (extras) {
+		    m = 2
+		    if (raw != NULL) {
+			buf = impl3r (out, k, m)
+	                call amovr (Memr[raw+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		    if (sky != NULL) {
+			buf = impl3r (out, k, m)
+	                call amovr (Memr[sky+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		    if (sig != NULL) {
+			buf = impl3r (out, k, m)
+	                call amovr (Memr[sig+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		}
+	    }
+	    if (iap == naps) {
+		call apmw_saveim (apmw, out, fmt)
+		call apmw_close (apmw)
+	        call imunmap (out)
+	    }
+
+	    if (Memc[name] != EOS) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s")
+		    call pargi (AP_ID(ap))
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		call ap_log (Memc[str], YES, YES, NO)
+		call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+	    }
+
+	case ECHELLE:
+	    do l = 1, nsubaps {
+		if (nsubaps == 1)
+		    call strcpy (Memc[name], Memc[name1], SZ_LINE)
+		else {
+		    call sprintf (Memc[name1], SZ_LINE, "%s%0*d")
+			call pargstr (Memc[name])
+			call pargi (int(log10(real(nsubaps)))+1)
+			call pargi (l)
+		}
+		if (streq (Memc[input], Memc[name])) {
+		    call strcat (".ec", Memc[name1], SZ_LINE)
+		    call strcat (Memc[ksection], Memc[name1], SZ_LINE)
+		}
+
+		if (iap == 1) {
+		    out = immap (Memc[name1], NEW_COPY, in)
+
+		    IM_PIXTYPE(out) = TY_REAL
+		    IM_NDIM(out) = 1
+		    IM_LEN(out, 1) = ny
+		    IM_LEN(out, 2) = naps
+		    IM_LEN(out, 3) = 1
+		    if (extras) {
+			if (sky != NULL)
+			    IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+			if (raw != NULL)
+			    IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+			if (sig != NULL)
+			    IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		    }
+		    if (IM_LEN(out, 2) > 1)
+			IM_NDIM(out) = 2
+		    if (IM_LEN(out, 3) > 1)
+			IM_NDIM(out) = 3
+
+		    apmw = apmw_open (in, out, dispaxis, naps, ny)
+
+		    # Write BAND IDs.
+		    k = 1
+		    call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			call pargi (k)
+		    call sprintf (Memc[str], SZ_LINE,
+			"spectrum - background %s, weights %s, clean %s")
+			call pargstr (bkg)
+			call pargstr (wt)
+			call pargstr (clean)
+		    call imastr (out, Memc[str1], Memc[str])
+		    k = k + 1
+		    if (extras) {
+			if (raw != NULL) {
+			    call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+				call pargi (k)
+			    call sprintf (Memc[str], SZ_LINE,
+				"raw - background %s, weights none, clean no")
+				call pargstr (bkg)
+			    call imastr (out, Memc[str1], Memc[str])
+			    k = k + 1
+			}
+			if (sky != NULL) {
+			    call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+				call pargi (k)
+			    call sprintf (Memc[str], SZ_LINE,
+				"background - background %s")
+				call pargstr (bkg)
+			    call imastr (out, Memc[str1], Memc[str])
+			    k = k + 1
+			}
+			if (sig != NULL) {
+			    call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+				call pargi (k)
+			    call sprintf (Memc[str], SZ_LINE,
+				"sigma - background %s, weights %s, clean %s")
+				call pargstr (bkg)
+				call pargstr (wt)
+				call pargstr (clean)
+			    call imastr (out, Memc[str1], Memc[str])
+			}
+		    }
+
+		    # Write keyword to allow matching by subaperture.
+		    if (nsubaps > 1)
+			call imaddi (out, "SUBAP", l)
+
+		    do k = 1, naps {
+			low = AP_CEN(aps[k],apaxis) + AP_LOW(aps[k],apaxis)
+			high = AP_CEN(aps[k],apaxis) + AP_HIGH(aps[k],apaxis)
+			step = (high - low) / nsubaps
+			call apmw_setap (apmw, k, AP_ID(aps[k]),
+			    AP_BEAM(aps[k]), low+(l-1)*step, low+l*step)
+		    }
+		    do k = 1, naps {
+			if (AP_TITLE(aps[k]) != NULL) {
+			    call sprintf (Memc[str], SZ_LINE, "APID%d")
+				call pargi (k)
+			    call imastr (out, Memc[str],
+				Memc[AP_TITLE(aps[k])])
+			}
+		    }
+		} else {
+		    if (l == 1)
+			out = outsave
+		    else
+			out = immap (Memc[name1], READ_WRITE, 0)
+		}
+	    
+		k = iap
+		buf = impl2r (out, k)
+		call amovr (Memr[spec+(l-1)*ny], Memr[buf], ny)
+		if (extras) {
+		    m = 2
+		    if (raw != NULL) {
+			buf = impl3r (out, k, m)
+			call amovr (Memr[raw+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		    if (sky != NULL) {
+			buf = impl3r (out, k, m)
+			call amovr (Memr[sky+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		    if (sig != NULL) {
+			buf = impl3r (out, k, m)
+			call amovr (Memr[sig+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		}
+
+		if (iap == 1) {
+		    call apmw_saveim (apmw, out, fmt)
+		    call apmw_close (apmw)
+		}
+		if (l != 1 || iap == naps)
+		    call imunmap (out)
+		if (l == 1)
+		    outsave = out
+
+		if (nsubaps == 1) {
+		    call sprintf (Memc[str], SZ_LINE,
+			"EXTRACT - Aperture %d from %s --> %s")
+			call pargi (AP_ID(ap))
+			call pargstr (image)
+			call pargstr (Memc[name1])
+		} else {
+		    call sprintf (Memc[str], SZ_LINE,
+			"EXTRACT - Aperture %d-%d from %s --> %s")
+			call pargi (AP_ID(ap))
+			call pargi (l)
+			call pargstr (image)
+			call pargstr (Memc[name1])
+		}
+		call ap_log (Memc[str], YES, YES, NO)
+	    }
+
+	    call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+			
+	case ONEDSPEC:
+	    do l = 1, nsubaps {
+		apid = AP_ID(ap) + (l - 1) * 1000
+		low = AP_CEN(ap,apaxis) + AP_LOW(ap,apaxis)
+		high = AP_CEN(ap,apaxis) + AP_HIGH(ap,apaxis)
+		step = (high - low) / nsubaps
+		low = low + (l - 1) * step
+		high = low + step
+
+		call sprintf (Memc[str], SZ_LINE, "%s.%0*d")
+		    call pargstr (Memc[name])
+		    call pargi (int(log10(real(nsubaps)))+4)
+		    call pargi (apid)
+		out = immap (Memc[str], NEW_COPY, in)
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s.%0*d")
+		    call pargi (apid)
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		    call pargi (int(log10(real(nsubaps)))+4)
+		    call pargi (apid)
+		call ap_log (Memc[str], YES, YES, NO)
+
+		apmw = apmw_open (in, out, dispaxis, 1, ny)
+		call apmw_setap (apmw, 1, apid, AP_BEAM(ap), low, high)
+		if (AP_TITLE(ap) != NULL)
+		    call imastr (out, "APID1", Memc[AP_TITLE(ap)])
+
+		IM_PIXTYPE(out) = TY_REAL
+		IM_NDIM(out) = 1
+		IM_LEN(out, 1) = ny
+		IM_LEN(out, 2) = 1
+		IM_LEN(out, 3) = 1
+		if (extras) {
+		    if (sky != NULL)
+			IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		    if (raw != NULL)
+			IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		    if (sig != NULL)
+			IM_LEN(out, 3) = IM_LEN(out, 3) + 1
+		}
+		if (IM_LEN(out, 2) > 1)
+		    IM_NDIM(out) = 2
+		if (IM_LEN(out, 3) > 1)
+		    IM_NDIM(out) = 3
+
+		# Write BAND IDs.
+		k = 1
+		call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+		    call pargi (k)
+		call sprintf (Memc[str], SZ_LINE,
+		    "spectrum: background %s, weights %s, clean %s")
+		    call pargstr (bkg)
+		    call pargstr (wt)
+		    call pargstr (clean)
+		call imastr (out, Memc[str1], Memc[str])
+		k = k + 1
+		if (extras) {
+		    if (raw != NULL) {
+			call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			    call pargi (k)
+			call sprintf (Memc[str], SZ_LINE,
+			    "spectrum: background %s, weights none, clean no")
+			    call pargstr (bkg)
+			call imastr (out, Memc[str1], Memc[str])
+			k = k + 1
+		    }
+		    if (sky != NULL) {
+			call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			    call pargi (k)
+			call sprintf (Memc[str], SZ_LINE,
+			    "background: background %s")
+			    call pargstr (bkg)
+			call imastr (out, Memc[str1], Memc[str])
+			k = k + 1
+		    }
+		    if (sig != NULL) {
+			call sprintf (Memc[str1], SZ_LINE, "BANDID%d")
+			    call pargi (k)
+			call sprintf (Memc[str], SZ_LINE,
+			    "sigma - background %s, weights %s, clean %s")
+			    call pargstr (bkg)
+			    call pargstr (wt)
+			    call pargstr (clean)
+			call imastr (out, Memc[str1], Memc[str])
+		    }
+		}
+
+		buf = impl2r (out, 1)
+	        call amovr (Memr[spec+(l-1)*ny], Memr[buf], ny)
+		if (extras) {
+		    m = 2
+		    if (raw != NULL) {
+			buf = impl3r (out, 1, m)
+	                call amovr (Memr[raw+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		    if (sky != NULL) {
+			buf = impl3r (out, 1, m)
+	                call amovr (Memr[sky+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		    if (sig != NULL) {
+			buf = impl3r (out, 1, m)
+	                call amovr (Memr[sig+(l-1)*ny], Memr[buf], ny)
+			m = m + 1
+		    }
+		}
+
+		call apmw_saveim (apmw, out, fmt)
+		call apmw_close (apmw)
+		call imunmap (out)
+
+	    }
+
+	    call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+
+	case STRIP:
+	    do l = 1, nsubaps {
+		apid = AP_ID(ap) + (l - 1) * 1000
+		low = AP_CEN(ap,apaxis) + AP_LOW(ap,apaxis)
+		high = AP_CEN(ap,apaxis) + AP_HIGH(ap,apaxis)
+		step = (high - low) / nsubaps
+		low = low + (l - 1) * step
+		high = low + step
+
+		call sprintf (Memc[str], SZ_LINE, "%s.%0*d")
+		    call pargstr (Memc[name])
+		    call pargi (int(log10(real(nsubaps)))+4)
+		    call pargi (apid)
+		out = immap (Memc[str], NEW_COPY, in)
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s.%0*d")
+		    call pargi (apid)
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		    call pargi (int(log10(real(nsubaps)))+4)
+		    call pargi (apid)
+		call ap_log (Memc[str], YES, YES, NO)
+
+		apmw = apmw_open (in, out, dispaxis, 1, ny)
+		call apmw_setap (apmw, 1, apid, AP_BEAM(ap), low, high)
+		call sprintf (Memc[str], SZ_LINE, "%s - Aperture %d")
+		    call pargstr (IM_TITLE(out))
+		    call pargi (AP_ID(ap))
+		call strcpy (Memc[str], IM_TITLE(out), SZ_IMTITLE)
+		if (AP_TITLE(ap) != NULL)
+		    call imastr (out, "APID1", Memc[AP_TITLE(ap)])
+
+		IM_PIXTYPE(out) = TY_REAL
+		IM_NDIM(out) = 2
+		IM_LEN(out, 1) = ny
+		IM_LEN(out, 2) = high - low + 1
+
+		if (profile == NULL)
+		    call ap_strip (ap, low, high, out, dbuf, nc, nl, c1, l1,
+			sbuf, nx, ny, xs, ys)
+		else
+		    call ap_pstrip (ap, low, high, out, gain, Memr[spec],
+			Memr[profile], nx, ny, xs, ys)
+
+		call apmw_saveim (apmw, out, fmt)
+		call apmw_close (apmw)
+		call imunmap (out)
+	    }
+
+	    call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+
+	case NORM, FLAT:
+	    if (iap == 1) {
+	        out = immap (Memc[name], NEW_COPY, in)
+	        IM_PIXTYPE(out) = TY_REAL
+		if (imaccf (out, "CCDMEAN") == YES)
+		    call imdelf (out, "CCDMEAN")
+	        call ap_fitspec (ap, in, Memr[spec], ny)
+		k = YES
+	    } else {
+	        call ap_fitspec (ap, in, Memr[spec], ny)
+		k = NO
+	    }
+	    if (apaxis == 1) {
+		if (fmt == NORM)
+		    call ap_lnorm (ap, out, gain, dbuf, nc, nl, c1, l1,
+			Memr[spec], ny, ys, k)
+		else
+		    call ap_lflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec],
+		        sbuf, Memr[profile], nx, ny, xs, ys, k)
+	    } else {
+		if (fmt == NORM)
+		    call ap_cnorm (ap, out, gain, dbuf, nc, nl, c1, l1,
+			Memr[spec], ny, ys, k)
+		else
+		    call ap_cflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec],
+		        sbuf, Memr[profile], nx, ny, xs, ys, k)
+	    }
+	    if (iap == naps)
+	        call imunmap (out)
+
+	    if (Memc[name] != EOS) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s")
+		    call pargi (AP_ID(ap))
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		call ap_log (Memc[str], YES, YES, NO)
+		call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+	    }
+
+	case RATIO, FIT:
+	    if (iap == 1) {
+		out = immap (Memc[name], NEW_COPY, in)
+	        IM_PIXTYPE(out) = TY_REAL
+		k = YES
+	    } else
+		k = NO
+	    if (apaxis == 1) {
+		switch (fmt) {
+		case RATIO:
+		    call ap_lflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec],
+		        sbuf, Memr[profile], nx, ny, xs, ys, k)
+		case FIT:
+		    call ap_lfit (ap, out, gain, Memr[spec], Memr[profile],
+			nx, ny, xs, ys, k)
+		}
+	    } else {
+		switch (fmt) {
+		case RATIO:
+		    call ap_cflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec],
+		        sbuf, Memr[profile], nx, ny, xs, ys, k)
+		case FIT:
+		    call ap_cfit (ap, out, gain, Memr[spec], Memr[profile],
+			nx, ny, xs, ys, k)
+		}
+	    }
+	    if (iap == naps)
+	        call imunmap (out)
+
+	    if (Memc[name] != EOS) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s")
+		    call pargi (AP_ID(ap))
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		call ap_log (Memc[str], YES, YES, NO)
+		call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+	    }
+
+	case DIFF:
+	    if (iap == 1) {
+	        out = immap (Memc[name], NEW_COPY, in)
+	        IM_PIXTYPE(out) = TY_REAL
+		do k = 1, IM_LEN(in,2) {
+		    buf = impl2r (out, k)
+		    call amovr (Memr[imgl2r(in,k)], Memr[buf], IM_LEN(out,1))
+		}
+		k = NO
+	    } else
+		k = NO
+	    if (apaxis == 1)
+		call ap_ldiff (ap, out, gain, dbuf, nc, nl, c1, l1, Memr[spec],
+		    Memr[profile], nx, ny, xs, ys, k)
+	    else
+		call ap_cdiff (ap, out, gain, dbuf, nc, nl, c1, l1, Memr[spec],
+		    Memr[profile], nx, ny, xs, ys, k)
+	    if (iap == naps)
+	        call imunmap (out)
+
+	    if (Memc[name] != EOS) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s")
+		    call pargi (AP_ID(ap))
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		call ap_log (Memc[str], YES, YES, NO)
+		call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+	    }
+
+	case NOISE:
+	    if (iap == 1) {
+		low = clgetr ("dmin")
+		high = clgetr ("dmax")
+		l = clgetr ("nbins")
+		if (high < low) {
+		    step = low; low = high; high = step
+		}
+		step = (high - low) / l
+		call malloc (sum2, l, TY_REAL)
+		call malloc (sum4, l, TY_REAL)
+		call malloc (nsum, l, TY_INT)
+		call aclrr (Memr[sum2], l)
+		call aclrr (Memr[sum4], l)
+		call aclri (Memi[nsum], l)
+	    }
+	    call ap_noise (ap, gain, dbuf, nc, nl, c1, l1, sbuf, Memr[spec],
+		Memr[profile], nx, ny, xs, ys, Memr[sum2], Memr[sum4],
+		Memi[nsum], l, low, high)
+	    if (iap == naps) {
+		do k = 0, l-1 {
+		    m = Memi[nsum+k]
+		    if (m > 10) {
+			Memr[sum2+k] = sqrt (Memr[sum2+k] / (m - 1))
+			step = max (0., Memr[sum4+k] / m - Memr[sum2+k]**2)
+			Memr[sum4+k] = sqrt (sqrt (step / m))
+		    } else {
+			Memr[sum2+k] = 0.
+			Memr[sum4+k] = 0.
+		    }
+		}
+		call ap_nplot (image, in, Memr[sum2], Memr[sum4], l,
+		    low, high)
+		call mfree (sum2, TY_REAL)
+		call mfree (sum4, TY_REAL)
+		call mfree (nsum, TY_INT)
+	    }
+
+	    if (Memc[name] != EOS) {
+		call sprintf (Memc[str], SZ_LINE,
+		    "EXTRACT - Aperture %d from %s --> %s")
+		    call pargi (AP_ID(ap))
+		    call pargstr (image)
+		    call pargstr (Memc[name])
+		call ap_log (Memc[str], YES, YES, NO)
+		call ap_plot1 (gt, Memr[spec], ny, nsubaps)
+	    }
+	}
+
+	call gt_free (gt)
+	call sfree (sp)
+end
+
+
+# AP_SUM -- Simple, unweighted aperture sum.
+
+procedure ap_sum (ap, dbuf, nc, nl, c1, l1, sbuf, nx, ny, xs, ys, spec,
+	nsubaps, asi)
+
+pointer	ap			# Aperture structure
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+pointer	sbuf			# Sky values (NULL if none)
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Origin of sky array
+real	spec[ny, nsubaps]	# Spectrum
+int	nsubaps			# Number of subapertures
+pointer	asi			# Interpolator for edge pixel weighting
+
+int	i, ix, iy, ix1, ix2
+real	low, high, step, x1, x2, wt1, wt2, s, sval, skyval
+real	ap_cveval()
+pointer	cv, data, sky
+errchk	asifit
+
+begin
+	i = AP_AXIS(ap)
+	low = AP_CEN(ap,i) + AP_LOW(ap,i)
+	high = AP_CEN(ap,i) + AP_HIGH(ap,i)
+	step = (high - low) / nsubaps
+	cv = AP_CV(ap)
+	do iy = 1, ny {
+	    s = ap_cveval (cv, real (iy + ys - 1)) - c1 + 1
+	    call ap_asifit (dbuf+(iy+ys-1-l1)*nc, nc, xs[iy]-c1+1,
+		low+s, high+s, data, asi)
+#	    data = dbuf + (iy + ys - 1 - l1) * nc + xs[iy] - c1 - 1
+#	    if (asi != NULL)
+#		call asifit (asi, Memr[data], nc-xs[iy]+c1)
+	    do i = 1, nsubaps {
+	        x1 = max (0.5, low + (i - 1) * step + s) + c1 - xs[iy]
+	        x2 = min (nc + 0.49, low + i * step + s) + c1 - xs[iy]
+	        if (x2 <= x1) {
+		    spec[iy,i] = 0.
+		    next
+	        }
+	        ix1 = nint (x1)
+	        ix2 = nint (x2)
+
+		# Compute end pixel weights.  Remember asi is offset by 1.
+		call ap_edge (asi, x1+1, x2+1, wt1, wt2)
+
+		# Sum pixels.
+		sval = wt1 * Memr[data+ix1] + wt2 * Memr[data+ix2]
+		do ix = ix1+1, ix2-1
+		    sval = sval + Memr[data+ix]
+
+		# Subtract sky if desired.
+	        if (sbuf != NULL) {
+		    sky = sbuf + (iy - 1) * nx - 1
+		    skyval = wt1 * Memr[sky+ix1] + wt2 * Memr[sky+ix2]
+	            do ix = ix1+1, ix2-1
+		        skyval = skyval + Memr[sky+ix]
+		    sval = sval - skyval
+	        }
+
+		# Save extracted pixel value.
+		spec[iy,i] = sval
+	    }
+	}
+end
+
+
+# AP_EDGE -- Compute edge weights.
+
+procedure ap_edge (asi, x1, x2, wt1, wt2)
+
+pointer	asi			#I Image interpolator pointer
+real	x1, x2			#I Aperture edges
+real	wt1, wt2		#I Weights
+
+int	ix1, ix2
+real	a, b
+real	asieval(), asigrl()
+
+begin
+	    # Edge pixel centers.
+	    ix1 = nint (x1)
+	    ix2 = nint (x2)
+
+	    # Default weights are fractions of pixel.
+	    if (ix1 == ix2) {
+		wt1 = (x2 - x1)
+		wt2 = 0
+	    } else {
+		wt1 = (ix1 - x1 + 0.5)
+		wt2 = (x2 - ix2 + 0.5)
+	    }
+
+	    # If there is an interpolator compute fraction of integral.
+	    # We require that data and integrals be positive.
+	    if (asi != NULL) {
+		if (asieval (asi, real(ix1)) > 0) {
+		    b = asigrl (asi, ix1-0.5, ix1+0.5)
+		    if (b > 0) {
+			if (ix1 == ix2)
+			    a = asigrl (asi, x1, x2)
+			else
+			    a = asigrl (asi, x1, ix1+0.5)
+			if (a > 0 && a < b)
+			    wt1 = a / b
+		    }
+		}
+		if (ix1 != ix2 && asieval (asi, real(ix2)) > 0) {
+		    b = asigrl (asi, ix2-0.5, ix2+0.5)
+		    if (b > 0) {
+			a = asigrl (asi, ix2-0.5, x2)
+			if (a > 0 && a < b)
+			    wt2 = a / b
+		    }
+		}
+	    }
+end
+
+
+# AP_STRIP -- Simple, unweighted aperture strip.
+# Interpolate so that the lower edge of the aperture is the first pixel.
+
+procedure ap_strip (ap, aplow, aphigh, out, dbuf, nc, nl, c1, l1, sbuf, nx, ny,
+	xs, ys)
+
+pointer	ap			# Aperture structure
+real	aplow, aphigh		# Aperture limits
+pointer	out			# Output IMIO pointer
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+pointer	sbuf			# Sky values (NULL if none)
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Origin of sky array
+
+int	i, na, iy, ix1, ix2, nasi
+real	low, high, s, x, ap_cveval(), asieval()
+pointer	obuf, cv, asi, data, sky, ptr, imps2r()
+
+begin
+	i = AP_AXIS(ap)
+	low = aplow - c1 + 1
+	high = aphigh - c1 + 1
+	cv = AP_CV(ap)
+	call asiinit (asi, II_LINEAR)
+
+	na = IM_LEN(out,2)
+	obuf = imps2r (out, 1, ny, 1, na)
+	call aclrr (Memr[obuf], na * ny)
+
+	do iy = 1, ny {
+	    i = iy + ys - 1
+	    s = ap_cveval (cv, real (i))
+	    ix1 = max (1, nint (low + s) - 1)
+	    ix2 = min (nc, nint (high + s) + 1)
+	    nasi = ix2 - ix1 + 1
+	    if (nasi < 3)
+		next
+	    data = dbuf + (i - l1) * nc + ix1 - 1
+	    iferr (call asifit (asi, Memr[data], nasi))
+		next
+	    
+	    x = low + s - ix1 + 1
+	    ptr = obuf + iy - 1
+	    if (sbuf == NULL) {
+	        do i = 1, na {
+		    if (x >= 1 && x <= nasi)
+		        Memr[ptr] = asieval (asi, x)
+		    x = x + 1.
+		    ptr = ptr + ny
+	        }
+	    } else {
+		sky = sbuf + (iy - 1) * nx + nint (low + s) - xs[iy] + c1 - 2
+	        do i = 1, na {
+		    if (x >= 1 && x <= nasi)
+		        Memr[ptr] = asieval (asi, x) - Memr[sky+i]
+		    x = x + 1.
+		    ptr = ptr + ny
+	        }
+	    }
+	}
+
+	call asifree (asi)
+end
+
+
+# AP_PSTRIP -- Profile based strip.
+# Interpolate the profile spectrum so that the lower aperture edge is the
+# first pixel.
+
+procedure ap_pstrip (ap, aplow, aphigh, out, gain, spec, profile, nx, ny,
+	xs, ys)
+
+pointer	ap			# Aperture structure
+real	aplow, aphigh		# Aperture limits
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+real	spec[ny]		# Spectrum
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Origin of profile array
+
+int	na, ix, iy
+real	low, high, s, x, ap_cveval(), asieval()
+pointer	sp, cv, asi, data, impl2r()
+
+begin
+	call smark (sp)
+	call salloc (data, nx, TY_REAL)
+
+	ix = AP_AXIS(ap)
+	low = aplow
+	high = aphigh
+	cv = AP_CV(ap)
+	na = IM_LEN(out,2)
+	call asiinit (asi, II_LINEAR)
+
+	do iy = 1, ny {
+	    s = spec[iy] / gain
+	    do ix = 1, nx
+		Memr[data+ix-1] = s * profile[iy,ix]
+	    call asifit (asi, Memr[data], nx)
+	    s = ap_cveval (cv, real (iy+ys-1)) - xs[iy] + 1
+	    x = low + s
+	    do ix = 1, na {
+		profile[iy,ix] = asieval (asi, x)
+		x = x + 1
+	    }
+	}
+
+	do ix = 1, na
+	    call amovr (profile[1,ix], Memr[impl2r(out,ix)], ny)
+
+	call asifree (asi)
+end
+
+
+# AP_ASIFIT -- Return interpolation pointer and data pointer.
+#
+# The main reason for this routine is to shift the origin of the data by
+# one pixel so that the interpolator may be called to evaluate across
+# the extent of the first and last pixels.  This means the calling program
+# will reference asi fit between 1.5 and N+1.5.  It also means the returned
+# data pointer may start before the first point but will never be
+# dereferenced outside of the data range.
+
+procedure ap_asifit (dbuf, nc, xs, low, high, data, asi)
+
+pointer	dbuf			#I Data buffer pointer
+int	nc			#I Size of data buffer
+int	xs			#I Start of aperture array (in dbuf coords)
+real	low			#I Low aperture edge (in dbuf coords)
+real	high			#I High aperture edge (in dbuf coords)
+pointer	data			#O Data pointer
+pointer	asi			#I ASI pointer
+
+int	i, ix1, ix2, n
+real	x1, x2
+pointer	fit
+
+begin
+	# Check for in bounds data.
+	x1 = max (0.5, low)
+	x2 = min (nc + 0.49, high)
+	if (x1 >= x2)
+	    return
+
+	# Set data pointer relative to the aperture start with an offset for
+	# one indexing; i.e. pixel i is referenced as Memr[data+i].  The
+	# aperture start may put this outside the data buffer but we expect
+	# routines using the pointer to never index outside of the buffer.
+
+	data = (dbuf + xs - 1) - 1
+
+	# If not using an interpolator we are done.
+
+	if (asi == NULL)
+	    return
+
+	# If the aperture, with one extra pixel on each end for integration
+	# across the end pixel, is within the data buffer then fit an
+	# interpolator directly.  Otherwise we need to use a temporary
+	# padded buffer.  The origin of the fitted buffer is relative
+	# to the data pointer.  Note that this means that evaluating the
+	# fit requires the aperture start coordinates to be incremented
+	# by 1.
+
+	ix1 = 0
+	ix2 = nint (x2) + 1 - (xs - 1)
+	n = ix2 + ix1 + 1
+	if (data + ix1 >= dbuf && data + ix2 <= dbuf + nc - 1) {
+	    call asifit (asi, Memr[data+ix1], n)
+	    return
+	}
+
+	# One or the other end point is out of bounds so to avoid potential
+	# NAN and segmentation errors use an internal array to pad.
+
+	call malloc (fit, n, TY_REAL)
+	do i = 0, n-1 {
+	    if (data + i < dbuf)
+		Memr[fit+i] = Memr[dbuf]
+	    else if (data + i > dbuf + nc - 1)
+		Memr[fit+i] = Memr[dbuf+nc-1]
+	    else
+		Memr[fit+i] = Memr[data+i]
+	}
+	call asifit (asi, Memr[fit], n)
+	call mfree (fit, TY_REAL)
+end