Initial commit

10 files changed, 2423 insertions, 0 deletions

diff --git a/noao/twodspec/longslit/doc/extinction.hlp b/noao/twodspec/longslit/doc/extinction.hlp
new file mode 100644
index 00000000..39579a07
--- /dev/null
+++ b/noao/twodspec/longslit/doc/extinction.hlp
@@ -0,0 +1,87 @@
+.help extinction May87 noao.twodspec.longslit
+.ih
+NAME
+extinction -- Apply atmospheric extinction corrections
+.ih
+USAGE
+extinction images
+.ih
+PARAMETERS
+.ls input
+List of input images to be extinction corrected.
+.le
+.ls output
+List of output extinction corrected images.  Output images may be the
+same as the input images.
+.le
+.ls extinction = "onedstds$kpnoextinct.dat"
+Extinction file to be used.  The standard extinction files:
+
+.nf
+	onedstds$kpnoextinct.dat - KPNO standard extinction
+	onedstds$ctioextinct.dat - CTIO standard extinction
+.fi
+.le
+.ih
+DESCRIPTION
+The specified images are corrected for atmospheric extinction according
+to the formula
+
+    correction factor = 10 ** (0.4 * airmass * extinction)
+
+where the extinction is a tabulated function of the wavelength.  The
+extinction file contains lines of wavelength and extinction at that
+wavelength.  The units of the wavelength must be the same as those of
+the dispersion corrected images; i.e. Angstroms.   If the image is
+dispersion corrected in logarithmic wavelength intervals (DC-FLAG = 1)
+the task will convert to wavelength and so the extinction file must
+still be wavelength.  The table values are interpolated
+to the wavelengths of the image pixels and the correction applied to
+the pixel values.  Note that the image pixel values are modifed.
+
+The airmass is sought in the image header under the name AIRMASS.  If the
+airmass is not found then it is computed from the zenith distance (ZD in hours)
+using the approximation formula from Allen's "Astrophysical Quantities", 1973,
+page125 and page 133
+
+	AIRMASS = sqrt (cos (ZD) ** 2 + 2 * scale + 1)
+
+where the atmospheric scale height is set to be 750.  If the parameter ZD
+is not found then it must be computed from the hour angle (HA in hours),
+the declination (DEC in degrees), and the observation latitude (LATITUDE
+in degress).   The hour angle may be computed from the right ascension
+(RA in hours) and siderial time (ST in hours).  Computed quantities are
+recorded in the image header.  Flags indicating extinction correction are
+also set in the image header.
+
+The image header keyword DISPAXIS must be present with a value of 1 for
+dispersion parallel to the lines (varying with the column coordinate) or 2
+for dispersion parallel to the columns (varying with line coordinate).
+This parameter may be added using \fBhedit\fR.  Note that if the image has
+been transposed (\fBimtranspose\fR) the dispersion axis should still refer
+to the original dispersion axis unless the physical world coordinate system
+is first reset (see \fBwcsreset\R).  This is done in order to allow images
+which have DISPAXIS defined prior to transposing to still work correctly
+without requiring this keyword to be changed.
+.ih
+EXAMPLES
+1. A set of dispersion corrected images is extinction corrected in-place as
+follows:
+
+.nf
+	cl> extinction img* img*
+.fi
+
+2. To keep the uncorrected image:
+
+.nf
+	cl> extinction nite1.004 nite1ext.004
+.fi
+
+3.  If the DISPAXIS keyword is missing and the dispersion is running
+vertically (varying with the image lines):
+
+.nf
+	cl> hedit *.imh dispaxis 2 add+
+.fi
+.endhelp
diff --git a/noao/twodspec/longslit/doc/fccoeffs b/noao/twodspec/longslit/doc/fccoeffs
new file mode 100644
index 00000000..ab8de92f
--- /dev/null
+++ b/noao/twodspec/longslit/doc/fccoeffs
@@ -0,0 +1,210 @@
+From davis Tue May 18 15:09:59 1993
+Received: by tucana.tuc.noao.edu (4.1/SAG.tucana.12)
+        id AA26431; Tue, 18 May 93 15:09:56 MST; for sites
+Date: Tue, 18 May 93 15:09:56 MST
+From: davis (Lindsey Davis)
+Message-Id: <9305182209.AA26431@tucana.tuc.noao.edu>
+To: belkine@mesiob.obspm.circe.fr
+Subject: RE: geomap
+Cc: sites
+
+
+
+Igor,
+
+    The following is a copy of a mail message I sent to another user who made
+the same request regarding geomap. I hope this is of use to you.
+
+
+                                           Lindsey Davis
+
+###############################################################################
+
+
+    Jeannette forwarded your request for a detailed description of the
+geomap output format to me. This format was originally intended to be
+for the internal use of geomap only, but the following should help you
+decode it.
+
+    1. For simple linear geometric transformations you will see the
+following two entries in the fit record. Surface1 describes the linear
+portion of the fit; surface2 describes the residual distortion map
+which is always 0 for linear fits.
+
+        surface1        11
+            surface(xfit) surface(yfit)    (surface type 1=cheb, 2=leg, 3=poly)
+            xxorder(xfit) yxorder(yfit)    (always 2)
+            xyorder(xfit) yyorder(yfit)    (always 2)
+            xxterms(xfit) yxterms(yfit)    (always 0)
+            xmin(xfit)   xmin(yfit)        (geomap input or data)
+            xmax(xfit)   xmax(yfit)        (geomap input or data)
+            ymin(xfit)   ymin(yfit)        (geomap input or data)
+            ymax(xfit)   ymax(yfit)        (geomap input or data)
+            a            d
+            b            e
+            c            f
+        surface2        0
+
+This above describes the following linear surfaces.
+
+        xfit = a + b * x + c * y  (polynomial)
+        yfit = d + e * x + f * y
+
+        xfit = a + b * xnorm + c * ynorm  (chebyshev)
+        yfit = d + e * xnorm + f * ynorm
+
+        xfit = a + b * xnorm + c * ynorm  (legendre)
+        yfit = d + e * xnorm + f * ynorm
+
+        xnorm = (2 * x - (xmax + xmin)) / (xmax - xmin)
+        ynorm = (2 * y - (ymax + ymin)) / (ymax - ymin)
+
+Xnorm and ynorm are the input x and y values normalized between -1.0
+and 1.0.
+
+
+
+
+      2. For a higher order fit, say xorder=4 yorder=4 and xterms=yes,
+the format is more complicated.  The second surface is computed by fitting
+the higher order surface to the residuals of the first fit. The geomap
+output will look something like the following.
+
+        surface1        11
+            surface(xfit) surface(yfit)    (surface type 1=cheb, 2=leg, 3=poly)
+            xxorder(xfit) yxorder(yfit)    (always 2)
+            xyorder(xfit) yyorder(yfit)    (always 2)
+            xxterms(xfit) yxterms(yfit)    (always 0)
+            xmin(xfit)   xmin(yfit)        (geomap input or data)
+            xmax(xfit)   xmax(yfit)        (geomap input or data)
+            ymin(xfit)   ymin(yfit)        (geomap input or data)
+            ymax(xfit)   ymax(yfit)        (geomap input or data)
+            a            d
+            b            e
+            c            f
+        surface2        24
+            surface(xfit) surface(yfit)    (surface type 1=cheb, 2=leg, 3=poly)
+            xxorder(xfit) yxorder(yfit)    (4)
+            xyorder(xfit) yyorder(yfit)    (4)
+            xxterms(xfit) yxterms(yfit)    (1 in this case)
+            xmin(xfit)   xmin(yfit)        (geomap input or data)
+            xmax(xfit)   xmax(yfit)        (geomap input or data)
+            ymin(xfit)   ymin(yfit)        (geomap input or data)
+            ymax(xfit)   ymax(yfit)        (geomap input or data)
+            C00(xfit)    C00(yfit)
+            C10(xfit)    C10(yfit)
+            C20(xfit)    C20(yfit)
+            C30(xfit)    C30(yfit)
+            C01(xfit)    C01(yfit)
+            C11(xfit)    C11(yfit)
+            C21(xfit)    C21(yfit)
+            C31(xfit)    C31(yfit)
+            C02(xfit)    C02(yfit)
+            C12(xfit)    C12(yfit)
+            C22(xfit)    C22(yfit)
+            C32(xfit)    C32(yfit)
+            C03(xfit)    C03(yfit)
+            C13(xfit)    C13(yfit)
+            C23(xfit)    C23(yfit)
+            C33(xfit)    C33(yfit)
+
+
+where the Cmn are the coefficients of the polynomials Pmn, and the Pmn
+are defined as follows
+
+            Pmn = x ** m * y ** n    (polynomial)
+
+            Pmn = Pm(xnorm) * Pn(ynorm)     (chebyshev)
+
+                P0(xnorm) = 1.0
+                P1(xnorm) = xnorm
+                Pm+1(xnorm) = 2.0 * xnorm * Pm(xnorm) - Pm-1(xnorm)
+                xnorm = (2 * x - (xmax + xmin)) / (xmax - xmin)
+
+                P0(ynorm) = 1.0
+                P1(ynorm) = ynorm
+                Pn+1(ynorm) = 2.0 * ynorm * Pn(ynorm) - Pn-1(ynorm)
+                ynorm = (2 * y - (ymax + ymin)) / (ymax - ymin)
+
+            Pmn = Pm(xnorm) * Pn(ynorm)     (legendgre)
+
+                P0(xnorm) = 1.0
+                P1(xnorm) = xnorm
+                Pm+1(xnorm) = ((2m + 1) * xnorm * Pm(xnorm) - m * Pm-1(xnorm))/
+                              (m + 1)
+                xnorm = (2 * x - (xmax + xmin)) / (xmax - xmin)
+
+                P0(ynorm) = 1.0
+                P1(ynorm) = ynorm
+                Pn+1(ynorm) = ((2n + 1) * ynorm * Pn(ynorm) - n * Pn-1(ynorm))/
+                              (n + 1)
+                ynorm = (2 * y - (ymax + ymin)) / (ymax - ymin)
+
+
+Hopefully I have copied this all down correctly. The main points to remember
+is that the mangitudes of the coefficients reflect both the function type
+(polynomial, chebyshev, or legendre) and the normalization (xmin, xmax,
+ymin, ymax).
+
+    Hope this helps you out and write back if you have more questions.
+
+                                             Lindsey Davis
+
+=======================================
+
+# <Date>
+begin	<name>
+	task fitcoords
+	axis	1	# Axis of fitted value
+	surface	24	# The number of following parameters/coefficients
+            surface	# surface type 1=chebyshev, 2=legendre
+            xorder	# X order
+            yorder	# Y order
+            xterms	# Cross terms?  0=no, 1=yes (always 1 for fitcoords)
+            xmin	# Minimum x value in fit - usually 1
+            xmax	# Maximum x value in fit - usually image dimension
+            ymin	# Minimum y value in fit - usually 1
+            ymax	# Maximum y value in fit - usually image dimension
+            C00		# Coefficients (shown for xorder=4 and yorder=4)
+            C10
+            C20
+            C30
+            C01
+            C11
+            C21
+            C31
+            C02
+            C12
+            C22
+            C32
+            C03
+            C13
+            C23
+            C33
+
+
+The fit is a sum of the form:
+
+	fit = sum(m=0 to xorder-1) sum(n=0 to yorder-1) {Cmn*Pm(x')*Pn(y')}
+
+where the cross-terms may or may not be included depending on the xterms
+parameter.  Cross-terms are always used in FITCOORDS.
+
+The coefficients are defined in terms of normalized independent variables
+in the range -1 to 1.  If x and y are actual values then the normalized
+variables, x' and y', are defined using the data range parameters as:
+
+	x' = (2 * x - (xmax + xmin)) / (xmax - xmin)
+	y' = (2 * y - (ymax + ymin)) / (ymax - ymin)
+
+The Pi(z), where z is either x' or y', are defined iteratively as follows:
+
+	# Chebyshev
+	P0(z) = 1.0
+	P1(z) = z
+	Pi+1(z) = 2.0 * z * Pi(z) - Pi-1(z)
+
+	# Legendre
+	P0(z) = 1.0
+	P1(z) = z
+	Pi+1(z) = ((2i + 1) * z * Pi(z) - i * Pi-1(z)) / (i + 1)
diff --git a/noao/twodspec/longslit/doc/fceval.hlp b/noao/twodspec/longslit/doc/fceval.hlp
new file mode 100644
index 00000000..87d258c0
--- /dev/null
+++ b/noao/twodspec/longslit/doc/fceval.hlp
@@ -0,0 +1,87 @@
+.help fceval Aug03 noao.twodspec.longslit
+.ih
+NAME
+fceval -- Evaluate coordinates using the FITCOORDS solutions
+.ih
+USAGE
+fceval input output fitnames
+.ih
+PARAMETERS
+.ls input
+Input text file of pixel coordinates.  This may be "STDIN" to read
+coordinates from the terminal or pipe.
+.le
+.ls output
+Output text file of pixel coordinates and fitted coordinates.  This may
+be "STDOUT" to write coordinates to the terminal or pipe.
+.le
+.ls fitnames  
+Names of the user coordinate maps to evaluate.
+.le
+.ls database = "database"
+Database containing the coordinate maps.
+.le
+.ih
+DESCRIPTION
+This task transforms pixel coordinates to the world coordinates fit with
+FITCOORDS.  When there is no map for an axis the identify transform is
+used.  If there are more the one map for an axis the average of the mapped
+coordinates is output.  This is the same behavior as TRANSFORM.
+
+The input file consists of two columns giving the x and y pixel values
+in the frame of the untransformed image data.  The output is a file
+with four columns giving the input x any y pixel values and the
+user coordinates fit by FITCOORDS.
+
+Two typical uses for this task are to look up world coordinates for
+points in the untransformed data and to generate transformations using
+GEOMAP and GEOTRAN.
+.ih
+EXAMPLES
+1. Evaluate a wavelength and slit position fit where the input pixel coordinates
+are entered interactively and the output is written to the terminal.
+
+.nf
+    cl> fceval STDIN STDOUT arcfit,std
+    1 1
+    1. 1. 20.60425149463117 4202.47202514205
+    60 1
+    60. 1. 79.60425149463118 4203.316616448186
+    1 512
+    1. 512. 19.15606081299484 7356.089801036373
+    60 512
+    60. 512. 78.15606081299485 7355.042495319318
+.fi
+
+In this case the first axis corresponds to the spatial dimension and
+the second to the dispersion dimension.  The arcfit was created using
+Angstroms and so the units of the last column is Angstroms.
+
+2. One use of this task is to generate the inverse transformation from
+that produced by TRANSFORM.  The steps are: 1) produce a grid of
+coordinates using LISTPIX and FCEVAL, 2) convert the user coordinates to
+pixel coordinates in the transformed data using WCSCTRAN, 3) fit a
+transformation using GEOMAP, and 4) transform the data with GEOTRAN.
+
+.nf
+    cl> listpix orig[*:5,*:5] wcs=physical verb- |
+    >>> fceval STDIN STDOUT arcfit,std |
+    >>> wcsctran STDIN coords trans world logical columns="3 4"
+    cl> geomap coords geomap.db 1 61 1 512
+    cl> geotran trans origNEW geomap.db coords flux+
+.fi
+
+This example uses pipes to eliminate intermediate files.  But these
+files can be useful for understanding the process.  LIXTPIX is used to
+generate a grid of points with some subsampling.  Be sure to use "physical"
+for the coordinate system otherwise the grid of x and y values will be
+for the subsection.  The order of the columns will be appropriate for
+GEOMAP to compute the inverse transformation.  By reversing the order
+of the columns one could generate a transformation similar to that
+produced by TRANSFORM in order to use features in GEOTRAN not provided
+by TRANSFORM.  However, the world coordinate system information will
+not be automatically set.
+.ih
+SEE ALSO
+fitcoords, transform, geomap, geotran
+.endhelp
diff --git a/noao/twodspec/longslit/doc/fitcoords.hlp b/noao/twodspec/longslit/doc/fitcoords.hlp
new file mode 100644
index 00000000..a376ee74
--- /dev/null
+++ b/noao/twodspec/longslit/doc/fitcoords.hlp
@@ -0,0 +1,287 @@
+.help fitcoords Apr00 noao.twodspec.longslit
+.ih
+NAME
+fitcoords -- Fit user coordinates to the image coordinates
+.ih
+USAGE
+fitcoords images fitname
+.ih
+PARAMETERS
+.ls images
+List of images containing the feature coordinates to be fit.  If the
+parameter \fIcombine\fR is yes then feature coordinates from all the images
+are combined and fit by a single function.  Otherwise the feature coordinates
+from each image are fit separately.
+.le
+.ls fitname = "" 
+If the input images are combined and fit by a single function then the fit
+is stored under this name.  If the images are not combined then the
+fit for each image is stored under the name formed by appending the image
+name to this name.  A null prefix is acceptable when not combining but it
+is an error if combining a list of images.
+.le
+.ls interactive = yes
+Determine coordinate fits interactively?
+.le
+.ls combine = no
+Combine the coordinates from all the input images and fit them by a single
+function?  If 'no' then fit the coordinates from each image separately.
+.le
+.ls database = "database"
+Database containing the feature coordinate information used in fitting the
+coordinates and in which the coordinate fit is recorded.
+.le
+.ls deletions = "deletions.db"
+Deletion list file.  If not null then points whose coordinates match those in
+this file (if it exists) are initially deleted from the fit.
+If the fitting is done interactively then the coordinates of
+any deleted points (after exiting from the interactive fitting) are recorded
+in this file.
+.le
+.ls function = "chebyshev"
+Type of two dimensional function to use in fitting the user coordinates.
+The choices are "chebyshev" polynomial and "legendre" polynomial.
+The function may be abbreviated.  If the task is interactive then
+the user may change the function later.
+.le
+.ls xorder = 6
+Order of the mapping function along the first image axis.
+The order is the number of polynomial terms.  If the task is interactive
+then the user may change the order later.
+.le
+.ls yorder = 6
+Order of the mapping function along the second image axis.
+The order is the number of polynomial terms.  If the task is interactive
+then the user may change the order later.
+.le
+.ls logfiles = "STDOUT,logfile"
+List of files in which to keep logs containing information about
+the coordinate fit.  If null then no log is kept.
+.le
+.ls plotfile = "plotfile"
+Name of file to contain metacode for log plots.  If null then no log plots
+are kept.  When the fitting is interactive the last graph is recorded in
+the plot file and when not interactive a default plot is recorded.
+.le
+.ls graphics = "stdgraph"
+Graphics output device.
+.le
+.ls cursor = ""
+Graphics cursor input.  If null the standard graphics cursor is used.
+.le
+.bp
+.ih
+CURSOR COMMANDS
+
+.nf
+?  List commands
+c  Print data values for point nearest the cursor
+d  Delete the point or set of points with constant x, y, or z
+	nearest the cursor (p, x, y, z,)
+f  Fit surface
+l  Graph the last set of points (in zoom mode)
+n  Graph the next set of points (in zoom mode)
+p  Graph all features
+q  Quit
+r  Redraw a graph
+u  Undelete the point or set of points with constant x, y, or z
+	nearest the cursor (p, x, y, z,)
+w  Window the graph.  Type '?' to the "window:" prompt for more help.
+x  Select data for the x axis (x, y, z, s, r)
+y  Select data for the y axis (x, y, z, s, r)
+z  Zoom on the set of points with constant x, y, or z (x, y, z)
+   Unzoom with p
+
+:corners	Show the fitted values for the corners of the image
+:function type	Set the function for the fitted surface
+		(chebyshev, legendre)
+:show		Show the fitting parameters
+:xorder value	Set the x order  for the fitted surface
+:yorder value	Set the y order  for the fitted surface
+.fi
+.ih
+DESCRIPTION
+A two dimensional function of the image coordinates is fitted to the user
+coordinates from the specified images;
+
+.nf
+	user coordinate = function (column, line)
+
+			or
+
+		      z = s (x, y)
+.fi
+
+The coordinates from all the input images may be combined in a single fit or
+the coordinates from each image may be fit separately.  If the
+coordinates from the input images are combined then the fitted function
+is recorded in the database under the specified name.  If
+the coordinates are fit separately the fitted function is recorded under
+a name formed by appending the image name to the specified root name.
+
+When the task is interactive the user is first queried whether to perform
+the fitting interactively.  The user may answer "yes", "no", "YES", or "NO"
+to the query.  The lowercase responses apply only to the current fit
+and the uppercase responses apply to all remaining fits.  When the
+fitting is done interactively the user may change the fitted function and
+orders iteratively, delete individual coordinates or entire features,
+and graph the fit and residuals in a number ways.
+The CURSOR COMMANDS section describes the graphics cursor keystrokes
+which are available.  When selecting data for the graph axes the
+follow definitions apply:
+
+.nf
+	x	Input image column positions
+	y	Input image line positions
+	z	Input user coordinates
+	s	Fitted user coordinates
+	r	Residuals (s - z)
+.fi
+
+A very useful feature is zooming, deleting, or undeleting a subset of data
+points.  The subsets
+are defined as points with the same x, y, or z value as the point indicated
+by the cursor when typing (z)oom, (d)elete, or (u)ndelete.
+
+When a satisfactory coordinate fit has been determined exit with the (q)uit
+key.  The user is asked if the fit is to be recorded in the database.
+
+If a deletion list file is specified then the coordinates of any
+points deleted interactively are recorded in this file.  This file then can
+be read by subsequent fits to initially delete points with matching
+coordinates.  This is generally used when fitting a series of images
+non-interactively.
+
+Information about the fitted function may be recorded.  Textual information
+is written to the specified log files (which may include the standard
+output STDOUT).  The last interactive plot or a default non-interactive
+plot is written the specified plot file which may be examined and spooled
+at a later time.
+
+
+FITCOORDS DATABASE
+
+The FITCOORDS fits are stored in text files in the subdirectory given by
+the "database" parameter.  The name of the file is fc<fitname> where
+<fitname> is the specified fit name.  The database text file contains
+blocks of lines beginning with a time stamp followed by line with the
+"begin" keyword.  The value following "begin" is the fit name, which is
+often the name of the image used for the fit.  If there is more than one
+block with the same fit name then the last one is used.
+
+The "task" keyword will has the value "fitcoords" and the "axis" keyword
+identifies the axis to which the surface fit applies.  An axis of 1 refers
+to the first or x axis (the first dimension of the image) and 2 refers to
+the second or y axis.
+
+The "surface" keyword specifies the number of coefficients for the surface
+fit given in the following lines .  The surface fit is produced by an IRAF
+math package called "gsurfit".  The coefficients recorded in the database
+are intented to be internal to that package.  However the following
+describes how to interpret the coefficients.
+
+The first 8 lines specify:
+
+.nf
+   function - Function type (1=chebyshev, 2=legendre)
+     xorder - X "order" (highest power of x)
+     yorder - Y "order" (highest power of y)
+     xterms - Cross-term type (always 1 for FITCOORDS)
+       xmin - Minimum x over which the fit is defined
+       xmax - Maximum x over which the fit is defined
+       ymin - Minimum y over which the fit is defined
+       ymax - Maximum y over which the fit is defined
+.fi
+
+The polynomial coefficients follow in array order with the x index
+varying fastest:
+
+.nf
+	C00
+	C10
+	C20
+	...
+	C<xorder-1>0
+	C01
+	C11
+	C21
+	...
+	C<xorder-1>1
+	...
+	C<xorder-1><yorder-1>
+.fi
+
+The surface fitting functions have the form
+
+.nf
+	fit(x,y) = Cmn * Pmn
+.fi
+
+where the Cmn are the coefficients of the polynomials terms Pmn, and the Pmn
+are defined as follows:
+
+.nf
+Chebyshev: Pmn = Pm(xnorm) * Pn(ynorm)
+
+	   xnorm = (2 * x - (xmax + xmin)) / (xmax - xmin)
+	   ynorm = (2 * y - (ymax + ymin)) / (ymax - ymin)
+
+	   P0(xnorm) = 1.0
+	   P1(xnorm) = xnorm
+	   Pm+1(xnorm) = 2.0 * xnorm * Pm(xnorm) - Pm-1(xnorm) 
+
+	   P0(ynorm) = 1.0
+	   P1(ynorm) = ynorm
+	   Pn+1(ynorm) = 2.0 * ynorm * Pn(ynorm) - Pn-1(ynorm) 
+
+Legendre:  Pmn = Pm(xnorm) * Pn(ynorm)
+
+	   xnorm = (2 * x - (xmax + xmin)) / (xmax - xmin)
+	   ynorm = (2 * y - (ymax + ymin)) / (ymax - ymin)
+
+	   P0(xnorm) = 1.0
+	   P1(xnorm) = xnorm
+	   Pm+1(xnorm) = ((2m+1)*xnorm*Pm(xnorm)-m*Pm-1(xnorm))/(m+1)   
+
+	   P0(ynorm) = 1.0
+	   P1(ynorm) = ynorm
+	   Pn+1(ynorm) = ((2n+1)*ynorm*Pn(ynorm)-n*Pn-1(ynorm))/(n+1)   
+.fi
+
+Notice that the x and y values are first normalized to the interval -1 to 1
+over the range of the surface as given by the xmin, xmax, ymin, and ymax
+elements of the database description.
+.ih
+EXAMPLES
+A number of strong arc lines are identified along one column of an arc
+calibration image "arc001".  The arc lines are then reidentified at every
+20th column.  A two dimensional dispersion solution is determined as follows:
+
+	cl> fitcoords arc001 fit.
+
+The fitting is done interactively and deleted points are recorded.
+The fit is recorded under the name fit.arc001.  A set of similar arc
+calibrations are fit non-interactively, with the same points deleted,
+as follows:
+
+	cl> fitcoords arc* interactive=no
+
+Several stellar spectra are identified at different positions along the slit
+and traced to other lines.  A fit to the geometric distortion is determined
+with the command:
+
+	cl> fitcoords star001,star003,star005 fitname=distortion combine=yes
+
+In this case the coordinates from all the tracings are combined in a single
+fit called distortion.
+
+The plots in the plot file are spooled to the standard plotting device as
+follows:
+
+	cl> gkimosaic plotfile
+
+\fBGkimosaic\fR is in the \fBplot\fR package.
+.ih
+SEE ALSO
+transform
+.endhelp
diff --git a/noao/twodspec/longslit/doc/fluxcalib.hlp b/noao/twodspec/longslit/doc/fluxcalib.hlp
new file mode 100644
index 00000000..ee38cee5
--- /dev/null
+++ b/noao/twodspec/longslit/doc/fluxcalib.hlp
@@ -0,0 +1,106 @@
+.help fluxcalib Oct86 noao.twodspec.longslit
+.ih
+NAME
+fluxcalib -- Apply flux calibration
+.ih
+USAGE
+fluxcalib images fluxfile
+.ih
+PARAMETERS
+.ls input
+List of input images to be flux calibrated.
+.le
+.ls output
+List of output flux calibrated images.  The output images may be the same
+as the input images.  The output image will be of type real regardless
+of the input pixel type.
+.le
+.ls fluxfile
+Flux calibration file from \fBonedspec.sensfunc\fR.
+.le
+.ls fnu = no
+Convert the flux calibration to flux per unit frequency (F-nu)?
+.le
+.ls exposure = "otime"
+Exposure time keyword in image headers.
+.le
+.ih
+DESCRIPTION
+The specified images are flux calibrated using a flux calibration image
+file derived from the \fBonedspec\fR package using standard stars.
+The flux calibration pixel values are in magnitudes and the pixel coordinates
+are in wavelength.  The multiplicative calibration factor is given by the
+formula
+
+    factor = 10 ** (-0.4 * calibration) / exposure / dispersion.
+
+Since the calibration data has units of (instrumental intensity) /
+(ergs/cm**2), the exposure time for the image must be in seconds and the
+pixel dispersion in wavelength/pixel to yield units of
+ergs/cm**2/sec/wavelength.
+
+The calibration wavelengths are interpolated to the wavelengths
+of the image pixels and the correction applied to the pixel values.
+Note that the image pixel values are modified.
+
+If flux per unit frequency is requested then the flux values are multiplied
+by
+
+	wavelength ** 2 / velocity of light (in Angstroms/sec)
+
+to yield units of ergs/cm**2/Hz/sec/(wavelength/Angstrom).  Note that normally
+the wavelength units should be Angstroms.
+
+It is possible to flux calibrate images which are binned in logarithmic
+wavelength intervals.  The point to note is that the units of the flux
+calibrated image will be the same.  Therefore, rebinning to linear
+wavelength coordinates requires only interpolation and not flux conservation.
+When extracting standard stars from logarithmicaly bin spectra for determination
+of a flux calibration it is necessary to rebin the extracted one dimensional
+spectra to linear wavelength (required by \fBonedspec\fR) conserving
+flux so that the instrumental counts are preserved.
+
+The image header keyword DISPAXIS must be present with a value of 1 for
+dispersion parallel to the lines (varying with the column coordinate) or 2
+for dispersion parallel to the columns (varying with line coordinate).
+This parameter may be added using \fBhedit\fR.  Note that if the image has
+been transposed (\fBimtranspose\fR) the dispersion axis should still refer
+to the original dispersion axis unless the physical world coordinate system
+is first reset (see \fBwcsreset\R).  This is done in order to allow images
+which have DISPAXIS defined prior to transposing to still work correctly
+without requiring this keyword to be changed.
+.ih
+EXAMPLES
+Standard stars were observed and extracted to one dimensional spectra.
+The standard stars are then used to determine a flux calibration using
+the \fBonedspec\fR package.  A set of dispersion and extinction corrected
+images is flux calibrated in-place with the command
+
+.nf
+	cl> fluxcalib img* img* sens.0000
+.fi
+
+where "sens.0000" is the calibration file produced by the task
+\fBonedspec.sensfunc\fR.
+
+To keep the uncalibrated image:
+
+.nf
+	cl> fluxcalib n1ext.004 n1extf.004 sens.0000
+.fi
+
+3.  If the DISPAXIS keyword is missing and the dispersion is running
+vertically (varying with the image lines):
+
+.nf
+	cl> hedit *.imh dispaxis 2 add+
+.fi
+.ih
+REVISIONS
+.ls FLUXCALIB V2.10
+The output pixel type is now forced to be real.
+.le
+.ih
+SEE ALSO
+onedspec.standard onedspec.sensfunc
+.endhelp
diff --git a/noao/twodspec/longslit/doc/illumination.hlp b/noao/twodspec/longslit/doc/illumination.hlp
new file mode 100644
index 00000000..5697bfad
--- /dev/null
+++ b/noao/twodspec/longslit/doc/illumination.hlp
@@ -0,0 +1,220 @@
+.help iillumination Jul86 noao.twodspec.longslit
+.ih
+NAME
+iillumination -- Determine iillumination calibrations
+.ih
+USAGE
+iillumination images iilluminations
+.ih
+PARAMETERS
+.ls images
+Images to use in determining iillumination calibrations.  These are
+generally sky spectra.  An image section may be used to select only a
+portion of the image.
+.le
+.ls iilluminations
+Iillumination calibration images to be created.  Each iillumination image is
+paired with a calibration image.  If the image exists then it will be modified
+otherwise it is created.
+.le
+.ls interactive = yes
+Graph the average spectrum and select the dispersion bins
+and graph and fit the slit profile for each dispersion bin interactively?
+.le
+.ls bins = ""
+Range string defining the dispersions bins within which the slit profiles
+are determined.  If the range string is null then the dispersion
+bins are determined by the parameter \fInbins\fR.
+.le
+.ls nbins = 5
+If the dispersion bins are not specified explicitly by the parameter
+\fIbins\fR then the dispersion range is divided into this number of
+nearly equal bins.
+.le
+.ls sample = "*"
+Sample of points to use in fitting each slit profile.
+The sample is selected with a range string.
+.le
+.ls naverage = 1
+Number of sample points to average or median before fitting a function.
+If the number is positive the average of each set of naverage sample
+points is formed while if the number is negative then the median of each set
+of points (in absolute value) is formed.  This subsample of points is
+used in fitting the slit profile.
+.le
+.ls function = "spline3"
+Function to fit to each dispersion bin to form the iillumination function.
+The options are "spline1", "spline3", "legendre", and "chebyshev".
+.le
+.ls order = 1
+Order of the fitting function or the number of spline pieces.
+.le
+.ls low_reject = 0., high_reject = 0.
+Rejection limits below and above the fit in units of the residual sigma.
+.le
+.ls niterate = 1
+Number of rejection iterations.
+.le
+.ls grow = 0
+Reject additional points within this distance of points exceeding the
+rejection threshold.
+.le
+.ls interpolator = "poly3"
+Interpolation type.  One of "nearest", "linear", "poly3", "poly5", or
+"spline3".
+.le
+.ls graphics = "stdgraph"
+Graphics output device.  May be one of the standard devices "stdgraph",
+"stdplot", or "stdvdm" or an explicit device.
+.le
+.ls cursor = ""
+Graphics input device.  May be either null for the standard graphics cursor
+or a file containing cursor commands.
+.le
+.ih
+CURSOR KEYS
+The interactive curve fitting package \fBicfit\fR is used to fit a function
+to the average calibration spectrum.  Additional help on using this package
+and the cursor keys is available under the name "icfit".
+
+When the dispersion bins are set graphically the following cursor keys are
+defined.
+
+.ls ?
+Clear the screen and print a menu of the cursor options.
+.le
+.ls i
+Initialize the sample ranges.
+.le
+.ls q
+Exit interactive dispersion bin selection.
+.le
+.ls s
+Set a bin with the cursor.  This may be repeated any number of times.
+Two keystrokes are required to mark the two ends of the bin.
+.le
+
+The parameters are listed or set with the following commands which may be
+abbreviated.  To list the value of a parameter type the command alone.
+
+.nf
+:bins value		Iillumination bins
+:show			Show the values of all the parameters
+.fi
+.ih
+DESCRIPTION
+An iillumination calibration, in the form of an image, is created for each
+longslit calibration image, normally a sky spectrum.  The iillumination
+calibration is determined by fitting functions across the slit (the slit
+profiles) at a number of points along the dispersion, normalizing each fitted
+function to unity at the center of the slit, and interpolating the iillumination
+between the dispersion points.  The fitted data is formed by dividing the
+dispersion points into a set of bins and averaging the slit profiles within
+each bin.  The interpolation type is a user parameter.
+
+The image header keyword DISPAXIS must be present with a value of 1 for
+dispersion parallel to the lines (varying with the column coordinate) or 2
+for dispersion parallel to the columns (varying with line coordinate).
+This parameter may be added using \fBhedit\fR.  Note that if the image has
+been transposed (\fBimtranspose\fR) the dispersion axis should still refer
+to the original dispersion axis unless the physical world coordinate system
+is first reset (see \fBwcsreset\fR).  This is done in order to allow images
+which have DISPAXIS defined prior to transposing to still work correctly
+without requiring this keyword to be changed.
+
+If the output image does not exist it is first created with unit iillumination
+everywhere.  Subsequently the iillumination is only modified in those regions
+occupied by the input image.  Thus, an image section in the input image may
+be used to select the data to be used and for which an iillumination calibration
+will be determined.  This ability is particularly userful when dealing with
+multiple slits or to exclude regions outside the slit.
+
+The dispersion bins may be selected by a range string (\fIbins\fR) or,
+if no range string is given, by the number of bins into which the dispersion
+range is to be divided (\fInbins\fR).  When the interactive parameter
+is set (\fIinteractive\fR) then the average spectrum is graphed and the
+bins may be set using the cursor or with a colon command.  Once the bins
+have been selected exit with (q)uit to continue to the slit profile fitting.
+
+Fitting of the slit profiles is done using the interactive curve fitting
+package (\fBicfit\fR).  The parameters determining the fit are the
+sample points, the averaging bin size, the fitting function,
+the order of the function, the rejection sigmas, the number of
+rejection iterations, and the rejection width.
+The sample points for the average slit profile are selected by a range string.  
+Points in the slit profile not in the sample are not used in determining
+the fitted function.  The selected sample points may be binned into a
+set of averages or medians which are used in the function fit instead of the
+sample points with the averaging bin size parameter
+\fInaverage\fR.  This parameter selects the number of sample points to be
+averaged if its value is positive or the number of points to be medianed
+if its value is negative (naturally, the absolute value is used for the
+number of points).  A value of one uses all sample points without binning.
+The fitted function may be used to reject points from the fit using the
+parameters \fIlow_reject, high_reject, niterate\fR and \fIgrow\fR.  If
+one or both of the rejection limits are greater than zero then the sigma
+of the residuals is computed and points with residuals less than
+\fI-low_reject\fR times the sigma and greater than \fIhigh_reject\fR times
+the sigma are removed and the function fitted again.  In addition points
+within a distance given by the parameter \fIgrow\fR of the a rejected point
+are also rejected.  A value of zero for this parameter rejects only the
+points exceeding the rejection threshold.  Finally, the rejection procedure
+may be iterated the number of times given by the parameter \fIniterate\fR.
+
+The fitted functions may be examined and modified interactively when the
+parameter \fIinteractive\fR is set.  The user is asked before each dispersion
+bin whether to perform the fit interactively.  The possible response are
+"no", "yes", "NO", and "YES".  The lower case responses only affect the
+specified dispersion bin while the upper case responses affect all following
+dispersion bins for the current image.  Thus, if the response is "NO" then
+no further prompts or interactive curve fitting need be performed while if
+the response is "YES" there are no further prompts but the slit profile
+for each dispersion bin must be graphed and exited with (q)uit.
+Changes to the fitting parameters remain in effect until they are next
+changed.  This allows the fitting parameters to be selected from only the first
+dispersion bin without requiring each dispersion bin to be graphed and
+confirmed.
+
+When a dispersion bin is to be fitted interactively the average slit profile
+and the fitted function or the residuals of the fit are graphed.
+Deleted points are marked with an x and rejected points by a diamond.
+The sample regions are indicated along the bottom of the graph.
+The cursor keys and colon commands are used to change the values
+of the fitting parameters, delete points, and window and expand the
+graph.  When the fitted function is satisfactory exit with
+with a carriage return or 'q'.  The prompt for the next dispersion bin will
+then be given until the last dispersion bin has been fit.  The iillumination
+calibration image is then created.
+.ih
+EXAMPLES
+1. To create an iillumination image non-interactively:
+
+.nf
+	cl> iillumination sky illum nbins=8 order=20 interactive=no
+.fi
+
+2. To determine independent iilluminations for a multislit image determine the
+image sections defining each slit.  Then the iillumination functions are
+computed as follows:
+
+.nf
+	cl> iillumination sky[10:20,*],sky[35:45,*] illum,illum
+.fi
+
+3. Generally the slit image sections are prepared in a file which is then
+used to define the lists of input images and iilluminations.
+
+.nf
+	cl> iillumination @slits @illums
+.fi
+
+3.  If the DISPAXIS keyword is missing and the dispersion is running
+vertically (varying with the image lines):
+
+.nf
+	cl> hedit *.imh dispaxis 2 add+
+.fi
+.ih
+SEE ALSO
+icfit, response
+.endhelp
diff --git a/noao/twodspec/longslit/doc/lscombine.hlp b/noao/twodspec/longslit/doc/lscombine.hlp
new file mode 100644
index 00000000..764c3b1b
--- /dev/null
+++ b/noao/twodspec/longslit/doc/lscombine.hlp
@@ -0,0 +1,296 @@
+.help lscombine Jun04 noao.twodspec.longslit
+.ih
+NAME
+lscombine -- Combine longslit images
+.ih
+USAGE
+lscombine input output
+.ih
+PARAMETERS
+.ls input
+List of input two-dimensional images to combine.  This task is typically
+used with dispersion calibrated longslit images though it will work with
+any 2D images.
+.le
+.ls output
+Output combined image.
+.le
+.ls headers = "" (optional)
+Optional output multiextension FITS file where each extension is a dataless
+headers from each input image.
+.le
+.ls bpmasks = "" (optional)
+Optional output bad pixel mask with good values of 0 and bad values of 1.
+Output pixels are marked as bad when no input pixels contributed to the
+output pixel.  The file name is also added to the output image header under
+the keyword BPM.
+.le
+.ls rejmask = "" (optional)
+Optional output mask file identifying rejected or excluded pixels.  The
+pixel mask is the size of the output image but there is one extra dimension
+with length equal to the number of input images.  Each element of the
+highest dimension is a mask corresponding to an input image with values of
+1 for rejected or excluded pixels and values of 0 for pixels which were
+used.  The order of the masks is the order of the input images and image
+header keywords, indexed by the pixel coordinate of the highest dimension
+identify the input images.  Note that the pixel positions are in the output
+pixel coordinate system.
+.le
+.ls nrejmasks = "" (optional)
+Optional output pixel mask giving the number of input pixels rejected or
+excluded from the input images.
+.le
+.ls expmasks = "" (optional)
+Optional output exposure mask giving the sum of the exposure values of
+the input images with non-zero weights that contributed to that pixel.
+Since masks are integer, the exposure values may be scaled to preserve
+dynamic range and fractional significance.  The scaling values are given in
+the header under the keywords MASKSCAL and MASKZERO.  Exposure values are
+computed from the mask values by scale * value + zero where scale is the
+value of the MASKSCAL keyword and zero is the value of the MASKZERO
+keyword.
+.le
+.ls sigma = "" (optional)
+Optional output sigma image.  The sigma is the standard deviation,
+corrected for a finite population, of the input pixel values (excluding
+rejected pixels) about the output combined pixel values.
+.le
+
+.ls logfile = "STDOUT" (optional)
+Optional output log file.  If no file is specified then no log information is
+produced.  The special filename "STDOUT" prints log information to the
+terminal.
+.le
+
+.ls interptype = "spline3"
+Image interpolation type for any resampling prior to combining.
+The allowed types are "nearest" (nearest neighbor), "linear" (bilinear),
+"poly3" (bicubic polynomial), "poly5" (biquintic polynomial), and "spline3"
+(bicubic polynomial).
+.le
+.ls x1 = INDEF, y1 = INDEF
+User coordinates of the first output column and line.  If INDEF then it
+is based on the smallest value over all the images.
+.le
+.ls x2 = INDEF, y2 = INDEF
+User coordinates of the last output column and line.  If INDEF then it
+is based on the largest value over all the images.
+.le
+.ls dx = INDEF, dy = INDEF
+User coordinate pixel interval of the output.  If INDEF then the it
+is based on smallest interval (i.e. highest dispersion) over all the images.
+.le
+.ls nx = INDEF, ny = INDEF
+Number of output pixels.  If INDEF then it is based on the values of the
+other coordinate parameters.
+.le
+
+.ls combine = "average" (average|median|sum)
+Type of combining operation performed on the final set of pixels (after
+offsetting, masking, thresholding, and rejection).  The choices are
+"average", "median", or "sum".  The median uses the average of the two central
+values when the number of pixels is even.  For the average and sum, the
+pixel values are multiplied by the weights (1 if no weighting is used)
+and summed.  The average is computed by dividing by the sum of the weights.
+If the sum of the weights is zero then the unweighted average is used.
+.le
+.ls reject = "none" (none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip)
+Type of rejection operation performed on the pixels remaining after offsetting,
+masking and thresholding.  The algorithms are described in the
+DESCRIPTION section.  The rejection choices are:
+
+.nf
+      none - No rejection
+    minmax - Reject the nlow and nhigh pixels
+   ccdclip - Reject pixels using CCD noise parameters
+  crreject - Reject only positive pixels using CCD noise parameters
+   sigclip - Reject pixels using a sigma clipping algorithm
+ avsigclip - Reject pixels using an averaged sigma clipping algorithm
+     pclip - Reject pixels using sigma based on percentiles
+.fi
+
+.le
+.ls outtype = "real" (none|short|ushort|integer|long|real|double)
+Output image pixel datatype.  The pixel datatypes are "double", "real",
+"long", "integer", unsigned short "ushort", and "short" with highest
+precedence first.  If "none" is specified then the highest precedence
+datatype of the input images is used.  When there is a mixture of
+short and unsigned short images the highest precedence become integer.
+The datatypes may be abbreviated to a single character.
+.le
+.ls outlimits = ""
+Output region limits in pixels specified as pairs of whitespace separated
+values.  The first two numbers are the limits along the first output image
+dimension, the next two numbers are the limits along the second dimension,
+and so on.  If the higher dimension limits are not specified they default
+to the full range.  Therefore, if no limits are specified then the full
+output is created.  Note that the output size is computed from all the
+input images including offsets if specified and the coordinates are
+relative to that size.
+.le
+.ls masktype = "none" (none|goodvalue)
+Type of pixel masking to use.  If "none" then no pixel masking is done
+even if an image has an associated  pixel mask.  Otherwise the
+value "goodvalue" will use any mask specified for the image under
+the BPM keyword.  The values of the mask will be interpreted as
+zero for good pixels and non-zero for bad pixels.  The mask pixels
+are assumed to be registered with the image pixels.
+.le
+.ls blank = 0.
+Output value to be used when there are no pixels.
+.le
+
+.ls scale = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>)
+Multiplicative image scaling to be applied.  The choices are none, multiply
+by the reciprocal of the mode, median, or mean of the specified statistics
+section, multiply by the reciprocal of the exposure time in the image header,
+multiply by the values in a specified file, or multiply by a specified
+image header keyword.  When specified in a file the scales must be one per
+line in the order of the input images.
+.le
+.ls zero = "none" (none|mode|median|mean|@<file>|!<keyword>)
+Additive zero level image shifts to be applied.  The choices are none, add
+the negative of the mode, median, or mean of the specified statistics
+section, add the values given in a file, or add the values given by an
+image header keyword.  When specified in a file the zero values must be one
+per line in the order of the input images.  File or keyword zero offset
+values do not allow a correction to the weights.
+.le
+.ls weight = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>)
+Weights to be applied during the final averaging.  The choices are none,
+the mode, median, or mean of the specified statistics section, the exposure
+time, values given in a file, or values given by an image header keyword.
+When specified in a file the weights must be one per line in the order of
+the input images and the only adjustment made by the task is for the number of
+images previously combined.   In this case the weights should be those
+appropriate for the scaled images which would normally be the inverse
+of the variance in the scaled image.
+.le
+.ls statsec = ""
+Section of images to use in computing image statistics for scaling and
+weighting.  If no section is given then the entire region of the input is
+sampled (for efficiency the images are sampled if they are big enough).
+When the images are offset relative to each other one can precede the image
+section with one of the modifiers "input", "output", "overlap".  The first
+interprets the section relative to the input image (which is equivalent to
+not specifying a modifier), the second interprets the section relative to
+the output image, and the last selects the common overlap and any following
+section is ignored.
+.le
+.ls  expname = ""
+Image header keyword to be used with the exposure scaling and weighting
+options.  Also if an exposure keyword is specified that keyword will be
+added to the output image using a weighted average of the input exposure
+values.
+.le
+
+.ce
+Algorithm Parameters
+.ls lthreshold = INDEF, hthreshold = INDEF
+Low and high thresholds to be applied to the input pixels.  This is done
+before any scaling, rejection, and combining.  If INDEF the thresholds
+are not used.
+.le
+.ls nlow = 1,  nhigh = 1 (minmax)
+The number of low and high pixels to be rejected by the "minmax" algorithm.
+These numbers are converted to fractions of the total number of input images
+so that if no rejections have taken place the specified number of pixels
+are rejected while if pixels have been rejected by masking, thresholding,
+or nonoverlap, then the fraction of the remaining pixels, truncated
+to an integer, is used.
+.le
+.ls nkeep = 1
+The minimum number of pixels to retain or the maximum number to reject
+when using the clipping algorithms (ccdclip, crreject, sigclip,
+avsigclip, or pclip).  When given as a positive value this is the minimum
+number to keep.  When given as a negative value the absolute value is
+the maximum number to reject.  The latter is in addition to pixels
+missing due to non-overlapping offsets, bad pixel masks, or thresholds.
+.le
+.ls mclip = yes (ccdclip, crreject, sigclip, avsigcliip)
+Use the median as the estimate for the true intensity rather than the
+average with high and low values excluded in the "ccdclip", "crreject",
+"sigclip", and "avsigclip" algorithms?  The median is a better estimator
+in the presence of data which one wants to reject than the average.
+However, computing the median is slower than the average.
+.le
+.ls lsigma = 3., hsigma = 3. (ccdclip, crreject, sigclip, avsigclip, pclip)
+Low and high sigma clipping factors for the "ccdclip", "crreject", "sigclip",
+"avsigclip", and "pclip" algorithms.  They multiply a "sigma" factor
+produced by the algorithm to select a point below and above the average or
+median value for rejecting pixels.  The lower sigma is ignored for the
+"crreject" algorithm.
+.le
+.ls rdnoise = "0.", gain = "1.", snoise = "0." (ccdclip, crreject)
+CCD readout noise in electrons, gain in electrons/DN, and sensitivity noise
+as a fraction.  These parameters are used with the "ccdclip" and "crreject"
+algorithms.  The values may be either numeric or an image header keyword
+which contains the value.  The noise model for a pixel is:
+
+.nf
+    variance in DN = (rdnoise/gain)^2 + DN/gain + (snoise*DN)^2
+    variance in e- = (rdnoise)^2 + (gain*DN) + (snoise*(gain*DN))^2
+		   = rdnoise^2 + Ne + (snoise * Ne)^2
+.fi
+
+where DN is the data number and Ne is the number of electrons.  Sensitivity
+noise typically comes from noise introduced during flat fielding.
+.le
+.ls sigscale = 0.1 (ccdclip, crreject, sigclip, avsigclip)
+This parameter determines when poisson corrections are made to the
+computation of a sigma for images with different scale factors.  If all
+relative scales are within this value of unity and all relative zero level
+offsets are within this fraction of the mean then no correction is made.
+The idea is that if the images are all similarly though not identically
+scaled, the extra computations involved in making poisson corrections for
+variations in the sigmas can be skipped.  A value of zero will apply the
+corrections except in the case of equal images and a large value can be
+used if the sigmas of pixels in the images are independent of scale and
+zero level.
+.le
+.ls pclip = -0.5 (pclip)
+Percentile clipping algorithm parameter.  If greater than
+one in absolute value then it specifies a number of pixels above or
+below the median to use for computing the clipping sigma.  If less
+than one in absolute value then it specifies the fraction of the pixels
+above or below the median to use.  A positive value selects a point
+above the median and a negative value selects a point below the median.
+The default of -0.5 selects approximately the quartile point.
+.le
+.ls grow = 0.
+Radius in pixels for additional pixel to be rejected in an image with a
+rejected pixel from one of the rejection algorithms.  This applies only to
+pixels rejected by one of the rejection algorithms and not the masked or
+threshold rejected pixels.
+.le
+.ih
+DESCRIPTION
+\fBLSCOMBINE\fR combines two-dimensional longslit images by first
+resampling them to a common world coordinate system, if not already on
+the same system, and then combining the matching pixels.  The final world
+coordinate system is specified by parameters or by looking at the maximum
+ranges and minimum intervals over the input data.
+
+Algorithmically it is a combination of the tasks \fBTRANSFORM\fR (using
+the WCS) and \fBIMCOMBINE\fR.  When executing it will generate temporary
+images ("lsc*") and masks ("mlsc*") if the images are not already on a
+common world coordinate system.  The user only need be aware of this
+in case of an unexpected abort leaving these files behind.
+
+Rather than repeat the details the user should consult the descriptions
+for \fBTRANSFORM\fR and \fBIMCOMBINE\fR ignoring parameters which are
+not part of this task.
+.ih
+EXAMPLES
+.nf
+    cl> lscombine obj* lscomb
+.fi
+.ih
+NOTES
+.ls LSCOMBINE: V2.12.3
+This is a new task in this relese.
+.le
+.ih
+SEE ALSO
+transform, imcombine. odcombine 
+.endhelp
diff --git a/noao/twodspec/longslit/doc/lslit.ms b/noao/twodspec/longslit/doc/lslit.ms
new file mode 100644
index 00000000..de35424f
--- /dev/null
+++ b/noao/twodspec/longslit/doc/lslit.ms
@@ -0,0 +1,712 @@
+.nr PS 9
+.nr VS 10
+.ps 9
+.vs 10
+.po 0.50i
+.nr PO 0.50i
+.ll 7.0i
+.nr LL 7.0i
+.nr PD 1v
+.EQ
+delim	$$
+.EN
+.TL
+Reduction of long slit spectra with IRAF
+.AU
+Francisco Valdes
+.AI
+IRAF Group, Central Computer Services, National Optical Astronomy Observatories
+P.O. Box 26732, Tucson, Arizona, 85726
+March 1986
+.AB
+Tools for the reduction of long slit spectra within the Interactive
+Data Reduction and Analysis Facility (IRAF) at the National Optical
+Astronomy Observatory (NOAO) are described.  The user interface
+(commands and special features) and the algorithms are discussed.
+Application of the reduction package to multi-slit images is briefly
+outlined.  The author developed and supports the package at NOAO.
+.AE
+.LP
+
+.ce
+\fB1. Introduction\fR
+.PP
+This paper describes the tools currently available within the Interactive Data
+Reduction and Analysis Facility (IRAF) at the National Optical
+Astronomy Observatories (NOAO) for the reduction of long slit spectra.
+The reduction tools, called tasks, are organized as an IRAF package
+called \fBlongslit\fR.  The tasks in the package are summarized below.
+
+.TS
+center;
+n n.
+apdefine \&- Define apertures for 1D aperture extraction	identify \&- Identify features
+apextract \&- Extract 1D aperture spectra	illumination \&- Determine illumination calibration
+background \&- Fit and subtract a line or column background	reidentify \&- Reidentify features
+extinction \&- Apply atmospheric extinction corrections to images	response \&- Determine response calibration
+fitcoords \&- Fit user coordinates to image coordinates	setimhdr \&- Set longslit image header parameters
+fluxcalib \&- Apply flux calibration to images	transform \&- Transform longslit images to user coordinates
+.TE
+
+.PP
+Since there are many types of long slit spectra, detectors, and
+astronomical goals we do not describe a reduction procedure or path.
+Reduction manuals giving cookbook instructions for the reduction of
+certain types of data at NOAO are available from the Central Computer
+Services Division.  Instead, each task is discussed separately.  The
+primary emphasis is on the algorithms.
+.PP
+The following terminology is used in this paper.  A \fIlong slit
+spectrum\fR is a two dimensional image.  The two image axes are
+called \fIaxis 1\fR and \fIaxis 2\fR and the pixel coordinates are
+given in terms of \fIcolumns\fR and \fIlines\fR.  The long slit
+axes are called the \fIdispersion axis\fR and the \fIslit
+axis\fR.  The reduction tasks do not require a particular orientation
+of the dispersion and slit axes, however, these axes should be
+fairly closely aligned with the image axes.  \fBIn the remainder of
+this paper the slit axis will correspond to image axis 1 and
+the dispersion axis with image axis 2\fR.
+.PP
+There are five types of operations performed by the tasks in the
+\fBlongslit\fR package: (1) detector response calibration, (2) geometric
+distortion and coordinate rectification, (3) background sky subtraction,
+(4) flux calibration, and (5) aperture extraction of one dimensional spectra.
+These are listed in the order in which they are usually performed and in
+which they are discussed in this paper.  There is also an initialization
+task, \fBsetimhdr\fR, and a general routine, \fBicfit\fR, used in may of the
+long slit tasks.  These are described first.
+.SH
+SETIMHDR - Set long slit image header parameters
+.PP
+The tasks in the \fBlongslit\fR package use information contained in the IRAF
+image header.  The task \fBsetimhdr\fR sets a required parameter in the image
+header advising the long slit tasks which image axis corresponds to the
+dispersion axis; the tasks work equally well with the dispersion axis
+aligned with the image lines or the image columns.  This is generally
+the first task executed when reducing long slit spectra.
+.SH
+ICFIT - The IRAF Interactive Curve Fitting routine
+.PP
+Many of the tasks in the IRAF which fit a one dimensional function
+utilize the same powerful interactive curve fitting routine called
+\fBicfit\fR.  This routine allows the user to perform sophisticated
+function fitting interactively and graphically or to specify the
+function fitting parameters in advance and run the task
+non-interactively.  That this routine is used in many tasks also has
+the advantage that the user need not learn a new set of commands and
+features for each task requiring function fitting.
+.PP
+The features of the this curve fitting tool include:
+.IP (1)
+A choice of four fitting functions; Chebyshev polynomial, Legendre polynomial,
+a linear spline, and a cubic spline.
+.nr PD 0v
+.IP (2)
+A choice of the polynomial order or the number of spline pieces.
+.IP (3)
+Deletion of individual points from the fit.
+.IP (4)
+Selection of a sample or subset of points to be fit (excluding the rest).
+.IP (5)
+Iterative deletion of points with large residuals from the fitted function.
+.IP (6)
+Binning sets of neighboring points into averages or medians which are then
+fit instead of the individual points.
+.nr PD 1v
+.LP
+In addition to the above features the interactive graphics mode allows
+the user to:
+.IP (1)
+Iterate any number of times on the fitting parameters.
+.nr PD 0v
+.IP (2)
+Display the fit in several different ways; residuals, ratios, and the fit
+overplotted on the data points.
+.IP (3)
+Manipulate the graphs using a large set of commands for formating and
+expanding any part of a graph for detailed examination.
+.IP (4)
+Produce copies of the graphs with a snap-shot command.
+.nr PD 1v
+.PP
+For the applications described in this paper the most important features
+are the ability to adjust the function order, exclude bad points, and
+select subsets of points to be fit.  Other useful features are taking the
+median or average of a set of points before fitting and iteratively
+rejecting deviant points.  When used non-interactively the user
+selects the function and the order.  The \fBlongslit\fR tasks using the
+interactive curve fitting routine are \fBbackground\fR, \fBidentify\fR,
+\fBillumination\fR, and \fBresponse\fR.
+
+
+.ce
+\fB2. Detector Response Calibrations\fR
+.PP
+The relative response of the pixels in the detector and the transmission
+of the spectrograph along the slit are generally not uniform.  Outside
+of the \fBlongslit\fR package are IRAF tasks for creating \fIflat fields\fR
+from quartz lamp calibration images which correct for small scale response
+variations.  Flat fields, however, do not correct for spectrograph
+transmission variations or any large scale response patterns.  The tasks
+\fBresponse\fR and \fBillumination\fR are specially designed for long slit
+spectra to correct both the small scale variations as well as
+larger scale response patterns and slit illumination and transmission effects.
+.PP
+These algorithms make the assumption that the wavelength and slit axis
+are very nearly aligned with the image lines and columns.  If this is
+not true then the images must be aligned first or alternate response
+calibration methods used.
+.SH
+RESPONSE - Determine response calibration
+.PP
+The task \fBresponse\fR is used with calibration images which (1)
+do not have any intrinsic structure along the slit dimension and (2)
+have a smooth spectrum without emission or absorption features.
+Typically the calibration images consist of quartz lamp exposures.
+The idea is to determine a response correction that turns an observed
+calibration image into one which is identical at all points along the
+slit.
+.PP
+From (1) a one dimensional spectrum is obtained by averaging along the
+slit; i.e. averaging the columns.  Based on (2) a smoothing function is
+fit to the one dimensional spectrum to reduce noise and eliminate
+response effects which are coherent in wavelength such as fringing.
+The response correction for each pixel is then obtained by dividing
+each point along the slit (the columns) by the smoothed one dimensional
+spectrum.
+.PP
+The purpose of fitting a function to the one dimensional spectrum is to
+reduce noise and to remove coherent response effects which are not part
+of the true quartz spectrum.  Examples of coherent response effects are
+fringing and regions of low or high response running along the slit
+dimension which are, therefore, not averaged out in the one dimensional
+spectrum.  The choice of smoothing function is dictated by the behavior
+of the particular detector.  Difficult cases are treated with the
+interactive graphical function fitting routine \fBicfit\fR.  For the
+automated case the user specifies the smoothing function and order.
+.PP
+This calibration algorithm has the advantage of removing spatial
+frequencies at almost all scales; in particular, there is no modeling
+of the response pattern along the slit dimension.  The only modeling is
+the fit to the \fBaverage\fR spectrum of the calibration source.  In
+tests at NOAO this algorithm was able to reduce the response variations
+to less 0.2%, to correct for a broad diagonal region of low response in
+one of the CCD detectors (the CRYOCAM), and to remove strong fringing
+in spectra taken in the red portion of the spectrum where the detector
+is particularly subject to fringing.
+.PP
+One feature common to \fBresponse\fR and \fBillumination\fR is that
+the algorithm can be restricted to a section of the calibration image.
+The response corrections are then determined only within that section.
+If a response image does not exist initially then the response values outside
+the section are set to unity.  If the response image does exist then
+the points outside the section are not changed.  This feature is used
+with data containing several slits on one image such as produced by
+the multi-slit masks at Kitt Peak National Observatory.
+.PP
+When there are many calibration images this algorithm may be applied to
+each image separately or to an average of the images.  If applied
+separately the response images may be averaged or applied to the
+appropriate long slit spectra; typically the one nearest the object
+exposure in time or telescope position.  The task allows a list of
+calibration images from which a set of response corrections is
+determined.
+.PP
+Figure 1 shows a portion of an average quartz spectrum ratioed with the
+smooth fit to the spectrum.  It is one of the graphs which can be
+produced with the \fBicfit\fR routine and, with the other figures in
+this paper, illustrates the formating,
+zooming, and snap-shot capabilities in IRAF.  The figure shows considerable
+structure of periodic high response lines and fringing which, because
+they are primarily aligned with the image lines, are still present in
+the average quartz spectrum.  Note that this is not the response
+since it is the average of all the columns; an actual response column
+would have much larger variations including pixel-to-pixel response
+differences as well as large scale response patterns such as the diagonal
+structure mentioned previously.
+.SH
+ILLUMINATION - Determine illumination calibration
+.PP
+The task \fBillumination\fR corrects for large scale variations along
+the slit and dispersion dimensions due to illumination or spectrograph
+transmission variations (often called the \fIslit profile\fR).  When
+the detector response function is determined from quartz calibration
+images, using \fBresponse\fR, an illumination error may be introduced
+due to differences in the way the spectrograph is illuminated by the
+quartz lamp compared to that of an astronomical exposure.  This
+violates the the assumption that the calibration spectrum has no
+intrinsic structure along the slit.  \fBIllumination\fR is also used
+when only the small scale response variations have been removed using a
+flat field correction.
+.PP
+The approach to determining the response correction is similar to that
+described for \fBresponse\fR.  Namely, the response correction is the
+ratio of a calibration image to the expected calibration image.  Again,
+the expected calibration image is that which has no structure along the
+slit.  Calibration images may be quartz lamp exposures, assuming there
+is no illumination problem, and blank sky exposures.  In the worst
+case, object exposures also may be used if the extent of the object in
+the slit is small.
+.PP
+There are several important differences between this algorithm and that
+of \fBresponse\fR:
+.IP (1)
+The spectra are not required to be smooth in wavelength and may contain
+strong emission and absorption lines.
+.nr PD 0v
+.IP (2)
+The response correction is a smooth, large scale function only.
+.IP (3)
+Since the signal-to-noise of spectra from blank sky and object images is
+lower than quartz calibration images, steps must be taken to minimize noise.
+.IP (4)
+Care must be taken that the spectral features do not affect the
+response determination.
+.nr PD 1v
+.PP
+The algorithm which satisfies these requirements is as follows.  First the
+calibration spectrum is binned in wavelength.  This addresses the
+signal-to-noise consideration (3) and is permitted because only large
+scale response variations are being determined (2).  Next a smoothing
+function is fit along the slit dimension in each bin; i.e. each
+wavelength bin is smoothed to reduce noise and determine the large
+scale slit profile.  Then each bin is normalized to the central point
+in the slit to remove the spectral signature of the calibration image.
+Finally, the binned response is interpolated back to the
+original image size.
+.PP
+The normalization to the central point in the slit is an assumption
+which limits the ability of the illumination algorithm to correct
+for all wavelength dependent response effects.  There is a wavelength
+dependence, however, in that the slit profile is a function of the
+wavelength though normalized to unity at the central point of the
+slit.
+.PP
+The wavelength bins and bin widths need not be constant.  The bins are
+chosen to sample the large scale variations in the slit profile as a
+function of wavelength, to obtain good signal statistics, and to avoid
+effects due to variations in the positions and widths of strong
+emission lines.  This last point means that bin boundaries should not
+intersect strong emission lines though the bin itself may and should
+contain strong lines.  Another way to put this criterion is that
+changes in the data in the wavelength bins should be small when the
+bin boundaries are changed slightly.
+.PP
+The bins may be set interactively using a graph of the average
+spectrum or automatically by dividing the dispersion axis into a
+specified number of equal width bins.  When the number of bins is small
+(and the number of wavelength points in each bin is large) bin
+boundary effects are likely to be insignificant.
+A single bin consisting of all wavelengths, i.e. the sum of all the image
+lines, may be used if no wavelength dependence is expected in the
+response.  Illumination effects introduced with \fBresponse\fR,
+however, appear as wavelength dependent variations in the slit
+profile.
+.PP
+Smoothing of each bin along the slit dimension is done with the
+interactive curve fitting routine.  The curve fitting may be done
+graphically and interactively on any set of bins or automatically by
+specifying the function and order initially.  The fitting should be
+done interactively (at least on the first bin) in order to exclude
+objects when the sky is not truly blank and contains faint objects or
+when object exposures must be used to determine the slit profile.
+.PP
+As with \fBresponse\fR, several blank sky images may be available
+(though this is less often true in practice).  An illumination
+correction may be determined for each calibration image or one
+illumination correction may be computed from the average of the
+calibration images.  Also the illumination response correction may be
+determined for only a section of the calibration image so as to be
+applicable to multi-slit data.
+.PP
+Figure 2 shows the fit to one of the wavelength bins; lines 1 to 150 have been
+summed and the sum is plotted as a function of slit position (column).
+The data is from a response image produced by \fBresponse\fR.  This
+figure illustrates a number of things.  \fBIllumination\fR may be run
+on a response image to remove the large scale illumination and slit
+transmission effects.  This creates a flat field in a manner different than
+normal surface fitting.  The figure shows that response effects occur
+at all scales (keeping in mind that the pixel-to-pixel response has
+been largely averaged out by summing 150 columns).  It also illustrates
+how the illumination algorithm works for a typical slit profile.  In
+this example about half the large scale variation in the slit profile
+is due to illumination effects and half is real slit transmission
+variations.  For a blank sky or object image the main differences
+would be larger data values (hundreds to thousands) and possibly
+objects present in the slit to be excluded from the fit.
+
+
+.ce
+\fB3. Distortion Corrections and Coordinate Transformations\fR
+.PP
+The removal of geometric distortions and the application of coordinate
+transformations are closely related.  Both involve applying a
+transformation to the observed image to form the desired final image.
+Generally, both steps are combined into a single image transformation
+producing distortion corrected images with linear wavelength
+coordinates (though the pixel interval may be logarithmic).
+This differs from other systems (for example, the Kitt Peak IPPS) which
+perform distortion corrections on each axis independently and then
+apply a dispersion correction on the distortion corrected image.
+While this approach is modular it requires several transformations of
+the images and does not couple the distortions in each dimension into
+a single two dimensional distortion.
+.PP
+To transform long slit images requires (1) identifying spectral
+features and measuring their positions in arc lamp or sky
+exposures at a number of points in the image, (2) determining the
+distortions in the slit positions at a number of points along the
+dispersion axis using either calibration images taken with special
+masks or narrow objects such as stars,
+(3) determining a transformation function between the image
+coordinates and the user coordinates for the measured wavelength and
+slit positions, (4) and interpolating the images to a uniform grid in
+the user coordinates according to the transformation function.  The
+coordinate feature information and the transformation functions are
+stored in a database.  If needed, the database may be examined and
+edited.
+.PP
+An important part of this task is the feature center determination.  This
+algorithm is described in a separate section below.
+.SH
+IDENTIFY - Identify features
+.PP
+The tasks \fBidentify\fR and \fBreidentify\fR are general tools used
+for one dimensional, multi-aperture, multi-slit, echelle, and long slit
+spectra.  The tasks are also general in the sense that they are used to
+identify features in any one dimensional vector.  For long slit
+reductions they are used to identify and trace objects in the slit and
+to identify, trace, and determine wavelength solutions for spectral
+features from arc calibration images and from sky and object
+exposures.
+.PP
+\fBIdentify\fR is used to identify emission or absorption features in a
+one dimensional projection of an image.  This projection consists of an
+image line or column or the
+average of many lines or columns.  Averaging is used to increase the
+signal in weak features and provide better accuracy in determining the
+one dimensional positions of the features.  The identified features are
+assigned user coordinates.  The user coordinates will ultimately define
+the final coordinates of the rectified images.
+.PP
+For determining the distortions along the slit, the positions of object
+profiles or profiles obtained with multi-aperture masks in the slit
+are measured at a reference line.  The user coordinates are then taken to be
+the positions at this reference line.  The
+coordinate rectification will then correct for the distortion to bring the
+object positions at the other lines to the same position.
+(Note that it is feasible to make an actual coordinate transformation of
+the spatial axis to arc seconds or some other units).
+.PP
+For wavelength features arc calibration images are generally used,
+though sky and object exposures can also be used if necessary.  After
+marking a number of spectral features and assigning them wavelength
+coordinates a \fIdispersion solution\fR can be computed relating the
+image coordinate to the wavelength; $lambda~=~f(l)$, where $lambda$ is
+wavelength and $l$ is the image line.  The dispersion
+solution is determined using the \fBicfit\fR routines described
+earlier.  This dispersion solution is used in the long slit package
+only as an aid in finding misidentified lines and to automatically add
+new features from a wavelength list.  The dispersion solution actually
+used in transforming the images is a two dimensional function
+determined with the task \fBfitcoords\fR.
+.PP
+Figure 3 shows a graph from \fBidentify\fR used on a Helium-Neon-Argon
+arc calibration image.  Only three lines were identified interactively
+and the reminder were added automatically from a standard line list.
+Note that the abscissa is in wavelength units and the ordinate is
+displayed logarithmically.  The latter again illustrates the flexibility
+the user has to modify the graph formats.  Each marked feature is
+stored in a database and is automatically reidentified at other columns
+in the image with \fBreidentify\fR.
+.SH
+REIDENTIFY - Reidentify features
+.PP
+The task \fBreidentify\fR automatically reidentifies the spectral and
+object features and measures their positions at a number of other
+columns and lines starting from those identified interactively with
+\fBidentify\fR.  The algorithms and the feature information produced is
+the same as that of \fBidentify\fR including averaging a number of
+lines or columns to enhance weak features.  The automatic tracing can
+be set to stop or continue when a feature fails to be found in a new
+column or line; failure is defined by the position either becoming
+indeterminate or shifting by more than a specified amount
+(\fIcradius\fR defined in the next section).
+.SH
+CENTER1D - One dimensional feature centering
+.PP
+The one dimensional position of a feature is determined by solving the equation
+
+.EQ
+define I0 'I sub 0'
+define XC 'X sub c'
+.EN
+.EQ (1)
+int ( I - I0 ) f( X - XC ) dX~=~0
+.EN
+
+where $I$ is the intensity at position $X$, $I0$ is the continuum
+intensity, $X$ is the vector coordinate, and $XC$ is the desired
+feature position.  The convolution function $f(X- XC )$ is a
+sawtooth as shown in figure 4.  For absorption features the negative of this
+function is used.  The figure defines the parameter \fIfwidth\fR which
+is set to be approximately the width of the feature.  If it is too
+large the centering may be affected by neighboring features and if it
+is too small the accuracy is worse.
+.PP
+For emission features the continuum, $I0$, is assumed to be zero.
+For absorption features the continuum
+is the maximum value in the region around the initial guess
+for $XC$.  The size of the region on each side of the initial guess is
+the sum of \fIfwidth\fR/2, to allow for the feature itself, \fIcradius\fR,
+to allow for the uncertainty in the feature position, and \fIfwidth\fR, for a
+buffer.  Admittedly this is
+not the best continuum but it contains the fewest assumptions and is
+tolerant of nearby contaminating features.
+.PP
+Equation (1) is solved iteratively starting with the initial position.
+When successive positions agree within 0.1% of a pixel the position is
+returned.  If the position wanders further than the user defined
+distance \fIcradius\fR from the initial guess or outside of the data
+vector then the position is considered to be indefinite.
+.SH
+FITCOORDS - Fit user coordinates to image coordinates
+.PP
+Let us denote the image coordinates of a point in the two dimensional
+image as $(c,~l)$ where $c$ is the column coordinate
+and $l$ is the line coordinate.  Similarly, denote the
+long slit coordinates as $(s,~lambda )$ where $s$ is
+the slit position and $lambda$ is the wavelength.
+The results of \fBidentify\fR and \fBreidentify\fR is a set of points
+$(c,~l,~s)$ and $(c,~l,~lambda )$ recorded in the database.
+.PP
+Two dimensional functions of the image coordinates are fit to the user
+coordinates for each set of slit and wavelength features,
+$s~=~t sub s (c, l)$ and $lambda~=~t sub lambda (c, l)$, which are
+stored in the database.
+Note that the second function is a two dimensional dispersion solution.
+It is this function which is used to transform the long slit images to
+linear wavelength coordinates.  Many images may be used to create a
+single transformation or each calibration images may be used separately
+to create a set of transformations.
+.PP
+This task has both an interactive and non-interactive mode.  For the
+non-interactive mode the user specifies the transformation function,
+either a two dimensional Chebyshev or Legendre polynomial, and separate
+orders for the column and line axes.  When run interactively the
+user can try different functions and orders, delete bad points, and
+examine the data and the transformation in a variety of graphical formats.
+The interactive option is quite useful in initially setting the
+transformation function parameters and deleting bad points.
+The two dimensional function fitting routine is similar in spirit to the
+\fBicfit\fR one dimensional function fitting routine.  It is possible
+that this routine may find uses in other IRAF tasks.
+.PP
+Figure 5 shows a graph from \fBfitcoords\fR.  The feature image coordinates
+of four objects in the slit (the first of which is very weak)
+from \fBidentify\fR and \fBreidentify\fR are plotted.  This information
+is used to measure the distortion of the spectrograph in the slit axis.
+This example shows particularly gross distortions; often the distortions
+would not be visible in such a graph, though expanding it would make
+the distortion visible.  The transformation surface fit to this data
+removes this distortion almost entirely as seen in the residual plot
+of figure 6.  Figure 7 shows the equivalent residual plot for the
+wavelength coordinates; a two dimensional dispersion solution.
+.SH
+TRANSFORM - Transform long slit images to user coordinates
+.PP
+The coordinate transformations determined with the task \fBfitcoords\fR are
+read from the database.  The transformations are evaluated on a grid of
+columns and lines, $s sub i~=~t sub s (c sub i , l sub i )$ and
+$lambda sub i~=~t sub lambda (c sub i , l sub i )$.
+If no transformation is defined for a particular dimension then a unit
+transformation is used.  If more than one transformation for a dimension
+is given then a set of points is computed for each transformation.
+The inverse transformations are obtained by fitting transformation
+functions of the same type and orders to the set of slit position and
+wavelength points.  Note how this allows combining separate
+transformations into one inverse transformation.
+.PP
+The inverse transformations, $c~=~t sub c (s, lambda )$ and
+$l~=~t sub l (s, lambda )$, are used to rectify a set of input images.
+The user specifies a linear grid for the transformed images by defining some
+subset of the starting and ending coordinates, the pixel interval, and the
+number of points.  In addition the pixel interval can be specified to be
+logarithmic; used primarily on the wavelength axis for radial
+velocity studies.  The inverse transformations define the image column
+and line to be interpolated in the input image.  The user has the choice
+of several types of image interpolation; bilinear, bicubic, and biquintic
+polynomials and bicubic spline.  In addition the interpolation
+can be specified to conserve flux by multiplying the interpolated value
+by the Jacobian of the transformation.
+.PP
+The wavelength of the first pixel and the pixel wavelength interval are
+recorded in image headers for later use in making plots and in the
+\fBonedspec\fR package.  In addition a flag is set in the header indicating
+that the image has been dispersion corrected.
+
+
+.ce
+\fB4. Background Subtraction\fR
+.SH
+BACKGROUND - Fit and subtract a line or column background
+.PP
+If required, the background sky at each wavelength is subtracted from
+the objects using regions of the slit not occupied by the object.
+This must be done on coordinate rectified images since the lines or
+columns of the image must correspond exactly to the same wavelength.
+A set of points along the slit dimension, which are representative of the
+background, are chosen interactively.  Generally this will consist of two
+strips on either side of the object spectrum.
+At each wavelength a low order function is fit to the sky points and then
+subtracted from the entire line or column.
+.PP
+Ideally the response corrections and coordinate rectification will make
+the background sky constant at all points on the slit at each
+wavelength and the subtracted background is just a constant.  However, if
+desired a higher order function may be used to correct for
+deficiencies in the data.  A possible problem is focus variations which
+cause the width of the sky emission lines to vary along the slit.  One
+may partially compensate for the focus variations by using a higher
+order background fitting function.
+.PP
+The background fitting uses the
+interactive curve fitting routine \fBicfit\fR described earlier.
+Figure 8 shows a graph from \fBbackground\fR illustrating how the user
+sets two sample regions defining the sky (indicated a the bottom of
+the graph).
+
+
+.ce
+\fB5. Flux Calibration\fR
+.SH
+EXTINCTION - Apply atmospheric extinction corrections to images
+.PP
+A set of coordinate rectified images is corrected for atmospheric
+extinction with the task \fBextinction\fR.  The extinction correction
+is given by the formula
+
+.EQ
+    roman {correction~factor}~=~10 sup {0.4~E sub lambda~A}
+.EN
+
+where $E sub lambda$ are tabulated extinctions values and $A$ is the air
+mass of the observation (determined from information in the image
+header).  The tabulated extinctions are interpolated to the wavelength of
+each pixel and the correction applied to the input pixel value to form
+the output pixel value.  The user may supply the extinction table but
+generally a standard extinction table is used.
+.PP
+The air mass is sought in the image header under the keyword AIRMASS.
+If the air mass is not found then it is computed from the zenith
+distance, ZD, using the approximation formula from Allen's
+"Astrophysical Quantities", 1973, pages 125 and 133
+
+.EQ
+	A = ( cos ( roman ZD ) sup 2~+~2 s~+~1) sup half
+.EN
+
+where $s$, the atmospheric scale height, is set to be 750.  If the
+zenith distance is not found then it must be computed from the
+hour angle, the declination, and the observation latitude.   The
+hour angle may be computed from the right ascension and the siderial time.
+Computed quantities are recorded in the image header.
+Flags indicating extinction correction are also set in the image
+header.
+.SH
+FLUXCALIB - Apply flux calibration to images
+.PP
+The specified images are flux calibrated using a flux calibration file
+derived with the \fBonedspec\fR package using standard stars.  The
+standard stars are extracted from response corrected, coordinate
+rectified, and background subtracted long slit images using the tasks
+\fBapdefine\fR and \fBapextract\fR.  The standard stars must not be
+extinction corrected because this is done by the \fBonedspec\fR flux
+calibration algorithms.  The user may specify flux per unit wavelength,
+$roman F sub lambda$, or flux per unit frequency, $roman F sub nu$.
+The flux is computed using the exposure time and dispersion from the
+image headers and a flux calibration flag is set.
+
+
+.ce
+\fB6. Extraction of One Dimensional Spectra\fR
+.PP
+The user may wish to extract one dimensional spectra at various points
+along the slit.  As mentioned earlier, this is necessary if observations
+of standard stars are to be used to calibrate the fluxes.  The flux
+calibration values are determined from one dimensional spectra of standard
+stars using the \fBonedspec\fR package.  The tools to extract
+one dimensional aperture spectra from long slit spectra are \fBapdefine\fR and
+\fBapextract\fR.
+.SH
+APDEFINE - Define apertures for 1D aperture extraction
+.PP
+Extraction apertures are defined as a list consisting of an
+aperture number and lower and upper limits for the aperture.  The aperture
+limits are specified as column or line positions which need not be
+integers.  The user may create a file containing these
+aperture definitions with an editor or use the interactive
+graphics task \fBapdefine\fR.
+.PP
+\fBApdefine\fR graphs the sum of a number of lines or columns (depending
+on the dispersion axis) and allows the user to interactively define and
+adjust apertures either with the cursor or using explicit commands.
+If an aperture definition file exists the apertures are indicated on
+the graph initially.  When the user is done a new aperture definition
+file is written.
+.SH
+APEXTRACT - Extract 1D aperture spectra
+.PP
+One dimensional aperture spectra are extracted from a list of
+long slit images using an aperture definition file.  The extraction
+consists of the sum of the pixels, including partial pixels, at
+each column or line along the dispersion axis between the aperture limits.
+.PP
+More sophisticated algorithms than simple strip extraction are available
+in IRAF and will soon be incorporated in the long slit package.  The
+other extraction tasks trace the positions of features, i.e. the aperture
+is not fixed at certain columns or lines, and allow weighted extractions
+and detecting and removing bad pixels such as cosmic rays.  The
+weighted extractions can be chosen to be optimal in a statistical sense.
+
+
+.ce
+\fBConclusion\fR
+.PP
+The IRAF long slit reduction tasks have been used at NOAO for about six
+months and have yielded good results.  The package does not contain specific
+analysis tasks.  Some analysis task will be added in time.  The package
+is part of the software distributed with release of the IRAF.  The
+author of this paper wrote and supports the tasks described here.
+Any comments are welcome.
+.sp5
+.ll 4.2i
+.nr LL 4.2i
+.LP
+\fBCaptions for Figures:\fP
+.sp 1
+Figure 1.  Ratio of average quartz spectrum to fit of a 20 piece cubic spline
+for determination of response correction using \fBresponse\fR.
+
+Figure 2.  Fit of 4 piece cubic spline to the slit profile from the average
+of the first 150 lines in a response image using \fBillumination\fR.
+
+Figure 3.  Identification of emission lines from the central column of a
+Helium-Neon-Argon spectrum using task \fBidentify\fR.
+
+Figure 4.  Sawtooth convolution function of width \fIfwidth\fR used in the
+profile centering algorithm.
+
+Figure 5.  Graph of stellar object positions identified with \fBidentify\fR,
+traced with \fBreidentify\fR, and graphed by \fBfitcoords\fR showing the
+spectrograph distortions.
+
+Figure 6.  Residuals of the fit of a two dimensional 6th order Chebyshev
+polynomial to the data of figure 5 using \fBfitcoords\fR.
+
+Figure 7.  Residuals of the fit of a two dimensional 6th order Chebyshev
+polynomial to the image positions of wavelength features using \fBfitcoords\fR.
+
+Figure 8.  Constant background fit to a line of an object spectrum using
+\fBbackground\fR.  The marks at the bottom of the graph indicate the
+set of points used in the fit.
diff --git a/noao/twodspec/longslit/doc/response.hlp b/noao/twodspec/longslit/doc/response.hlp
new file mode 100644
index 00000000..61a7b34a
--- /dev/null
+++ b/noao/twodspec/longslit/doc/response.hlp
@@ -0,0 +1,178 @@
+.help response Aug86 noao.twodspec.longslit
+.ih
+NAME
+response -- Determine response calibrations
+.ih
+USAGE
+response calibration normalization response
+.ih
+PARAMETERS
+.ls calibration
+Images to use in determining response calibrations.  These are
+generally quartz continuum spectra.  An image section may be used to select
+only a portion of the image.
+.le
+.ls normalization
+Images to use determining the normalization spectrum.  In almost all cases
+the normalization images are the same as the calibration images or a
+subsection of the calibration images.
+.le
+.ls responses
+Response calibration images to be created.  Each response image is paired
+with a calibration image.  If the image exists then it will be modified
+otherwise it is created.
+.le
+.ls interactive = yes
+Graph the average calibration spectrum and fit the normalization spectrum
+interactively?
+.le
+.ls threshold = INDEF
+Set the response to 1 when the normalization spectrum or input image data
+fall below this value.  If INDEF then no threshold is applied.
+.le
+.ls sample = "*"
+Sample of points to use in fitting the average calibration spectrum.
+The sample is selected with a range string.
+.le
+.ls naverage = 1
+Number of sample points to average or median before fitting the function.
+If the number is positive the average of each set of naverage sample
+points is formed while if the number is negative then the median of each set
+of points (in absolute value) is formed.  This subsample of points is
+used in fitting the normalization spectrum.
+.le
+.ls function = "spline3"
+Function to fit to the average image spectrum to form the normalization
+spectrum.  The options are "spline1", "spline3", "legendre", and "chebyshev".
+.le
+.ls order = 1
+Order of the fitting function or the number of spline pieces.
+.le
+.ls low_reject = 0., high_reject = 0.
+Rejection limits below and above the fit in units of the residual sigma.
+.le
+.ls niterate = 1
+Number of rejection iterations.
+.le
+.ls grow = 0
+Reject additional points within this distance of points exceeding the
+rejection threshold.
+.le
+.ih
+CURSOR KEYS
+The interactive curve fitting package \fBicfit\fR is used to fit a function
+to the average calibration spectrum.  Help for this package is found
+under the name "icfit".
+.ih
+DESCRIPTION
+A response calibration, in the form of an image, is created for each input
+image, normally a quartz spectrum.  The response calibration is formed by
+dividing the calibration image by a normalization spectrum which is the
+same at all points along the spatial axis.  The normalization spectrum is
+obtained by averaging the normalization image across the dispersion to form
+a one dimensional spectrum and smoothing the spectrum by fitting a
+function.  The threshold value does not apply to creating or fitting of
+the normalization spectrum but only the final creation of the response
+values.  When normalizing (that is dividing the data values by the
+fit to the normalization spectrum) only pixels in which both the fitted
+normalization value and the data value are above the threshold are
+computed.  If either the normalization value or the data value is below
+the threshold the output response value is one.
+
+The image header keyword DISPAXIS must be present with a value of 1 for
+dispersion parallel to the lines (varying with the column coordinate) or 2
+for dispersion parallel to the columns (varying with line coordinate).
+This parameter may be added using \fBhedit\fR.  Note that if the image has
+been transposed (\fBimtranspose\fR) the dispersion axis should still refer
+to the original dispersion axis unless the physical world coordinate system
+is first reset (see \fBwcsreset\fR).  This is done in order to allow images
+which have DISPAXIS defined prior to transposing to still work correctly
+without requiring this keyword to be changed.
+
+If the output image does not exist it is first created with unit response
+everywhere.  Subsequently the response is only modified in those regions
+occupied by the input calibration image.  Thus, image sections may be used
+to select regions in which the response is desired.  This ability is
+particularly useful when dealing with multiple slits within an image or to
+exclude regions outside the slit.
+
+Normally the normalization images are the same as the calibration images.
+In other words the calibration image is normalized by the average spectrum
+of the calibration image itself.  Sometimes, however, the normalization
+image may be a smaller image section of the calibration image to avoid
+contaminating the normalization spectrum by effects at the edge of the
+slit.  Again, this may be quite useful in multi-slit images.
+
+The normalization spectrum is smoothed by fitting a function
+using the interactive curve fitting package (\fBicfit\fR).  The
+parameters determining the fitted normalization spectrum are the sample
+points, the averaging bin size, the fitting function, the order of the
+function, the rejection sigmas, the number of rejection iterations, and
+the rejection width.  The sample points for the average spectrum are
+selected by a range string.  Points in the normalization spectrum not in the
+sample are not used in determining the fitted function.  The selected
+sample points may be binned into a set of averages or medians which are
+used in the function fit instead of the sample points with the
+averaging bin size parameter \fInaverage\fR.  This parameter selects
+the number of sample points to be averaged if its value is positive or
+the number of points to be medianed if its value is negative
+(naturally, the absolute value is used for the number of points).  A
+value of one uses all sample points without binning.  The fitted
+function may be used to reject points from the fit using the parameters
+\fIlow_reject, high_reject, niterate\fR and \fIgrow\fR.  If one or both
+of the rejection limits are greater than zero then the sigma of the
+residuals is computed and points with residuals less than
+\fI-low_reject\fR times the sigma and greater than \fIhigh_reject\fR
+times the sigma are removed and the function fitted again.  In addition
+points within a distance given by the parameter \fIgrow\fR of the a
+rejected point are also rejected.  A value of zero for this parameter
+rejects only the points exceeding the rejection threshold.  Finally,
+the rejection procedure may be iterated the number of times given by
+the parameter \fIniterate\fR.
+
+The fitted function may be examined and modified interactively when the
+parameter \fIinteractive\fR is set.  In this case the normalization spectrum
+and the fitted function or the residuals of the fit are graphed.
+Deleted points are marked with an x and rejected points by a diamond.
+The sample regions are indicated along the bottom of the graph.
+The cursor keys and colon commands are used to change the values
+of the fitting parameters, delete points, and window and expand the
+graph.  When the fitted function is satisfactory exit with a carriage
+return or 'q' and the calibration image will be created.  Changes in
+the fitted parameters are remembered from image to image within the
+task but not outside the task.
+
+When the task finishes creating a response image the fitting parameters
+are updated in the parameter file.
+.ih
+EXAMPLES
+1. To create a response image non-interactively:
+
+	cl> response quartz quartz response order=20 interactive=no
+
+2. To determine independent responses for a multislit image determine the
+image sections defining each slit.  Then the responses are computed as
+follows:
+
+.nf
+	cl> response quartz[10:20,*],quartz[35:45,*] \
+	>>> quartz[12:18,*],quartz[12:18,*] resp,resp
+.fi
+
+Generally the slit image sections are prepared in a file which is then
+used to define the lists of input images and response.
+
+.nf
+	cl> response @slits @slits @responses
+.fi
+
+3.  If the DISPAXIS keyword is missing and the dispersion is running
+vertically (varying with the image lines):
+
+.nf
+	cl> hedit *.imh dispaxis 2 add+
+.fi
+.ih
+SEE ALSO
+icfit, iillumination
+.endhelp
diff --git a/noao/twodspec/longslit/doc/transform.hlp b/noao/twodspec/longslit/doc/transform.hlp
new file mode 100644
index 00000000..6955b51e
--- /dev/null
+++ b/noao/twodspec/longslit/doc/transform.hlp
@@ -0,0 +1,240 @@
+.help transform Sep87 noao.twodspec.longslit
+.ih
+NAME
+transform -- Transform longslit images to user coordinates
+.ih
+USAGE
+transform input output fitnames
+.ih
+PARAMETERS
+.ls input
+List of input images to be transformed.
+.le
+.ls output
+List of output images.  The number of output images in the list must
+match the number of input images.
+.le
+.ls minput = ""
+List of input masks or references.  This mask is used to create an output
+mask and is currently not used in the calculation of the output pixel
+values.  The list may be empty, a single element to apply to all input
+images, or a list that matches the input list.  A element in the list
+may be "BPM" to use the mask referenced by the standard bad pixel mask
+keyword "BPM", "!<keyword>" to use another header keyword pointing to a
+mask, or a mask filename.  The mask file is typically a pixel list file
+but it may also be an image.  The mask values are interpreted as zero and
+greater than zero with the actual values ignored.  The mask is assumed to
+be registered with the input and no coordinate system matching is used.
+The mask maybe smaller or larger than the input image with non-overlapping
+pixels ignored and missing pixels assumed to be zero valued.  The mask
+.le
+.ls moutput = ""
+List of output masks to be created.  The list may be empty or must match
+the input list.  Output masks may be specified even if no input mask is
+specified, in which case the output mask will identify pixels which map
+to regions outside the input images (also see the \fIblank\fR parameter).
+If an explicit extension is not specified a FITS mask is extension is
+created unless the environment variable "masktype" is set to "pl".
+.le
+.ls fitnames  
+Names of the user coordinate maps in the database to be used in the transform.
+If no names are specified, using the null string "", the world coordinate
+system (WCS) of the image is used.  This latter case may be used to
+resample previously WCS calibrated images to a different linear range
+or sampling.
+.le
+.ls database = "database"
+Database containing the coordinate map to be used in transforming the images.
+.le
+.ls interptype = "spline3"
+Image interpolation type.  The allowed types are "nearest" (nearest neighbor),
+"linear" (bilinear), "poly3" (bicubic polynomial), "poly5" (biquintic
+polynomial), and "spline3" (bicubic polynomial).
+.le
+.ls flux = yes
+Conserve flux per pixel?  If "no" then each output pixel is simply interpolated
+from the input image.  If "yes" the interpolated output pixel value is
+multiplied by the Jacobean of the transformation (essentially the ratio of
+pixel areas between the output and input images).
+.le
+.ls x1 = INDEF, y1 = INDEF
+User coordinates of the first output column and line.  If INDEF then the
+smallest value corresponding to a pixel from the image used to create the
+coordinate map is used.  These values are in user units regardless of whether
+logarithmic intervals are specified or not.
+.le
+.ls x2 = INDEF, y2 = INDEF
+User coordinates of the last output column and line.  If INDEF then the
+largest value corresponding to a pixel from the image used to create the
+coordinate map is used.  These values are in user units regardless of whether
+logarithmic intervals are specified or not.
+.le
+.ls dx = INDEF, dy = INDEF
+Output pixel intervals.  If INDEF then the interval is set to yield the
+specified number of pixels.  Note that for logarithmic intervals the
+interval must be specified as a base 10 logarithm (base 10) and not in
+user units.
+.le
+.ls nx = INDEF, ny = INDEF
+Number of output pixels.  If INDEF and if the pixel interval is also INDEF then
+the number of output pixels is equal to the number of input pixels.
+.le
+.ls xlog = no, ylog = no
+Convert to logarithmic intervals?  If "yes" the output pixel intervals
+are logarithmic.
+.le
+.ls blank = INDEF
+Value to put in the output transformed image when it transforms to regions
+outside the input image.  The special value INDEF will use the nearest
+input pixel which is the behavior before the addition of this parameter.
+Using special blank values allows other software to identify such out
+of input pixels.  See also the \fImoutput\fR parameter to identify
+out of input pixels in pixel masks.
+.le
+.ls logfiles = "STDOUT,logfile"
+List of files in which to keep a log.  If null, "", then no log is kept.
+.le
+.ih
+DESCRIPTION
+The coordinate maps U(X,Y) and V(X,Y), created by the task \fBfitcoords\fR,
+are read from the specified database coordinate fits or from the
+world coordinate system (WCS) of the image.  X and Y are the original
+untransformed pixel coordinates and U and V are the desired output user or
+world coordinates (i.e. slit position and wavelength).  If a coordinate map
+for only one of the user coordinates is given then a one-to-one mapping
+is assumed for the other such that U=X or V=Y.  The coordinate maps are
+inverted to obtain X(U,V) and Y(U,V) on an even subsampled grid of U and
+V over the desired output image coordinates.  The X and Y at each output
+U and V used to interpolate from the input image are found by linear
+interpolation over this grid.  X(U,V) and Y(U,V) are not determined at
+every output point because this is quite slow and is not necessary since
+the coordinate surfaces are relatively slowly varying over the subsampling
+(every 10th output point).
+
+The type of image interpolation is
+selected by the user.  Note that the more accurate the interpolator the
+longer the transformation time required.  The parameter \fIflux\fR selects
+between direct image interpolation and a flux conserving interpolation.
+Flux conservation consists of multiplying the interpolated pixel value by
+the Jacobean of the transformation at that point.  This is essentially
+the ratio of the pixel areas between the output and input images.  Note
+that this is not exact since it is not an integral over the output pixel.
+However, it will be very close except when the output pixel size is much
+greater than the input pixel size.  A log describing the image transformations
+may be kept or printed on the standard output.
+
+The output coordinate grid may be defined by the user or allowed to
+default to an image of the same size as the input image spanning the
+full range of user coordinates in the coordinate transformation maps.
+When the coordinate maps are created by the task \fBfitcoords\fR the
+user coordinates at the corners of the image are recorded in the
+database.  By default these values are used to set the limits of the
+output grid.  If a pixel interval is not specified then an interval
+yielding the specified number of pixels is used.  The default number of
+pixels is that of the input image.  Note that if a pixel interval is
+specified then it takes precedence over the number of pixels.
+
+The pixel intervals may also be logarithmic if the parameter \fIxlog\fR or
+\fIylog\fR is "yes".  Generally, the number of output pixels is specified
+in this case .  However, if the interval is specified it must be a base
+10 logarithmic interval and not in units of the x and y limits which are
+specified in user units.
+
+The transformation from the desired output pixel to the input image may
+fall outside of the input image.  In this case the output pixel may be
+set to the nearest pixel value in the input image or to a particular value
+using the \fIblank\fR parameter.  Also if an output mask is created this
+pixels will have a value of one in the mask.
+
+The parameters \fIminput\fR and \fImoutput\fR provide for input and output
+pixel masks.  An input mask is not used in calculating the transformed
+pixel value but is used to identify the output pixels in the output mask
+which make a significant contribution to the interpolated value.  The
+significance is determined as follows.  The input mask values above zero
+are converted to one hundred.  The mask is then interpolated in the same
+way as the input image.  Any interpolated value of ten or greater is then
+given the value one in the output mask.  This means if all the input pixels
+had mask values of zero a result of zero means no bad pixels were used.
+If all the input pixels had values of 100 then the result will be 100 and
+the output mask will flag this as a bad pixel.  Other values are produced
+by a mixture of good and bad pixels weighted by the interpolation kernel.
+The choice of 10% is purely empirical and gives an approximate identification
+of significant affected pixels.
+zero and
+is created with values of 100
+
+.ih
+EXAMPLES
+Arc calibration images were used to determine a two dimensional dispersion
+map called dispmap.  Stellar spectra were used to determine a two dimensional
+distortion map call distort.  These maps where made using the task
+\fBfitcoords\fR. To transform a set of input images into linear wavelength
+between 3800 and 6400 Angstroms (the user coordinate units) with a dispersion
+of 3 Angstroms per pixel:
+
+.nf
+	cl> transform obj001,obj002 out001,out002 dispmap,distort \
+	>>> y1=3800 y2=6400 dy=3
+.fi
+
+To use logarithmic intervals in the wavelength to yield the same number of
+pixels in the output images as in the input images:
+
+.nf
+	cl> transform obj001,obj002 out001,out002 dispmap,distort \
+	>>> y1=3800 y2=6400 ylog=yes
+.fi
+.ih
+TIMINGS
+The following timings were obtained for transforming a 511x512 real
+image to another 511x512 real image using two Chebyshev transformation
+surface functions (one for the dispersion axis, "henear", and one in
+spatial axis, "object") of order 6 in both dimensions created with the
+task \fBfitcoords\fR.  The times are for a UNIX/VAX 11/750.
+
+.nf
+cl> $transform input output henear,object interp=linear
+TIME (transform)  173.73  5:13  55%
+cl> $transform input output henear,object interp=poly3
+TIME (transform)  266.63  9:17  42%
+cl> $transform input output henear,object interp=spline3
+TIME (transform)  309.05  6:11  83%
+cl> $transform input output henear,object interp=spline3
+TIME (transform)  444.13  9:44  76%
+cl> $transform input output henear interp=linear
+TIME (transform)  171.32  7:24  38%
+cl> $transform input output henear interp=spline3
+TIME (transform)  303.40  12:17  41%
+cl> $transform input output henear,object interp=spline3 flux=no
+TIME (transform)  262.42  10:42  40%
+.fi
+
+The majority of the time is due to the image interpolation and not evaluating
+the transformation functions as indicated by the last three examples.
+.ih
+NOTES
+.ls TRANSFORM: V2.12.2
+The use of bad pixel masks, a specified "blank" value, and use of a WCS
+to resample a WCS calibrated image was added.
+.le
+.ls TRANSFORM: V2.6
+With Version 2.6 of IRAF the algorithm used to invert the user
+coordinate surfaces, U(X,Y) and V(X,Y) to X(U,V) and Y(U,V), has been
+changed.  Previously surfaces of comparable order to the original
+surfaces were fit to a grid of points, i.e. (U(X,Y), V(X,Y), X) and
+(U(X,Y), V(X,Y), Y), with the same surface fitting routines used in
+\fBfitcoords\fR to obtain the input user coordinate surfaces.  This
+method of inversion worked well in all cases in which reasonable
+distortions and dispersions were used.  It was selected because it was
+relatively fast.  However, it cannot be proved to work in all cases; in
+one instance in which an invalid surface was used the inversion was
+actually much poorer than expected.  Therefore a more direct iterative
+inversion algorithm is now used.  This is guaranteed to give the
+correct inversion to within a set error (0.05 of a pixel in X and Y).
+It is slightly slower than the previous algorithm but it is still not
+as major a factor as the image interpolation itself.
+.le
+.ih
+SEE ALSO
+fitcoords
+.endhelp