Initial commit

17 files changed, 5084 insertions, 0 deletions

diff --git a/noao/twodspec/multispec/doc/MSalgo.ms b/noao/twodspec/multispec/doc/MSalgo.ms
new file mode 100644
index 00000000..0c64e2b3
--- /dev/null
+++ b/noao/twodspec/multispec/doc/MSalgo.ms
@@ -0,0 +1,1032 @@
+.de FX
+.nf
+.ps -2
+.ss 25
+.cs R 25
+..
+.de EX
+.ps +2
+.ss
+.cs R
+.fi
+..
+.EQ
+delim $$
+.EN
+.RP
+.TL
+Algorithms for the Multi-Spectra Extraction Package
+.AU
+Francisco Valdes
+.K2
+.TU
+.AB
+The algorithms for the Multi-Spectra Extraction Package (\fBmultispec\fR)
+in the Image Reduction and Analysis Facility (\fBIRAF\fR) is described.
+The basic aspects of the general two dimensional aperture spectra extraction
+problem are first discussed.
+The specific algorithms for extraction of multi-aperture plate and
+Echelle digital data are presented.  Results of the authors experiments
+with this type of data are included.
+The detailed specification of the package is given in a second document,
+\fIDetailed Specifications for the Multi-Spectra Extraction Package\fB.
+.AE
+.NH
+Introduction
+.PP
+There are an increasing number of astronomical instruments which produce
+multiple spectra on a two dimensional detector.
+The basic concept is to use one dimension for wavelength,
+the dispersion dimension, and the other, the cross dimension, for
+packing additional information during a single exposure.
+For example, the cross dimension can be different objects or
+different spectral orders.  The classic multi-spectra instrument is
+the Echelle spectrograph.  New instruments are the aperture plate and
+Medusa spectrographs.
+.PP
+There is an additional aspect of the multi-spectra format; namely,
+the individual spectra can contain spatial data.  An example of
+this would be multiple slit spectra in which each slit spectra contains
+sky signal and object signal.  In the following
+discussion we limit the spectra to be simple aperture spectra in
+which we only desire to sum the intensities at each wavelength to form
+a one dimensional spectrum.
+.PP
+The analysis of multi-spectra aperture data consists of two steps; the
+separation and extraction into individual aperture spectra
+and the calibration and measurement of the spectra.  These steps can
+either be incorporated into one analysis package or two separate
+packages.  There are advantages to the first approach since some
+aspects of the individual spectra are directly related by the physical
+geometry of the multi-spectra format.  However, because packages for
+the analysis of individual spectra exist we begin by dividing the
+reduction into separate extraction and analysis tasks.  It is
+important to realize, however, that the existing analysis tools are not well
+suited to reducing the larger number of spectra and treating sets of
+spectra together.
+.PP
+The latter part of this paper describes the algorithms for the
+extraction of two types of data; the multi-aperture plate (MAP)
+and Echelle used with digital detectors.  However,
+it is important to keep the more general problem in mind
+and the remainder of this introduction considers the different
+conceptual aspects of the multi-spectra extraction task.
+Table 1 lists many of the general properties of multi-spectra aperture data.
+The other two columns give possible alternatives that each property may take.
+
+.TS
+center box;
+c s s
+c s s
+c c s
+= = =
+c || c | c.
+Table 1:  Aspects of Multi-Spectral Data
+
+Property	Alternatives
+detector	digital	photographic
+alignment	aligned	skewed
+blending	blended	unblended
+aperture	holes	slits
+spectral resolution	low	high
+.TE
+
+.PP
+The detector determines what kind of calibration procedures are
+required to produce intensity from the measurements.
+A digital detector requires sensitivity calibrations on all scales.
+This is the "flat field" problem.  There are also corrections for
+bias and dark current.  Photographic detectors require
+intensity calibration.  Data which are not aligned with the natural
+dimensions of the digital image require extra procedures.  Two types
+of non-alignment are a rotation of the dispersion dimension relative
+to the pixel dimension and a "wiggling" or "snaking" of the dispersion
+dimension.  Blending refers to the degree of contamination along the
+cross dimension.  Blended data requires extra effort to correct for
+the overlap between different spectra and to determine the background.
+The aperture defines the extent of the spectra in the cross dimension.
+The two most relevant choices are holes and slits.  In some
+instruments, like the Echelle, the size of the aperture can be varied
+at the time of the observations.  Finally, the spectral resolution is
+important in conjunction with digital detectors.  If the resolution is
+high then quartz flat field calibrations are relatively easy because
+the spectral
+signature need not be considered.  Otherwise, the flat field problem
+is more difficult because the gain variations of the detector
+must be separated from the natural spectral intensity variation of the
+quartz.
+.PP
+There is always some confusion of terms when talking about multi-spectra
+data; in particular, the terms x, y, line, and band.
+The image pixel dimensions are refered to as x and y.  We assume
+for the moment that the alignment of the multi-spectra format is such
+that x corresponds to the cross dimension and y to the dispersion
+dimension.  If this is not the case a rotation or interpolation
+program can be used to approximately orient the data in this way.
+A line is the set of intensity values as a function of x at constant y.
+In other words, a line is a cut across the dispersion dimension.
+A band is the average of more than one line.
+The image residing on disk will generally be organized
+such that x varies more rapidly and a line of the image is easily
+obtained.  In this form a display of the image will have the spectra
+running vertically.  In the Cyber extraction package the data is
+organized with x corresponding to the dispersion dimension.
+.NH
+Multi-Spectra Image Formats
+.PP
+The remainder of this paper will refer to two specfic and very
+different multi-spectra formats; the Kitt Peak Multi-Aperture Plate
+System and the Kitt Peak Echelle Spectrograph.
+.NH 2
+Kitt Peak Multi-Aperture Plate System
+.PP
+The reduction of data from multi-aperture plate observations is the 
+driving force for the development of a multi-spectra extraction
+package.  This application turns out to have most of the worse aspects
+of the properties listed in Table 1.  The multi-aperture plate spectrograph uses
+digital dectectors with low resolution, the spectra are blended and
+change alignment along the pixel dimension.  Furthermore, the camera
+has a variable point-spread function and focus,
+suffers from flexture problems, has a different illumination for
+the quartz than object exposures, and unexplained background level
+variations (in the CRYOCAM).  There are two detectors which have been
+used with the multi-aperture plate system, the Cryogenic Camera
+(CRYOCAM) and the High Gain Video Spectrometer (HGVS).
+.NH 2
+Echelle
+.PP
+As some of the algorithms were developed the Echelle data was brought
+to my attention.  It is considerably simpler than the MAP data because
+it is unblended and of high spectral resolution.
+Furthermore, the quartz exposure
+can be made wider than the object exposures and flexture is not a
+problem.  The principle problem in this data was the
+prevalence of cosmic rays.  It pointed to the need to maintain generality
+in dealing with both the MAP data and other types of
+multi-spectra data which have different profiles, may or may not be
+merged, and may or may not have different widths in quartz and object.
+Dealing with the cosmic ray problem lead to a very effective solution
+usable in both the Echelle and multi-aperture plate data.
+.NH
+User Level Extraction Logic
+.PP
+The user should generally only be concerned with the logical steps of
+extracting the individual spectra from the multi-spectra image.  This
+means that apart from specifying the detector system and the format
+he should not deal with details of the detector and the format.
+In the paper, 
+\fIDetailed Specifications for the Multi-Spectra Extraction Package\fB,
+the \fBIRAF\fR extraction package design and program specifications
+are described.
+.NH
+Flat Fields
+.PP
+There are two types of flat field situations depending on the spectral
+resolution.  When the resolution is high then the spectral signature of
+the continum calibration source, a quartz exposure, will be unimportant
+and variations in the signal will be due to detector sensitivity variations.
+In this case the quartz frame, or average of several frames, is the flat
+field and division of the object frames by the quartz frame is all that
+is required.  However, a special
+image division program is desirable to handle the region of low or absent
+signal between the the spectra.  This is described in section 4.2.
+.PP
+In the alternate case of lower resolution the quartz spectral signature is
+larger than the detector response variations.  A flat
+field in which the intrinsic quartz spectrum is removed is produced by
+assuming that the true value of a pixel is given by the smoothed average
+of the pixels near that point in position and wavelength and taking
+the ratio of the data value to the smoothed value.
+This requires a special smoothing program described in section 4.1.
+After the flat field is generated then the same image division
+program used for the Echelle data is applied.
+The image division and smoothing programs are general image operators and
+not specific to the Multi-Spectra Extraction Package.
+.NH 2
+MULTISPEC_FLAT
+.PP
+The multi-aperture plate data varies in both dimensions.  Thus, any averaging
+to smooth the image must take this variation into account.  In the Cyber
+a flat field for the multi-aperture plate data smooths across the dispersion
+by modeling the spectra.  This is a difficult task to do accurately because
+the true shape of the spectra is not known and the counts vary greatly
+and rapidly in this dimension.  This approach has the further difficulty
+that it is not possible to average several quartz exposures directly.
+.PP
+The alternate approach to modeling is statistical averaging.
+Averaging across the dispersion requires very high order polynomials
+because of the rapid variations;
+the spectra are typically spaced about 8 pixels apart and there are on
+the order of 50 spectra.  On the other hand, variations along the dispersion
+are much slower even if the spectra are slightly skewed; a bad case would
+have two peaks in 800 pixels along the y dimension.  This kind
+of variation is tractable with relatively simple averaging polynomials
+and is the one adopted for the multi-aperture plate data.
+.PP
+The flat fields are produced by a quadratic moving average along the
+y direction.  This means that the region centered at a given pixel
+is fitted by a least-squares quadratic polynomial and the value of the
+polynomial at that point is the appropriate statistical average.
+The width of the moving average is an adjustable parameter.
+At the edges of the frame where it is not possible to center a region of
+the specified width about the desired pixel the polynomial fit is used to
+extrapolate the average value to the edge.
+.PP
+Because the quadratic fit will
+be influenced by bad pixels an attempt is made to detect and smooth over
+the bad pixels.  This is accomplished by comparing the smoothed values to
+the observed values and ignoring pixels with a value of
+
+.EQ (1)
+	chi = | observed - smoothed | / sqrt smoothed
+.EN
+
+greater than a specified value.  Then the smoothing is recalculated and tested
+again for bad pixels.  This iteration continues until no further pixels
+are rejected.
+.PP
+Following the smoothing the flat field is produced by ratioing the raw quartz
+to the smoothed quartz.  Pixels of low signal (specified by the
+parameter \fIconstant\fR )
+are treated by the equation
+
+.EQ
+	r = (data + (constant - smoothed) ) / constant  .
+.EN
+
+The resultant flat field image is then divided into the object frames in
+the manner described in the next section.
+.PP
+Experience with data from the Cryogenic Camera has proved very good.
+The flat field which is produced can be examined on a display.  It
+shows fringing at red wavelengths and is not too strongly affected
+by bad pixels.  Some further effort, however, could go into smoothing
+over the bad pixels.
+.PP
+The smoothing operation on data from the Cryogenic Camera actually
+consists of four steps.  The quartz exposures are first averaged.
+The average quartz is rotated so that the dispersion
+direction is the most rapidly varying or x dimension.  Then the
+smoothing is performed along x followed by another rotation to return
+the flat field image to its original orientation.  The reason for the
+rotations is that they can be done quickly and efficiently whereas
+smoothing along the y dimension is very slow and coding an efficient
+version is much more complicated than doing a single line at a time.
+.NH 2
+FLAT_DIVIDE
+.PP
+The Echelle data has quartz frames which can be used directly as flat fields.
+One just has to divide the object frames by the quartz or average of several
+quartz.  However, in image division consideration has to be given the
+problem of division by zero or very small numbers.  In direct imaging this
+may not be much of a problem but in multi-spectra data the region between
+the spectra and near the edges of the spectra will have very low counts.
+Another aspect of image division for making flat field corrections is the
+scaling of the result.  The flat field integer image data must be large
+to give accurate relative response values.  However, one wants to divide
+an object frame by values near unity.
+This section describes a special image division operator allowing the user
+to specify how to handle these cases.
+.PP
+The parameters are a \fIdivision threshold\fR
+(default of zero) and a \fIthreshold violation value\fR.  Values of the
+denominator above the \fIthreshold\fR are treated separatedly from those
+below the \fIthreshold\fR.  The denominator image is scaled to have an
+average of one for pixels above the \fIthreshold\fR.  The pixel by pixel
+division is then performed for those points for which the denominator
+is above the \fIthreshold\fR.  Pixels for which the denominator is below the
+\fIthreshold\fR are set to the \fIthreshold violation value\fR in the resultant
+image if the \fIviolation value\fR is specified.  If the value is not
+specified then the numerator value is taken as the resultant value.
+The divisions can be done in place or the result put into a new image file.
+.PP
+For the multi-spectra situation where the object spectra have a
+smaller width than the quartz, as in the Echelle, one can either
+set the \fIthreshold
+violation value\fR to zero or not set it at all resulting in either
+exactly zero or background values between the spectra while still flattening
+the spectra.  This allows looking at the flattened spectra without the
+annoying "grass" between the spectra caused by dividing by small
+values.
+.NH
+Extraction
+.NH 2
+MULTIAP_EXTRACT
+.PP
+The extraction of spectra from multi-aperture plate images consists of
+a series of steps.  The steps are executed from a script.
+The command
+
+.FX
+ms> multiap_extract "ap165.*", "", 165, 50
+.EX
+
+will take the flattened images, ap165.*, from aperture plate 165 with 50
+spectra and automatically locate the spectra, model the profiles, and
+extract the one dimensional spectra.  The script consists of
+a number of steps as described below.
+.PP
+\fBFind_spectra\fR (section 6) initializes the \fBmultispec\fR data file
+and does a peak search to determine the initial positions of the
+spectra.
+\fBFind_bckgrnd\fR fits a polynomial of order 1 (or more) for the pixels which
+are not near the spectra as defined by \fBfind_spectra\fR.
+.PP
+The spectra are then modeled in bands of 32 lines by the model profiles
+described in section 8.1.  The first \fBmodel_fit\fR uses three Gaussian
+parameters for
+each spectra measuring the peak intensity, peak position, and width.
+The second \fBmodel_fit\fR adds a fourth parameter to modify the wings of the
+profile.
+.PP
+The \fBmodel_extract\fR program extracts the spectra line by line and also
+detects and removes cosmic ray events which do not fit the model
+profiles (see section 9).
+In outline, the extraction of blended spectral data uses the
+model profiles to determine the fraction of light
+from each of the neighboring spectra at the pixel in question.  The
+appropriate fraction of the
+.ul
+observed
+pixel intensity (minus the background) is
+assigned to the luminosities of the spectra.  There are two versions
+of the \fBmodel_extract\fR extraction.  The first simultaneously fits the
+peak intensity of all the spectra and the second uses the
+data value at the peak of each spectra to normalize the model.  The
+first method is slow and accurate and the second is fast and approximate.  
+Because the models are used in extraction only to define the relative
+contributions of neighboring spectra to the total observed pixel luminosity
+the speed of the approximate method far outweighs the need for
+accuracy.  However, cleaning the spectra of cosmic rays is a different
+matter and is discussed further in section 12.
+.PP
+After the extraction the spectra are correlated with the aperture plate
+description using \fBap_plate\fR (see Section 10) to determine the
+relative wavelength offsets and assign identification information to
+the spectra.
+.PP
+For successive frames it is not necessary to resort to the initial
+steps of finding the spectra and fitting from scratch.  The \fBcopy_params\fR
+routine makes a new copy of the fitting database.  Small shifts in positions
+of the spectra and the peak intensities are determined by doing a two
+parameter fit for the peak and position using the previously determined
+shape parameters.
+Changes in the shape of the spectra are then determined by the three
+and four parmater fits.  Because the solution is likely to be close to
+the previously determined shape the transfering of one solution from a
+previously solved image is faster than starting from scratch.
+Note that the shapes as well as the positions and peak intensities
+of all frames including the object exposures are allowed to change.
+.PP
+The spectra are then extracted from the image by \fBmodel_extract\fR and the
+process repeats for the succeeding images.
+.PP
+One useful feature is the ability to specify the bands or lines to be
+modeled or extracted.
+This feature is useful for diagnosising the programs quickly.
+The default is all bands or lines.
+.NH 2
+ECHELLE_EXTRACT
+.PP
+The extraction of the unblended Echelle spectra is performed
+begins in a similar way with \fBfind_spectra\fR and \fBfind_bckgrnd\fR.
+The extraction and cleaning, however, uses \fBstrip_extract\fR which
+adds up the instrumental counts for each unblended spectra at each
+wavelength to get the total luminosity.
+.NH
+FIND_SPECTRA -- Finding the Spectra
+.PP
+The first step in the extraction and processing of multi-spectra data is
+to locate the spectra.  This can be done interactively by
+the user but it is far preferable to automate the process;
+particularly since multi-spectra data can have a large number of
+spectra and frames.  The approach is to find the peaks in a line, or
+average of lines, sort the peaks found in some manner, such as by
+strength, and select the expected number of peaks from the top of the
+list.
+Beyond this simple outline are several algorithmic details such as how
+to define and locate valid peaks and how to sort the list of peaks.
+Peak finding is a general problem and a subroutine for peak finding is
+described below.  The \fBfind_spectra\fR program provides an
+interface between the \fBmultispec\fR data file and the
+general peak finding algorithm.
+.PP
+The \fBpeaks\fR function takes arrays of x (position) and y (value) points
+and the number of
+points in the arrays and returns the number of peaks found.  It also
+returns the estimated positions of the peaks in the x array and the
+extimated peak values in the y array in order of peak likelihood.
+There is one user parameter, the smoothing \fIwidth\fR.
+The choice of the \fIwidth\fR parameter is dicatated by how closely and how
+wide the peaks are to be expected.
+The algorithm takes a region of \fIwidth\fR points
+centered on each x point and fits a quadratic;
+
+.EQ
+y sub fit = a + b x + c x sup 2~.
+.EN
+
+A peak is defined
+when the slopes, $b sub 1$ and $b sub 2$, of two neighboring points
+$x sub 1$ and $x sub 2$ change
+sign from positive to negative and the curvatures, $c sub 1$ and $c
+sub 2$, are less than -0.001 for both points.
+The quadratic polynomials define two estimated peak positions
+
+.EQ
+x sub 1 sub peak = x sub 1 - b sub 1 / (2 * c sub 1 ),~~
+x sub 2 sub peak = x sub 2 - b sub 2 / (2 * c sub 2 )~.
+.EN
+
+The offsets are then normalized to give a linear interpolation
+fraction
+$f = ( x sub 1 sub peak - x sub 1 ) / ( x sub 2 sub peak - x sub 1 sub
+peak )$ in the interval between the two points.
+The estimated position of the peak is then
+
+.EQ
+x sub peak = f * ( x sub 1 - x sub 2 )
+.EN
+
+and the estimated peak value is the average value of the two quadratic
+polynomials at $x sub peak$.  The curvature at the peak is
+estimated by $c sub peak = c sub 1 + f * (c sub 1 - c sub 2 )$.
+Finally, the peaks are sorted by the magnitude of the peak curvature.
+.PP
+This peak finding algorithm works quite well.  I have also used it to
+automatically locate peaks in the extracted one dimensional spectra
+and then do peak correlations between spectra to find a relative
+wavelength solution.  Some such use of this program may be implemented
+in either future additions to the Multi-Spectra Extraction Package or
+the Spectral Reduction Package.
+.PP
+In \fBfind_spectra\fR the number of spectra to be found is specified by
+the user.  The user should have previously looked at an image
+on a display or done a profile plot across the
+dispersion to count the observed spectra.
+Additional parameters specify the columns in which the spectra
+are to be found and the minimum separation and width of the spectra.
+The column specification allows the elimination of problems with defective
+areas of the detector such as the LED in the Cryogenic Camera.  The minimum
+width and separation provide algorithmic constraints to the spectra finding
+procedure.
+.PP
+The peaks are found at two or more points in the
+multi-spectra image for a band of 32 lines using a
+\fBpeaks\fR \fIwidth\fR parameter of 5.  After the peaks are found
+at a number of bands in the image a linear fit is made to determine any small
+slope of the spectra relative to the columns.
+The reason for specifying only a few bands is that the process of
+finding the peaks is moderately slow and only two bands are needed for
+determining the initial position angle of the spectra to the y
+dimension.  Furthermore, some bands do not give a satisfactory result
+because of extraneous data such as the LED in the CRYOCAM or bad
+focus.  Another possibility is that a spectrum may go off the edge
+and in order to "find" the spectrum for partial extraction bands that
+include the on image part of the spectrum must be specified.
+.NH
+FIND_BCKGRND -- Background
+.PP
+The background on a multi-spectra image is the result of very broad
+scattering as opposed to the narrower scattering which produces
+distinguishable wings on individual spectra.
+Modeling of the background in a Cryogenic Camera multi-aperture plate
+image shows that the background is well explained by a broad
+scattering function.
+It is not reasonable, however, to model the scattering to this detail
+in actual extractions.
+Instead a smooth polynomial is fitted to the pixels not covered by
+spectra.  The order of the polynomial is a specified parameter.
+For the CRYOCAM MAP data a quadratic is appropriate.
+.PP
+The algorithm is the same for all multi-spectra data except for the
+choice of parameters.  First, the location of the spectra must be
+determined.  This is done by the \fBfind_spectra\fR program.  There
+are two principle parameters, a buffer distance and the order of the
+fitting polynomial.  Each line, or average of several lines, is fitted
+by least-squares for the points lying farther than the buffer
+distance from any spectra.  If there are no points which completely
+stradle the spectra, i.e. located on each side of the spectra, then
+the order of the fitting polynomial is ignored and a constant, or
+first order polynomial, is determined.
+A hidden parameter specifying the columns allowed for searching for
+background points is available so that bad parts of the image can be
+ignored.
+.PP
+A difference in philosophy with the Cyber programs is that the
+determined background is not subtracted from the image data.  It is
+instead kept in the database for the image.  Generally, it is better to
+modify the basic image data as little as possible.  This is the approach
+used in the Multi-Spectra Extraction Package.
+.NH
+Spectra Profiles
+.NH 2
+MODEL_FIT -- Models for Multi-Spectra Images
+.PP
+The object of modeling is to separate blended spectra for extraction
+and to remove artifacts, such as cosmic rays, which do not fit
+the model.  The models should have the minimum number of parameters
+which give residuals approaching the detector statistics, they
+should incorporate constraints  based on the physics of the
+detector/camera system,  and the models must be ammenable to a
+statistical fitting algorithm which is stable.
+There are a large number of possibilities.
+.PP
+An important point to bear in mind during the following discussion is
+the necessary accuracy of the model fitting.  In the design proposed
+here the model fitting is not used for determining the smooth quartz.
+Use of a model for making a flat field would require a very accurate
+model and using an average quartz is not possible.  However, for
+extraction the model is used only to indicate the
+relative fraction of light for each spectrum when the spectra are
+blended.  The cleaning application is more critical but not nearly so
+much as in the flat field modeling.  Thus, though we do a good job of
+model fitting (better the the Cyber modeling) some additional features
+such as smoothing along the spectra are not included.
+Also, though some improvement can be gained by the additional shape
+parameters in the fit, they are not necessary for the required purpose
+and can be left out leading to a faster extraction.
+.PP
+During the course of my investigation I tried more than one hundred
+models and combinations of constraints.  Some general results of this
+study follow.
+The model which I find gives the best results has six parameters not
+counting the background.  The model is defined by the following
+equations where x is the cross dimension.
+
+.EQ (1)
+I = I sub 0 exp (-s * ( DELTA x sup 2 ))
+.EN
+.EQ
+DELTA x = (x - x sub 0 )
+.EN
+.EQ
+s = s sub 0 + s sub 1 y sup 2 + s sub 2 y sup 3
+.EN
+.EQ
+y = DELTA x / sqrt { DELTA x sup 2 + x sub 1 sup 2 }
+.EN
+
+The model consists of a intensity scale parameter, $I sub 0$,
+and a profile which is
+written in a Gaussian form.  The center of the profile is given by
+the parameter $x sub 0$.  The profile is not exactly Gaussian because the
+scale, $s$, is not a constant but depends on $DELTA x$.  The scale
+function has three terms; a constant term, $s sub 0$, which is the scale
+near the center of the profile, and even and odd terms, $s sub 1$
+and $s sub 2$,
+which change the scale in the wings of the profile.
+.PP
+The characteristic of the profile which must be satisfied is that at
+large distances from the profile center the scale is positive.  This
+requirement means that the profile will be monotonically decreasing at
+large distances and will have a finite luminosity.  This point was
+crucial in determining the form of the scale function.  A straight
+power series in $DELTA x$ does not work because power series diverge.
+Instead, the scale function is defined in terms of a separation
+variable $y$ which is bounded by -1 and 1 at infinite separation and is
+zero at zero separation.  The parameter $x sub 1$ defines a characteristic
+distance where the character of $y$ changes from going as $DELTA x$ to
+asymptotic to 1.  The parameters are, thus, $I sub 0$, $x sub 0$, $s sub 0$,
+$s sub 1$, $s sub 2$, $x sub 1$.
+.PP
+An important property of this model is that the terms have a physical
+interpretation.  The profile scale and profile center are obvious and
+any model must include them.  It is the remaining terms, $s sub 0$,
+$s sub 1$, $s sub 2$,
+and $x sub 1$, which are called the shape parameters, which are interesting.
+In an ideal aperture plate system the shape of a profile would be
+given by the projection of the circular aperture into the cross dimension:
+
+.EQ
+P( DELTA x ) = sqrt {1 - a DELTA x sup 2}
+.EN
+
+where the constant a is related to the size of the hole by 
+
+.EQ
+a = 1 / r sup 2
+.EN
+
+For small $DELTA x$ the profile can be expressed in the Gaussian form with
+a scale
+
+.EQ
+s = a( 1/2  + a DELTA x sup 2 + ...)
+.EN
+
+Thus, even in a perfect aperture plate system a Gaussian form shows the
+scale increasing from a central value determined by the size of the hole.
+In other words, the profile decreases more sharply than a Gaussian.
+.PP
+However, no aperture plate system is ideal because the thickness of
+the aperture plate is finite and there is scattering and changes in
+the focus of the system.  One must
+convolve the profile above with a scattering/focus function.  One can show
+that for reasonable functions, exponential and Gaussian,
+with some scale b the final profile is a function of the ratio b/a.
+If the ratio is less than 1 then the profile will be more like that of
+the hole and the profile will be sharper than a Gaussian in the wings.
+If the ratio is much greater than 1 then the profile will become the
+scattering profile at large separations.  Simulations using Gaussian
+and exponential scattering profiles show behaviors very much like the
+profile (1) with $s sub 1$ greater than zero when b/a < 1
+meaning the profile becomes sharper (than a Gaussian) in the wings
+and $s sub 1$ < 0 when b/a > 1.
+Thus, $s sub 1$ defines the scale of the scattering profile relative
+to the hole size.
+The size of the hole is incorporated into the parameter $x sub 1$.
+The parameter $s sub 2$ allows an asymmetry in the profile.
+.PP
+An interesting property of the scale function is that it is all right
+for it to be negative at small distances from the profile center.  This
+occurs when $s sub 0$ is negative.  The effect of this, provided $s$
+becomes positive at large distances, is to give a two horned profile.
+This, in fact, is observed when the focus of the system becomes very
+poor.
+.PP
+The best fits (least chi-square or rms residual) are
+obtained when each spectrum at each wavelength has independent
+parameters.  However, this sometimes gives rise to unphysical results.
+If left entirely unconstrained the parameter fitting algorithm can
+make one line broad and dominant and a neighboring line weak and
+sharp.
+This is not, of course, a property of the camera or detector.
+Thus, constraints based on the physics of the
+camera/detector are used.  This means that the shape
+parameters $s sub 0$, $s sub 1$, $s sub 2$, and $x sub 1$
+are coupled locally by making them vary as a polynomial of position
+across the dispersion.  One might also
+constrain the variation of the shape along the spectrum as is done in
+the Cyber.  This is not needed because there are no drastic differences
+between the fitted parameters at neighboring points along the spectra.
+.PP
+My experience with the Cyrogenic Camera system has shown the
+following.  The focus ripples twice across the CCD with the
+propagation angle being approximately 30 degrees from the long dimension.
+The change in focus is partly just a scale change.  This is seen in
+the correlation of $s sub 0$ with the image scale found by \fBap_plate\fR.
+The shape parameter $s sub 1$ changes sign from positive to
+negative indicating that when the focus is good the profile
+decreases faster than a Gaussian and when the focus is bad it decreases
+slower.  Occassionally the focus is very bad and $s sub
+0$ is negative and $s sub 1$ is small and positive causing a broad two
+horned profile.  The
+assymetry parameter, $s sub 2$, is useful only when the signal is strong near
+the peak of a quartz exposure.  It is not really necessary to include
+it in the model fits.  The assymetry parameter was dominant, however,
+in some Echelle data which were clearly asymmetric.  The value of
+$x sub 1$ is
+not highly sensitive and can be fixed for a given hole size.  Large
+changes in the hole size would require resetting $x sub 1$.
+The use of the four parameters, $I sub 0$, $x sub 0$, $s sub 0$,
+and $s sub 1$, allow good fits
+to all the data I've examined including those in which the peak/valley
+intensity ratio across the spectra is about 1.1.  It is the importance
+of the parameter $s sub 1$ which improves the fitting dramatically over the
+Cyber three parameter fitting (in addition to a different fitting
+algorithm).
+.PP
+The heart of profile fitting is the solution of the multi-parameter
+least-squares problem.  In a blended multi-spectra image the profile
+parameters of one spectra are affected by its neighbors which are,
+in turn, affected by their neighbors and so forth.  The key to this
+type of problem is to realize that only nearest neighbors affect the
+profile parameters and this leads to a "banded" least-squares matrix.
+A banded matrix is one in which cross terms far from the diagonal are
+zero.  Solution of banded matrices is much more efficient than solving
+the entire matrix.  This allows solution for more than 100 parameters
+simultaneously in a short time.
+.PP
+Use of the banded multi-parameter solution has the restriction, however,
+that there can be no parameters in the model which are not local to
+the profiles.  This affects the way
+global constraints are applied to the parameters.  In particular,
+the way the shape parameters are constrained to vary smoothly across the
+detector.
+The shape parameters are first found as independent parameters by the
+banded matrix solution and then smoothed by a polynomial in x.
+.PP
+An area which was extensively investigated was the appropriate
+weighting to use for the model fitting.  The most likely choices are
+weighting by $1 / sqrt data$ and unit weight corresponding to
+$chi sup 2$
+and least squares fitting.  It was found that the two methods
+agreed fairly closely but that the least squares fitting was more
+appropriate because the blending correction depends largely on the
+value of the peak intensity and less on the exact fit of the wings.
+With $chi sup 2$ the peak is fit with less accuracy in order to improve
+the fit in the wings of the profile.  In some cases this gave clear
+errors in estimating the peak intensity and, hence, the proper contributions
+between the blended spectra were not made.
+.PP
+Now follows the details of the fitting algorithm.
+The algorithm is a series of script steps in \fBmultiap_extract\fR
+which call the model fitting program \fBmodel_fit\fR with different
+parameters.  In the script all bands are fit, $x sub 1$ is fixed,
+and the asymmetry shape parameter $s sub 2$ is ignored.
+The four parameter fit is applied to bands of 32 lines.  The band
+solutions are linearly interpolated to the full image and then only
+the intensity scale parameter is calculated for each line during the
+extraction of the spectra with \fBmodel_extract\fR.
+.PP
+The general fitting scheme proceeds as follows:
+.LP
+1.  Fit the three parameters $I sub 0$, $x sub 0$, $s sub 0$ with
+$x sub 1$ fixed and $s sub 1$ and $s sub 2$
+zero.  This is precisely a Gaussian fit.  The three parameters are
+determined simultaneously for all the lines at once using the banded
+matrix method.  Thus for 50 lines the solution has 150 variables.
+After each fit the scale
+parameter $s sub 0$ is smoothed by a polynomial in x.  The polynomial is
+taken with seven terms.
+.LP
+2.  Once the improvement in each iteration becomes less than a
+specified amount (2% in rms residual) the next parameter $s sub 1$ is added.
+The solution has two steps: fit for $s sub 0$ and $s sub 1$ with $I sub 0$
+and $x sub 0$ fixed and
+then fit $I sub 0$ and $x sub 0$ with $s sub 0$ and $s sub 1$ fixed.  As before the scale terms
+are smoothed by a seventh order polynomial.  Attempts to solve for all
+four parameters a once gave unstable results for reasons I don't
+understand.
+.LP
+3.  If desired, the last shape parameter $s sub 2$ can be added by solving
+for $s sub 0$, $s sub 1$, and $s sub 2$ while holding $I sub 0$ and
+$x sub 0$ fixed and then solving for
+$I sub 0$ and $x sub 0$.  This provides some improvement when the signal is very
+strong but is generally not needed in the multi-aperture plate data.
+It can be an important term in the Echelle data.
+.LP
+4.  It is possible to then adjust $x sub 1$ followed by steps 2 or 3.
+However, this gives very little improvement and $x sub 1$ should be fixed for
+each type of data.
+.LP
+5.  During the final extraction when individual lines are evaluated a one
+parameter fit is used to find $I sub 0$ for each spectra.  This is
+rather slow, however, on the order of 3 hours per frame.  By using
+the pixel value near $x sub 0$ as the value for $I sub 0$ the extraction is reduced
+to 13 minutes per frame (see section 12).
+.PP
+In addition to the preceeding steps the fitting algorithm applies some
+heuristic constraints.  These constraints limit how far the peak positions can
+shift in each iteration, require the peak intensity to remain positive, and
+limit the scale function to be positive at large values of y.
+.NH 2
+STRIP_EXTRACT -- Unblended Profiles
+.PP
+For unblended multi-spectra data the profiles can be anything.  The profiles
+are obtained by averaging a number of lines (say 32) and normalizing
+at some point like the peak value.  These profiles are then used for
+detecting bad pixels, such as cosmic rays, and correcting for them as
+discussed in the section on cleaning.  Modeling using the \fBmodel_fit\fR
+program is only used on Echelle data to find peak positions
+accurately in order to follow any curvature of the spectra.
+.NH
+Extraction and Cleaning
+.PP
+The extraction of spectra are done separately from the modeling.  It is
+possible to extract spectra without any modeling at all using
+\fBstrip_extract\fR.  The extraction step also allows the user to specify
+if cleaning of the spectra for cosmic rays is desired.  Also modifying
+the image is an option.
+.NH 2
+MODEL_EXTRACT
+.PP
+Extraction and cleaning using a model fit is described here.
+First the $I sub 0$ values for the model profiles are determined for
+all the spectra in a line either by multi-parameter fitting or by
+taking the peak value.  The pixel values are then compared to the
+model in a chi-squared way:
+
+.EQ
+r = (data - model) / sqrt model
+.EN
+
+If the value of r is larger than a set amount, say 5, then the pixel
+value is set to that of the model.  Since the "bad" pixel may affect
+the intensity scale $I sub 0$ the cleaning is iterated until no further
+pixels are changed.
+.PP
+The fitting of the data from an individual line of data to the model profiles
+is the key element in this scheme.  The best method is to use all the
+data in a multi-parameter least square fit.  This minimizes the effect
+of bad pixels on the estimated profile which is the essence of this
+cleaning method.  However, while the time required to do this for one
+line is not great, it adds up to nearly three hours for the 800 lines
+in a CRYOCAM frame.  A quick alternative is to scale the model profile
+by the data value at the peak position.  This is
+quite fast.  However, if the peak has a cosmic ray event or is
+otherwise bad then the estimated profile will not correspond closely
+to the data profile and the cleaning procedure will make gross errors.
+The limited experience I've had with the Echelle and MAP data
+has worked well with using the peak estimate.  However, the possible
+problems make me nervous and some compromise based on using more than
+the peak to estimate the intensity scale of the profile to the data
+needs to be found.  This is important because much of the feedback on
+the \fBmultispec\fR package from Paul Hintzen and Caty Pilachowski
+have dealt with
+the particular usefulness of a good cosmic ray cleaning algorithm in
+extracting multi-spectra data.
+.NH 2
+STRIP_EXTRACT
+.PP
+Removing cosmic rays is the major part of Echelle extraction.
+Because these are unblended spectra of arbitrary shape a strip
+extraction is all that is needed.
+The cleaning is done by the same algorithm used for the multi-aperture
+plates except that the profiles are found, as described earlier, by
+averaging a number of lines.
+The intensity scaling is determined from either a least-square fit
+or the peak value.
+The same question about the appropriate way to
+determine the fit of the profiles to the data discussed previously
+applies here except since the spectra are not blended the spectra
+can be treated separately in any least square fitting.
+.NH
+AP_PLATE -- Aperture Plate Correlation
+.PP
+The final step in the extraction of a multi-aperture plate image is to
+correlate the spectra with the on-line database description of the 
+drilled hole positions.  This allows for estimates of relative wavelength
+offsets and the identification of the spectra with the ID, RA, and DEC
+parameters.
+.PP
+The spectra are fitted to the known aperture plate drilled positions, given in
+millimeters, to find an \fIangle\fR for the aperture plate relative to the
+detector x dimension and the image \fIscale\fR in pixels / millimeter,
+
+.EQ
+x sub fit = a + scale (x sub drill cos (angle) + y sub drill sin (angle))~.
+.EN
+
+If the number of spectra is less than that given by the aperture plate drilled
+positions then a correlation is done leaving out sequences of
+consecutive holes until the fit residual is minimized.  If the number of
+spectra is greater than that supposedly drilled then sequences of
+consecutive peaks are left out of the fit to minimize the residual.
+The missing holes or extra peaks are printed out and, if allowed, the aperture
+plate description file is modified, otherwise the program terminates.
+In all cases if the final fit residual is greater than 1
+pixel the program will terminate.
+The program prints out the \fIangle\fR of the aperture plate and the \fIscale\fR
+which is also stored in the database.
+.PP
+An indication that a large part of the focus variation is purely a
+scale change is that the derived image \fIscale\fR correlates very well with
+the width of the spectra as derived from the profile fitting.  I
+estimate that at least 50% of the focus variation is a scale
+variation.  This is good news in the sense that a scale variation will
+be taken out in the dispersion solution and lines in different parts
+of the detector will become more similiar after the solution.
+However, the scale variation does not account for all the profile
+shape changes and there is indeed a change in the point spread function
+across the detector.
+.NH
+Problems
+.PP
+There a few problems which I have not been able to resolve or have not
+had the time to consider.  The problems which are largely intractable
+without a great deal of effort are the unexplained background
+variations and deconvolving the spectra for the variation in the
+point-spread-function.  The background variations are abrupt increases
+in the background in parts of the CRYOCAM detector.  The step edge sometimes
+occurs under the spectra and so any smooth polynomial fit to the
+background will not be very good.  The modeling of the multi-aperture
+plate profiles provides information about the point-spread function
+but a deconvolution of the variation in the PSF is a difficult problem
+and not warrented at this time.
+.PP
+I had expected that the large scale response of the CRYOCAM could be
+corrected by determining an overall average quartz spectrum from all the
+extracted quartz spectra and then dividing the object spectra in each
+hole by the ratio of the average quartz spectra from that hole to the
+overall average quartz spectrum.  This was attempted and it was found
+to work only partially.  Specifically, while there might be a 20%
+difference between a spectrum on the edge and one near the center of
+the detector the quartz correction left a 10% difference in the object
+spectra.  This is apparently due to a poor illumination by the quartz
+light source which does not correspond to the telescope illumination.
+This quartz correction technique may be made available to users if
+desired.
+.NH
+Comparison with the Cyber Extraction Package
+.PP
+The discussion of this section must be qualified by the fact that I
+have not used the Cyber Extraction Package and I base my understanding on the
+algorithms from the Multi-Aperture Plate Data Reduction Manual and
+conversations with knowledgable people.  There are many differences
+both major and minor and this section only seeks to mention the
+some of the important differences.  In the Cyber package:
+
+The background is subtracted from the images as a preliminary process.
+
+The background is either constant or linear across the spectra.
+
+The flat fields are produced by modeling the quartz and data from
+several quartz exposures cannot be easily combined.
+
+The initial peak finding and aperture plate correlation algorithm is less
+automated in determining missing or additional holes.
+
+The model fitting uses only a three parameter Gaussian model
+and the algorithms do not yield results when the focus becomes poor.
+
+The fitting algorithm is neighbor subtraction rather than full
+simultaneous solution for all the profiles.
+
+The model fitting is applied only to a quartz and the model is transfered to
+object exposures.  This does not allow the shape of the profiles to
+change with time as the telescope moves.
+
+The modelling does not couple solutions for neighboring spectra
+across the dispersion as is suggested in this design and it does smooth
+along the spectra which is not done in this proposal.
+
+The extraction is only to some specified sigma in the model profile and
+there is no attempt to correct for blending.
+
+There is no cleaning of the spectra.
+.NH
+Discussion
+.PP
+The only data which has passed beyond the extraction phase using the
+algorithms described here was that of Paul Hintzen.
+Comparison of data reduced by the TV package for
+spectra extracted by both the Cyber package and the techniques of the
+suggested \fBmultispec\fR package were quite comparable.  To the level he
+examined
+the spectra there was no clear increase in accuracy though the \fBmultispec\fR
+extractions generally had higher counts due to the full extraction of
+the blended spectra.  The big advantages found were
+the ability to extract all the data even when the focus
+became very poor and the success of the cosmic ray cleaning
+algorithm.  Thus, Hintzen feels that the need for speed in the extraction
+(primarily dependent on the cleaning algorithm)
+is modified significantly by the difficulty of dealing with cosmic
+rays in the TV spectral analysis programs.  More exploration
+of techniques for determining the profile intensity scale from the
+model without the full multi-parameter solution is warrented for this
+reason.
+.PP
+I have extracted some Echelle data including field flattening.  The
+data had a considerable number of cosmic rays which were removed
+quite well.  The extracted spectra were put into a CAMERA format
+for further analysis.
+.PP
+The programs were recently applied to a long slit analysis problem
+being studied by Vesa Junkkarinen.  The image was already flat fielded.
+The data had two closely spaced and very faint diffuse objects and scattered
+light from a nearby QSO.
+The three spectra were so weak and closely spaced
+that the automatic finding was not used.  However, the rest of the modeling
+and extraction were applied directly.
+The \fBfind_bckgrnd\fR program, whose original purpose was to correct for
+scattered light, worked well to extrapolate the sky across the
+image.  The model fitting accurately followed
+the peaks of the spectra but the shape fitting was only moderately accurate
+since the model shape parameters are not suited to modeling galaxies.
+It successfully extracted spectra with a minimum of effort on my part.
+Analysis of the extracted spectra and comparison with other techniques
+must still be done.  The conclusions to be drawn from this experiment are
+that with sufficiently general multi-spectra tools multiple objects in
+long slit spectroscopy can be handled.
+.PP
+One area in which I do not have practical experience is
+the extraction of HGVS data.  I believe
+the proposed design will work on this type of data.
+.PP
+A point which needs to be considered in the final design are the
+formats of the data files.  The currently used one dimensional spectra
+formats are an IIDS format and a CAMERA image format.
+The formating of data files for the current spectral analysis packages by
+\fBto_iids\fR starts from the \fBmultispec\fR database and throws away a lot 
+of information about the spectra.  
+Some refinement of this database should focus on the format
+to be used by a new \fBIRAF\fR spectral analysis package.
+.PP
+It should be pointed out that many of the operations can have
+alternate algorithms substituted.  In particular, the smoothing
+algorithm for the multi-aperture plate flat fields can be replaced by
+some other scheme.  The links between the multi-parameter fitting
+program and the model have been made very general for investigating
+a broad range of models.  Thus, it is also possible to substitute
+additional model profiles with relative ease.
+.PP
+Estimates of excution time are taken from the experimental C programs
+implementing the algorithms of this design and they are only
+approximate estimates.  The steps corresponding
+to \fBdebias\fR, \fBmultispec_flat\fR, and \fBflat_divide\fR for
+the multi-aperture data from the CRYOCAM take
+about 1 hour for a typical set of frames, say 5 to 15.  This includes
+debiasing, triming, computing a flat field from several quartz frames
+and dividing the quartz into the object frames.
+.PP
+The CRYOCAM \fBmultiap_extract\fR phase takes about 40 minutes for the modeling of a frame using 32 lines per band and either 3 hours for an extraction
+using the profile fitting
+method or 14 minutes for extraction using the peak profile scaling
+method.
+.PP
+Finally, the \fBto_iids\fR takes about 3 minutes per frame.  It takes
+this long because it has to convert the \fBmultispec\fR database organized across
+the dispersion into formats in which the data is stored as consecutive
+spectra; i.e. a type of rotation operation.
diff --git a/noao/twodspec/multispec/doc/MSalgo_c.doc b/noao/twodspec/multispec/doc/MSalgo_c.doc
new file mode 100644
index 00000000..b3322dff
--- /dev/null
+++ b/noao/twodspec/multispec/doc/MSalgo_c.doc
@@ -0,0 +1,522 @@
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+         Algorithms for the Multi-Spectra Extraction Package
+                       Analysis and Discussion
+                           December 2, 1983
+
+
+
+1. Disclaimer
+
+    This  should  not  be  taken as a statement of how the algorithms of
+the final package should  function;  this  is  merely  an  analysis  and
+discussion  of  the  algorithms,  and  should  be  followed  by  further 
+discussion  before  we  decide  what  course  to  follow  in  the  final 
+package.   We  may very well decide that the level of effort required to
+implement  rigorously  correct  nonlinear  fitting  algorithms  is   not 
+justified  by  the  expected scientific usage of the package.  Before we
+can decide that, though, we need an accurate estimate of  the  level  of
+effort required.
+
+In  attacking  nonlinear  surface  fitting  problems  it is important to
+recognize that almost any techniques can  be  made  to  yield  a  result
+without  the  program crashing.  Production of a result (extraction of a
+spectrum)  does  not  mean  that  the  algorithm  converged,  that   the 
+solution   is   unique,   that  the  model  is  accurate,  or  that  the 
+uncertainties in the computed coefficients have been minimized.
+
+
+
+2. Multispec Flat (pg. 4)
+
+    This sounds like a classical high pass  filter  and  might  be  best
+implemented  via  convolution.   Using  a  convolution  operator  with a
+numerical kernel has  the  advantage  that  the  filter  can  be  easily
+modifed  by  resampling  the  kernel  or  by  changing  the  size of the
+kernel.  It is also quite efficient.   The  boundary  extension  feature
+of  IMIO  makes  it  easy  to  deal  with  the  problem  of  the  kernel 
+overlapping  the  edge  of  the  image  in  a  convolution.   Since  the 
+convolution  is  one-dimensional  (the  image is only filtered in Y), it
+will always be desirable to transpose the image.
+
+The method used to detect and reject bad pixels (eqn 1) is not  correct.
+The  rejection  criteria  should  be invariant with respect to a scaling
+of  the  pixel  values.   If  the  data  has  gone  through  very   much 
+processing  (i.e.,  dtoi  on  photographic  data),  the relation between
+photon  counts  and  pixel  value  may  be  linear,  but  the  scale  is 
+unknown.   Rejection  by  comparison  of  a  data  value to a "smoothed"
+value is more commonly done as follows:
+
+        reject if:  abs (observed - smoothed) > (K * sigma)
+
+where sigma is the noise sigma of the  data,  generally  a  function  of
+the signal.
+
+It  is  often  desirable  in rejection algorithms to be able to specify,
+
+
+                                  -1-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+as an option, that all pixels within a specified radius of a  bad  pixel
+be  rejected,  rather  than  just the pixel which was detected.  This is
+only  unnecessary  if  the  bad  pixels  are  single  pixel  events  (no 
+wings).   Region  rejection makes an iterative rejection scheme converge
+faster, as well  as  rejecting  the  faint  wings  of  the  contaminated
+region.
+
+
+
+2.1 Dividing by the Flat (pg. 5)
+
+    There  is  no  mention of any need for registering the flat with the
+data field.  Is it safe  to  assume  that  the  quartz  and  the  object
+frames  are  precisely  registered?   What  if  the  user  does  in fact
+average several quartz frames taken  over  a  period  of  time?   (Image
+registration  is  a  general  problem  that  is probably best left until
+solved in IMAGES).
+
+
+
+3. Multiap Extraction (pg. 5-6, 8-13)
+
+    The thing that bothers me most about  the  modeling  and  extraction
+process  is  that  the  high  signal  to noize quartz information is not
+used  to  full  advantage,  and  the  background  is  not  fitted   very 
+accurately.   The  present  algorithms will work well for high signal to
+noise data, but will result  in  large  (percentage)  errors  for  faint
+spectra.
+
+Basically,  it  seems to me that the high signal to noise quartz spectra
+should, in many cases, be used to determine the position  and  shape  of
+the  spectral  lines.   This  is  especially attractive since the quartz
+and spectra appear  to  be  closely  registered.   Furthermore,  if  the
+position-shape   solution   and   extraction   procedures  are  separate 
+procedures, there is nothing to prevent one from applying  both  to  the
+object  spectum  if  necessary for some reason (i.e., poor registration,
+better signal  to  noise  in  the  object  spectrum  in  the  region  of
+interest,  signal  dependent  distortions, lack of a quartz image, etc.,
+would all justify use of the object frame).  It should  be  possible  to
+model  either  the  quartz or the object frame, and to reuse a model for
+more than one extraction.
+
+Let  us  divide  the  process  up  into  two  steps,   "modeling",   and 
+"extraction"  (as  it  is  now).   The  "calibration  frame"  may be the
+quartz, an averaged quartz, or the object frame.  Ideally it  will  have
+a  high  signal  to  noise ratio and any errors in the background should
+be negligible compared to the signal.
+
+We do not solve  for  the  background  while  modeling  the  calibration
+frame;  we  assume  that  the  background  has  been  fitted by any of a
+variety  of  techniques  and  a  background  frame  written  before  the 
+calibration  frame  is  modeled.   A  "swath"  is the average of several
+image lines, where an image line  runs  across  the  dispersion,  and  a
+
+
+                                  -2-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+column along the dispersion.
+
+
+
+3.1 Modeling
+
+    I  would  set  the thing up to start fitting at any arbitrary swath,
+rather than the first swath, because it not much harder,  and  there  is
+no  guarantee  that  the  calibration frame will have adequate signal to
+noise in the first swath (indeed often the lowest signal to  noise  will
+be  found  there).   We  define the "center" swath as the first swath to
+be fitted, corresponding to the highest signal to noise  region  of  the
+calibration  frame.   By  default  the  center swath should be the swath
+used by find_spectra, especially if there is  significant  curvature  in
+the spectra.
+
+algorithm model_calibration_frame
+
+begin
+        extract center swath
+        initialize coeff using centers from find_spectra
+        model center swath (nonlinear)
+
+        for (successive swaths upward to top of frame) {
+            extract swath
+            initialize coeff to values from last fit
+            model swath (nonlinear)
+            save coeff in datafile
+        }
+
+        set last-fit coeff to values for center swath
+        for (successive swaths downward to bottom of frame) {
+            extract swath
+            initialize coeff to values from last fit
+            model swath (nonlinear)
+            save coeff in datafile
+        }
+
+        smooth model coeff (excluding intensity) along the dispersion 
+            [high freq variations in spectra center and shape from line]
+            [to line are nonphysical]
+        variance of a coeff at line-Y from the smoothed model value is
+            a measure of the uncertainty in that coeff.
+end
+
+
+I   would  have  the  background  fitting  routine  write  as  output  a 
+background frame, the name of which would  be  saved  in  the  datafile,
+rather  than  saving  the  coeff  of  the  bkg fit in the datafile.  The
+background  frame  may  then  be  produced  by  any  of  a   number   of 
+techniques;  storing  the  coefficients  of  the bkg fit in the datafile
+limits the technique used to a particular model.  For  similar  reasons,
+the  standard  bkg  fitting  routine  should  be broken up into a module
+
+
+                                  -3-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+which determines the region to be fitted, and a module  which  fits  the
+bkg pixels and writes the bkg image.
+
+For  example,  if  the  default  background fitting routine is a line by
+line  routine,  the  output  frame  could  be  smoothed  to  remove  the 
+(nonphysical)  fluctuations  in  the  background  from  line to line.  A
+true two dimensional background  fitting  routine  may  be  added  later
+without   requiring   modifications   to  the  datafile.   Second  order 
+corrections could be made to the background by  repeating  the  solution
+using the background fitted by the extraction procedure.
+
+
+procedure extract_swath
+
+begin
+        extract raw swath from calibration frame
+        extract raw swath from background frame
+        return (calib swath minus bkg swath)
+end
+
+
+The  algorithm  used to simultaneously model all spectra in a swath from
+across the dispersion is probably the most difficult and time  consuming
+part  of  the  problem.   The problem is nonlinear in all but one of the
+four or more parameters for each spectra.   You  have  spent  a  lot  of
+time  on  this  and  we  are probably not going to be able to improve on
+your algorithms significantly, though the generation of  the  matrix  in
+each step can probably be optimized significantly.
+
+The  analytic  line-profile  model  is  the most general and should work
+for  most  instruments  with  small  circular  apertures,  even  in  the 
+presence  of  severe  distortions.   It  should be possible, however, to
+fit a simpler model given by  a  lookup  table,  solving  only  for  the
+position  and  height  of  each spectra.  This model may be adequate for
+many instruments, should be must faster to fit,  and  may  produce  more
+accurate  results  since  there  are  fewer  parameters in the fit.  The
+disadvantage of an  empirical  model  is  that  it  must  be  accurately
+interpolated  (including  the  derivatives),  requiring  use  of  spline 
+interpolation or a similar technique (I have  tried  linear  and  it  is
+not  good  enough).  Vesa has implemented procedures for fitting splines
+and evaluating their derivatives.
+
+Fitting the empirical model simultaneously  to  any  number  of  spectra
+should  be  straightforward  provided the signal to noise is reasonable,
+since there are few parameters in the  fit  and  the  matrix  is  banded
+(the  Marquardt  algorithm  would work fine).  However, if you ever have
+to deal with data where a very faint or nonexistent spectra is  next  to
+a  bright  one,  it  may  be difficult to constrain the fit.  I doubt if
+the present approach of smoothing the coeff across  the  dispersion  and
+iterating  would  work  in  such  a case.  The best approach might be to
+fix the center of the faint spectra relative  to  the  bright  one  once
+the  signal  drops  below  a  certain  level, or to drop it from the fit
+entirely.  This requires that the matrix be able to change  size  during
+
+
+                                  -4-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+the fit.
+
+algorithm fit_empirical_model
+
+begin
+        [upon entry, we already have an initial estimate of the coeff]
+
+        # Marquardt (gradient expansion) algorithm.  Make 2nd order
+        # Taylor's expansion to chisquare near minimum and solve for
+        # correction vector which puts us at minimum (subject to
+        # Taylor's approx).  Taylor's approximation rapidly becomes
+        # better as we near the minimum of the multidimensional
+        # chisquare, hence convergence is extremely rapid given a good
+        # starting estimate.
+
+        repeat {
+            evaluate curvature matrix using current coeff.
+            solve banded curvature matrix 
+
+            compute error matrix
+            for (each spectra)
+                if (uncertainty in center coeff > tol) {
+                    fix center by estimation given relative spacing
+                        in higher signal region
+                    remove spectra center from solution
+                }
+
+            if (no center coeff were rejected)
+                tweak correction vector to accelerate convergence
+            new coeff vector = old coeff vector + correction vector
+            compute norm of correction vector
+        } until (no more center coeff rejected and norm < tolerance)
+
+        compute final uncertainties
+end
+
+
+The  following  is  close  to what is currently done to fit the analytic
+model, as far as  I  can  remember  (I  have  modified  it  slightly  to
+stimulate  discussion).   The  solution  is  broken up into two parts to
+reduce the number of free parameters and  increase  stability.   If  the
+uncertainty  in  a  free  parameter  becomes large it is best to fix the
+parameter (it is particularly easy for this data  to  estimate  all  but
+the  intensity  parameter).   A fixed parameter is used in the model and
+affects the solution but is not solved for (i.e., like the background).
+
+The analytic fit will  be  rather  slow,  even  if  the  outer  loop  is
+constrained  to  one  iteration.   If it takes (very rough estimates) .5
+sec to set up the banded matrix and .3 sec to  solve  it,  3  iterations
+to  convergence,  we  have  5  sec  per  swath.  If we have an 800 lines
+broken into swaths of 32 lines, the total  is  125  sec  per  image  (to
+within a factor of 5).
+
+
+
+                                  -5-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+algorithm fit_analytic_model
+
+begin
+        [upon entry, we already have an initial estimate of the coeff]
+
+        repeat {
+            save coeff
+            solve for center,height,width of each line with second
+                order terms fixed (but not necessarily zero)
+            apply constraints on line centers and widths
+            repeat solution adding second order coeff (shape terms)
+
+            compute error matrix
+            for (each coeff)
+                if (uncertainty in coeff > tol) {
+                    fix coeff value to reasonable estimate
+                    remove coeff from solution
+                }
+
+            compute total correction vector given saved coeff
+            if (no coeff were rejected)
+                tweak correction vector to accelerate convergence
+            compute norm of correction vector
+        } until (no additional coeff rejected and norm < tolerance)
+
+        compute final uncertainties
+end
+
+
+
+3.2 Extraction
+
+    The  purpose of extraction is to compute the integral of the spectra
+across the dispersion, producing I(y) for each spectra.  An estimate  of
+the  uncertainty  U(y)  should  also  be produced.  The basic extraction
+techniques  are  summarized  below.   The  number  of  spectra,  spectra 
+centers,  spectra  width  and  shape parameters are taken from the model
+fitted  to  the  calibration  frame  as  outlined  above.   We  make   a 
+simultaneous  solution  for  the  profile  heights and the background (a
+linear problem), repeating the solution independently for each  line  in
+the  image.   For  a  faint  spectrum,  it is essential to determine the
+background accurately, and we can do that safely here since  the  matrix
+will be very well behaved.
+
+    (1) Aperture  sum.   All  of the pixels within a specified radius of
+        the spectra  are  summed  to  produce  the  raw  integral.   The
+        background  image  is  also  summed  and subtracted to yield the
+        final integral.  The radius may be a constant or a  function  of
+        the  width  of  the profile.  Fractional pixel techniques should
+        be used to minimize sampling effects at the  boundaries  of  the
+        aperture.   Pixel  rejection  is  not possible since there is no
+        fitted surface.  The model is  used  only  to  get  the  spectra
+        center  and  width.   This  technique is fastest and may be best
+
+
+                                  -6-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+        if the profile is difficult to model, provided the  spectra  are
+        not crowded.
+
+    (2) Weighted  aperture  sum.   Like  (1),  except  that  a weighting
+        function is applied to correct  for  the  effects  of  crowding.
+        The  model  is  fitted  to  each  object  line,  solving  for  I 
+        (height) and B (background) with  all  other  parameters  fixed.
+        This  is  a  linear  solution  of  a banded matrix and should be
+        quite fast provided the model  can  be  sampled  efficiently  to
+        produce  the  matrix.   It  is possible to iterate to reject bad
+        pixels.  The weight for  a  spectra  at  a  data  pixel  is  the
+        fractional  contribution  of that spectra to the integral of the
+        contributions of all spectra.
+
+    (3) Fit and integrate the model.  The model is fitted as in  (2)  to
+        the   data   pixels  but  the  final  integral  is  produced  by 
+        integrating  the  model.   This   technique   should   be   more 
+        resistant  to  noise  in  the  data  than is (2), because we are
+        using the high signal to  noise  information  in  the  model  to
+        constrain  the  integral.   More  accurate  photometric  results 
+        should therefore be possible.
+
+
+Method (2) has  the  advantage  that  the  integral  is  invariant  with
+respect  to  scale  errors in the fitted models, provided the same error
+is made in each model.  Of course, the same  error  is  unlikely  to  be
+made  in  all  models contributing to a point; it is more likely that an
+error will put more energy into  one  spectra  at  the  expense  of  its
+neighbors.   In  the  limit as the spectra become less crowded, however,
+the effects of errors  in  neighboring  spectra  become  small  and  the
+weighted  average  technique looks good; it becomes quite insensitive to
+errors in the model and in the fit.  For  crowded  spectra  there  seems
+no  alternative  to a good multiparameter fit.  For faint spectra method
+(3) is probably best, and  fitting  the  background  accurately  becomes
+crucial.
+
+In  both (2) and (3), subtraction of the scaled models yields a residual
+image which can be used to evaluate at a glance the quality of the  fit.
+Since  most  all  of  the  effort in (2) and (3) is in the least squares
+solution and the pixel rejection, it might be desirable to  produce  two
+integrals  (output  spectra),  one  for  each  algorithm, as well as the
+uncertainty  vector  (computed  from  the  covariance  matrix,  not  the 
+residual).
+
+
+
+3.3 Smoothing Coefficient Arrays
+
+    In  several  places  we have seen the need for smoothing coefficient
+arrays.  The use of polynomials for  smoothing  is  questionable  unless
+the   order   of  the  polynomial  is  low  (3  or  less).   High  order 
+polynomials are  notoriously  bad  near  the  endpoints  of  the  fitted
+array,   unless  the  data  curve  happens  to  be  a  noisy  low  order 
+
+
+                                  -7-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+polynomial  (rare,  to  say  the  least).   Convolution   or   piecewise 
+polynomial  functions  (i.e., the natural cubic smoothing spline) should
+be considered if there is any reason to  believe  that  the  coefficient
+array  being  smoothed  may  have  high  frequency  components which are
+physical and must be followed (i.e., a bend or kink).
+
+
+
+3.4 Weighting (pg. 11)
+
+    The first weighting scheme (1 / sqrt (data)) seems inverted  to  me.
+The  noise  goes  up  as  with the signal, to be sure, but the signal to
+noise usually goes up faster.  Seems to me the  weight  estimate  should
+be  sqrt(data).   It  also  make  more sense to weight the least blended
+(peak) areas most.
+
+
+
+3.5 Rejection criteria (pg. 13)
+
+    The same comments apply to this rejection criterium  as  in  section
+2.   I  assume  that  "(data  -  model)"  is supposed to be "abs (data -
+model").
+
+
+
+3.6 Uncertainties and Convergence Criteria
+
+    I got the impression that you were using the residual  of  the  data
+minus  the  fitted  surface  both  as the convergence criterium and as a
+measure of the errors in the fit.  It is  neither;  assuming  a  perfect
+model, the residual gives only a measure of the noise in the data.
+
+Using   the   residual   to  establish  a  convergence  criterium  seems 
+reasonable except that it is hard to reliably  say  what  the  criterium
+should  be.   Assuming  that  the  algorithm converges, the value of the
+residual when convergence is acheived is in  general  hard  to  predict,
+so   it  seems  to  me  to  be  difficult  to  establish  a  convergence 
+criterium.  The conventional way  to  establish  when  a  nonlinear  fit
+converges  is  by measuring the norm of the correction vector.  When the
+norm becomes less than some small number the algorithm is said  to  have
+converged.   The  multidimensional  chisquare  surface is parabolic near
+the minimum and a good nonlinear algorithm will  converge  very  rapidly
+once it gets near the minimum.
+
+The  residual  is a measure of the overall goodness of fit, but tells us
+nothing about the uncertainties in the individual  coefficients  of  the
+model.    The  uncertainties  in  the  coefficients  are  given  by  the 
+covariance or error matrix (see Bevington pg. 242).  It is  ok  to  push
+forward  and  produce  an extraction if the algorithm fails to converge,
+but  ONLY  provided  the  code  gives  a  reliable   estimate   of   the 
+uncertainty.
+
+
+
+                                  -8-
+MULTISPEC (Dec83)        Multispec Algorithms          MULTISPEC (Dec83)
+
+
+
+3.6 Evaluating the Curvature Matrix Efficiently
+
+    The  most  expensive  part  of  the  reduction  process  is probably
+evaluating the model to form the curvature matrix at each  iteration  in
+the  nonlinear  solution.   The  most efficient way to do this is to use
+lookup tables.  If the profile shape does not change,  the  profile  can
+be  sampled,  fitted  with a spline, and the spline evaluated to get the
+zero through second derivatives for the curvature matrix.  This  can  be
+done  even  if  the  width  of  the  profile  changes by adding a linear
+term.  If the shape of the profile has to change, it is  still  possible
+to  sample  either  the  gaussian or the exponential function with major
+savings in computation time.
+
+
+
+3.7 Efficient Extraction (pg. 12)
+
+    The reported time of 3 cpu hours to  extract  the  spectra  from  an
+800  line  image  is  excessive for a linear solution.  I would estimate
+the time required for the 800 linear  banded  matrix  solutions  at  4-8
+minutes,  with  a  comparable  time  required  for matrix setup if it is
+done efficiently.  I suspect that the present code  is  not  setting  up
+the   linear   banded   matrix   efficiently  (not  sampling  the  model 
+efficiently).  Pixel rejection should not seriously affect  the  timings
+assuming that bad pixels are not detected in most image lines.
+
+
+
+4. Correcting for Variations in the PSF
+
+    For  all  low  signal  to  noise data it is desirable to correct for
+variations in the point  spread  function,  caused  by  variable  focus,
+scattering,  or  whatever.   This does not seem such a difficult problem
+since the width of the line profile  is  directly  correlated  with  the
+width  of  the  PSF and the information is provided by the current model
+at each point in each extracted spectrum.   The  extracted  spectra  can
+be  corrected  for the variation in the PSF by convolution with a spread
+function the width of which varies along the spectrum.
diff --git a/noao/twodspec/multispec/doc/MSalgo_c.hlp b/noao/twodspec/multispec/doc/MSalgo_c.hlp
new file mode 100644
index 00000000..4b9c3356
--- /dev/null
+++ b/noao/twodspec/multispec/doc/MSalgo_c.hlp
@@ -0,0 +1,449 @@
+.help multispec Dec83 "Multispec Algorithms"
+.ce
+Algorithms for the Multi-Spectra Extraction Package
+.ce
+Analysis and Discussion
+.ce
+December 2, 1983
+
+.sh
+1. Disclaimer
+
+    This should not be taken as a statement of how the algorithms of the
+final package should function; this is merely an analysis and discussion
+of the algorithms, and should be followed by further discussion before we
+decide what course to follow in the final package.  We may very well decide
+that the level of effort required to implement rigorously correct nonlinear
+fitting algorithms is not justified by the expected scientific usage of
+the package.  Before we can decide that, though, we need an accurate estimate
+of the level of effort required.
+
+In attacking nonlinear surface fitting problems it is important to recognize
+that almost any techniques can be made to yield a result without the program
+crashing.  Production of a result (extraction of a spectrum) does not mean
+that the algorithm converged, that the solution is unique, that the model
+is accurate, or that the uncertainties in the computed coefficients have
+been minimized.
+
+.sh
+2. Multispec Flat (pg. 4)
+
+    This sounds like a classical high pass filter and might be best implemented
+via convolution.  Using a convolution operator with a numerical kernel has
+the advantage that the filter can be easily modifed by resampling the kernel
+or by changing the size of the kernel.  It is also quite efficient.  The
+boundary extension feature of IMIO makes it easy to deal with the problem of
+the kernel overlapping the edge of the image in a convolution.  Since the
+convolution is one-dimensional (the image is only filtered in Y), it will
+always be desirable to transpose the image.
+
+The method used to detect and reject bad pixels (eqn 1) is not correct.
+The rejection criteria should be invariant with respect to a scaling of the
+pixel values.  If the data has gone through very much processing (i.e.,
+dtoi on photographic data), the relation between photon counts and pixel value
+may be linear, but the scale is unknown.  Rejection by comparison of a data
+value to a "smoothed" value is more commonly done as follows:
+
+	reject if:  abs (observed - smoothed) > (K * sigma)
+
+where sigma is the noise sigma of the data, generally a function of the signal.
+
+It is often desirable in rejection algorithms to be able to specify,
+as an option, that all pixels within a specified radius of a bad pixel
+be rejected, rather than just the pixel which was detected.  This is only
+unnecessary if the bad pixels are single pixel events (no wings).  Region
+rejection makes an iterative rejection scheme converge faster, as well as
+rejecting the faint wings of the contaminated region.
+
+.sh
+2.1 Dividing by the Flat (pg. 5)
+
+    There is no mention of any need for registering the flat with the data
+field.  Is it safe to assume that the quartz and the object frames are
+precisely registered?  What if the user does in fact average several quartz
+frames taken over a period of time?  (Image registration is a general
+problem that is probably best left until solved in IMAGES).
+
+.sh
+3. Multiap Extraction (pg. 5-6, 8-13)
+
+    The thing that bothers me most about the modeling and extraction
+process is that the high signal to noize quartz information is not used to
+full advantage, and the background is not fitted very accurately.  The
+present algorithms will work well for high signal to noise data, but 
+will result in large (percentage) errors for faint spectra.
+
+Basically, it seems to me that the high signal to noise quartz spectra
+should, in many cases, be used to determine the position and shape of the
+spectral lines.  This is especially attractive since the quartz and spectra
+appear to be closely registered.  Furthermore, if the position-shape solution
+and extraction procedures are separate procedures, there is nothing to prevent
+one from applying both to the object spectum if necessary for some reason
+(i.e., poor registration, better signal to noise in the object spectrum in
+the region of interest, signal dependent distortions, lack of a quartz image,
+etc., would all justify use of the object frame).  It should be possible to
+model either the quartz or the object frame, and to reuse a model for more
+than one extraction.
+
+Let us divide the process up into two steps, "modeling", and "extraction"
+(as it is now).  The "calibration frame" may be the quartz, an averaged
+quartz, or the object frame.  Ideally it will have a high signal to noise
+ratio and any errors in the background should be negligible compared to
+the signal.
+
+We do not solve for the background while modeling the calibration frame;
+we assume that the background has been fitted by any of a variety of
+techniques and a background frame written before the calibration frame
+is modeled.  A "swath" is the average of several image lines, where an
+image line runs across the dispersion, and a column along the dispersion.
+
+.sh
+3.1 Modeling
+
+    I would set the thing up to start fitting at any arbitrary swath, rather
+than the first swath, because it not much harder, and there is no guarantee
+that the calibration frame will have adequate signal to noise in the first
+swath (indeed often the lowest signal to noise will be found there).
+We define the "center" swath as the first swath to be fitted, corresponding
+to the highest signal to noise region of the calibration frame.  By default
+the center swath should be the swath used by find_spectra, especially if
+there is significant curvature in the spectra.
+
+.ks
+.nf
+algorithm model_calibration_frame
+
+begin
+	extract center swath
+	initialize coeff using centers from find_spectra
+	model center swath (nonlinear)
+
+	for (successive swaths upward to top of frame) {
+	    extract swath
+	    initialize coeff to values from last fit
+	    model swath (nonlinear)
+	    save coeff in datafile
+	}
+
+	set last-fit coeff to values for center swath
+	for (successive swaths downward to bottom of frame) {
+	    extract swath
+	    initialize coeff to values from last fit
+	    model swath (nonlinear)
+	    save coeff in datafile
+	}
+
+	smooth model coeff (excluding intensity) along the dispersion 
+	    [high freq variations in spectra center and shape from line]
+	    [to line are nonphysical]
+	variance of a coeff at line-Y from the smoothed model value is
+	    a measure of the uncertainty in that coeff.
+end
+.fi
+.ke
+
+
+I would have the background fitting routine write as output a background
+frame, the name of which would be saved in the datafile, rather than saving
+the coeff of the bkg fit in the datafile.  The background frame may then
+be produced by any of a number of techniques; storing the coefficients of
+the bkg fit in the datafile limits the technique used to a particular model.
+For similar reasons, the standard bkg fitting routine should be broken up
+into a module which determines the region to be fitted, and a module which
+fits the bkg pixels and writes the bkg image.
+
+For example, if the default background fitting routine is a line by line
+routine, the output frame could be smoothed to remove the (nonphysical)
+fluctuations in the background from line to line.  A true two dimensional
+background fitting routine may be added later without requiring modifications
+to the datafile.  Second order corrections could be made to the background
+by repeating the solution using the background fitted by the extraction
+procedure.
+
+
+.ks
+.nf
+procedure extract_swath
+
+begin
+	extract raw swath from calibration frame
+	extract raw swath from background frame
+	return (calib swath minus bkg swath)
+end
+.fi
+.ke
+
+
+The algorithm used to simultaneously model all spectra in a swath from
+across the dispersion is probably the most difficult and time consuming
+part of the problem.  The problem is nonlinear in all but one of the four
+or more parameters for each spectra.  You have spent a lot of time on this
+and we are probably not going to be able to improve on your algorithms
+significantly, though the generation of the matrix in each step can
+probably be optimized significantly.
+
+The analytic line-profile model is the most general and should work for most
+instruments with small circular apertures, even in the presence of severe
+distortions.  It should be possible, however, to fit a simpler model given
+by a lookup table, solving only for the position and height of each spectra.
+This model may be adequate for many instruments, should be must faster to
+fit, and may produce more accurate results since there are fewer parameters
+in the fit.  The disadvantage of an empirical model is that it must be
+accurately interpolated (including the derivatives), requiring use of spline
+interpolation or a similar technique (I have tried linear and it is not good
+enough).  Vesa has implemented procedures for fitting splines and evaluating
+their derivatives.
+
+Fitting the empirical model simultaneously to any number of spectra should
+be straightforward provided the signal to noise is reasonable, since there
+are few parameters in the fit and the matrix is banded (the Marquardt
+algorithm would work fine).  However, if you ever have to deal with data
+where a very faint or nonexistent spectra is next to a bright one, it may
+be difficult to constrain the fit.  I doubt if the present approach of
+smoothing the coeff across the dispersion and iterating would work in such
+a case.  The best approach might be to fix the center of the faint spectra
+relative to the bright one once the signal drops below a certain level,
+or to drop it from the fit entirely.  This requires that the matrix be able
+to change size during the fit.
+
+.ks
+.nf
+algorithm fit_empirical_model
+
+begin
+	[upon entry, we already have an initial estimate of the coeff]
+
+	# Marquardt (gradient expansion) algorithm.  Make 2nd order
+	# Taylor's expansion to chisquare near minimum and solve for
+	# correction vector which puts us at minimum (subject to
+	# Taylor's approx).  Taylor's approximation rapidly becomes
+	# better as we near the minimum of the multidimensional
+	# chisquare, hence convergence is extremely rapid given a good
+	# starting estimate.
+
+	repeat {
+	    evaluate curvature matrix using current coeff.
+	    solve banded curvature matrix 
+
+	    compute error matrix
+	    for (each spectra)
+		if (uncertainty in center coeff > tol) {
+		    fix center by estimation given relative spacing
+			in higher signal region
+		    remove spectra center from solution
+		}
+
+	    if (no center coeff were rejected)
+		tweak correction vector to accelerate convergence
+	    new coeff vector = old coeff vector + correction vector
+	    compute norm of correction vector
+	} until (no more center coeff rejected and norm < tolerance)
+
+	compute final uncertainties
+end
+.fi
+.ke
+
+
+The following is close to what is currently done to fit the analytic
+model, as far as I can remember (I have modified it slightly to stimulate
+discussion).  The solution is broken up into two parts to reduce the number
+of free parameters and increase stability.  If the uncertainty in a free
+parameter becomes large it is best to fix the parameter (it is particularly
+easy for this data to estimate all but the intensity parameter).  A fixed
+parameter is used in the model and affects the solution but is not solved
+for (i.e., like the background).
+
+The analytic fit will be rather slow, even if the outer loop is constrained
+to one iteration.  If it takes (very rough estimates) .5 sec to set up the
+banded matrix and .3 sec to solve it, 3 iterations to convergence, we have
+5 sec per swath.  If we have an 800 lines broken into swaths of 32 lines,
+the total is 125 sec per image (to within a factor of 5). 
+
+
+.ks
+.nf
+algorithm fit_analytic_model
+
+begin
+	[upon entry, we already have an initial estimate of the coeff]
+
+	repeat {
+	    save coeff
+	    solve for center,height,width of each line with second
+		order terms fixed (but not necessarily zero)
+	    apply constraints on line centers and widths
+	    repeat solution adding second order coeff (shape terms)
+
+	    compute error matrix
+	    for (each coeff)
+		if (uncertainty in coeff > tol) {
+		    fix coeff value to reasonable estimate
+		    remove coeff from solution
+		}
+
+	    compute total correction vector given saved coeff
+	    if (no coeff were rejected)
+		tweak correction vector to accelerate convergence
+	    compute norm of correction vector
+	} until (no additional coeff rejected and norm < tolerance)
+
+	compute final uncertainties
+end
+.fi
+.ke
+
+.sh
+3.2 Extraction
+
+    The purpose of extraction is to compute the integral of the spectra
+across the dispersion, producing I(y) for each spectra.  An estimate of
+the uncertainty U(y) should also be produced.  The basic extraction techniques
+are summarized below.  The number of spectra, spectra centers, spectra width
+and shape parameters are taken from the model fitted to the calibration
+frame as outlined above.  We make a simultaneous solution for the profile
+heights and the background (a linear problem), repeating the solution
+independently for each line in the image.  For a faint spectrum, it is
+essential to determine the background accurately, and we can do that safely
+here since the matrix will be very well behaved.
+.ls 4
+.ls (1)
+Aperture sum.  All of the pixels within a specified radius of the spectra
+are summed to produce the raw integral.  The background image is also summed
+and subtracted to yield the final integral.  The radius may be a constant or a
+function of the width of the profile.  Fractional pixel techniques should
+be used to minimize sampling effects at the boundaries of the aperture.
+Pixel rejection is not possible since there is no fitted surface.  The model
+is used only to get the spectra center and width.  This technique is fastest
+and may be best if the profile is difficult to model, provided the spectra
+are not crowded.
+.le
+.ls (2)
+Weighted aperture sum.  Like (1), except that a weighting function is
+applied to correct for the effects of crowding.  The model is fitted
+to each object line, solving for I (height) and B (background) with all
+other parameters fixed.  This is a linear solution of a banded matrix and
+should be quite fast provided the model can be sampled efficiently to
+produce the matrix.  It is possible to iterate to reject bad pixels.
+The weight for a spectra at a data pixel is the fractional contribution
+of that spectra to the integral of the contributions of all spectra.
+.le
+.ls (3)
+Fit and integrate the model.  The model is fitted as in (2) to the data
+pixels but the final integral is produced by integrating the model.
+This technique should be more resistant to noise in the data than is (2),
+because we are using the high signal to noise information in the model to
+constrain the integral.  More accurate photometric results should therefore
+be possible.
+.le
+.le
+
+
+Method (2) has the advantage that the integral is invariant with respect
+to scale errors in the fitted models, provided the same error is made in
+each model.  Of course, the same error is unlikely to be made in all
+models contributing to a point; it is more likely that an error will put
+more energy into one spectra at the expense of its neighbors.  In the limit
+as the spectra become less crowded, however, the effects of errors in
+neighboring spectra become small and the weighted average technique looks
+good; it becomes quite insensitive to errors in the model and in the fit.
+For crowded spectra there seems no alternative to a good multiparameter
+fit.  For faint spectra method (3) is probably best, and fitting the
+background accurately becomes crucial.
+
+In both (2) and (3), subtraction of the scaled models yields a residual
+image which can be used to evaluate at a glance the quality of the fit.
+Since most all of the effort in (2) and (3) is in the least squares solution
+and the pixel rejection, it might be desirable to produce two integrals
+(output spectra), one for each algorithm, as well as the uncertainty vector
+(computed from the covariance matrix, not the residual).
+
+.sh
+3.3 Smoothing Coefficient Arrays
+
+    In several places we have seen the need for smoothing coefficient arrays.
+The use of polynomials for smoothing is questionable unless the order of
+the polynomial is low (3 or less).  High order polynomials are notoriously
+bad near the endpoints of the fitted array, unless the data curve happens
+to be a noisy low order polynomial (rare, to say the least).  Convolution or
+piecewise polynomial functions (i.e., the natural cubic smoothing spline)
+should be considered if there is any reason to believe that the coefficient
+array being smoothed may have high frequency components which are physical and
+must be followed (i.e., a bend or kink).
+
+.sh
+3.4 Weighting (pg. 11)
+
+    The first weighting scheme (1 / sqrt (data)) seems inverted to me.
+The noise goes up as with the signal, to be sure, but the signal to noise
+usually goes up faster.  Seems to me the weight estimate should be sqrt(data).
+It also make more sense to weight the least blended (peak) areas most.
+
+.sh
+3.5 Rejection criteria (pg. 13)
+
+    The same comments apply to this rejection criterium as in section 2.
+I assume that "(data - model)" is supposed to be "abs (data - model").
+
+.sh
+3.6 Uncertainties and Convergence Criteria
+
+    I got the impression that you were using the residual of the data minus
+the fitted surface both as the convergence criterium and as a measure of the
+errors in the fit.  It is neither; assuming a perfect model, the residual gives
+only a measure of the noise in the data.
+
+Using the residual to establish a convergence criterium seems reasonable
+except that it is hard to reliably say what the criterium should be.
+Assuming that the algorithm converges, the value of the residual when
+convergence is achieved is in general hard to predict, so it seems to me to
+be difficult to establish a convergence criterium.  The conventional way
+to establish when a nonlinear fit converges is by measuring the norm of
+the correction vector.  When the norm becomes less than some small number
+the algorithm is said to have converged.  The multidimensional chisquare
+surface is parabolic near the minimum and a good nonlinear algorithm will
+converge very rapidly once it gets near the minimum.
+
+The residual is a measure of the overall goodness of fit, but tells us
+nothing about the uncertainties in the individual coefficients of the model.
+The uncertainties in the coefficients are given by the covariance or error
+matrix (see Bevington pg. 242).  It is ok to push forward and produce an
+extraction if the algorithm fails to converge, but ONLY provided the code 
+gives a reliable estimate of the uncertainty.
+
+.sh
+3.6 Evaluating the Curvature Matrix Efficiently
+
+    The most expensive part of the reduction process is probably evaluating
+the model to form the curvature matrix at each iteration in the nonlinear
+solution.  The most efficient way to do this is to use lookup tables.
+If the profile shape does not change, the profile can be sampled, fitted
+with a spline, and the spline evaluated to get the zero through second
+derivatives for the curvature matrix.  This can be done even if the width
+of the profile changes by adding a linear term.  If the shape of the profile
+has to change, it is still possible to sample either the gaussian or the
+exponential function with major savings in computation time.
+
+.sh
+3.7 Efficient Extraction (pg. 12)
+
+    The reported time of 3 cpu hours to extract the spectra from an 800 line
+image is excessive for a linear solution.  I would estimate the time required
+for the 800 linear banded matrix solutions at 4-8 minutes, with a comparable
+time required for matrix setup if it is done efficiently.  I suspect that the
+present code is not setting up the linear banded matrix efficiently (not
+sampling the model efficiently).  Pixel rejection should not seriously affect
+the timings assuming that bad pixels are not detected in most image lines.
+
+.sh
+4. Correcting for Variations in the PSF
+
+    For all low signal to noise data it is desirable to correct for variations
+in the point spread function, caused by variable focus, scattering, or
+whatever.  This does not seem such a difficult problem since the width of
+the line profile is directly correlated with the width of the PSF and the
+information is provided by the current model at each point in each extracted
+spectrum.  The extracted spectra can be corrected for the variation in the
+PSF by convolution with a spread function the width of which varies along
+the spectrum.
+.endhelp
diff --git a/noao/twodspec/multispec/doc/MSspecs.doc b/noao/twodspec/multispec/doc/MSspecs.doc
new file mode 100644
index 00000000..09955e9c
--- /dev/null
+++ b/noao/twodspec/multispec/doc/MSspecs.doc
@@ -0,0 +1,698 @@
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+   Detailed Specifications for the Multi-Spectra Extraction Package
+                              F. Valdes
+                           December 8, 1983
+
+
+1.  Introduction
+
+     The  multi-spectra  extraction  package  (MULTISPEC)  provides  the 
+basic tools for modeling, cleaning, and extracting spectra  from  images
+containing  multiple  aperture  spectra running roughly parallel.  These
+tools will generally be combined  in  reduction  script  tasks  but  may
+also be used directly for non-standard analysis.
+
+     This  design  presents  the requirements and specifications for the
+MULTISPEC package.  Details concerning the algorithms  are  given  in  a
+separate document, Algorithms for the Multi-Spectra Extraction Package.
+
+
+2.  Input Data Requirements
+
+    The  input  data  for  the MULTISPEC package consists of image files
+containing one or more aperture spectra.  The spectra  are  required  to
+run   roughly  parallel  to  each  other  and  parallel  to  the  second 
+digitization  axis.   The  latter  requirement  may  require  a  general 
+rotation  and  interpolation  image operator.  The images are assumed to
+be  corrected  to  linear  relative  intensity.   Thus,  the  steps   of 
+correcting   digital   detector  images  for  dark  current,  bias,  and 
+pixel-to-pixel sensitivity variations must  be  performed  before  using
+the MULTISPEC tasks.
+
+     Because  the  the MULTISPEC package is being developed concurrently
+with the IRAF standard image processing tools  this  document  specifies
+the   requirements  for  the  preliminary  image  processing  needed  to 
+prepare digital detector images for the MULTISPEC package.
+
+
+2.1  Basic Digital Detector Reduction Tasks
+
+     The prelimary reduction of multi-spectra  images  uses  CL  scripts
+based  on  general  image  operators.   Some  of  the  scripts  are  for 
+specific instruments or specific reduction  applications  and  some  are
+generally   useful  image  processing  tasks.   The  scripts  allow  the 
+specification  of  many  images  for  which  the  operations   will   be 
+repetitively applied.
+
+     The  following  CL  scripts  are  required  to reduce multi-spectra
+images from digital detectors.
+
+
+        debias          multispec_flat          flat_divide
+
+
+
+
+
+                                  -1-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+debias
+    The files in a list of  filenames  are  automatically  debiased  and
+    trimmed.   This  routine  will  be  instrument  specific but used by
+    other reduction tasks beyond MULTISPEC.
+
+multispec_flat
+    The files in a list of quartz  multi-spectra  filenames  are  added,
+    the  result  is  smoothed  along  the dispersion dimension, and then
+    the original image is divided by the smoothed  image  to  produce  a
+    flat  field  image.   The  unsmoothed  to smoothed ratio is computed
+    only  if  the  value  of  the  smoothed  pixel  is  greater  than  a 
+    specified  amount.   Otherwise,  the  ratio  is  set to unity.  This
+    routine is not instrument specific but is used  only  for  MULTISPEC
+    reductions.
+
+flat_divide
+    The  files  in  a  list of filenames are each divided by a specified
+    flat field image.  This routine is  not  instrument  or  application
+    specific.
+
+    The  required  general image processing programs needed to implement
+these scripts are specified below.
+
+
+(1) A routine to compute the average value from a specified area of  the
+    image.  Used to determine the average bias value from a bias strip.
+
+(2) A  routine  to  trim  a specified portion of an image.  Used to trim
+    the bias strip.
+
+(3) Routines to multiply and subtract images by  a  constant.   Used  to
+    scale  images  such as dark exposures and to remove the average bias
+    value as obtained by (1) above.
+
+(4) Routines to subtract, add, and  divide  images.   Used  to  subtract
+    dark  current  and  bias  exposures,  to  add  several  exposures to
+    increase the signal-to-noise, and to divide by a flat  field  image.
+    The  divide  routine  must  give the user the option to substitute a
+    constant or ignore any divisions in which the  denominator  is  less
+    than a specified value.
+
+(5) A  routine  to  rotate  or  transpose  an  image.  Used to align the
+    spectra along lines or columns.
+
+(6) A  routine  to  apply  a  filter  to  lines  of  the   image.    For 
+    multi-spectra  images  a  smooth  quartz  is  produced  by  using  a 
+    running  quadratic  filter  along  each  line  of   the   dispersion 
+    dimension.   The  filter  must  be  able  to  recognize  bad  pixels 
+    (specified by a user defined threshold) and  remove  them  from  the
+    filtering operation.
+
+
+
+
+
+                                  -2-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+3.  Requirements for the MULTSPEC Package
+
+     The MULTISPEC package shall satisfy the following requirements.
+
+(1) The component programs shall be CL callable.
+
+(2) The  programs  shall interact only through image files and MULTISPEC
+    data files.
+
+(3) It shall be possible to extract spectra without modeling.
+
+(4) The entire image shall be extracted and not limited by  failures  in
+    the algorithms.
+
+(5) It  shall  be  possible  to  specify specific lines or swaths in the
+    image on which to operate.
+
+(6) CL scripts shall be provided for the  common  data  sources.   These
+    scripts will work automatically.
+
+The follow functions shall be provided:
+
+o   Make  an  initial  rough but automated identification of the spectra
+    locations.
+
+o   Provide for a user identification list for  the  spectra  locations.
+    This  list  shall  be  of  the  standard  image cursor type to allow
+    generation of the list with the standard image cursor programs.
+
+o   Determine and correct for a slowly varying background.
+
+o   Reliably and accurately trace spectra in the presence  of  geometric
+    distortions (pincushion, s, shear, etc.).
+
+o   Extract spectra by one of:
+
+    a.  Strips  of  constant width about the located spectra.  The width
+        may be specified to fractions of  a  pixel  and  the  extraction
+        will use fractional pixel interpolation.  l
+
+    b.  Strips of width proportional to a Gaussian width parameter.
+
+    c.  Modeling  to  obtain  estimates  of  the  total luminosity.  The
+        estimate will be the integral of the model.
+
+    d.  Summation   of   the   data   pixel   values   with   fractional  
+        contributions  of  the  pixel  value  to  the  spectra  based on
+        modeling.
+
+o   An option shall be available to  specify  whether  to  ignore  blank
+    pixels or use interpolated values.
+
+    o   Programs  shall  be  provided to produce data files which can be
+
+
+                                  -3-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+        accessed by one dimensional  spectroscopic  reduction  routines.
+        At a minimum these formats shall include:
+
+        a.  Reduction  to  an  image  file  consisting  of  one line per
+            extracted spectrum
+
+        b.  The  standard  IIDS  format   available   with   the   CYBER 
+            Multi-Aperture Plate programs
+
+
+3.2  Modeling Requirements
+
+     The  modeling  of multi-spectra images, particularly in the case of
+blended spectra, shall:
+
+(1) Model blended spectra with  sufficient  reliability  and  robustness
+    that  a  reasonable  solution is always obtained, though of possibly
+    limited usefulness.
+
+(2) The modeling shall provide estimates for the  uncertainties  in  the
+    fitted parameters as a function of position along the spectrum.
+
+(3) Remove  cosmic  rays  and other defective pixels by reference to the
+    model.
+
+(4) Allow the transfer of a model solution  for  one  image  to  another
+    image.
+
+(5) Display  numerically and graphically the data, the fitted model, and
+    the residuals.
+
+
+4.  Program Specifications
+
+
+4.1  Basic Programs
+
+     The basic programs of the package are general purpose  tools  which
+initialize  a  MULTISPEC  data  file  and  perform  a single fundamental
+operation on the data in the MULTISPEC data file.   There  is  one  data
+file  associated  with  each  image.   The  data file is hidden from the
+user and so the user need not be aware  of  the  data  file.   The  data
+files  are  referenced  only the image filename specified in the program
+parameters.  The data files contain such  information  as  a  processing
+history,  the  spectra  positions  and extracted luminosities, the model
+parameters (one set for each spectra for each modelled  image  line  (or
+swath),  etc.   The  programs  generally are allowed to specify specific
+lines, columns, and/or spectra on which to operate.   The  line,  column
+and  spectra  specifications  are given as strings which contain numbers
+separated by whitespace, commas,  and  the  range  indicator  "-".   The
+script  tasks  of  section  4.2  will  combine  these  basic programs to
+perform a general multi-spectra extraction.
+
+
+
+                                  -4-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+    ap_plate            copy_params     find_spectra    convolve
+    fit_bckgrnd         find_bckgrnd    line_list       model_extrac
+    model_fit           model_image     model_list      sigma_extract
+    strip_extract       to_iids         to_image        to_onedspec
+
+ap_plate
+    The information from an on-line data  file  containing  descriptions
+    of  all  the aperture plates prepared at Kitt Peak is read to find a
+    specified  aperture  plate.   The  drilled  aperture  positions  are 
+    correlated  with  the  spectra  in  the  image  to  deduce  relative 
+    wavelength offsets.  The identifications for  the  spectra  as  well
+    as  other  auxiliary  information  is recorded in the data file.  If
+    no  image  file  is  specified  then   only   the   aperture   plate 
+    information   is   printed.    This   program   is   used   in   the  
+    MULTIAP_EXTRACT program.  This  program  is  not  essential  to  the
+    operation of the MULTISPEC package.
+
+                    Multi-Spectra image  image = 
+                         Aperture plate  plate = 
+                                        (mode = ql)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+                                  -5-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+The Background
+    The  are  two possibilities for dealing with the background.  In the
+    first  case,  FIT_BCKGRND,  the  background  will   be   fitted   by 
+    polynomials  and  the  coefficients  stored  in  the  MULTISPEC data
+    file.  These coefficients are then used by  the  other  programs  to
+    estimate   the  background  at  the  spectra.   The  second  option, 
+    FIND_BCKGRND, generates a background  image  in  which  the  spectra
+    and  other  selected  areas are set to blank pixels.  Then a general
+    image interpolator is used fill in the blank pixels with  background
+    estimates.   The  other  MULTISPEC  programs  will  then access this
+    background frame.  The background frame image name  will  be  stored
+    in the MULTISPEC data file and the image header.
+
+
+    fit_bckgrnd
+        Fit  a  background  in  a  MULTISPEC image by a polynomial using
+        pixels not near the spectra and in  the  user  specified  swaths
+        and  columns.   The buffer distance is in pixels and refers to a
+        minimum distance from the center of any  spectrum  beyond  which
+        the  background  pixels  are found.  Blank pixels are ignored in
+        the background fit.  Deviant pixels will be rejected.
+
+                        Multi-Spectra image  image = 
+                        Buffer from spectra  buffer = 12
+                           Polynomial order  order = 3
+                            Lines per swath (lines_per_swath = 32)
+                              Swaths to fit (swaths = 1-1000)
+                             Columns to fit (columns = 1-1000)
+                        Rejection threshold (threshold = 5)
+                  Print general diagnostics (verbose = no)
+                                            (mode = ql)
+
+    find_bckgrnd
+        The spectra within a buffer distance  and  specified  areas  are
+        set  to  blank  pixels  and  the  remaining  pixels  copied to a
+        background image file.
+
+                        Multi-Spectra image  image = 
+                           Background image  background =
+                        Buffer from spectra  buffer = 12
+                            Lines to ignore (lines = )
+                          Columns to ignore (columns = )
+                                            (mode = ql)
+
+convolve
+    A program will be provided to reduce either the  extracted  spectrum
+    or the modeled image to a common point-spread function.
+
+
+
+
+
+
+
+
+                                  -6-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+copy_params
+    Create  a  MULTISPEC  data  file  for  a new image using appropriate
+    MULTISPEC parameters from an old image.  The  old  image  must  have
+    been  processed  to find the spectra using FIND_SPECTRA and possibly
+    model fit.
+
+                Old Multi-Spectra image  old_image = 
+                New Multi-Spectra image  new_image = 
+                                        (mode = ql)
+
+find_spectra
+    Initially locate the spectra in a MULTISPEC  image.   The  positions
+    of  the  spectra  within the range of columns are determined for the
+    starting line and then the spectra are tracked within the  range  of
+    lines.   The  minimum  separation  and minimum width would generally
+    be set for a particular instrument.   If  the  automatic  search  is
+    not  used  then a list of cursor positions is read from the standard
+    input.
+
+                    Multi-Spectra image  image = 
+                       Automatic search  auto = yes
+                          Starting line  start_line =
+                     Minimum separation (min_sep = 1)
+                          Minimum width (min_width = 1)
+                        Averaging width (average = 32)
+                        Lines to search (lines = 1-1000)
+                      Columns to search (columns = 1-1000)
+              Print general diagnostics (verbose = no)
+                                        (mode = ql)
+
+line_list
+    For the  specified  lines  in  the  image  print  the  image  column
+    number,  data  value  (possibly as a swath average), the model value
+    at that point (i.e. the sum of  the  model  contributions  from  all
+    the  spectra),  the  background  value,  and the residual.  Plotting
+    scripts may be written using this routine to show the quality  of  a
+    model fit.
+
+                    Multi-Spectra image  image = 
+                          Lines to list (lines = 1-1000)
+                                        (mode = ql)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+                                  -7-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+model_extract
+    A  previously  fitted  model  is  used  to extract the spectra total
+    luminosity by apportioning the data values to spectra in  the  ratio
+    indicated  by  the model.  If the clean option is specified then the
+    model is used to detect pixels which deviate from  the  model  by  a
+    specified  amount.   The  model  value replaces the deviant pixel in
+    the extraction and, if specified, also in the image file.
+
+                    Multi-Spectra image  image = 
+                       Lines to extract (lines = 1-1000)
+                          Clean spectra (clean = yes)
+                     Cleaning threshold (threshold = 5)
+                           Modify image (modify = yes)
+              Print general diagnostics (verbose = no)
+                                        (mode = ql)
+
+model_fit
+    A specified model is iteratively fitted to the data in each  of  the
+    specified  lines  (or  swaths)  until  the  RMS  residual  fails  to 
+    decrease.  The models  are  selected  by  a  string.   The  possible
+    values are
+
+        (null string) - initialize the model
+        i - fit only the intensity scale
+        ip - fit the intensity scale and the position
+        ips1 - fit the intensity scale, position, and one parameter shape
+        ips2 - fit the intensity scale, position, and two parameter shape
+        ips3 - fit the intensity scale, position, and three parameter shape
+        ips4 - fit the intensity scale, position, and four parameter shape
+    These  models  will  be  combined in a script to search for the best
+    fit.
+
+    The initial shape parameters will generally be set by scripts for  a
+    particular data reduction.
+
+                    Multi-Spectra image  image = 
+                             Model type  model = 
+                        Lines per swath (lines_per_swath = 32)
+                        Swaths to model (swaths = 1-1000)
+                         Initial shape1 (shape1 = .1 )
+                         Initial shape2 (shape2 = 0 )
+                         Initial shape3 (shape3 = 0 )
+                         Initial shape4 (shape4 = 5 )
+              Print general diagnostics (verbose = no)
+                                        (mode = ql)
+
+
+
+
+
+
+
+
+
+
+                                  -8-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+model_image
+    An  image  file of the fitted model is created.  This image may then
+    be displayed or a residual image may be calculated and displayed.
+
+                    Multi-Spectra image  image = 
+                            Model image  model =
+                                        (mode = ql)
+    .nf
+    .le
+    .ls model_list
+    For the specified lines and spectra the model is listed.
+    The listing gives, for each spectra,
+    the spectrum number, the line number, the fitted position,
+    the estimated wavelength, the
+    extracted luminosity, the intensity scale, model width parameters, and
+    the background polynomial coefficients.  This routine can be used in scripts
+    to plot the extracted spectra, the trend of width with wavelength, and so
+    forth.
+    
+    .nf
+                    Multi-Spectra image  image = 
+                          Lines to list (lines = 1-1000)
+                        Spectra to list (spectra = 1-1000)
+                                        (mode = ql)
+
+sigma_extract
+    A previously fitted model is used to extract the spectra  luminosity
+    within  a  specified  sigma  of  the peak.  Because the model is not
+    necessarily a Gaussian the sigma is used to  compute  the  intensity
+    ratio  of  the  cutoff to the peak assumining a Gaussian profile and
+    then the data is extracted to the point the  model  intensity  falls
+    below  that  cutoff.   If  the  clean  option  is specified then the
+    model is used to detect pixels which deviate from  the  model  by  a
+    specified  amount.   The  model  value replaces the deviant pixel in
+    the extraction and, if specified, also in the image file.
+
+                    Multi-Spectra image  image = 
+                 Sigma extraction width  width = 1.
+                       Lines to extract (lines = 1-1000)
+                          Clean spectra (clean = yes)
+                     Cleaning threshold (threshold = 5)
+                           Modify image (modify = yes)
+              Print general diagnostics (verbose = no)
+                                        (mode = ql)
+
+
+
+
+
+
+
+
+
+
+
+                                  -9-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+strip_extract
+    A strip of constant width about the spectra positions is  extracted.
+    If  cleanning  is  desired  a  smoothed  estimate  of the profile is
+    obtained by averaging a  number  of  lines  about  the  line  to  be
+    cleaned.   After  fitting  for  the intensity scale pixels are found
+    which deviate from the profile by a specified amount.   The  profile
+    value   replaces  the  deviant  pixel  in  the  extraction  and,  if 
+    specified, also in the image file.  No prior  modeling  is  required
+    to use this extraction routine.
+
+                    Multi-Spectra image  image = 
+                 Strip extraction width  width = 1.
+                       Lines to extract (lines = 1-1000)
+                          Clean spectra (clean = yes)
+                     Cleaning threshold (threshold = 5)
+              Lines per profile average (averge_lines = 32)
+                           Modify image (modify = yes)
+              Print general diagnostics (verbose = no)
+                                        (mode = ql)
+
+to_iids
+    For  a  specified prefix, files of the form prefix.nn, where nn is a
+    specified spectra  number,  are  created  containing  the  extracted
+    spectra  for  all  the  specified  image  files.   The format of the
+    files is the IIDS format  developed  for  the  CYBER  Multi-Aperture
+    Plate Extractions.
+
+                    Multi-Spectra image  images = 
+                   IIDS filename prefix  iids_file = 
+                      Spectra to format (spectra = 1-1000)
+                                        (mode = ql)
+
+to_image
+    An  image  file  containing  one  line of the extracted luminosities
+    for each specified spectra in the specified MULTISPEC image.
+
+                    Multi-Spectra image  in_image = 
+                Extracted spectra image  out_image = 
+                                Spectra (spectra = 1-1000)
+                                        (mode = ql)
+
+to_onedspec
+    The extractions are converted to an as yet to  be  specified  format
+    for use in the ONEDSPEC reduction package.
+
+                   Multi-Spectra images  images = 
+                     ONEDSPEC data file  onedspec_file = 
+                                Spectra (spectra = 1-1000)
+                                        (mode = ql)
+
+
+
+
+
+
+                                 -10-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+4.2  General MULTISPEC CL Scripts
+
+     The  general  MULTISPEC CL scripts perform a series of steps needed
+to extract the spectra from a specified  list  of  image  files.   These
+steps  have  been found to generally perform the desired extraction task
+fully.
+
+
+        multiap_extract         echelle_extract
+
+multiap_extract
+    The specified multi-aperture plate  images  are  extracted.   If  no
+    starting  solution  image,  one which has previously been extracted,
+    is specified then the script performs an automatic  search  for  the
+    specified  number  of  spectra.   Otherwise  the  solution  from the
+    starting image is used as the  initial  model.   The  background  is
+    then  determined.   This is followed by a series of fitting steps on
+    swaths of data.  (For further details on the fitting steps  see  the
+    Algorithms   paper).    A   MODEL_EXTRACT   and   cleaning  follows. 
+    Finally, the extraction is correlated with  the  specified  aperture
+    plate  using  AP_PLATE.   If  there  was no starting image then this
+    extraction becomes the initial solution  image.   Subsequent  images
+    are extracted starting from the initial solution image.
+
+                  Multi-Aperture images  images = 
+                 Initial solution image  initial = 
+                  Aperture plate number  plate = 
+                      Number of spectra  nspectra = 
+                                        (mode = ql)
+
+echelle_extract
+    The   specified  echelle  images  are  extracted.   If  no  starting 
+    solution  image,  one  which  has  previously  been  extracted,   is 
+    specified  then  the  script  performs  an  automatic search for the
+    specified  number  of  orders.   Otherwise  the  solution  from  the 
+    starting   image  is  used  as  the  initial  starting  point.   The 
+    background  is  then  determined.   Finally  a   STRIP_EXTRACT   and 
+    cleaning  is  performed.   If  there was no starting image then this
+    extraction becomes the initial solution  image.   Subsequent  images
+    are extracted starting from the initial solution image.
+
+                         Echelle images  images = 
+                 Initial solution image  initial = 
+                       Number of orders  norders = 
+                       Extraction width  width =
+                                        (mode = ql)
+
+
+5.  Outline of a MULTISPEC Reduction
+
+     The  following  outline  is for the reduction of a cryogenic camera
+multi-aperture plate.  All the programmer supplied  default  values  are
+used.
+
+
+                                 -11-
+MULTISPEC (Oct83)  Multi-Spectra Extraction Package    MULTISPEC (Oct83)
+
+
+
+        (1)  rcamera mtb, "ap165.", "s", "3-9"
+        (2)  debias "ap165.*"
+        (3)  multispec_flat "ap165.[36]", "ap165.flat"
+        (4)  flat_divide "ap165.*", "ap165.flat"
+        (5)  multiap_extract "ap165.*", "", 165, 50
+        (6)  to_onedspec "ap165.*", oned165
+
+
+(1) The  data  is  read  from  the observing tape(s) using RCAMERA.  The
+    image files created are ap165.3, ap165.4,  ...,  ap165.9.   This  is
+    easily  accomplished  by  using  the filename prefix "ap165." in the
+    RCAMERA program.  The raw images may be examined at  this  point  on
+    a display.
+
+(2) The  images  are  debiased  using DEBIAS with all the "ap165." files
+    specified.  The debias program  knows  about  the  location  of  the
+    bias strip for the cryogenic camera.
+
+(3) A  a  flat  field  is  created  using  MULTISPEC_FLAT  in  which the
+    desired  quartz  frames  are  specified  and  a  flat  field   image 
+    filename  is  defined.  The created flat field image may be examined
+    on an image display if desired.
+
+(4) All the  debiased  images  are  divided  by  the  flat  field  using
+    FLAT_DIVIDE.
+
+(5) The  script  MULTIAP_EXTRACT  is  run  in  which  the aperture plate
+    number, the number of spectra, and the image files to  be  extracted
+    are  specified.   The  number  of  spectra  is found by examining an
+    image on an image display or by plotting a cut  across  the  spectra
+    using a general image profile program.
+
+(6) Finally,  the  extracted  spectra  are  formatted  for  the ONEDSPEC
+    package using TO_ONEDSPEC with the extracted images specified.
diff --git a/noao/twodspec/multispec/doc/MSspecs.hlp b/noao/twodspec/multispec/doc/MSspecs.hlp
new file mode 100644
index 00000000..92f285ed
--- /dev/null
+++ b/noao/twodspec/multispec/doc/MSspecs.hlp
@@ -0,0 +1,659 @@
+.help multispec Oct83 "Multi-Spectra Extraction Package"
+.sp 3
+.ce
+Detailed Specifications for the Multi-Spectra Extraction Package
+.ce
+F. Valdes
+.ce
+December 8, 1983
+.sh
+1.  Introduction
+
+     The multi-spectra extraction package (MULTISPEC) provides the basic tools
+for modeling, cleaning, and extracting spectra from images
+containing multiple aperture spectra running roughly parallel.
+These tools will generally be combined in reduction script tasks
+but may also be used directly for non-standard analysis.
+
+     This design presents the requirements and specifications
+for the MULTISPEC package.  Details concerning the 
+algorithms are given in a separate document, Algorithms for the
+Multi-Spectra Extraction Package.
+.sh
+2.  Input Data Requirements
+
+    The input data for the MULTISPEC package consists of image
+files containing one or more aperture spectra.  The spectra are
+required to run roughly parallel to each other and parallel to the
+second digitization axis.  The latter requirement may require a
+general rotation and interpolation image operator.  The images are
+assumed to be corrected to linear relative intensity.  Thus, the
+steps of correcting digital detector images for dark current, bias, and
+pixel-to-pixel sensitivity variations must be performed before using
+the MULTISPEC tasks.
+
+     Because the MULTISPEC package is being developed
+concurrently with the IRAF standard image processing
+tools this document specifies the requirements for the preliminary
+image processing needed to prepare digital detector images for the MULTISPEC
+package.
+.sh
+2.1  Basic Digital Detector Reduction Tasks
+
+     The prelimary reduction of multi-spectra images uses CL scripts
+based on general image operators.
+Some of the scripts are for specific instruments or specific
+reduction applications and some are generally useful image processing
+tasks.  The scripts allow the specification of many images for which
+the operations will be repetitively applied.
+
+     The following CL scripts are required to reduce multi-spectra images
+from digital detectors.
+
+
+.nf
+	debias		multispec_flat		flat_divide
+.fi
+.ke
+.ks
+.ls 4 debias
+The files in a list of filenames are automatically debiased and trimmed.
+This routine will be instrument specific but used by other reduction
+tasks beyond MULTISPEC.
+.le
+.ke
+.ks
+.ls multispec_flat
+The files in a list of quartz multi-spectra filenames are added,
+the result is smoothed
+along the dispersion dimension, and then the original image is divided
+by the smoothed image to produce a flat field image.  The unsmoothed
+to smoothed ratio is computed only if the value of the smoothed
+pixel is greater than a specified amount.  Otherwise, the ratio is set
+to unity.  This routine is not instrument specific but is used only
+for MULTISPEC reductions.
+.le
+.ke
+.ks
+.ls flat_divide
+The files in a list of filenames are each divided by a specified flat
+field image.  This routine is not instrument or application specific.
+.le
+.ke
+
+    The required general image processing programs needed to implement these
+scripts are specified below.
+
+.ls (1)
+A routine to compute the average value from a specified area of the
+image.  Used to determine the average bias value from a bias strip.
+.le
+.ls (2)
+A routine to trim a specified portion of an image.  Used to trim the
+bias strip.
+.le
+.ls (3)
+Routines to multiply and subtract images by a constant.  Used to scale
+images such as dark exposures and to remove the average bias value as
+obtained by (1) above.
+.le
+.ls (4)
+Routines to subtract, add, and divide images.  Used to subtract dark
+current and bias exposures, to add several exposures to increase the
+signal-to-noise, and to divide by a flat field image.
+The divide routine must give the user the option to substitute a constant or
+ignore any divisions in which the denominator is less than a specified value.
+.le
+.ls (5)
+A routine to rotate or transpose an image.  Used to align the spectra
+along lines or columns.
+.le
+.ls (6)
+A routine to apply a filter to lines of the image.  For multi-spectra images
+a smooth quartz is produced by using a running quadratic filter along each
+line of the dispersion dimension.  The filter must be able to recognize
+bad pixels (specified by a user defined threshold) and remove them from the
+filtering operation.
+.le
+.sh
+3.  Requirements for the MULTSPEC Package
+
+     The MULTISPEC package shall satisfy the following requirements.
+.ls (1)
+The component programs shall be CL callable.
+.le
+.ls (2)
+The programs shall interact only through image files and MULTISPEC data files.
+.le
+.ls (3)
+It shall be possible to extract spectra without modeling.
+.le
+.ls (4)
+The entire image shall be extracted and not limited by failures in the
+algorithms.
+.le
+.ls (5)
+It shall be possible to specify specific lines or swaths in the image
+on which to operate.
+.le
+.ls (6)
+CL scripts shall be provided for the common data sources.  These scripts
+will work automatically.
+.le
+
+The follow functions shall be provided:
+.ls o
+Make an initial rough but automated identification of the spectra
+locations.
+.le
+.ls o
+Provide for a user identification list for the spectra locations.
+This list shall be of the standard image cursor type to allow generation
+of the list with the standard image cursor programs.
+.le
+.ls o
+Determine and correct for a slowly varying background.
+.le
+.ls o
+Reliably and accurately trace spectra in the presence of geometric
+distortions (pincushion, s, shear, etc.).
+.le
+
+Extract spectra by one of:
+.ls a.
+Strips of constant width about the located spectra.  The width may be specified
+to fractions of a pixel and the extraction will use fractional pixel
+interpolation.
+l
+.le
+.ls b.
+Strips of width proportional to a Gaussian width parameter.
+.le
+.ls c.
+Modeling to obtain estimates of the total luminosity.  The estimate will
+be the integral of the model.
+.le
+.ls d.
+Summation of the data pixel values with fractional contributions of the
+pixel value to the spectra based on modeling.
+.le
+.le
+.ls o
+An option shall be available to specify whether to ignore blank pixels
+or use interpolated values.
+.ls o
+Programs shall be provided to produce data files which can be accessed
+by one dimensional spectroscopic reduction routines.  At a minimum
+these formats shall include:
+.ls a.
+Reduction to an image file consisting of one line per extracted
+spectrum
+.le
+.ls b.
+The standard IIDS format available with the CYBER Multi-Aperture Plate
+programs
+.le
+.le
+.sh
+3.2  Modeling Requirements
+
+     The modeling of multi-spectra images, particularly in the case of
+blended spectra, shall:
+.ls (1)
+Model blended spectra with sufficient reliability and robustness that
+a reasonable solution is always obtained, though of possibly limited
+usefulness.
+.le
+.ls (2)
+The modeling shall provide estimates for the uncertainties in the fitted
+parameters as a function of position along the spectrum.
+.le
+.ls (3)
+Remove cosmic rays and other defective pixels by reference to the model.
+.le
+.ls (4)
+Allow the transfer of a model solution for one image to another image.
+.le
+.ls (5)
+Display numerically and graphically the data, the fitted model, and
+the residuals.
+.le
+.sh
+4.  Program Specifications
+.sh
+4.1  Basic Programs
+
+     The basic programs of the package are general purpose tools which
+initialize a MULTISPEC data file and perform a single fundamental operation
+on the data in the MULTISPEC data file.  There is one data file associated
+with each image.  The data file is hidden from the user and so the user
+need not be aware of the data file.
+The data files are referenced only the image filename specified in the
+program parameters.
+The data files contain such information as a processing history, the
+spectra positions and extracted luminosities, the model parameters (one
+set for each spectra for each modelled image line (or swath), etc.
+The programs generally are allowed to specify specific
+lines, columns, and/or spectra on which to operate.
+The line, column and spectra specifications are given as strings which
+contain numbers separated by whitespace, commas, and the range indicator
+"-".  The script tasks
+of section 4.2 will combine these basic programs to perform a general
+multi-spectra extraction.
+
+.ks
+.nf
+    ap_plate		copy_params	find_spectra	convolve
+    fit_bckgrnd		find_bckgrnd	line_list	model_extrac
+    model_fit		model_image	model_list	sigma_extract
+    strip_extract	to_iids		to_image	to_onedspec
+.fi
+.ke
+.ks
+.ls ap_plate
+The information from an on-line data file containing descriptions of all
+the aperture plates prepared at Kitt Peak is read to find a specified
+aperture plate.  The drilled aperture positions are correlated with the
+spectra in the image to deduce relative wavelength offsets.  The
+identifications for the spectra as well as other auxiliary information
+is recorded in the data file.
+If no image file is specified then only the aperture
+plate information is printed.  This program is used in the MULTIAP_EXTRACT
+program.  This program is not essential to the operation of the MULTISPEC
+package.
+
+.nf
+                Multi-Spectra image  image = 
+                     Aperture plate  plate = 
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls The Background
+The are two possibilities for dealing with the background.  In the first
+case, FIT_BCKGRND, the background will be fitted by polynomials and
+the coefficients stored in the MULTISPEC data file.  These coefficients
+are then used by the other programs to estimate the background at the
+spectra.  The second option, FIND_BCKGRND, generates a background image in which
+the spectra and other selected areas are set to blank pixels.  Then a
+general image interpolator is used fill in the blank pixels with background
+estimates.  The other MULTISPEC programs will then access this background
+frame.  The background frame image name will be stored in the MULTISPEC
+data file and the image header.
+
+.ls fit_bckgrnd
+Fit a background in a MULTISPEC image by a polynomial using pixels
+not near the spectra and in the user specified swaths and columns.
+The buffer distance is in pixels and refers to a minimum distance from
+the center of any spectrum beyond which the background pixels are found.
+Blank pixels are ignored in the background fit.  Deviant pixels will be
+rejected.
+
+.nf
+                Multi-Spectra image  image = 
+                Buffer from spectra  buffer = 12
+                   Polynomial order  order = 3
+                    Lines per swath (lines_per_swath = 32)
+                      Swaths to fit (swaths = 1-1000)
+                     Columns to fit (columns = 1-1000)
+		Rejection threshold (threshold = 5)
+	  Print general diagnostics (verbose = no)
+                                    (mode = ql)
+.fi
+.le
+.ls find_bckgrnd
+The spectra within a buffer distance and specified areas are set to blank
+pixels and the remaining pixels copied to a background image file.
+
+.nf
+                Multi-Spectra image  image = 
+		   Background image  background =
+                Buffer from spectra  buffer = 12
+                    Lines to ignore (lines = )
+                  Columns to ignore (columns = )
+                                    (mode = ql)
+.fi
+.le
+.le
+.ke
+.ks
+.ls convolve
+A program will be provided to reduce either the extracted spectrum or
+the modeled image to a common point-spread function.
+.le
+.ke
+.ks
+.ls copy_params
+Create a MULTISPEC data file for a new image using
+appropriate MULTISPEC parameters from an old image.
+The old image must have been processed to find the spectra using FIND_SPECTRA
+and possibly model fit.
+
+.nf
+            Old Multi-Spectra image  old_image = 
+            New Multi-Spectra image  new_image = 
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls find_spectra
+Initially locate the spectra in a MULTISPEC image.
+The positions of the spectra within the range of columns are determined
+for the starting line and then the spectra are tracked within the
+range of lines.  The minimum separation
+and minimum width would generally be set for a particular instrument.
+If the automatic search is not used then a list of cursor positions is
+read from the standard input.
+
+.nf
+                Multi-Spectra image  image = 
+		   Automatic search  auto = yes
+		      Starting line  start_line =
+		 Minimum separation (min_sep = 1)
+		      Minimum width (min_width = 1)
+                    Averaging width (average = 32)
+                    Lines to search (lines = 1-1000)
+                  Columns to search (columns = 1-1000)
+	  Print general diagnostics (verbose = no)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls line_list
+For the specified lines in the image print the image column
+number, data value (possibly as a swath average), the model value at that
+point (i.e. the sum of the model contributions from all the spectra),
+the background value, and the residual.
+Plotting scripts may be written using this routine to
+show the quality of a model fit.
+
+.nf
+                Multi-Spectra image  image = 
+                      Lines to list (lines = 1-1000)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls model_extract
+A previously fitted model is used to extract the spectra total luminosity
+by apportioning the data values to spectra in the ratio indicated by the
+model.  If the clean option is specified then the model is used to detect
+pixels which deviate from the model by a specified amount.
+The model value replaces the deviant pixel in the extraction and, if specified,
+also in the image file.
+
+.nf
+                Multi-Spectra image  image = 
+                   Lines to extract (lines = 1-1000)
+                      Clean spectra (clean = yes)
+                 Cleaning threshold (threshold = 5)
+                       Modify image (modify = yes)
+	  Print general diagnostics (verbose = no)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls model_fit
+A specified model is iteratively fitted to the data in each of the specified
+lines (or swaths) until the RMS residual fails to decrease.  The models
+are selected by a string.  The possible values are
+
+.nf
+    (null string) - initialize the model
+    i - fit only the intensity scale
+    ip - fit the intensity scale and the position
+    ips1 - fit the intensity scale, position, and one parameter shape
+    ips2 - fit the intensity scale, position, and two parameter shape
+    ips3 - fit the intensity scale, position, and three parameter shape
+    ips4 - fit the intensity scale, position, and four parameter shape
+.fi
+These models will be combined in a script to search for the best fit.
+
+The initial shape parameters will generally be set by scripts for a
+particular data reduction.
+
+.nf
+                Multi-Spectra image  image = 
+                         Model type  model = 
+                    Lines per swath (lines_per_swath = 32)
+                    Swaths to model (swaths = 1-1000)
+		     Initial shape1 (shape1 = .1 )
+		     Initial shape2 (shape2 = 0 )
+		     Initial shape3 (shape3 = 0 )
+		     Initial shape4 (shape4 = 5 )
+	  Print general diagnostics (verbose = no)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls model_image
+An image file of the fitted model is created.  This image may then be displayed
+or a residual image may be calculated and displayed.
+
+.nf
+                Multi-Spectra image  image = 
+		        Model image  model =
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls model_list
+For the specified lines and spectra the model is listed.
+The listing gives, for each spectra,
+the spectrum number, the line number, the fitted position,
+the estimated wavelength, the
+extracted luminosity, the intensity scale, model width parameters, and
+the background polynomial coefficients.  This routine can be used in scripts
+to plot the extracted spectra, the trend of width with wavelength, and so
+forth.
+
+.nf
+                Multi-Spectra image  image = 
+                      Lines to list (lines = 1-1000)
+                    Spectra to list (spectra = 1-1000)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls sigma_extract
+A previously fitted model is used to extract the spectra luminosity
+within a specified sigma of the peak.  Because the model is not necessarily
+a Gaussian the sigma is used to compute
+the intensity ratio of the cutoff to the peak assuming a Gaussian profile
+and then the data is extracted to the point the model intensity falls below that
+cutoff.  If the clean option is specified then the model is used to detect
+pixels which deviate from the model by a specified amount.
+The model value replaces the deviant pixel in the extraction and, if specified,
+also in the image file.
+
+.nf
+                Multi-Spectra image  image = 
+             Sigma extraction width  width = 1.
+                   Lines to extract (lines = 1-1000)
+                      Clean spectra (clean = yes)
+                 Cleaning threshold (threshold = 5)
+                       Modify image (modify = yes)
+	  Print general diagnostics (verbose = no)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls strip_extract
+A strip of constant width about the spectra positions is extracted.
+If cleanning is desired a smoothed estimate of the profile is obtained
+by averaging a number of lines about the line to be cleaned.  After fitting
+for the intensity scale pixels are found which deviate from the profile by
+a specified amount.
+The profile value replaces the deviant pixel in the extraction and,
+if specified, also in the image file.  No prior modeling is required
+to use this extraction routine.
+
+.nf
+                Multi-Spectra image  image = 
+             Strip extraction width  width = 1.
+                   Lines to extract (lines = 1-1000)
+                      Clean spectra (clean = yes)
+                 Cleaning threshold (threshold = 5)
+          Lines per profile average (averge_lines = 32)
+                       Modify image (modify = yes)
+	  Print general diagnostics (verbose = no)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls to_iids
+For a specified prefix, files of the form prefix.nn, where nn is a specified
+spectra number, are created containing the extracted spectra for all
+the specified image files.  The format of the files is the IIDS format
+developed for the CYBER Multi-Aperture Plate Extractions.
+
+.nf
+                Multi-Spectra image  images = 
+               IIDS filename prefix  iids_file = 
+                  Spectra to format (spectra = 1-1000)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls to_image
+An image file containing one line of the extracted luminosities for each
+specified spectra in the specified MULTISPEC image.
+
+.nf
+                Multi-Spectra image  in_image = 
+            Extracted spectra image  out_image = 
+                            Spectra (spectra = 1-1000)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls to_onedspec
+The extractions are converted to an as yet to be specified format for
+use in the ONEDSPEC reduction package.
+
+.nf
+               Multi-Spectra images  images = 
+                 ONEDSPEC data file  onedspec_file = 
+                            Spectra (spectra = 1-1000)
+                                    (mode = ql)
+.fi
+.le
+.ke
+.sh
+4.2  General MULTISPEC CL Scripts
+
+     The general MULTISPEC CL scripts perform a series of steps needed to
+extract the spectra from a specified list of image files.  These steps have
+been found to generally perform the desired extraction task fully.
+
+
+.nf
+	multiap_extract		echelle_extract
+.fi
+.ks
+.ls multiap_extract
+The specified multi-aperture plate images are extracted.
+If no starting solution image, one which has previously been extracted,
+is specified then the script performs an automatic search for the
+specified number of spectra.
+Otherwise the solution from the starting image is used as the initial
+model.  The background is then determined.
+This is followed by a series of fitting steps on swaths of data.
+(For further details on the fitting steps see the Algorithms paper).
+A MODEL_EXTRACT and cleaning follows.
+Finally, the extraction is correlated with the specified aperture plate
+using AP_PLATE.
+If there was no starting image then this extraction becomes the
+initial solution image.
+Subsequent images are extracted starting from the initial solution image.
+
+.nf
+              Multi-Aperture images  images = 
+             Initial solution image  initial = 
+              Aperture plate number  plate = 
+                  Number of spectra  nspectra = 
+                                    (mode = ql)
+.fi
+.le
+.ke
+.ks
+.ls echelle_extract
+The specified echelle images are extracted.
+If no starting solution image, one which has previously been extracted,
+is specified then the script performs an automatic search for the
+specified number of orders.
+Otherwise the solution from the starting image is used as the initial
+starting point.  The background is then determined.
+Finally a STRIP_EXTRACT and cleaning is performed.
+If there was no starting image then this extraction becomes the
+initial solution image.
+Subsequent images are extracted starting from the initial solution image.
+
+.nf
+                     Echelle images  images = 
+             Initial solution image  initial = 
+                   Number of orders  norders = 
+		   Extraction width  width =
+                                    (mode = ql)
+.fi
+.le
+.sh
+5.  Outline of a MULTISPEC Reduction
+
+     The following outline is for the reduction of a cryogenic camera
+multi-aperture plate.  All the programmer supplied default values are
+used.
+
+.nf
+	(1)  rcamera mtb, "ap165.", "s", "3-9"
+	(2)  debias "ap165.*"
+	(3)  multispec_flat "ap165.[36]", "ap165.flat"
+	(4)  flat_divide "ap165.*", "ap165.flat"
+	(5)  multiap_extract "ap165.*", "", 165, 50
+	(6)  to_onedspec "ap165.*", oned165
+.fi
+
+.ls (1)
+The data is read from the observing tape(s) using RCAMERA.
+The image files created are ap165.3, ap165.4, ..., ap165.9.  This is
+easily accomplished by using the filename prefix "ap165." in the RCAMERA
+program.  The raw images may be examined at this point on a display.
+.le
+.ls (2)
+The images are debiased using DEBIAS with all the "ap165." files specified.
+The debias program knows about the location of the bias strip for the
+cryogenic camera.
+.le
+.ls (3)
+A a flat field is created
+using MULTISPEC_FLAT in which the desired quartz frames are specified
+and a flat field image filename is defined.  The created flat field
+image may be examined on an image display if desired.
+.le
+.ls (4)
+All the debiased images are divided by the flat field using FLAT_DIVIDE.
+.le
+.ls (5)
+The script MULTIAP_EXTRACT is run in which the aperture plate number,
+the number of spectra, and the image files to be extracted are specified.
+The number of spectra is found by examining an image on an image display
+or by plotting a cut across the spectra using a general image profile
+program.
+.le
+.ls (6)
+Finally, the extracted spectra are formatted for the ONEDSPEC package
+using TO_ONEDSPEC with the extracted images specified.
+.le
+.endhelp
diff --git a/noao/twodspec/multispec/doc/MSspecs_c.hlp b/noao/twodspec/multispec/doc/MSspecs_c.hlp
new file mode 100644
index 00000000..848d589d
--- /dev/null
+++ b/noao/twodspec/multispec/doc/MSspecs_c.hlp
@@ -0,0 +1,243 @@
+
+.help multispec Nov82 "Multispec Specifications"
+.ce
+Comments on Multispec Package Specifications
+.ce
+November 8, 1983
+
+
+
+    The basic package structure and the decomposition of the package into
+tasks looks good.  The requirements for both general operators and canned
+procedures are addressed well.  I got the impression that you have a pretty
+clear idea of what you want to do (which is the thing I am most looking for
+when I read a specs document), but I confess to having to reread the document
+several times to figure out what you have in mind.  Your writing style is
+very terse and leaves much up to the reader!
+
+Most of my comments have to do with details.  These are presented in the
+order in which they occurred while reading the document.  These comments
+apply only to the specs document.  I have started going over the algorithms
+paper, mostly when I could not understand a section of the specs document,
+but I have not finished it yet.
+
+.sh
+General Comments
+.ls 4
+.ls (1)
+When eventually we write the user documentation, the nomenclature
+should be carefully explained up front.  Users will tend to confuse
+image lines and spectral lines, but there is little we can do about
+that other than to make the distinction clear.  The term "band" is
+confusing because it normally refers to the third dimension of an
+image and that is not how it is used here.  A better term might be
+"swath".  In what follows I will continue to use the term band, but
+it is definitely not too late to change.
+.le
+.ls (2)
+It seems to me that the concept of a band or swath is a detail of how
+the algorithm works and should not have such a prominent place in the
+user interface to the package.  Several of the routines require that
+image coordinates be entered in units of band number and column.
+This introduces an unnecessary coupling between two input parameters
+and forces the user to convert from line number to band number.  The
+result will be that the user will be reluctant to change the number
+of lines per band (I'll bet that you have kept this a constant in
+using the prototype).  My inclination would be to have the user enter
+all coordinates in units of lines and columns, and have the program
+select the nearest band depending on the band width parameter.
+The band width could then be easily changed depending on the data,
+without need to respecify the region of the image to be processed.
+.le
+.ls (3)
+Routines all over the system will have an option for printing extra
+information, i.e., a verbose mode of execution.  I think we should
+standardize on the name of this parameter.  "Verbose" seems to me
+more descriptive than "print", and is consistent with UNIX terminology.
+.le
+.le
+
+.sh
+Pages 3,4
+.ls
+.ls (1)
+Functions for extracting spectra.  I assume "strips of constant
+width" means aperture sum out to a specified x-radius from the
+center of a spectra.  Can the radius be specified in fractional
+pixels, and if so, does the routine do fractional pixel interpolation.
+What happens if there are blank pixels in the aperture?
+
+If extraction is based on the model, I gather that you are still
+summing data pixel values, using a weight for each spectra based
+on the modeled contribution of each spectra to the data pixel.  In
+other words we are still taking an aperture sum, but with allowances
+for crowding.  This has the disadvantage that if we sum way out into
+the wings, we will be adding noise to the aperture sum, degrading signal
+to noise.
+
+Extraction based on integration of the model rather than
+the data should be available as another extraction procedure; this may
+yield better photometric results.  I would eventually like to compare
+the two approaches with artificial data.  Also by integrating the model
+there is no need to "clean" (I assume that deviant pixels are detected
+and rejected when the model is fitted, or the model will not be
+accurate).  Blank pixels should be recognized and ignored when fitting
+the model.
+.le
+
+.ls (2)
+I gather that all extracted spectra for an image are put into a single
+imagefile.  This is fine, even desirable, as long as it is ok if all
+spectra share the same header, and as long as all we want to output
+is intensity versus wavelength.  If it is desired to also output the
+signal to noise or whatever than another scheme may be needed.
+.le
+.ls (3)
+The text file output form ('c'pg.4) should be engineered with the idea
+that the user will take the data away in cardimage form.  From the
+description it sounds like there is one pixel (wavelength bin) per
+line in the text file.  This has its advantages, but is not what one
+wants for a cardimage file, which always writes 80 chars per line.
+Also, the detailed technical specs should give some details about
+such a format; it is a major part of the user interface and people
+will want to know what this format is going to look like.  In a way
+it is more important to specify formats like this than the calling
+sequences of the tasks, because it is harder to change after the
+package is released, and other program are written to read the
+text format spectra.
+.le
+.ls (4)
+To item 3.2 (2) (on uncertainty estimates) I would add "as a function
+of position along the spectrum".
+.le
+.le
+
+.sh
+4.1 Basic Programs
+.ls
+.ls (1)
+Evidently there is a datafile associated with each image.  What is
+the function of the datafile?  Is it transparent to the user?  How
+much is stored in the image header and how much in the datafile?
+.le
+.ls (2)
+The distinction between "line_list" and "model_list" is confusing.
+Does "line_list" print the sum of the models for all the spectra
+a each column?  Please specify the form of the output for this
+procedure in more detail.  The line_list and model_list procedures
+are natural candidates for use with the "lists" utilities for
+extracting columns, plotting one column against another, etc. I
+could not tell whether or not this would work well from the info
+given.
+.le
+.ls (3)
+"ap_plate":  "The identifications for the spectra ... is recorded."
+Is recorded where?  In the datafile?  Is this information essential
+to the operation of multispec, or is it merely passed on through
+multispec?
+.le
+.ls (4)
+"find_background":  Might be more aptly named "fit_background".
+I would expect "find" to mean find which regions of the image are
+background and which are spectra.  Find is spatial, fit is grayscale.
+
+We need to decide whether we want to specify polynomials in IRAF by
+the order (0,1,2, etc.) or by the number of coefficients or terms.
+It seems to me that people are most used to talking about second,
+third, fifth etc. order polynomials and that we might better specify
+polynomials with an "order" parameter rather than a "terms" param.
+
+Buffer radius or diameter?  I would assume radius, but it is not
+clear from the docs.  What is being "searched"?  Shouldn't that read
+"bands to be fitted".  The "colummns" parameter should permit a list
+of ranges of columns; I couldn't tell whether this was the case
+from the specs.  Cursor input may be desirable here.
+
+Blank pixels should be detected and ignored when fitting the
+background.  Are deviant pixels detected and rejected?  This is
+generally a desirable option in a bkg fit.  You may be able to
+decompose this routine (internally) into a find_background and
+a fit_background, making use of the Images background fitting
+routines, though these generate an image as output rather than the
+coeff of the fitted functions.  I wuld guess that you are storing
+the bkg coeff for each band in the datafile from the description,
+and that the fit is strictly one-dimensional.
+
+If only a limited number of bands are fitted, what do you do about
+the other bands if the bkg fit is one-dimensional?  Is the user
+req'd to use the same bands range when they do the extraction?
+.le
+
+.ls (5)
+"find_spectra".  It is not clear how this routine uses cursor input.
+Perhaps you should have a gcur type parameter.  Reading cursor
+coordinates from the standard input may be the way to go, but you
+should explain how this is going to work.
+.le
+.ls (6)
+"line_list".  One output line per image line?  One or more spectra
+per output line?  Output should be suitable for further processing
+with the LISTS package utilities (i.e., getcol, and the graphics
+utility which will plot or overplot lists).  The specs should
+specify the form of the output.
+.le
+.ls (7)
+I assume that the extraction procedures extract spectra which
+are put somewhere.  Where, in the datafile?  If the image is
+to be cleaned, it would be safer to write a new output image,
+or at least to rename the original.  It is strange to have these
+two quite different functions in the same module.
+.le
+.ls (8)
+"model_fit".  The range of modeling options is impressive, good
+stuff.  However, there must be something better than magic integer
+numbers for specifying the model to be fitted.  Perhaps the
+strings  "i, ip, ipw, ipw2, ipw3, ipw4", where 'i' is for intensity,
+'p' for position, and 'w' for width.
+
+How are the "initial parameters" specified?
+.le
+.ls (9)
+"model_list".  Again, I can only guess from the description what the
+output will look like.  It sounds like it might be best to have
+this routine print data for only one spectra at a time, particularly
+if the lists package is to be used for analysis.  It might be good
+to have the line number in the output somewhere, especially if the
+wavelength information is not available.
+.le
+.le
+
+.sh
+4.2 Scripts
+.ls
+.ls (1)
+It sounds like there is no easy alternative to an automatic search
+for the line centers.  This is best as long as it works, but the
+users will want easy way to use the cursor available as an option.
+A script such as this can easily use the line plot routine Images
+to make a plot and generate a list of line centers, without even
+requiring find_spectra to be able to access the cursor (and perhaps
+it should not if the script can do it).  The graphics cursor should
+be used here rather than the image cursor.
+.le
+.le
+
+.sh
+5. Example
+.ls
+.ls (1)
+The rcamera example is in error.  Rcamera, as implemented, has only
+three query mode params, while you show four in the example.
+I believe the ranges string should be quoted and should be the second
+argument.
+
+The last command should be "to_onedspec", not "onedspec".
+.le
+.ls (2)
+5.(5): It seems strange to make the user manually count 50 spectra
+by examining a plot.  If the program automatically finds centers,
+this should not be necessary; if the user interactively marks centers,
+it is not necessary.
+.le
+.le
+.endhelp
diff --git a/noao/twodspec/multispec/doc/findpeaks.hlp b/noao/twodspec/multispec/doc/findpeaks.hlp
new file mode 100644
index 00000000..f6118281
--- /dev/null
+++ b/noao/twodspec/multispec/doc/findpeaks.hlp
@@ -0,0 +1,88 @@
+.help findpeaks Jul84 noao.twodspec.multispec
+.ih
+NAME
+findpeaks -- Find peaks in a multi-spectra image
+.ih
+USAGE
+findpeaks image lines contrast
+.ih
+PARAMETERS
+.ls image
+Image to be searched.
+.le
+.ls lines
+Sample image lines in which the peaks are to be found.
+.le
+.ls contrast
+Maximum contrast between the highest peak and the lowest peak.
+.le
+.ls separation = 5
+Minimum separation in pixels between acceptable peaks.
+.le
+.ls edge = 0
+Minimum distance in pixels to the edge of the image for acceptable peaks.
+.le
+.ls threshold = 0.
+The minimum acceptable peak pixel value.
+.le
+.ls min_npeaks = 1
+Minimum number of peaks to be found.  It is an error for fewer than
+this number of peaks to be found.
+.le
+.ls max_npeaks = 1000
+Maximum number of peaks to be found.  If more than this number of peaks
+is found then only the those with the highest peak values are accepted.
+.le
+.ls columns = '*'
+Columns to be searched.
+.le
+.ls naverage = 20
+Number of image lines around the sample line to be averaged before
+finding the peaks.
+.le
+.ls debug = no
+Print detailed information on the progress of the peak finding algorithm.
+.le
+.ih
+DESCRIPTION
+For each specified sample image line the number of peaks and their column
+positions in the image are determined.
+The number of peaks and their positions are assumed to correspond to points
+along the spectra.  This information is entered in the MULTISPEC database.
+
+The \fInaverage\fR image lines about the specified sample line are first
+averaged.  The local maxima in the average line are then located
+in the specified columns more than the minimum distance from the edge of the
+image.  A minimum peak pixel value cutoff is determined as the maximum of
+the specified \fIthreshold\fR and \fIcontrast\fR times the largest peak pixel
+value.  All local maxima with pixel values below the cutoff are rejected.
+Next all peaks with separations less than \fIseparation\fR from a stronger
+peak are rejected.  Finally, if there are more than \fImax_npeaks\fR remaining
+only the \fImax_npeaks\fR strongest peaks are accepted.  If fewer
+than \fImin_npeaks\fR are found then the task quits with an error.
+
+If the number of spectra has been previously determined, such as by an earlier
+use of \fBfindpeaks\fR, then it is an error if a different number of
+peaks is found.
+.ih
+EXAMPLES
+The parameters of this task provide a great deal of flexibility in
+automatically determining the number and positions of the peaks.  
+The most automatic method just uses the contrast to limit the acceptable
+peaks:
+
+	cl> findpeaks image.db 1 .1
+
+However, if the number of spectra in the image is known:
+
+	cl> findpeaks image.db 1 0 min=10 max=10
+
+or if a threshold is known:
+
+	cl> findpeaks image.db 1 0 threshold = 1000
+
+For a noisy image the separation parameter can be set to eliminate spurious
+noise peaks near the peaks to be found:
+
+	cl> findpeaks image.db 1 .1 sep=20
+.endhelp
diff --git a/noao/twodspec/multispec/doc/fitfunc.hlp b/noao/twodspec/multispec/doc/fitfunc.hlp
new file mode 100644
index 00000000..09510bb7
--- /dev/null
+++ b/noao/twodspec/multispec/doc/fitfunc.hlp
@@ -0,0 +1,73 @@
+.help fitfunction Jul84 noao.twodspec.multispec
+.ih
+NAME
+fitfunction -- Fit a function to the spectra parameter values
+.ih
+USAGE
+fitfunction image
+.ih
+PARAMETERS
+.ls image
+Image in which the parameter values are to be fitted.
+.le
+.ls parameter = "x0"
+Parameter to be fit.  The legal minimum match abbreviated parameters
+are x0, s0, s1, s2.
+.le
+.ls lines = "*"
+Sample image lines to be used in the function fit.
+.le
+.ls spectra = "*"
+Spectra for which the parameters are to be fit.
+.le
+.ls function = "interpolation spline"
+Fitting function to be used.  The function is specified as a string
+which may be minimum match abbreviated.  The functions currently available
+are:
+.ls interpolation spline
+Interpolation spline of specified order.
+.le
+.ls smoothing spline
+Smoothing spline of specified order and number of polynomial pieces.
+.le
+.le
+.ls spline_order = 4
+Order of the fitting spline.  The order must be even.
+The minimum value is 2 and maximum value is determined from the number of
+sample lines in the fit.
+.le
+.ls spline_pieces = 1
+The number of polynomial pieces in a smoothing spline.
+The minimum value is 1 and the maximum value is determined from the number of
+sample lines in the fit.
+.le
+.ih
+DESCRIPTION
+A function is fit to the parameter values previously determined at the sample
+lines for each spectrum.  The function coefficients are stored in the
+database and the fitted values replace the original values at all the sample
+lines (not just the sample lines used in the fit).  The type of function,
+the parameter to be fitted, the sample lines used in the fit, and the
+spectra to be fitted are all selected by the user.  The function is
+extrapolated to cover all image lines.
+
+The values of the function fit at arbitrary image lines may be listed
+with \fBmslist\fR.
+.ih
+EXAMPLES
+The extraction of the spectra requires that a fitting function be
+determined for the spectra positions.  This is done by:
+
+	cl> fitfunction image
+
+To smooth the parameter "s0" in model \fIgauss5\fR with a cubic spline
+and leave out a bad point at sample line 7:
+
+.nf
+	cl> fitfunction image parmeter=s0 function=smooth \
+	>>> lines="1-6,8-"
+.fi
+.ih
+SEE ALSO
+mslist
+.endhelp
diff --git a/noao/twodspec/multispec/doc/fitgauss5.hlp b/noao/twodspec/multispec/doc/fitgauss5.hlp
new file mode 100644
index 00000000..bcb37276
--- /dev/null
+++ b/noao/twodspec/multispec/doc/fitgauss5.hlp
@@ -0,0 +1,148 @@
+.help fitgauss5 Jul84 noao.twodspec.multispec
+.ih
+NAME
+fitgauss5 -- Fit spectra profiles with five parameter Gaussian model
+.ih
+USAGE
+fitgauss5 image start
+.ih
+PARAMETERS
+.ls image
+Image to be modeled.
+.le
+.ls start
+Starting sample line containing the initial model parameters.
+.le
+.ls lower = -10
+Lower limit for the profile fit relative to each spectrum position.
+.le
+.ls upper = 10
+Upper limit for the profile fit relative to each spectrum position.
+.le
+.ls lines = "*"
+Sample image lines to be fit.
+.le
+.ls spectra = "*"
+Spectra to be fit.
+.le
+.ls naverage = 20
+Number of data lines to be averaged about each sample image line before
+model fitting.
+.le
+.ls factor = 0.05
+The model fit to each line is iterated until the RMS error between the
+model line and the data line improves by less than this factor.
+.le
+.ls track = yes
+Track the model solution from the starting line to the other sample lines?
+.le
+.ls algorithm = 1
+Parameter fitting algorithm to use.  Legal values are 1 and 2.
+.le
+.ls fit_i0 = yes
+Fit the profile scale parameters i0?
+.le
+.ls fit_x0 = yes
+Fit the spectra position parameters x0?
+.le
+.ls fit_s0 = yes
+Fit the spectra shape parameters s0?
+.le
+.ls fit_s1 = no
+Fit the spectra shape parameters s1?
+.le
+.ls fit_s2 = no
+Fit the spectra shape parameters s2?
+.le
+.ls smooth_s0 = yes
+Fit a smoothing spline to the shape parameters s0 after each iteration?
+.le
+.ls smooth_s1 = yes
+Fit a smoothing spline to the shape parameters s1 after each iteration?
+.le
+.ls smooth_s2 = yes
+Fit a smoothing spline to the shape parameters s2 after each iteration?
+.le
+.ls spline_order = 4
+Order of the smoothing spline to be fit to the shape parameters.
+.le
+.ls spline_pieces = 3
+Number of polynomial pieces for the smoothing spline.
+.le
+.ls verbose = no
+Print general information about the progress of the model fitting.
+.le
+.ih
+DESCRIPTION
+The spectra profiles in the interval (\fIlower, upper\fR) about each
+spectrum position are fit with a five parameter Gaussian model for
+the specified sample lines of the image.  For a description of
+the model see \fBgauss5\fR.  The model fitting is performed using 
+simultaneous linearized least squares on the selected model profile
+parameters as determined by the \fIalgorithm\fR for the specified
+\fIspectra\fR.  The parameter fitting technique computes correction
+vectors for the parameters until the RMS error of the model image line
+to the data image line, which is an average of \fInaverage\fR lines
+about the sample line, improves by less than \fIfactor\fR.
+A solution which increases the RMS error of the model is not allowed.
+
+If the parameter \fItrack\fR is yes then the initial model parameters are
+those given in the database for the sample line \fIstart_line\fR.  From
+this starting point the model parameters are iterated to a best fit at
+each specified sample line and then the best fit is used as the starting
+point at the next line.  The tracking sequence is from the starting line
+to the last line and then, starting again from the starting line, to
+the first line.  Note that the model parameters, including the starting
+spectra positions, need be set only at the starting line.
+
+If \fItrack\fR is no then each specified sample line is fitted independently
+from the initial model parameters previously set for that line.  This option
+is used to add additional parameters to the model after an
+initial solution has been obtained or to refit a new image whose database
+was created as a copy of the database of a previously fit image.
+
+The shape parameters s0, s1, and s2 can be smoothed by fitting a spline of
+specified \fIorder\fR and number of spline pieces, \fInpp\fR to the
+parameters as a function of spectra position.
+The smoothing is performed after each iteration and before
+computing the next RMS error.  The smoothing is a form of local constraint
+to keep neighboring spectra from having greatly different shapes.
+The possibility of such erroneous solutions being obtained is present in
+very blended data.
+
+In \fIverbose\fR mode the RMS errors of each iteration are printed on the
+standard output.
+
+The selection of the parameters to be fit and the order in which they are
+fit is determined by \fIalgorithm\fR.  These algorithms are:
+
+.ls 4 1
+This algorithm fits the selected parameters (\fIfit_i0, fit_x0,
+fit_s0, fit_s1, fit_s2\fR) for the selected \fIspectra\fR simultaneously.
+.le
+.ls 4 2
+This algorithm begins by fitting the parameters i0, x0, and s0
+simultaneously.  Note that the values of s1 and s2 are used but are
+kept fixed.  Next the parameters s0 and s1 (the shape) are fit simultaneously
+keeping i0, x0, and s2 fixed followed by fitting i0 and x0 while
+keeping s0, s1, and s2 (the shape) fixed.  If either of these fits
+fails to improve the RMS then the algorithm terminates.
+Also, if after the two steps (the fit of s0 and s1 followed by the fit
+of i0 and x0), the RMS of the fit has not improved by more than the
+user specified factor the algorithm also terminates.  This algorithm has been
+found to be the best way to fit highly blended spectra.
+.le
+.ih
+EXAMPLES
+The default action is to fit Gaussian profiles to the spectra and trace
+the fit from the starting line.  An example of this is:
+
+	cl> fitgauss5 image 1
+
+To fit heavily blended spectra with the four parameter model (i0, x0, s0, s1):
+
+	cl> fitgauss5 image 1 algorithm=2
+.ih
+SEE ALSO
+findspectra
+.endhelp
diff --git a/noao/twodspec/multispec/doc/modellist.hlp b/noao/twodspec/multispec/doc/modellist.hlp
new file mode 100644
index 00000000..70e95ce4
--- /dev/null
+++ b/noao/twodspec/multispec/doc/modellist.hlp
@@ -0,0 +1,52 @@
+.help modellist Jul84 noao.twodspec.multispec
+.ih
+NAME
+modellist -- List data and model pixel values
+.ih
+USAGE
+modellist image lines
+.ih
+PARAMETERS
+.ls image
+Image whose model is to be listed.
+.le
+.ls lines
+Sample lines to be listed.
+.le
+.ls model = "gauss5"
+Profile model to be used to create the model line.
+The only model currently defined is \fIgauss5\fR.
+.le
+.ls columns = "*"
+Image columns to be listed.
+.le
+.ls naverage = 20
+The number of image lines to be averaged to form the data values.
+.le
+.ls lower = -10
+Lower limit of model profiles measured in pixels from the spectra centers.
+.le
+.ls upper = 10
+Upper limit of model profiles measured in pixels from the spectra centers.
+.le
+.ih
+DESCRIPTION
+The model of the image for the selected sample \fIlines\fR
+are used to generate model image lines.  Only the model \fIgauss5\fR is
+currently available.  The output format is column, sample line, image pixel
+value, and model pixel value.  The image pixel data are formed by averaging
+\fInaverage\fR lines about the sample lines.
+.ih
+EXAMPLES
+To list the image and model pixel values for the first sample line after
+fitting the \fIgauss5\fR model with \fBfitgauss5\fR:
+
+	cl> modellist image 1 >outputlist
+
+The list file \fIoutputlist\fR can be used with the \fBlists\fR and
+\fBplot\fR packages to graph the image and model lines or to compute
+and graph residuals.
+.ih
+SEE ALSO
+newimage
+.endhelp
diff --git a/noao/twodspec/multispec/doc/msextract.hlp b/noao/twodspec/multispec/doc/msextract.hlp
new file mode 100644
index 00000000..fa361b38
--- /dev/null
+++ b/noao/twodspec/multispec/doc/msextract.hlp
@@ -0,0 +1,172 @@
+.help msextract Jul84 noao.twodspec.multispec
+.ih
+NAME
+msextract -- Extract spectra from a multi-spectra image
+.ih
+USAGE
+msextract image output
+.ih
+PARAMETERS
+.ls image
+Image to be extracted.
+.le
+.ls output
+Filename for the three dimensional image to be created containing the
+extracted spectra.
+.le
+.ls lower = -10
+Lower limit of the integral for integrated spectra or the first column of the
+strip spectra.  It is measured in pixels from the spectrum center
+defined by the position function in the MULTISPEC database.
+.le
+.ls upper = 10
+Upper limit of the integral for integrated spectra or (approximately) the
+last column of the strip spectra.  It is measured in pixels from the
+spectrum center defined by the position function in the MULTISPEC database.
+.le
+.ls spectra = "*"
+Spectra to be extracted.
+.le
+.ls lines = "*"
+Image lines to be extracted.
+.le
+.ls ex_model = no
+Extract model spectra fit to the image spectra?
+.le
+.ls integrated = yes
+Extract integrated spectra?
+.le
+.ls unblend = no
+Correct for blending in the extracted spectra?
+.le
+.ls clean = yes
+Replace bad pixels with model values?  The following parameters are used:
+.ls  nreplace = 1000.
+Maximum number of pixels to be replaced per image line when cleaning with
+model \fIgauss5\fR or maximum number of pixels to be replaced per spectrum when
+cleaning with model \fIsmooth\fR.
+.le
+.ls sigma_cut = 4.
+Cleaning threshold in terms of sigma of the fit.
+.le
+.ls niterate = 1
+Maximum number of cleaning iterations per line when cleaning with model
+\fIgauss5\fR.
+.le
+.le
+.ls model = "smooth"
+Choice of \fIgauss5\fR or \fIsmooth\fR.  Minimum match abbreviation is
+allowed.  This parameter is required only if \fIex_model\fR = yes
+or \fIclean\fR = yes.
+.le
+.ls naverage = 20
+Number of lines to be averaged in model \fIsmooth\fR.
+.le
+.ls fit_type = 2
+Model fitting algorithm for model \fIgauss5\fR.
+.le
+.ls interpolator = "spline3"
+Type of image interpolation function to be used.
+The choices are "nearest", "linear", "poly3", "poly5", and "spline3".
+Minimum match abbreviation is allowed.
+.le
+.ls verbose = no
+Print verbose output?
+.le
+.ih
+DESCRIPTION
+The MULTISPEC database describing the spectra positions and shapes
+is used to guide the extraction of the spectra in the multi-spectra image.
+The user selects the \fIspectra\fR and image
+\fIlines\fR to be extracted and whether to extract integrated or strip spectra.
+In addition options are available to extract model spectra, replace bad
+pixels by model spectra values, and correct for blending of the spectra.
+The \fIoutput_file\fR three dimensional
+image consists of one band (the third dimension) per extracted spectrum,
+the extracted lines (the second dimension) and either one column for
+the integrated luminosity or the number of columns in the extracted strip.
+
+Integrated spectra (\fIintegrated\fR = yes) are extracted by summing
+the pixel or model values over the specified limits \fIlower\fR and \fIupper\fR
+measured relative to the spectra centers defined by the position functions in
+the database.  Partial pixel sums are used at the endpoints.
+
+Strip spectra (\fIintegrated\fR = no) are extracted by image interpolation
+of the image line or model profiles to obtain a line of values for
+each spectrum and for each image line.  The length of the strip is the
+smallest integer containing the interval between \fIlower\fR and \fIupper\fR.
+The strips for each spectrum are aligned so that the first column is a distance
+\fIlower\fR from the spectrum center as given by the position function in the
+database.
+
+If \fIex_model\fR = yes, \fIunblend\fR = yes, or \fIclean\fR = yes model
+spectra are fit to the spectra in the image.  There are two models:
+a five parameter Gaussian profile called \fIgauss5\fR and profiles obtained
+by averaging \fInaverage\fR image lines surrounding the image line being
+modeled called \fIsmooth\fR.  The model is selected either when the parameter
+\fIunblend\fR = yes or with the parameter \fImodel\fR.  If \fIunblend\fR = yes
+then the model is \fIgauss5\fR regardless of the value of \fImodel\fR.
+
+When \fIex_model\fR = yes the effect is to substitute model spectra for the
+image spectra in the output extraction image.
+
+When \fIclean\fR = yes pixels with large residuals from the model are
+detected and removed from the model fit.  The selected model is
+fit to the pixels which are not in the bad pixel list (not yet implemented)
+and which have not been removed from the model fit.  The sigma of the fit
+is computed.  Deviant pixels are detected by comparing them to the model
+to determine if they differ by more than \fIsigma_cut\fR times the sigma.
+The model fit is iterated, removing deviant pixels at each iteration, until
+no more pixels are found deviant or \fInreplace\fR pixels have been found.
+The pixels removed or in the bad pixel list are then replaced with
+model values.  (To clean an image with this algorithm see \fBnewimage\fR.)
+
+There are some technical differences in the model fitting and cleaning
+algorithms for the two models.  In model \fIsmooth\fR
+the fit for the profile scale factors is done independently for each spectrum
+and automatically corrected when a bad pixel is detected.  This fitting process
+is fast and rigorous.  The parameter \fInreplace\fR in this model refers to
+the maximum number of pixels replaced \fIper spectrum\fR.
+
+In model \fIgauss5\fR, however, the profile scale factors are fit
+to the entire image line (hence its ability to fit blended spectra).
+There are two fitting algorithms; a rigorous simultaneous fit
+and an approximate method.  The simultaneous fit is selected when
+\fIfit_type\fR = 1.  This step is relatively slow. The
+alternative method of \fIfit_type\fR = 2 sets the scale factor for each
+spectrum by taking the median scale, where scale = data / model profile,
+for the three pixels nearest the center of the profile.  The median
+minimizes the chance of a large error due to a single bad pixel.  This
+scale may be greatly in error in the case of extreme blending but is also
+quite fast; the extraction time is reduced by at least 40%.
+The steps of profile fitting and deviant pixel detection are alternated
+and the maximum number of iterations through these two steps is
+set by \fIniterate\fR.  The default of 1 means that the model fitting is not
+repeated after detecting deviant pixels.
+
+When \fIunblend\fR = yes the \fIgauss5\fR model
+is fitted to the image spectra (including possible cleaning).
+The relative contributions to the total image pixel value from each of the
+blended spectra are determined from the model and applied toward either the
+integrated or strip spectra.  If \fIex_model\fR = yes then this option has
+no effect other than to force the selection of model \fIgauss5\fR.
+
+The option \fIverbose\fR is used to print the image lines being extracted
+and the number of pixels replaced by the cleaning process.
+.ih
+EXAMPLES
+To extract all the integrated spectra from all the image lines:
+
+	cl> msextract image image.ms
+
+To extract model strip spectra:
+
+	cl> msextract image image.ms ex_model=yes int=no
+
+To extract integrated spectra without any modeling:
+
+	cl> msextract image image.ms clean=no
+.ih
+SEE ALSO
+newimage
+.endhelp
diff --git a/noao/twodspec/multispec/doc/mslist.hlp b/noao/twodspec/multispec/doc/mslist.hlp
new file mode 100644
index 00000000..461b52b4
--- /dev/null
+++ b/noao/twodspec/multispec/doc/mslist.hlp
@@ -0,0 +1,77 @@
+.help mslist Jul84 noao.twodspec.multispec
+.ih
+NAME
+mslist -- List entries in a MULTISPEC database
+.ih
+USAGE
+mslist image keyword lines spectra
+.ih
+PARAMETERS
+.ls image
+Image whose MULTISPEC database entries are to be listed.
+.le
+.ls keyword
+Keyword for the database entry to be listed.  The keywords are:
+.ls header
+List general header information.
+.le
+.ls comments
+List the comments.
+.le
+.ls samples
+List the sample image lines.
+.le
+.ls x0
+List the spectra positions for the specified sample lines and spectra.
+.le
+.ls i0
+List the model profile scales for the specified sample lines and spectra.
+.le
+.ls s0, s1, or s2
+List the gauss5 model shape parameter s0, s1, or s2 for the specified sample
+lines and spectra.
+.le
+.ls gauss5
+List the gauss5 model parameters x0, i0, s0, s1, and s2 for the specified
+sample lines and spectra.
+.le
+.ls x0 spline
+List the spline evaluation of the spectra positions for the specified
+image lines and spectra.
+.le
+.ls s0 spline, s1 spline, or s2 spline
+List the spline evaluation of the gauss5 model shape parameters s0, s1, or s2
+for the specified image lines and spectra.
+.le
+.le
+.ls lines
+Lines to be listed.  For the entries x0, i0, s0, s1, s2, and gauss5 the
+lines refer only to the sample image lines. For the spline entries the
+lines refer to the image lines at which the spline is to be evaluated.
+.le
+.ls spectra
+Spectra to be listed.
+.le
+.ls titles = no
+Print additional titles?
+.le
+.ih
+DESCRIPTION
+This task is a general MULTISPEC database listing tool.  A keyword is selected
+and the referenced data is listed.  Some entries require the specification of
+the desired sample or image lines and the desired spectra.
+.ih
+EXAMPLES
+To list the spectra positions for spectrum 3 at all the sample lines:
+
+	cl> mslist image x0 "*" 3
+
+To list the model profile scale parameter for sample line 1:
+
+	cl> mslist image i0 1 "*"
+
+To list the gauss5 model parameters for spectra 2 and 3 and sample lines 5
+and 7:
+
+	cl> mslist image gauss5 "5,7" "2-3" titles+
+.endhelp
diff --git a/noao/twodspec/multispec/doc/msplot.hlp b/noao/twodspec/multispec/doc/msplot.hlp
new file mode 100644
index 00000000..f08eac1b
--- /dev/null
+++ b/noao/twodspec/multispec/doc/msplot.hlp
@@ -0,0 +1,44 @@
+.help msplot Oct85 noao.twodspec.multispec
+.ih
+NAME
+msplot -- Plot data and model image line
+.ih
+USAGE
+msplot image line
+.ih
+PARAMETERS
+.ls image
+Image to be plotted.
+.le
+.ls line
+The image line to be plotted.  Actually the nearest sample line will be
+plotted.
+.le
+.ls naverage = 20
+Number of image lines to average about the specified line.
+.le
+.ls lower = -10., upper = 10.
+Limits of the model profiles relative to the center of each profile.
+.le
+.ls graphics = "stdgraph"
+Graphics output device.
+.le
+.ls cursor = ""
+Graphics cursor input.  If a file is given then the cursor input is taken
+from the file.  If no file is given then the standard graphics cursor will
+be used.
+.le
+.ih
+DESCRIPTION
+A line of image data and the profile model for the line is graphed.
+The model is graphed with a dashed line.  The graph may be then expanded,
+manipulated, and printed with the standard cursor mode commands.
+.ih
+EXAMPLES
+To plot the model fit for image sample for image line 400:
+
+	cl> msplot sample 400
+.ih
+SEE ALSO
+modellist
+.endhelp
diff --git a/noao/twodspec/multispec/doc/msset.hlp b/noao/twodspec/multispec/doc/msset.hlp
new file mode 100644
index 00000000..689e525a
--- /dev/null
+++ b/noao/twodspec/multispec/doc/msset.hlp
@@ -0,0 +1,104 @@
+.help msset Jul84 noao.twodspec.multispec
+.ih
+NAME
+msset -- Set entries in a MULTISPEC database
+.ih
+USAGE
+msset image keyword value
+.ih
+PARAMETERS
+.ls image
+Image in which the MULTISPEC database entries are to be modified or initialized.
+.le
+.ls keyword
+Keyword for the database entry to be set.  The keywords are:
+.ls nspectra
+Set the number of spectra in the header.
+.le
+.ls comments
+Add comments lines to the database comment block.
+.le
+.ls x0
+Set the spectra positions for the specified sample lines and spectra.
+.le
+.ls i0
+Set the model profile central intensities for the specified sample lines
+and spectra.
+.le
+.ls s0, s1, or s2
+Set the gauss5 model shape parameter s0, s1, or s2 for the specified sample
+lines and spectra.
+.le
+.le
+.ls value
+Value to be used for value input.
+.le
+.ls lines = "*"
+Sample lines to be affected by value input.
+.le
+.ls spectra = "*"
+Spectra to be affected by value input.
+.le
+.ls read_list = no
+If yes use list input and if no use value input.
+.le
+.ls list = ""
+List for list input.  See the description below for the appropriate format.
+.le
+.ih
+DESCRIPTION
+The entries in a MULTISPEC database associated with the image
+are modified or initialized.
+The parameters \fIimage\fR, \fIkeyword\fR, and \fIread_list\fR 
+determine the database to be operated upon, the database entry to
+be set, and the input type.  There are two forms of input;
+list input and value input.
+The input type is selected by the boolean parameter
+\fIread_list\fR.  For list input the parameter \fIlist\fR
+is used and for value input the parameter \fIvalue\fR and
+possibly the parameters \fIlines\fR and \fIspectra\fR are used.
+The required parameters and input formats for the different keywords
+are outlined below.
+.ls nspectra
+For list input the list format is the number of spectra and
+for value input the \fIvalue\fR parameter is the number of spectra.
+.le
+.ls comments
+For list input the list format is lines of comments and for value
+input \fIvalue\fR parameter is a comment string.
+.le
+.ls x0, i0, s0, s1, s2
+For list input the list format is sample line, spectrum number, and
+parameter value
+and for value input \fIlines\fR is a range string selecting the
+sample lines to be affected, \fIspectra\fR is a range string selecting
+the spectra to be affected, and \fIvalue\fR is the value to be set for all
+the selected lines and spectra.
+.le
+.ih
+EXAMPLES
+To add several comments to the database by query:
+
+.nf
+	cl> msset image "comments" read_list+
+	Input list> First comment here.
+	Input list> Second comment here.
+	Input list> <eof>
+.fi
+
+where <eof> is the end of file character terminating the list.
+To set the value of s0 to 1 for all the spectra in sample line 1:
+
+	cl> msset image "s0" 1
+
+To set the spectra positions from a list:
+
+	cl> msset image "x0" read_list+ list=positionlist
+
+To add a single comment such as in a script:
+
+	cl> msset image "comments" "Comment here."
+.ih
+SEE ALSO
+findspectra mslist
+.endhelp
diff --git a/noao/twodspec/multispec/doc/multispec.ms b/noao/twodspec/multispec/doc/multispec.ms
new file mode 100644
index 00000000..cc17352e
--- /dev/null
+++ b/noao/twodspec/multispec/doc/multispec.ms
@@ -0,0 +1,532 @@
+.EQ
+delim	$$
+.EN
+.TL
+The Multi-Spectra Extraction Package (multispec)
+.AU
+Francisco Valdes
+.AI
+IRAF Group
+.K2
+October 1984
+.NH
+Introduction
+.PP
+This document provides an introduction and overview of the multi-spectra
+extraction package \fBmultispec\fR.  Detailed descriptions and usage
+information for the tasks of the package are available in the manual
+pages.  The tasks in the package are:
+
+.TS
+center;
+n.
+findpeaks \&- Find the peaks
+fitfunction \&- Fit a function to the spectra parameter values
+fitgauss5 \&- Fit spectra profiles with five parameter Gaussian model
+modellist \&- List data and model pixel values
+msextract \&- Extract spectra
+mslist \&- List entries in a MULTISPEC database
+msplot \&- Plot a line of image and model data
+msset \&- Set entries in a MULTISPEC database
+newextraction \&- Create a new MULTISPEC extraction database
+newimage \&- Create a new multi-spectra image
+.TE
+
+.PP
+The \fBmultispec\fR package is a subpackage of the \fBtwodspec\fR package.
+It provides tools to locate, model, clean, correct for blending,
+and extract integrated or strip spectra from two dimensional, multi-spectra
+images.  These tools may be used directly or combined in scripts to
+extract specific types of spectra or spectra from specific instruments.
+Examples of the latter usage are the tasks in the image reduction package
+\fBcryomap\fR.
+.PP
+The extraction of spectra consists of locating pixels along each
+image line which intersect the spectra and recording either the sum of
+the pixels, \fIintegrated spectra\fR (some times referred to as
+one-dimensional spectra), or the set of pixels,
+\fIstrip spectra\fR, for each line and for each spectrum as output.
+The size and limits of the intersection region are specified by the
+user relative to the centers of the spectra.
+The locations of the spectra in each image line are determined separately
+so that the spectra need not be aligned along the columns of the image nor
+be perfectly straight.  However, since the extraction is done by image line,
+if the spectra are not aligned with the columns then the spectral resolution
+will be decreased.  If the spectra are aligned with the image lines then
+the image should be rotated or transposed using \fBimtranspose\fR.
+.PP
+The \fBmultispec\fR extraction produces three dimensional images with
+one image band (the third dimension) for each extracted spectrum
+and one line (the second dimension) for each extracted image line.
+For integrated spectra there is only one column
+while for strip spectra, the number of columns is equal to the extraction
+strip width.  The strips are aligned to the same positions relative to the
+spectra centers by image interpolation.  If desired the output extractions can
+be reformated in a variety of ways.
+.PP
+In addition to direct extraction of the image spectra the \fBmultispec\fR
+package provides for modeling the spectrum profiles.  The model
+may be extracted instead of the image spectra as either integrated or
+strip spectra.  The model may be used to correct for blending of the spectra
+and to detect and replace bad pixels.  The cleaning replaces data pixels which
+are discrepant from the model by the model values.
+.PP
+The modeling and cleaning features of the \fBmultispec\fR package can also
+be used for creating new multi-spectra images.  In other words a new
+image is created containing cleaned or model spectra for selected
+lines.
+.PP
+Section 2 gives an overview of the \fBmultispec\fR package and the extraction
+process.  The next section briefly describes the tasks in the package.
+This is followed by a description of the extraction database.
+The final section defines the model profiles used in the \fBmultispec\fR
+package.
+.NH
+Overview of the Multispec Package and the Extraction Process
+.PP
+The \fBmultispec\fR package consists of general and flexible tools
+for creating and manipulating databases which describe multi-spectra
+images.  The contents of the databases are described in a later section.
+Each database is associated with a particular image and is referenced
+through the image name.  The first positional argument in all the
+\fBmultispec\fR tasks is an image.  In the current version of the package each
+database exists as a separate binary file with a filename formed by adding
+the extension '.db' to the image name.  Note, however, that this
+need not be the case in the future.
+.PP
+The organization of the package as a set of tools operating on a database
+allows room for the package to evolve.  Different algorithms may be
+designed for different types of multi-spectra images by using combinations
+of the existing tools and by adding new tools.  The discussion below
+points out areas where new tasks might be added as well as citing the
+applicable existing tasks.
+.PP
+The extraction of spectra from a multi-spectra image consists of two
+basic steps; determining the locations of the spectra in the image and
+extracting the spectra.  The positions of the spectra in a multi-spectra
+image are determined at a set of "sample" image lines.  These positions
+are used to fit an interpolation function defining the spectrum positions
+at all the image lines.  This function is then used in the extraction of
+the spectra.
+.PP
+The sample image lines are chosen by the user when the database is first
+created by the task \fBnewextraction\fR.  An exception to this is when
+a template image is used (discussed below).  However, in this case the
+sample image lines are still those chosen by the user when the template
+image database was created.  The sample image lines may consist of
+anywhere from one image line to all the image lines.  The purpose
+of the sample lines is to sample the image often enough to follow changes
+in the positions and shapes of the spectra but to still minimize the
+time spent in finding the spectra and the size of the database.  The choice
+of sample lines also depends on the algorithm used to determine the
+positions of the spectra; a large number of sample
+lines for a fast, approximate method and a smaller number of lines
+for a complex and accurate method.  For example, in order to deal with
+very blended spectra the task \fBfitgauss5\fR provides a sophisticated
+model fitting algorithm.  This technique is computationally slow and, so,
+the user should not choose too many sample lines.
+.PP
+After the database has been created the minimum information needed
+for extraction is the spectrum positions at the sample lines.  There
+are many ways in which the positions may be determined.  Some
+possibilities are listed below.
+
+.IP (1)
+Enter the spectrum positions from a list using \fBmsset\fR.  The
+list might be generated from a graphics/cursor task.
+This is method is very time consuming when the number of spectra and
+the number of images are large.
+.IP (2)
+Determine the spectrum positions automatically by finding the peaks in
+each sample image line.  The task \fBfindpeaks\fR performs this function.
+.IP (3)
+Determine the spectrum positions at just one sample image line
+using either (1) or (2) and trace the spectra by a fast and refined
+peak finding method.  Such a task is desirable but is not a part of the
+current package.
+.IP (4)
+Determine the spectrum positions at just one sample image line
+using either (1) or (2) and trace the spectra by fitting model
+spectrum profiles.  The task \fBfitgauss5\fR does this using
+the model gauss5 described in section 5.  Additional model fitting
+tasks can be added as needed.
+.IP (5)
+Use the positions determined for a previous image and, if necessary,
+refine the positions.  \fBFitgauss5\fR is used to
+refine the spectrum positions at each sample line independently.
+
+.PP
+Several position finding algorithms may be used in stages to achieve
+the degree of accuracy required by the user.
+Thus, the first position determinations may be relatively crude and
+then, if needed, more sophisticated methods may be applied to refine the
+positions.  The task \fBfindpeaks\fR is a crude peak finder.  The positions
+are only determined to the nearest pixel.  The task \fBfitgauss5\fR is
+a sophisticated model fitting techique which is used after \fBfindpeaks\fR
+first determines the approximate positions of the spectra.
+.PP
+The determination of the spectra locations may be performed independently
+at each sample line as in (1) and (2) above or the spectra locations may
+be traced starting from one sample line as in (3) and (4).  The second method
+is preferable.  Generally, \fBfindpeaks\fR is used at only one sample line
+to initially determine the number and approximate locations of the spectra.
+\fBFitgauss5\fR then fits model gauss5 to the spectrum profiles and
+the model solution is used at the next sample line as the starting
+point for the next model fit.  In this manner the positions of
+the spectra are determined at the other sample image lines.
+.PP
+The results of the peak finding and profile fitting are improved
+by using an average of many image lines about the sample image line rather
+than just the sample image line by itself.  Both \fBfindpeaks\fR and
+\fBfitgauss5\fR have this ablility.
+.PP
+It is often the case that several multi-spectra images have essentially
+the same format; i.e. the same image size, the same number of spectra,
+and the same positions (either approximately or identically).
+Commonly, one of the images is used for calibrations and has strong,
+high signal-to-noise spectra while the other images have weaker spectra.
+In this case it is not necessary to repeat the position determinations.
+The spectrum positions in one of the images, generally the one with
+the strong calibration spectra, are determined first.  This image is
+then used as a "template" to provide the initial position estimates for
+the other images.  If the positions are identical no further work is needed,
+otherwise, the positions can be refined to correct for small changes in the
+positions and shapes of the spectra.
+.PP
+The task \fBnewextraction\fR creates new databases.  If a template image
+is specified then a copy is made of the template image database.  This means
+that the number of spectra and the sample image lines remain the same.
+If the spectrum positions are slightly different from the template image
+then the task \fBfitgauss5\fR is used to determine the new positions.
+.PP
+The spectrum positons and possibly any model parameters are interpolated
+from the sample lines to the remaining image lines by fitting a function
+to values at the sample lines.  In addition, the function fits may
+leave out poorly determined points and also smooth the values at the
+sample lines.  The task \fBfitfunction\fR fits selected functions of
+specified order to the selected spectra and sample image lines.
+.PP
+The extraction of the spectra from multi-spectra images is performed by
+the task \fBmsextract\fR.  The task extracts either integrated or strip
+spectra, either data or model values, with or without blending corrections,
+and with or without replacing bad pixels by model values.
+The user specifies the limits of the extraction
+strip as well as the spectra and image lines to be extracted.
+.PP
+For the simplest type of data extractions (basically strip extraction)
+no modeling is required.  Other types of extractions, such as model
+extractions and/or with cleaning and blending corrections require some
+degree of modeling.  There are two models which may be used;
+"smooth" and "gauss5".  These models are described in section 5.
+The model parameters for model gauss5 must be set by \fBfitgauss5\fR
+before \fBmsextract\fR is used.  Additional models may added for
+extraction as well as for the spectrum position determinations.
+.PP
+The model based features of \fBmsextract\fR -- model extractions
+and cleaning -- are available in the related task \fBnewimage\fR.
+This task creates new images which consist of either model spectra
+or cleaned data spectra.
+.PP
+The models in the \fBmultispec\fR package assume that the profiles
+go to zero; i.e. there is no background light.  Background light
+may be removed using \fBbackground\fR.  In the future a task will
+be provided create a mask defining the locations of the spectra from
+the database which can be used with general surface fitting tasks
+to create a background surface to be subtracted from the image.
+.PP
+The final step in using the \fBmultispec\fR package is to convert the
+extraction output to the desired format.  This may include graphs,
+card image formats, and files for the \fBonedspec\fR and \fBlongslit\fR
+packages.  Currently, the available formats are images and IIDS
+card images.
+.NH
+The Tasks of the Multispec Package
+.PP
+Use of the \fBmultispec\fR package begins with \fBnewextraction\fR and
+ends, usually, with \fBmsextract\fR.  In between there are tasks which
+update, refine or change the database and tasks which provide diagnositic
+information.  The informational tasks can be combined with tasks from
+other packages to produce tabular or graphical output.  The task
+\fBmsplot\fR is an example.  In this section a brief description of
+each task is given.  Further information about the tasks, including usage,
+is available in the manual pages.
+.SH
+findpeaks
+.IP
+Selected sample image lines are examined to determine the number and
+column positions of data peaks in the line.  An average of a number of image
+lines surrounding the sample lines is formed in which the local maxima
+are located.  Various criteria are applied to cull the list of local
+maxima to the desired peaks.  These criteria include a peak threshold,
+a maximum peak-to-peak contrast, a minimum peak separation, and a
+maximum number of peaks.  This task is used to determine crude, initial
+estimates for the spectrum positions.  It could be used alone for
+simple extractions.
+.SH
+fitfunction
+.IP
+This task has two roles.  It's primary role is to define the
+interpolation/extrapolation function for the spectra
+positions between the sample lines.  The fitting function can be
+either purely interpolative or may also provide smoothing of the
+parameters from the sample lines.  The second role is to provide
+smoothing of the model parameters along the dispersion and the
+ability to replace bad values by the function fit to the remaining
+parameters.  In this second role the user may iterate between
+smoothing and model fittng.  The functions are always defined between
+the first and last image lines.
+.SH
+fitgauss5
+.IP
+The model profiles gauss5, described in section 5, are fit to the
+selected spectra and sample lines.  The parameters to be determined
+and the fitting algorithm may also be selected.
+The model parameters are recorded in the database.
+The model may be tracked from a starting line to other sample image
+lines or each sample line may be fitted independently.
+This task is used to accurately determine the spectrum positions
+and provide an extraction model for heavily blended spectra.
+.SH
+modellist
+.IP
+For the selected sample image lines and image columns data
+and model values are listed.  This task is used to check how well
+the model fitting tasks (currently just \fBfitgauss5\fR) have fit
+the sample image line.  The task \fBmsplot\fR is used to produce
+graphical output.
+.SH
+msextract
+.IP
+This task does the actual extraction of spectra.  It requires that
+the spectrum positions are defined by fitting functions in the
+database.  If model gauss5 is to be used then the database must
+also contain the model parameters for the sample image lines.  It
+extracts integrated or strip spectra, using data or model values,
+with or without blending corrections, and with or without cleaning
+of bad pixels.
+.SH
+mslist
+.IP
+Of the diagnositic or informational tasks \fBmslist\fR is the most
+general.  The user selects the type of information from the database
+which is desired and it is then printed.  The types of information
+include the database header, the database comments, the spectra
+positions and model parameter values for the sample lines, and the
+interpolation/smoothing function values for any desired set of
+image lines.
+.SH
+msplot
+.IP
+This task extracts data and models values and plots them superposed.
+This task is used as a diagnositic tool to inspect how well model fitting
+represents the image spectra.
+.SH
+msset
+.IP
+This task is a general tool for modifying or setting some of the quantities 
+in the database.  The quantity to be changed or set is
+selected by a keyword and the values are input in two ways;
+with a list structured parameter (such a file containing the list of
+values or the standard input) or as a parameter value.  This task
+is the way a user may enter comments in the database or manually
+set the number and positions of the spectra.  It is also used to
+set the initial values for the gauss5 model parameters s0, s1, and s2
+prior to using \fBfitgauss5\fR.
+.SH
+newextraction
+.IP
+This task has three important roles.  First it creates the database
+associated with the multi-spectra image.  Second, it defines the sample
+image lines to be used.  The user can specify as many or as few sample lines
+as desired.  It should be kept in mind that the more sample lines used
+the larger the database becomes and the longer the processing time when
+modeling the spectra.  Finally, \fBnewextraction\fR allows
+a database from another image (called a template image) to initialize the
+database for the new multi-spectra image.  The template image is generally
+a calibration image with strong, well-defined spectra.
+Initializing a database with a template image saves time, reduces problems
+with bad pixels, and is more accurate when an image with weak spectra is
+to be extracted.
+.SH
+newimage
+.IP
+This task is similar to \fBmsextract\fR; it uses the same algorithms
+and parameters.  It differs in the type output.
+Rather than producing extracted integrated or strip spectra this task
+produces new image lines.  It is particularly useful for extracting
+model images to be compared against the original image or to
+produce images which have been cleaned.
+.NH
+The Multispec Database
+.PP
+The tasks in the \fBmultispec\fR package create and manipulate a database.
+The database contains a description of the multi-spectra image which
+is modified, refined, examined, or otherwise used by the tasks in the package.
+In the current version the database is a separate binary file with a filename
+formed by appending ".db" to the image name described by the database.
+.PP
+The database contains four basic types of data; general information,
+comments and history, position parameters, and model parameters.
+The data in the database is examined with the task \fBmslist\fR.
+The general information section, called the database header, contains the
+the name of the image, the size of the image, and the number of spectra in
+the image.  Once the number of spectra in the image has
+been entered in the database it is an error to attempt to change this
+number.  The database must be deleted and a new database created in order
+to change the number of spectra.
+.PP
+The comment and history section of the database contains text
+strings.  Each task which modifies the contents of the database places
+a dated history line in this section.  The user may also add comments
+with \fBmsset\fR.  Currently this information is not passed on to
+the extraction output.
+.PP
+There are three types of position information in the database.  The
+first is a set of sample image lines.  The sample lines are set when
+the database is created by \fBnewextraction\fR.  The sample lines select
+which image lines from the multi-spectra image are to be examined and used
+during the extraction.  Information from these sample lines, and only
+these sample lines, is entered in the database.  The sample lines
+may be listed with \fBmslist\fR.
+.PP
+The second type of position information is the positions of the
+spectra (centers) at each sample line.  These positions are initially
+set by either \fBfindpeaks\fR or, manually, by \fBmsset\fR.  The
+position information is refined by fitting model profiles.
+.PP
+The third type of position information is a function fit to the
+positions from all the sample lines for each spectrum.
+These function fits are produced by \fBfitfunction\fR.
+The functions define the positions of the spectra at all the image
+lines.  The spectra positions at the sample lines or the function
+evaluation for any image line may be listed with \fBmslist\fR.
+.PP
+The finally type of basic data contained in the database are
+model parameter values.  A model need not be used in the extraction
+but if one is used then the parameters determining the model profiles
+are recorded in the database.  The specific parameters depend on the
+model.  Currently the only model is \fIgauss5\fR.  The model and its
+parameters are described in section 5.
+.PP
+As with the spectra positions the parameters are stored in the database
+in two forms; as values for each spectrum at each sample image line
+and as function fits to the values at the sample lines which interpolate
+them to any image line.  The sample line values are
+set by the model fitting tasks and the function fits are set by
+\fBfitfunction\fR.  The parameter values at the sample lines or the
+function evaluations for any image lines may be listed with \fBmslist\fR.
+.NH
+Multispec Spectrum Profile Models
+.PP
+The spectra profiles in the image are modeled for many reasons:
+To provide accurate, subpixel position determinations, to extract model
+spectra or model images, to detect and replace bad pixels, and
+to estimate and correct for blending between the spectra.
+There are currently two models used in the \fBmultispec\fR package, "gauss5"
+and "smooth".
+.NH 2
+Model Gauss5
+.PP
+The gauss5 model profiles are Gaussian but with a scale which varies
+smoothly between the center and the edge of the profile.  There
+are five parameters:
+
+.RS
+.IP x0
+The column position in the image line of the center of the profile.
+.IP i0
+The intensity scale of the profile.  It corresponds to the intensity
+of the center of the profile.
+.IP s0
+The zeroth order, constant, term in the Gaussian scale.
+.IP s1
+The even first order term in the Gaussian scale.
+.IP s2
+The odd first order term in the Gaussian scale.
+.RE
+
+.PP
+The mathematical form of the the model is shown in equation (1):
+.EQ (1)
+roman profile (x)~=~i0 exp~left { -s( DELTA x )~DELTA x sup 2 right }
+.EN
+where
+.EQ
+DELTA x ~=~x~-~x0~,
+.EN
+.EQ
+s( DELTA x)~=~s0~+~s1~|y| +~s2~y~,
+.EN
+and
+.EQ
+y~=~ DELTA x / ( DELTA x sup 2 + alpha ) sup half ~.
+.EN
+The profile is defined within the user specified limits \fIlower\fR and
+\fIupper\fR measured relative to the the profile center and
+$alpha~=~(upper-lower)/4$.  The quantity $y$ lies in the range
+-1 to 1 over the interval in which the profile is defined.  The odd
+and even terms, s1 and s2, allow for symmetric and antisymmetric profile
+changes relative to a simple Gaussian profile.
+.PP
+The task \fBfitgauss5\fR fits the gauss5 model to the spectrum profiles in
+the sample image lines to determine one or more of the model parameters for
+each spectrum.  The parameter values are stored in the database for the image.
+In \fBmsextract\fR the model profiles for each
+image line are obtained by interpolating the profile shapes from the sample
+lines (with the model parameters in the database determined by
+\fBfitgauss5\fR) and then fitting only the intensity scale "i0".
+There are a number of technical details associated with the model fitting
+in each of these tasks which are discussed in the manual pages.
+.PP
+The gauss5 model is used to accurately determine the positions of the
+spectrum centers at the sample image lines.  Fitting simultaneously
+for the model parameters allows the spectra to be blended.
+This is the chief advantage of this model.
+This model is also used during extraction to correct for blending of
+the spectra and to detect and replace bad pixels.
+.NH 2
+Model Smooth
+.PP
+The spectrum profiles from the lines immediately preceeding
+the image line in which the spectrum profile is to be fit are shifted
+to a common center and averaged to form the model profile.
+An intensity scale factor is then determined which best fits the model
+profile to the image profile.  This is done for each spectrum in the
+image.  The scale factors are determined by least squares with
+possible bad pixel rejection.  Rejected pixels are eliminated
+when the image line is later used in forming new average model profiles.
+.PP
+The advantages of this model are that the image spectrum profiles may
+have any shape and the least squares fitting with bad pixel rejection
+is fast and rigorous.  By passing through the image lines sequentially
+the image lines need be accessed only once and the profile averages
+can be quickly updated for the next image line.
+.PP
+The disadvantages of this model are that the spectrum profiles cannot
+be blended and the model does not measure profile positions.
+This means that the spectrum profile positions must be
+known.  This model is suitable for model extractions and cleaning of
+bad pixels in unblended multi-spectra images.  It is available in
+the task \fBmsextract\fR.
+.bp
+.SH
+Glossary
+.LP
+\fBmultispec\fR
+.IP
+Acronym for Multi-Spectra Extraction as in \fBmultispec\fR Package.
+.LP
+integrated spectra
+.IP
+The spectra are extracted by integrating the pixel values across the spectrum
+to produce a single aperture luminosity value.
+.LP
+sample image line
+.IP
+The spectra positions and model profile shapes are determined at a set
+of image lines selected when the database is created.
+.LP
+strip spectra
+.IP
+The spectra are extracted as a strip of fixed with the spectra shifted by
+image interpolation to a common center.
diff --git a/noao/twodspec/multispec/doc/newextract.hlp b/noao/twodspec/multispec/doc/newextract.hlp
new file mode 100644
index 00000000..37123f28
--- /dev/null
+++ b/noao/twodspec/multispec/doc/newextract.hlp
@@ -0,0 +1,61 @@
+.help newextraction Jul84 noao.twodspec.multispec
+.ih
+NAME
+newextraction -- Initialize a new MULTISPEC extraction
+.ih
+USAGE
+newextraction image template
+.ih
+PARAMETERS
+.ls image
+Image to be extracted.
+.le
+.ls template
+The previously created database for the template image is used to initialize
+the new database.  If the null string is given then the database is not
+initialized.
+.le
+.ls sample_lines = "10x50"
+Sample image lines in which the spectra positions are to be determined and,
+optionally, modeled.  This parameter is not used if a template image is given.
+.le
+.ih
+DESCRIPTION
+To extract the spectra from a multi-spectra image a database must be created
+and associated with the image.  This task creates the database with a name
+formed by adding the extension '.db' and initializes some of the database
+entries.
+
+The sample lines are used to track the spectra positions and, if an analytic
+profile model is to be fit to the spectra, to map profile shape changes.
+The image lines only need be sampled enough to track \fInon-linear\fR position
+distortions and significant profile shape changes since interpolation
+is used between the sample lines.  Though specifying just one sample
+line is allowed using at least two sample lines is recommended to allow for
+any slope in the position of the spectra.  Specifying all the image lines
+will greatly increase the processing time and is never justified.
+
+Using a previous database to initialize the new database is useful if the
+new image is only slightly different in the positions and profiles of the
+spectra.  In some cases extraction may proceed immediately without any
+further position determination and modeling.  Further modeling
+and spectra position determinations will refine the previously determined
+parameters with an increase in execution time.  Using a template image is
+particularly important if the first image extracted has strong spectra
+and subsequent images have much weaker spectra since the automatic spectra
+position location and profile modeling may yield poor results for very weak
+spectra.
+.ih
+EXAMPLES
+To initialize a MULTISPEC database for extracting the spectra in
+the image \fIimage1\fR:
+
+	cl> newextraction image1 ""
+
+To create a new MULTISPEC database for extracting the spectra in
+the image \fIimage2\fR  using \fIimage1\fR as a template image:
+
+.nf
+	cl> newextraction image2 image1
+.fi
+.endhelp
diff --git a/noao/twodspec/multispec/doc/newimage.hlp b/noao/twodspec/multispec/doc/newimage.hlp
new file mode 100644
index 00000000..1ef7fbe0
--- /dev/null
+++ b/noao/twodspec/multispec/doc/newimage.hlp
@@ -0,0 +1,130 @@
+.help newimage Jul84 noao.twodspec.multispec
+.ih
+NAME
+newimage -- Create a new multi-spectra image
+.ih
+USAGE
+newimage image output
+.ih
+PARAMETERS
+.ls image
+Image to be used to create the new image.
+.le
+.ls output
+Filename for the new multi-spectra image.
+.le
+.ls lower = -10
+Lower limit for model profiles.  It is measured in pixels from the
+spectra centers defined by the position functions in the database.
+.le
+.ls upper = -10
+Upper limit for model profiles.  It is measured in pixels from the
+spectra centers defined by the position functions in the database.
+.le
+.ls lines = "*"
+Image lines of the multi-spectra image to be in the new multi-spectra image.
+.le
+.ls ex_model = no
+Create a model image?
+.le
+.ls clean = yes
+Replace bad pixels with model values?  The following parameters are used:
+.ls  nreplace = 1000.
+Maximum number of pixels to be replaced per image line when cleaning with
+model \fIgauss5\fR or maximum number of pixels to be replaced per spectrum when
+cleaning with model \fIsmooth\fR.
+.le
+.ls sigma_cut = 4.
+The cleaning threshold in terms of the predicted pixel sigma.
+.le
+.ls niterate = 1
+Maximum number of cleaning iterations per line when cleaning with model
+\fIgauss5\fR.
+.le
+.le
+.ls model = "smooth"
+Choice of \fIgauss5\fR or \fIsmooth\fR.  Minimum match abbreviation is
+allowed.  This parameter is required only if \fIex_model\fR = yes
+or \fIclean\fR = yes.
+.le
+.ls fit_type = 2
+Model fitting algorithm for model \fIgauss5\fR.
+.le
+.ls naverage = 20
+Number of lines to be averaged in model \fIsmooth\fR.
+.le
+.ls interpolator = "spline3"
+Type of image interpolation function to be used.
+The choices are "nearest", "linear", "poly3", "poly5", and "spline3".
+Minimum match abbreviation is allowed.
+.le
+.ls verbose = no
+Print verbose output?
+.le
+.ih
+DESCRIPTION
+A new multi-spectra image is created using the description of the
+multi-spectra image in the MULTISPEC database associated with \fIimage\fR.
+The user selects the image \fIlines\fR from the original image to be in
+the new image.  The options allow the creation of model images or images in
+which the bad or deviant pixels are replaced by model profile values.
+
+If \fIex_model\fR = yes or \fIclean\fR = yes model
+spectra are fit to the spectra in the image.  There are two models:
+a five parameter Gaussian profile called \fIgauss5\fR and profiles obtained
+by averaging \fInaverage\fR image lines surrounding the image line being
+modeled called \fIsmooth\fR.  The model is selected with the parameter
+\fImodel\fR.
+
+When \fIex_model\fR = yes an image containing model spectra is produced.
+
+When \fIclean\fR = yes pixels with large residuals from the model are
+detected and removed from the model fit.  The selected model is
+fit to the pixels which are not in the bad pixel list (not yet implemented)
+and which have not been removed from the model fit.  The sigma of the fit
+is computed.  Deviant pixels are detected by comparing them to the model
+to determine if they differ by more than \fIsigma_cut\fR times the sigma.
+The model fit is iterated, removing deviant pixels at each iteration, until
+no more pixels are found deviant or \fInreplace\fR pixels have been found.
+The pixels removed or in the bad pixel list are then replaced with
+model values.  (To clean and extract the spectra with this algorithm see
+\fBmsextract\fR.)
+
+There are some technical differences in the model fitting and cleaning
+algorithms for the two models.  In model \fIsmooth\fR
+the fit for the profile scale factors is done independently for each spectrum
+and automatically corrected when a bad pixel is detected.  This fitting process
+is fast and rigorous.  The parameter \fInreplace\fR in this model refers to
+the maximum number of pixels replaced \fIper spectrum\fR.
+
+In model \fIgauss5\fR, however, the profile scale factors are fit
+to the entire image line (hence its ability to fit blended spectra).
+There are two fitting algorithms; a rigorous simultaneous fit
+and an approximate method.  The simultaneous fit is selected when
+\fIfit_type\fR = 1.  This step is relatively slow. The
+alternative method of \fIfit_type\fR = 2 sets the scale factor for each
+spectrum by taking the median scale, where scale = data / model profile,
+for the three pixels nearest the center of the profile.  The median
+minimizes the chance of a large error due to a single bad pixel.  This
+scale may be greatly in error in the case of extreme blending but is also
+quite fast; the extraction time is reduced by at least 40%.
+The steps of profile fitting and deviant pixel detection are alternated
+and the maximum number of iterations through these two steps is
+set by \fIniterate\fR.  The default of 1 means that the model fitting is not
+repeated after detecting deviant pixels.
+
+The option \fIverbose\fR can be used to print the image lines being extracted
+and any pixels replaced by the cleaning process.
+.ih
+EXAMPLES
+To create a cleaned version of the image using model \fIsmooth\fR for cleaning:
+
+	cl> newimage image newimage
+
+To create an model image using model \fIgauss5\fR:
+
+	cl> newimage image newimage ex_model=yes model="gauss5"
+.ih
+SEE ALSO
+msextract
+.endhelp