Initial commit

27 files changed, 3755 insertions, 0 deletions

diff --git a/noao/twodspec/longslit/transform/Notes b/noao/twodspec/longslit/transform/Notes
new file mode 100644
index 00000000..16f5a7a3
--- /dev/null
+++ b/noao/twodspec/longslit/transform/Notes
@@ -0,0 +1,6 @@
+May 29, 1987
+
+If a user accidentally leaves the user coordinate as INDEF in tracing
+the spatial distortion then FITCOORDS uses the fitted coordinate
+which is the same as the pixel coordinate.  This causes incorrect
+results.  Some thought should be given to this situation.
diff --git a/noao/twodspec/longslit/transform/fcdbio.x b/noao/twodspec/longslit/transform/fcdbio.x
new file mode 100644
index 00000000..caf4ac5d
--- /dev/null
+++ b/noao/twodspec/longslit/transform/fcdbio.x
@@ -0,0 +1,99 @@
+include	<error.h>
+include	<math/gsurfit.h>
+include	<pkg/dttext.h>
+include	<units.h>
+
+# FC_DBWRITE -- Write an fitcoords database entry.
+
+procedure fc_dbwrite (database, fitname, axis, un, sf)
+
+char	database[ARB]		# Database
+char	fitname[ARB]		# Database fit name
+int	axis			# Axis for surface
+pointer	un			# Units pointer
+pointer	sf			# Surface pointer
+
+int	i, nsave
+pointer	dt, coeffs, sp, dbfile
+
+int	xgsgeti()
+pointer	dtmap1()
+
+begin
+	if (sf == NULL)
+	    return
+
+	call smark (sp)
+	call salloc (dbfile, SZ_FNAME, TY_CHAR)
+	call strcpy ("fc", Memc[dbfile], SZ_FNAME)
+	call imgcluster (fitname, Memc[dbfile+2], SZ_FNAME-2)
+	dt = dtmap1 (database, Memc[dbfile], APPEND)
+
+	call dtptime (dt)
+	call dtput (dt, "begin\t%s\n")
+	    call pargstr (fitname)
+	call dtput (dt, "\ttask\tfitcoords\n")
+	call dtput (dt, "\taxis\t%d\n")
+	    call pargi (axis)
+	if (un != NULL) {
+	    call dtput (dt, "\tunits\t%s\n")
+		call pargstr (UN_UNITS(un))
+	}
+
+	nsave = xgsgeti (sf, GSNSAVE)
+	call salloc (coeffs, nsave, TY_DOUBLE)
+	call xgssave (sf, Memd[coeffs])
+	call dtput (dt, "\tsurface\t%d\n")
+	    call pargi (nsave)
+	do i = 1, nsave {
+	    call dtput (dt, "\t\t%g\n")
+		call pargd (Memd[coeffs+i-1])
+	}
+
+	call sfree (sp)
+	call dtunmap (dt)
+end
+
+
+# LM_DBREAD -- Read an lsmap database entry.
+
+procedure lm_dbread (database, fitname, axis, un, sf)
+
+char	database[ARB]		# Database
+char	fitname[ARB]		# Fit name
+int	axis			# Axis for surface
+pointer	un			# Units pointer
+pointer	sf			# Surface pointer
+
+int	rec, ncoeffs
+pointer	dt, coeffs, sp, dbfile, units
+
+int	dtlocate(), dtgeti()
+pointer	dtmap1(), un_open()
+
+errchk	dtlocate(), dtgeti(), dtgad(), un_open()
+
+begin
+	un = NULL
+	sf = NULL
+	coeffs = NULL
+
+	call smark (sp)
+	call salloc (dbfile, SZ_FNAME, TY_CHAR)
+	call salloc (units, SZ_FNAME, TY_CHAR)
+	call strcpy ("fc", Memc[dbfile], SZ_FNAME)
+	call imgcluster (fitname, Memc[dbfile+2], SZ_FNAME-2)
+	dt = dtmap1 (database, Memc[dbfile], READ_ONLY)
+
+	rec = dtlocate (dt, fitname)
+	axis = dtgeti (dt, rec, "axis")
+	ifnoerr (call dtgstr (dt, rec, "units", Memc[units], SZ_FNAME))
+	    un = un_open (Memc[units])
+	ncoeffs = dtgeti (dt, rec, "surface")
+	call salloc (coeffs, ncoeffs, TY_DOUBLE)
+	call dtgad (dt, rec, "surface", Memd[coeffs], ncoeffs, ncoeffs)
+	call xgsrestore (sf, Memd[coeffs])
+
+	call sfree (sp)
+	call dtunmap (dt)
+end
diff --git a/noao/twodspec/longslit/transform/fcdlist.x b/noao/twodspec/longslit/transform/fcdlist.x
new file mode 100644
index 00000000..7b9816a7
--- /dev/null
+++ b/noao/twodspec/longslit/transform/fcdlist.x
@@ -0,0 +1,91 @@
+include	<mach.h>
+include	<error.h>
+
+# FC_DLIST -- Fit Coordinates Deletion List Procedures.
+
+# FC_DLREAD -- Fit Coordinates Deletion List Read.
+# Read the deletion list file and match points in the list with the data
+# and delete them.
+
+procedure fc_dlread (x, y, w, npts)
+
+real	x[npts]			# First coordinate to match
+real	y[npts]			# Second coordinate to match
+real	w[npts]			# Weight of coordinate
+int	npts			# Number of coordinates
+
+int	i, fd
+real	r
+char	file[SZ_FNAME]
+real	xdel, ydel
+
+int	access(), open(), fscan(), nscan()
+
+begin
+	call clgstr ("deletions", file, SZ_FNAME)
+
+	if (access (file, READ_ONLY, TEXT_FILE) == NO)
+	    return
+
+	fd = open (file, READ_ONLY, TEXT_FILE)
+
+	while (fscan (fd) != EOF) {
+	    call gargr (xdel)
+	    call gargr (ydel)
+
+	    if (nscan() != 2)
+		next
+
+	    do i = 1, npts {
+		r = sqrt ((x[i]-xdel)**2 + (y[i]-ydel)**2)
+		if (r < 10*EPSILONR)
+		    w[i] = 0.
+#	        if (x[i] != xdel)
+#		    next
+#	        if (y[i] != ydel)
+#		    next
+#		w[i] = 0.
+	    }
+	}
+
+	call close (fd)
+end
+
+
+# FC_DLWRITE -- Fit Coordinates Deletion List Write.
+
+procedure fc_dlwrite (x, y, w, npts)
+
+real	x[npts]			# First coordinate to match
+real	y[npts]			# Second coordinate to match
+real	w[npts]			# Weight of coordinate
+int	npts			# Number of coordinates
+
+int	i, fd
+char	file[SZ_FNAME]
+
+int	open()
+
+begin
+	call clgstr ("deletions", file, SZ_FNAME)
+
+	if (file[1] == EOS)
+	    return
+
+	iferr (call delete (file))
+	    ;
+	iferr (fd = open (file, NEW_FILE, TEXT_FILE)) {
+	    call erract (EA_WARN)
+	    return
+	}
+
+	do i = 1, npts {
+	    if (w[i] == 0.) {
+		call fprintf (fd, "%g %g\n")
+		    call pargr (x[i])
+		    call pargr (y[i])
+	    }
+	}
+
+	call close (fd)
+end
diff --git a/noao/twodspec/longslit/transform/fcfitcoords.x b/noao/twodspec/longslit/transform/fcfitcoords.x
new file mode 100644
index 00000000..13943302
--- /dev/null
+++ b/noao/twodspec/longslit/transform/fcfitcoords.x
@@ -0,0 +1,211 @@
+include	<pkg/gtools.h>
+include	<pkg/igsfit.h>
+include	<pkg/xtanswer.h>
+
+# FC_FITCOORDS -- Fit a surface to the user coordinates.
+
+procedure fc_fitcoords (fitname, database, list, logfiles, interactive)
+
+char	fitname[SZ_FNAME]	# Fitname
+char	database[SZ_FNAME]	# Database
+int	list			# List of images
+int	logfiles		# List of log files
+int	interactive		# Interactive?
+
+int	axis			# Axis of surface fit
+pointer	sf			# Surface pointer
+char	logfile[SZ_FNAME], labels[SZ_LINE, IGSPARAMS]
+bool	answer
+int	ncoords, logfd, axes[2]
+real	xmin, xmax, ymin, ymax
+pointer	gp, gplog, gt, coords, title, un
+
+int	imtgetim(), fntgfntb(), open(), igs_geti(), scan()
+real	xgseval()
+pointer	gopen(), gt_init()
+
+errchk	fc_getcoords
+
+begin
+	# Print a header to the log files giving the inputs.  This is
+	# done first so that if one of the logfiles is STDOUT the user
+	# will see that something is happening.
+
+	axis = 0
+	while (fntgfntb (logfiles, logfile, SZ_FNAME) != EOF) {
+	    logfd = open (logfile, APPEND, TEXT_FILE)
+	    call sysid (logfile, SZ_FNAME)
+	    call fprintf (logfd, "\n%s\n")
+		call pargstr (logfile)
+	    call fprintf (logfd, "  Longslit coordinate fit name is %s.\n")
+		call pargstr (fitname)
+	    call fprintf (logfd, "  Longslit database is %s.\n")
+		call pargstr (database)
+	    call fprintf (logfd, "  Features from images:\n")
+	    while (imtgetim (list, logfile, SZ_FNAME) != EOF) {
+		call fprintf (logfd, "    %s\n")
+		    call pargstr (logfile)
+	    }
+	    call imtrew (list)
+	    call close (logfd)
+	}
+	call fntrewb (logfiles)
+
+	# Get the coordinates for the specified images and axis.  The
+	# coordinates are returned in a pointer which must be explicitly
+	# freed.
+
+	call fc_getcoords (database, list, axis, xmin, xmax, ymin, ymax,
+	    coords, ncoords, labels, un)
+
+	# Read points from the deletion list.
+
+	switch (axis) {
+	case 1:
+	    call fc_dlread (Memr[coords+(Z-1)*ncoords],
+		Memr[coords+(Y-1)*ncoords], Memr[coords+(W-1)*ncoords], ncoords)
+	case 2:
+	    call fc_dlread (Memr[coords+(Z-1)*ncoords],
+		Memr[coords+(X-1)*ncoords], Memr[coords+(W-1)*ncoords], ncoords)
+	}
+
+	# Initialize the graphics.
+
+	if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	    call clgstr ("graphics", logfile, SZ_FNAME)
+	    gp = gopen (logfile, NEW_FILE, STDGRAPH)
+	}
+
+	# Set plot log.
+
+	gplog = NULL
+	call clgstr ("plotfile", logfile, SZ_FNAME)
+	if (logfile[1] != EOS) {
+	    logfd = open (logfile, APPEND, BINARY_FILE)
+	    gplog = gopen ("stdplot", APPEND, logfd)
+	} else
+	    gplog = NULL
+
+	gt = gt_init ()
+	call malloc (title, SZ_LINE, TY_CHAR)
+	call sprintf (Memc[title], SZ_LINE,
+	    "Fit User Coordinates to Image Coordinates for %s")
+	    call pargstr (fitname)
+	call gt_sets (gt, GTTITLE, Memc[title])
+	call mfree (title, TY_CHAR)
+
+	# Fit the surface.  The surface is defined over the full range of
+	# image coordinates.
+
+	call igs_setr (IGS_XMIN, xmin)
+	call igs_setr (IGS_XMAX, xmax)
+	call igs_setr (IGS_YMIN, ymin)
+	call igs_setr (IGS_YMAX, ymax)
+
+	switch (axis) {
+	case 1:
+	    if (Memr[coords+ncoords-1] == 1) {
+	        axes[1] = Y
+	        axes[2] = R
+	        call igs_fit2 (sf, gp, gplog, gt, axes, Memr[coords], ncoords,
+		    labels, interactive)
+	    } else {
+	        axes[1] = X
+	        axes[2] = R
+	        call igs_fit1 (sf, gp, gplog, gt, axes, Memr[coords], ncoords,
+		    labels, interactive)
+	    }
+	case 2:
+	    if (Memr[coords+ncoords-1] == 1) {
+	        axes[1] = X
+	        axes[2] = R
+	        call igs_fit3 (sf, gp, gplog, gt, axes, Memr[coords], ncoords,
+		    labels, interactive)
+	    } else {
+	        axes[1] = Y
+	        axes[2] = R
+	        call igs_fit1 (sf, gp, gplog, gt, axes, Memr[coords], ncoords,
+		    labels, interactive)
+	    }
+	}
+
+	# Close graphics.
+
+	if (gp != NULL)
+	    call gclose (gp)
+	if (gplog != NULL) {
+	    call gclose (gplog)
+	    call close (logfd)
+	}
+	call gt_free (gt)
+
+	# Print logs.
+
+	while (fntgfntb (logfiles, logfile, SZ_FNAME) != EOF) {
+	    logfd = open (logfile, APPEND, TEXT_FILE)
+	    call fprintf (logfd,
+		"  Map %s coordinates for axis %d using image features:\n")
+		call pargstr (labels[1, Z])
+		call pargi (axis)
+	    call fprintf (logfd, "  Number of feature coordnates = %d\n")
+		call pargi (ncoords)
+	    call igs_gets (IGS_FUNCTION, logfile, SZ_FNAME)
+	    call fprintf (logfd, "  Mapping function = %s\n")
+		call pargstr (logfile)
+	    call fprintf (logfd, "    X order = %d\n    Y order = %d\n")
+		call pargi (igs_geti (IGS_XORDER))
+		call pargi (igs_geti (IGS_YORDER))
+	    call fprintf (logfd,
+		"  Fitted coordinates at the corners of the images:\n")
+	    call fprintf (logfd, "    (%d, %d) = %g  (%d, %d) = %g\n")
+		call pargr (xmin)
+		call pargr (ymin)
+		call pargr (xgseval (sf, xmin, ymin))
+		call pargr (xmax)
+		call pargr (ymin)
+		call pargr (xgseval (sf, xmax, xmin))
+	    call fprintf (logfd, "    (%d, %d) = %g  (%d, %d) = %g\n")
+		call pargr (xmin)
+		call pargr (ymax)
+		call pargr (xgseval (sf, xmin, ymax))
+		call pargr (xmax)
+		call pargr (ymax)
+		call pargr (xgseval (sf, xmax, ymax))
+	    call close (logfd)
+	}
+	call fntrewb (logfiles)
+
+	# Write the fit to the database.
+
+	answer = true
+	if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	    call printf ("Write coordinate map to the database (yes)? ")
+	    call flush (STDOUT)
+	    if (scan() != EOF)
+	        call gargb (answer)
+	}
+	if (answer)
+	    call fc_dbwrite (database, fitname, axis, un, sf)
+
+	# Write list of deleted points.
+
+	if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	    switch (axis) {
+	    case 1:
+		call fc_dlwrite (Memr[coords+(Z-1)*ncoords],
+		    Memr[coords+(Y-1)*ncoords],
+		    Memr[coords+(W-1)*ncoords], ncoords)
+	    case 2:
+		call fc_dlwrite (Memr[coords+(Z-1)*ncoords],
+		    Memr[coords+(X-1)*ncoords],
+		    Memr[coords+(W-1)*ncoords], ncoords)
+	    }
+	}
+
+	# Free memory.
+
+	call mfree (coords, TY_REAL)
+	if (un != NULL)
+	    call un_close (un)
+	call xgsfree (sf)
+end
diff --git a/noao/twodspec/longslit/transform/fcgetcoords.x b/noao/twodspec/longslit/transform/fcgetcoords.x
new file mode 100644
index 00000000..dda1c0f0
--- /dev/null
+++ b/noao/twodspec/longslit/transform/fcgetcoords.x
@@ -0,0 +1,212 @@
+include	<imio.h>
+include	<mach.h>
+include	<mwset.h>
+include	<pkg/dttext.h>
+include	<pkg/igsfit.h>
+
+# FC_GETCOORDS -- Get feature coordinates for the specified axis and list
+# of images.  Determine the image dimensions.
+
+procedure fc_getcoords (database, list, axis, xmin, xmax, ymin, ymax,
+    coords, ncoords, labels, un)
+
+char	database[ARB]			# Database
+int	list				# List of images
+int	axis				# Image axis
+real	xmin, xmax			# Image X limits
+real	ymin, ymax			# Image Y limits
+pointer	coords				# Coordinate data pointer
+pointer	ncoords				# Number of coordinate points
+char	labels[SZ_LINE,IGSPARAMS]	# Axis labels
+pointer	un				# Units pointer
+
+char	image1[SZ_FNAME], image2[SZ_FNAME], root[SZ_FNAME], units[SZ_FNAME]
+int	i, j, rec, index, imin, imax, nfeatures, ntotal
+real	value, wt, ltm[2,2], ltv[2]
+pointer	dt, im, mw, ct, x, y, user
+
+int	fc_getim(), dtgeti(), dtscan(), mw_stati()
+real	mw_c1tranr()
+bool	strne()
+pointer	dtmap1(), immap(), mw_openim(), mw_sctran(), un_open()
+
+errchk	dtmap1, dtgstr, immap
+
+begin
+	x = NULL
+	ncoords = 0
+	ntotal = 0
+	axis = 0
+	imin = MAX_INT
+	imax = -MAX_INT
+	un = NULL
+
+	while (fc_getim (list, image1, SZ_FNAME) != EOF) {
+	    call strcpy ("id", root, SZ_FNAME)
+	    call imgcluster (image1, root[3], SZ_FNAME-2)
+	    dt = dtmap1 (database, root, READ_ONLY)
+	    do rec = 1, DT_NRECS(dt) {
+
+		iferr (call dtgstr (dt, rec, "task", image2, SZ_FNAME))
+		    next
+		if (strne ("identify", image2))
+		    next
+
+	        call dtgstr (dt, rec, "image", image2, SZ_FNAME)
+		call get_root (image2, root, SZ_FNAME)
+		if (strne (image1, root))
+		    next
+
+		# Map the 1D image section and determine the axis, the
+		# line or column in the 2D image, and the 2D image size.
+
+		im = immap (image2, READ_ONLY, 0)
+		j = IM_VMAP(im, 1)
+		switch (j) {
+		case 1:
+		    index = IM_VOFF (im, 2) + 1
+		case 2:
+		    index = IM_VOFF (im, 1) + 1
+		}
+		imin = min (imin, index)
+		imax = max (imax, index)
+
+		xmin = 1.
+		xmax = IM_SVLEN (im, 1)
+		ymin = 1.
+		ymax = IM_SVLEN (im, 2)
+
+		if (axis == 0)
+		    axis = j
+
+		if (j != axis) {
+		    call imunmap (im)
+		    call eprintf (
+       "Warning: Fit axes don't agree for combine option.  Ignoring %s.\n")
+		       call pargstr (image1)
+		    break
+		}
+
+		# Set the WCS to convert the feature positions from
+		# IDENTIFY/REIDENTIFY which are in "physical" coordinates
+		# to "logical" coordinates currently used by TRANSFORM.
+
+		mw = mw_openim (im)
+		call mw_seti (mw, MW_USEAXMAP, NO)
+		i = mw_stati (mw, MW_NPHYSDIM)
+		call mw_gltermr (mw, ltm, ltv, i)
+		if (ltm[1,1] == 0. && ltm[2,2] == 0.) {
+		    ltm[1,1] = ltm[2,1]
+		    ltm[2,1] = 0.
+		    ltm[2,2] = ltm[1,2]
+		    ltm[1,2] = 0.
+		    call mw_sltermr (mw, ltm, ltv, i)
+		} else if (ltm[1,2] != 0. || ltm[2,1] != 0.) {
+		    ltv[1] = 0.
+		    ltv[2] = 0.
+		    ltm[1,1] = 1.
+		    ltm[2,1] = 0.
+		    ltm[2,2] = 1.
+		    ltm[1,2] = 0.
+		    call mw_sltermr (mw, ltm, ltv, i)
+		}
+		call mw_seti (mw, MW_USEAXMAP, YES)
+		ct = mw_sctran (mw, "physical", "logical", 1)
+
+		# Allocate memory for the feature information and read
+		# the database.
+
+		ifnoerr (call dtgstr (dt, rec, "units", units, SZ_FNAME))
+		    un = un_open (units)
+		nfeatures = dtgeti (dt, rec, "features")
+		if (x == NULL) {
+		    call malloc (x, nfeatures, TY_REAL)
+		    call malloc (y, nfeatures, TY_REAL)
+		    call malloc (user, nfeatures, TY_REAL)
+		} else {
+		    call realloc (x, ncoords+nfeatures, TY_REAL)
+		    call realloc (y, ncoords+nfeatures, TY_REAL)
+		    call realloc (user, ncoords+nfeatures, TY_REAL)
+		}
+
+		do i = 1, nfeatures {
+		    j = dtscan (dt)
+		    call gargr (value)
+		    switch (axis) {
+		    case 1:
+			Memr[x+ncoords] = mw_c1tranr (ct, value)
+			Memr[y+ncoords] = index
+		    case 2:
+			Memr[x+ncoords] = index
+			Memr[y+ncoords] = mw_c1tranr (ct, value)
+		    }
+		    call gargr (value)
+		    call gargr (value)
+		    call gargr (wt)
+		    call gargr (wt)
+		    call gargr (wt)
+		    if (!IS_INDEF (value) && wt > 0.) {
+			Memr[user+ncoords] = value
+			ncoords = ncoords + 1
+		    }
+		    ntotal = ntotal + 1
+		}
+		call mw_close (mw)
+		call imunmap (im)
+	    }
+
+	    # Finish up
+	    call dtunmap (dt)
+	}
+
+	# Set coordinates.  Take error action if no features are found.
+
+	if (ncoords > 0) {
+	    call xt_sort3 (Memr[user], Memr[x], Memr[y], ncoords)
+	    call malloc (coords, ncoords*IGSPARAMS, TY_REAL)
+	    call amovr (Memr[x], Memr[coords+(X-1)*ncoords], ncoords)
+	    call amovr (Memr[y], Memr[coords+(Y-1)*ncoords], ncoords)
+	    call amovr (Memr[user], Memr[coords+(Z-1)*ncoords], ncoords)
+	    call amovkr (1., Memr[coords+(W-1)*ncoords], ncoords)
+
+	    call fc_setfeatures (Memr[coords], Memr[coords+(Z-1)*ncoords],
+		ncoords)
+
+	    call strcpy ("X (pixels)", labels[1,X], SZ_LINE)
+	    call strcpy ("Y (pixels)", labels[1,Y], SZ_LINE)
+	    call strcpy ("User", labels[1,Z], SZ_LINE)
+	    call strcpy ("Surface", labels[1,S], SZ_LINE)
+	    call strcpy ("Residuals", labels[1,R], SZ_LINE)
+	}
+
+	call mfree (x, TY_REAL)
+	call mfree (y, TY_REAL)
+	call mfree (user, TY_REAL)
+
+	if (ncoords == 0) {
+	    if (ntotal == 0)
+		call error (1, "No coordinates found in database")
+	    else
+		call error (1, "Only INDEF coordinates found in database")
+	}
+end
+
+
+# FC_SETFEATURES -- Set the feature numbers.
+
+procedure fc_setfeatures (features, user, npts)
+
+real	features[npts]			# Feature numbers
+real	user[npts]			# User coordinates
+int	npts				# Number of points
+
+int	i
+
+begin
+	features[1] = 1
+	do i = 2, npts {
+	    features[i] = features[i-1]
+	    if (user[i] != user[i-1])
+		features[i] = features[i] + 1
+	}
+end
diff --git a/noao/twodspec/longslit/transform/fcgetim.x b/noao/twodspec/longslit/transform/fcgetim.x
new file mode 100644
index 00000000..e76ba25a
--- /dev/null
+++ b/noao/twodspec/longslit/transform/fcgetim.x
@@ -0,0 +1,32 @@
+# FC_GETIM -- Get next image name with standard image extensions removed.
+# This is necessary to avoid having two legal image names refering to the
+# same image.
+
+int procedure fc_getim (list, image, maxchar)
+
+int	list		# Image list
+char	image[maxchar]	# Image name
+int	maxchar		# Maximum number of chars in image name
+
+int	i, stat, imtgetim(), strmatch()
+
+begin
+	stat = imtgetim (list, image, maxchar)
+
+	if (stat == EOF)
+	    return (stat)
+
+	i = strmatch (image, ".imh")
+	if (i > 0) {
+	    call strcpy (image[i], image[i-4], maxchar)
+	    return (stat)
+	}
+
+	i = strmatch (image, ".hhh")
+	if (i > 0) {
+	    call strcpy (image[i], image[i-4], maxchar)
+	    return (stat)
+	}
+
+  	return (stat)
+end
diff --git a/noao/twodspec/longslit/transform/fitcoords.x b/noao/twodspec/longslit/transform/fitcoords.x
new file mode 100644
index 00000000..e849caf2
--- /dev/null
+++ b/noao/twodspec/longslit/transform/fitcoords.x
@@ -0,0 +1,83 @@
+include	<error.h>
+include	<pkg/igsfit.h>
+include	<pkg/xtanswer.h>
+
+# T_FITCOORDS -- Fit a surface to the coordinates of longslit images.
+#
+# This is the CL entry for this task.  All the real work is done by
+# fc_fitcoords.
+
+procedure t_fitcoords ()
+
+int	list1			# Image list
+char	fitname[SZ_FNAME]	# Database name for coordinate fit
+char	database[SZ_FNAME]	# Database
+int	logfiles		# List of log files
+bool	combine			# Combine input data?
+int	interactive		# Interactive?
+
+char	image[SZ_FNAME], prompt[SZ_LINE]
+int	list2
+
+int	clgeti(), clpopnu(), imtopen(), fc_getim()
+bool	clgetb()
+
+begin
+	# Get the task parameters.
+
+	call clgstr ("fitname", fitname, SZ_FNAME)
+	call xt_stripwhite (fitname)
+	combine = clgetb ("combine")
+
+	if (combine && (fitname[1] == EOS))
+	    call error (0, "Fit name not specified")
+
+	call clgstr ("images", prompt, SZ_LINE)
+	list1 = imtopen (prompt)
+	call clgstr ("database", database, SZ_FNAME)
+	logfiles = clpopnu ("logfiles")
+	if (clgetb ("interactive"))
+	    interactive = YES
+	else
+	    interactive = ALWAYSNO
+
+	# Set the initial surface in the igsfit package.
+
+	call clgstr ("function", prompt, SZ_LINE)
+	call igs_sets (IGS_FUNCTION, prompt)
+	call igs_seti (IGS_XORDER, clgeti ("xorder"))
+	call igs_seti (IGS_YORDER, clgeti ("yorder"))
+
+	# For each fit ask the user whether to do the fit interactively.
+	# If combining the coordinates from all the images in the
+	# input list then pass the list directly to fc_fitcoords.
+	# Otherwise for each image in the list create a second list
+	# containing just that image.  A second list is needed because
+	# fc_fitcoords expects a list.
+
+	if (combine) {
+	    call sprintf (prompt, SZ_LINE, "Fit interactively")
+	    call xt_answer (prompt, interactive)
+	    call fc_fitcoords (fitname, database, list1, logfiles, interactive)
+
+	} else {
+	    while (fc_getim (list1, image, SZ_FNAME) != EOF) {
+	        list2 = imtopen (image)
+	        call sprintf (prompt, SZ_LINE, "Fit %s interactively")
+		    call pargstr (image)
+	        call xt_answer (prompt, interactive)
+		call sprintf (prompt, SZ_LINE, "%s%s")
+		    call pargstr (fitname)
+		    call pargstr (image)
+	        iferr (call fc_fitcoords (prompt, database, list2, logfiles,
+		    interactive))
+		    call erract (EA_WARN)
+	        call imtclose (list2)
+	    }
+	}
+
+	# Finish up.
+
+	call clpcls (logfiles)
+	call imtclose (list1)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/Revisions b/noao/twodspec/longslit/transform/igsfit/Revisions
new file mode 100644
index 00000000..92b36cca
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/Revisions
@@ -0,0 +1,42 @@
+.help revisions Jun88 noao.twodspec.longslit.transform.igsfit
+.nf
+ igsfit.x
+ igsnearest.x
+     GSCUR was being called with DOUBLE precision values. (12/22/87)
+
+ igsfit.x
+ igscolon.x
+ igsget.x
+     Added colon options to print fit at corners of surface. (8/10/87 Valdes)
+
+ ====
+ V2.5
+ ====
+
+noao$twodspec/longslit/transform/igsfit/*.x
+    Valdes, February 17, 1987
+    1.  GIO changes.
+
+noao$twodspec/longslit/transform/igsfit/igsfit.x
+noao$twodspec/longslit/transform/igsfit/igscolon.x
+    Valdes, January 16, 1987
+    1.  '?' now uses system page facility.
+    2.  Colon command dictionary and switch modified to use macro definitions.
+
+noao$twodspec/longslit/transform/igsfit/igsdelete.x
+noao$twodspec/longslit/transform/igsfit/igsundelete.x
+    Valdes, October 16, 1986
+    1.  Real line type specified in gseti call changed to integer.
+	This caused a crash on AOS/IRAF.
+
+========================================================
+
+From Valdes on Feb 7, 1986:
+
+1.  Bug fixed in deleting and undeleting points.
+------
+From Valdes on Jan 3, 1986:
+
+1.  Modified IGSFIT to allow zooming on constant x, constant y, constant z,
+and constant feature.
+.endhelp
diff --git a/noao/twodspec/longslit/transform/igsfit/igscolon.x b/noao/twodspec/longslit/transform/igsfit/igscolon.x
new file mode 100644
index 00000000..6847974a
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igscolon.x
@@ -0,0 +1,115 @@
+include	<gset.h>
+
+# List of colon commands
+define	CMDS "|show|function|xorder|yorder|corners|"
+
+define	SHOW		1	# Show parameters
+define	FUNCTION	2	# Set or show function type
+define	XORDER		3	# Set or show x order of function
+define	YORDER		4	# Set or show y order of function
+define	CORNERS		5	# Show corners
+
+# IGS_COLON -- Processes colon commands.
+
+procedure igs_colon (cmdstr, gp, sf)
+
+char	cmdstr[ARB]			# Command string
+pointer	gp				# GIO pointer
+pointer	sf				# Surface pointer
+
+char	cmd[SZ_LINE]
+int	ncmd, ival
+
+int	nscan(), strdic()
+real	xgseval()
+
+string	funcs "|chebyshev|legendre|"
+
+include	"igsfit.com"
+
+begin
+	# Use formated scan to parse the command string.
+	# The first word is the command and it may be minimum match
+	# abbreviated with the list of commands.
+
+	call sscan (cmdstr)
+	call gargwrd (cmd, SZ_LINE)
+	ncmd = strdic (cmd, cmd, SZ_LINE, CMDS)
+
+	switch (ncmd) {
+	case SHOW: # :show - Show the values of the fitting parameters.
+	    call gdeactivate (gp, AW_CLEAR)
+	    call printf ("function %s\n")
+		call pargstr (function)
+	    call printf ("xorder %d\n")
+		call pargi (xorder)
+	    call printf ("yorder %d\n")
+		call pargi (yorder)
+	    call printf ("Fitted coordinates at the corners of the images:\n")
+	    call printf ("    (%d, %d) = %g  (%d, %d) = %g\n")
+		call pargr (xmin)
+		call pargr (ymin)
+		call pargr (xgseval (sf, xmin, ymin))
+		call pargr (xmax)
+		call pargr (ymin)
+		call pargr (xgseval (sf, xmax, xmin))
+	    call printf ("    (%d, %d) = %g  (%d, %d) = %g\n")
+		call pargr (xmin)
+		call pargr (ymax)
+		call pargr (xgseval (sf, xmin, ymax))
+		call pargr (xmax)
+		call pargr (ymax)
+		call pargr (xgseval (sf, xmax, ymax))
+	    call printf ("rms %g\n")
+		call pargr (rms)
+	    call greactivate (gp, AW_PAUSE)
+
+	case FUNCTION: # :function - List or set the fitting function.
+	    call gargwrd (cmd, SZ_LINE)
+	    if (nscan() == 1) {
+		call printf ("function = %s\n")
+		    call pargstr (function)
+	    } else {
+		if (strdic (cmd, cmd, SZ_LINE, funcs) > 0)
+		    call strcpy (cmd, function, SZ_LINE)
+		else
+		    call printf ("Unknown or ambiguous function\n")
+	    }
+
+	case XORDER: # xorder: List or set the function order.
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("xorder %d\n")
+		    call pargi (xorder)
+	    } else if (ival < 2)
+		call printf ("xorder must be at least 2\n")
+	    else
+		xorder = ival
+
+	case YORDER: # yorder: List or set the function order.
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("yorder %d\n")
+		    call pargi (yorder)
+	    } else if (ival < 2)
+		call printf ("yorder must be at least 2\n")
+	    else
+		yorder = ival
+	case CORNERS: # corners: List coordinates at corners.
+	    call printf ("(%d,%d)=%g (%d,%d)=%g (%d,%d)=%g (%d,%d)=%g\n")
+		call pargr (xmin)
+		call pargr (ymin)
+		call pargr (xgseval (sf, xmin, ymin))
+		call pargr (xmax)
+		call pargr (ymin)
+		call pargr (xgseval (sf, xmax, xmin))
+		call pargr (xmin)
+		call pargr (ymax)
+		call pargr (xgseval (sf, xmin, ymax))
+		call pargr (xmax)
+		call pargr (ymax)
+		call pargr (xgseval (sf, xmax, ymax))
+	default:
+	    call printf ("Unrecognized or ambiguous command\007")
+	}
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsdelete.x b/noao/twodspec/longslit/transform/igsfit/igsdelete.x
new file mode 100644
index 00000000..3de2fb25
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsdelete.x
@@ -0,0 +1,103 @@
+include	<mach.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	<pkg/igsfit.h>
+
+# IGS_NEARESTD -- Nearest point to delete.
+
+int procedure igs_nearestd (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+
+pointer	gp			# GIO pointer
+int	ztype			# Zoom type
+int	refpt			# Reference point
+int	axis[2]			# Axes
+real	pts[npts, ARB]		# Data points
+int	npts			# Number of data points
+real	wx, wy			# Cursor coordinates
+int	wcs			# WCS
+
+int	i, j, x, y
+real	r2, r2min, x0, y0
+
+begin
+	x = axis[1]
+	y = axis[2]
+
+	call gctran (gp, wx, wy, wx, wy, wcs, 0)
+	r2min = MAX_REAL
+	j = 0
+
+	if (IS_INDEFI (ztype)) {
+	    do i = 1, npts {
+	        if (pts[i,W] == 0.)
+		    next
+	        call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+	} else {
+	    do i = 1, npts {
+	        if ((pts[i,ztype] != pts[refpt,ztype]) || (pts[i,W] == 0.))
+		    next
+	        call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+	}
+
+	return (j)
+end
+
+# IGS_DELETE -- Delete points or subsets.
+
+procedure igs_delete (gp, gt, ztype, refpt, axis, pts, npts, dtype)
+
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+int	ztype			# Zoom type
+int	refpt			# Reference point for deletion
+int	axis[2]			# Axes
+real	pts[npts, ARB]		# Data points
+int	npts			# Number of data points
+int	dtype			# Deletion type
+
+int	i, x, y
+real	xsize, ysize
+
+real	gt_getr()
+
+begin
+	x = axis[1]
+	y = axis[2]
+
+	xsize = gt_getr (gt, GTXSIZE)
+	ysize = gt_getr (gt, GTYSIZE)
+
+	switch (dtype) {
+	case X, Y, Z:
+	    do i = 1, npts {
+		if (!IS_INDEFI (ztype))
+		    if (pts[i,ztype] != pts[refpt,ztype])
+		        next
+		if (pts[i,dtype] != pts[refpt,dtype])
+		    next
+	        call gseti (gp, G_PMLTYPE, 0)
+	        call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+	        call gseti (gp, G_PMLTYPE, 1)
+	        call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+		pts[i,W] = 0.
+	    }
+	default:
+	    call gseti (gp, G_PMLTYPE, 0)
+	    call gmark (gp, pts[refpt,x], pts[refpt,y], GM_PLUS, xsize, ysize)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, pts[refpt,x], pts[refpt,y], GM_CROSS, xsize, ysize)
+	    pts[refpt,W] = 0.
+	}
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsfit.com b/noao/twodspec/longslit/transform/igsfit/igsfit.com
new file mode 100644
index 00000000..90bf90aa
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsfit.com
@@ -0,0 +1,10 @@
+# Common parameters.
+
+char	function[SZ_LINE]	# Surface function
+int	xorder			# X order of surface function
+int	yorder			# Y order of surface function
+real	xmin, xmax		# X range
+real	ymin, ymax		# Y range
+real	mean, rms		# Mean and RMS of fit
+
+common	/igscom/ xmin, xmax, ymin, ymax, xorder, yorder, function, mean, rms
diff --git a/noao/twodspec/longslit/transform/igsfit/igsfit.x b/noao/twodspec/longslit/transform/igsfit/igsfit.x
new file mode 100644
index 00000000..14e8e51e
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsfit.x
@@ -0,0 +1,373 @@
+include	<mach.h>
+include	<pkg/gtools.h>
+include <pkg/igsfit.h>
+
+define	HELP	"noao$lib/scr/igsfit.key"
+define	PROMPT	"fitcoords surface fitting options"
+
+
+# IGS_FIT1 -- Fit z = f(x, y)
+
+procedure igs_fit1 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive)
+
+pointer	sf			# GSURFIT pointer
+pointer	gp			# GIO pointer
+pointer	gplog			# GIO pointer for plot log
+pointer	gt			# GTOOLS pointer
+int	axis[2]			# Axis definitions
+real	pts[npts, ARB]		# Data
+int	npts			# Number of pts points
+char	labels[SZ_LINE, ARB]	# Identification labels
+int	interactive		# Interactive?
+
+extern	igs_solve1()
+
+begin
+	call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive,
+	    igs_solve1)
+end
+
+
+# IGS_FIT2 -- Fit z = x + f(y)
+
+procedure igs_fit2 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive)
+
+pointer	sf			# GSURFIT pointer
+pointer	gp			# GIO pointer
+pointer	gplog			# GIO pointer for plot log
+pointer	gt			# GTOOLS pointer
+int	axis[2]			# Axis definitions
+real	pts[npts, ARB]		# Data
+int	npts			# Number of pts points
+char	labels[SZ_LINE, ARB]	# Identification labels
+int	interactive		# Interactive?
+
+extern	igs_solve2()
+
+begin
+	call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive,
+	    igs_solve2)
+end
+
+
+# IGS_FIT3 -- Fit z = y + f(x)
+
+procedure igs_fit3 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive)
+
+pointer	sf			# GSURFIT pointer
+pointer	gp			# GIO pointer
+pointer	gplog			# GIO pointer for plot log
+pointer	gt			# GTOOLS pointer
+int	axis[2]			# Axis definitions
+real	pts[npts, ARB]		# Data
+int	npts			# Number of pts points
+char	labels[SZ_LINE, ARB]	# Identification labels
+int	interactive		# Interactive?
+
+extern	igs_solve3()
+
+begin
+	call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive,
+	    igs_solve3)
+end
+
+
+# IGS_FIT -- Interactive surface fitting.
+
+procedure igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive, igs_solve)
+
+pointer	sf			# GSURFIT pointer
+pointer	gp			# GIO pointer
+pointer	gplog			# GIO pointer for plot log
+pointer	gt			# GTOOLS pointer
+int	axis[2]			# Axis definitions
+real	pts[npts, ARB]		# Data
+int	npts			# Number of pts points
+char	labels[SZ_LINE, ARB]	# Identification labels
+int	interactive		# Interactive?
+extern	igs_solve()		# Surface solution routine
+
+int	i, newgraph, ztype, dtype, refpt, refpt1
+real	zval, zval1
+pointer	wts
+
+real	wx, wy
+int	wcs, key
+char	cmd[SZ_LINE]
+
+int	clgcur(), gt_gcur(), igs_nearest(), igs_nearestd(), igs_nearestu()
+errchk	igs_solve
+
+include	"igsfit.com"
+
+begin
+	# Compute a solution and set the residuals.
+
+	call igs_solve (sf, pts[1,X], pts[1,Y], pts[1,Z], pts[1,W], npts)
+	call xgsvector (sf, pts[1,X], pts[1,Y], pts[1,S], npts)
+	call asubr (pts[1,Z], pts[1,S], pts[1,R], npts)
+	call aavgr (pts[1,R], npts, mean, rms)
+	call igs_params (gt)
+
+	# Return if not interactive.
+
+	ztype = INDEFI
+	if ((gp == NULL) || (interactive == NO))
+	    goto 30
+
+	call malloc (wts, npts, TY_REAL)
+	call amovr (pts[1,W],  Memr[wts], npts)
+
+	call igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels)
+	newgraph = NO
+
+	# Read cursor commands.
+
+10	while (gt_gcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) != EOF) {
+	    switch (key) {
+	    case '?':
+		# Print help text.
+
+		call gpagefile (gp, HELP, PROMPT)
+
+	    case ':':
+		# List or set parameters
+
+		if (cmd[1] == '/')
+		    call gt_colon (cmd, gp, gt, newgraph)
+		else
+		    call igs_colon (cmd, gp, sf)
+
+	    # Set abscissa
+
+	    case 'x':
+		call printf ("Select abscissa (x, y, z, s, r): ")
+		if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+		    goto 10
+		call printf ("\n")
+
+		switch (key) {
+		case 'x':
+		    i = X
+		case 'y':
+		    i = Y
+		case 'z':
+		    i = Z
+		case 's':
+		    i = S
+		case 'r':
+		    i = R
+		default:
+		    call printf ("\07\n")
+		    goto 10
+		}
+
+		if (axis[1] != i) {
+		    axis[1] = i
+		    call gt_setr (gt, GTXMIN, INDEF)
+		    call gt_setr (gt, GTXMAX, INDEF)
+		}
+
+	    # Set ordinate
+
+	    case 'y':
+		call printf ("Select ordinate (x, y, z, s, r): ")
+		if(clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+		    goto 10
+		call printf ("\n")
+
+		switch (key) {
+		case 'x':
+		    i = X
+		case 'y':
+		    i = Y
+		case 'z':
+		    i = Z
+		case 's':
+		    i = S
+		case 'r':
+		    i = R
+		default:
+		    call printf ("\07\n")
+		    goto 10
+		}
+
+		if (axis[2] != i) {
+		    axis[2] = i
+		    call gt_setr (gt, GTYMIN, INDEF)
+		    call gt_setr (gt, GTYMAX, INDEF)
+		}
+
+	    case 'r':
+		newgraph = YES
+
+	    case 'z':
+		if (IS_INDEFI (ztype)) {
+		    refpt = igs_nearest (gp, ztype, refpt, axis, pts, npts, wx,
+			wy, wcs)
+
+		    call printf ("Zoom type (x, y, z): ")
+		    if (clgcur ("cursor",wx,wy,wcs,key,cmd,SZ_LINE) == EOF)
+			goto 10
+		    call printf ("\n")
+
+		    switch (key) {
+		    case 'x':
+		        ztype = X
+		    case 'y':
+		        ztype = Y
+		    case 'z':
+		        ztype = Z
+		    default:
+		        call printf ("\07\n")
+		        goto 10
+		    }
+
+		    newgraph = YES
+		}
+
+	    case 'p':
+		if (!IS_INDEFI (ztype)) {
+		    ztype = INDEFI
+		    newgraph = YES
+		}
+
+	    case 'l':
+		if (!IS_INDEFI (ztype)) {
+		    refpt1 = 0
+		    zval = pts[refpt, ztype]
+		    zval1 = -MAX_REAL
+		    do i = 1, npts {
+			if ((pts[i,ztype] < zval) && (pts[i,ztype] > zval1)) {
+			    refpt1 = i
+			    zval1 = pts[refpt1,ztype]
+			}
+		    }
+
+		    if (refpt1 != 0) {
+			refpt = refpt1
+		        newgraph = YES
+		    }
+		}
+
+	    case 'n':
+		if (!IS_INDEFI (ztype)) {
+		    refpt1 = 0
+		    zval = pts[refpt, ztype]
+		    zval1 = MAX_REAL
+		    do i = 1, npts {
+			if ((pts[i,ztype] > zval) && (pts[i,ztype] < zval1)) {
+			    refpt1 = i
+			    zval1 = pts[refpt1,ztype]
+			}
+		    }
+
+		    if (refpt1 != 0) {
+			refpt = refpt1
+		        newgraph = YES
+		    }
+		}
+
+	    case 'c':
+		# cursor read
+		i = igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+		call printf ("%g %g %g %g %g %g\n")
+		    call pargr (pts[i, X])
+		    call pargr (pts[i, Y])
+		    call pargr (pts[i, Z])
+		    call pargr (pts[i, W])
+		    call pargr (pts[i, S])
+		    call pargr (pts[i, R])
+
+	    case 'd':
+		i = igs_nearestd (gp, ztype, refpt, axis, pts, npts, wx, wy,
+		    wcs)
+		if (i == 0)
+		    goto 10
+
+	        call gscur (gp, real (pts[i,axis[1]]), real (pts[i,axis[2]]))
+
+		call printf ( "Delete 'p'oint or constant 'x', 'y', or 'z': ")
+		if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+		    goto 10
+		call printf ("\n")
+
+		switch (key) {
+		case 'p':
+		    dtype = 0
+		case 'x':
+		    dtype = X
+		case 'y':
+		    dtype = Y
+		case 'z':
+		    dtype = Z
+		default:
+		    call printf ("\07\n")
+		    goto 10
+		}
+
+		call igs_delete (gp, gt, ztype, i, axis, pts, npts, dtype)
+
+	    case 'u':
+		i = igs_nearestu (gp, ztype, refpt, axis, pts, npts, wx, wy,
+		    wcs)
+		if (i == 0)
+		    goto 10
+
+	        call gscur (gp, real (pts[i,axis[1]]), real (pts[i,axis[2]]))
+
+		call printf ( "Undelete 'p'oint or constant 'x', 'y', or 'z': ")
+		if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+		    goto 10
+		call printf ("\n")
+
+		switch (key) {
+		case 'p':
+		    dtype = 0
+		case 'x':
+		    dtype = X
+		case 'y':
+		    dtype = Y
+		case 'z':
+		    dtype = Z
+		default:
+		    call printf ("\07\n")
+		    goto 10
+		}
+
+		call igs_undelete (gp, gt, ztype, i, axis, pts, Memr[wts],
+		    npts, dtype)
+
+	    case 'f':
+		#call printf ("Fitting ...")
+		#call flush (STDOUT)
+		call igs_solve (sf,pts[1,X],pts[1,Y],pts[1,Z],pts[1,W],npts)
+		call xgsvector (sf, pts[1,X], pts[1,Y], pts[1,S], npts)
+		call asubr (pts[1,Z], pts[1,S], pts[1,R], npts)
+		call aavgr (pts[1,R], npts, mean, rms)
+		call igs_params (gt)
+		newgraph = YES
+
+	    case 'w':
+		call gt_window (gt, gp, "cursor", newgraph)
+
+	    case 'I':
+		call fatal (0, "Interrupt")
+
+	    default:
+		# Ring the bell.
+
+		call printf ("\07\n")
+	    }
+
+	    if (newgraph == YES) {
+		call igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels)
+		newgraph = NO
+	    }
+	}
+
+	call mfree (wts, TY_REAL)
+
+30	call igs_graph (gplog, gt, ztype, refpt, axis, pts, npts, labels)
+
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsget.x b/noao/twodspec/longslit/transform/igsfit/igsget.x
new file mode 100644
index 00000000..ccd1fb6c
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsget.x
@@ -0,0 +1,62 @@
+include	<pkg/igsfit.h>
+
+# IGS_GETI -- Get the value of an integer parameter.
+
+int procedure igs_geti (param)
+
+int	param			# IGS parameter
+
+include	"igsfit.com"
+
+begin
+	switch (param) {
+	case IGS_XORDER:
+	    return (xorder)
+	case IGS_YORDER:
+	    return (yorder)
+	default:
+	    call error (0, "igs_geti: Unknown parameter")
+	}
+end
+
+
+# IGS_GETS -- Get the value of a string parameter.
+
+procedure igs_gets (param, str, maxchar)
+
+int	param			# IGS parameter
+char	str[maxchar]		# String
+int	maxchar			# Maximum number of characters
+
+include	"igsfit.com"
+
+begin
+	switch (param) {
+	case IGS_FUNCTION:
+	    call strcpy (function, str, maxchar)
+	default:
+	    call error (0, "igs_gets: Unknown parameter")
+	}
+end
+
+
+# IGS_GETR -- Get the values of real valued fitting parameters.
+
+real procedure igs_getr (param)
+
+int	param			# Parameter to be get
+
+include	"igsfit.com"
+
+begin
+	switch (param) {
+	case IGS_XMIN:
+	    return (xmin)
+	case IGS_XMAX:
+	    return (xmax)
+	case IGS_YMIN:
+	    return (ymin)
+	case IGS_YMAX:
+	    return (ymax)
+	}
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsgraph.x b/noao/twodspec/longslit/transform/igsfit/igsgraph.x
new file mode 100644
index 00000000..83eba7e1
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsgraph.x
@@ -0,0 +1,73 @@
+include	<mach.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	<pkg/igsfit.h>
+
+procedure igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels)
+
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+int	ztype			# Zoom type
+int	refpt			# Reference point
+int	axis[2]			# Axis definitions
+real	pts[npts, ARB]		# Data
+int	npts			# Number of pts points
+char	labels[SZ_LINE, ARB]	# Data labels
+
+int	i, x, y
+real	xmin, xmax, ymin, ymax, xsize, ysize, gt_getr()
+
+begin
+	if (gp == NULL)
+	    return
+
+	x = axis[1]
+	y = axis[2]
+
+	call gt_sets (gt, GTXLABEL, labels[1, x])
+	call gt_sets (gt, GTYLABEL, labels[1, y])
+	xsize = gt_getr (gt, GTXSIZE)
+	ysize = gt_getr (gt, GTYSIZE)
+
+	call gclear (gp)
+
+	if (IS_INDEFI (ztype)) {
+	    call gascale (gp, pts[1, x], npts, 1)
+	    call gascale (gp, pts[1, y], npts, 2)
+	} else {
+	    xmin = MAX_REAL
+	    xmax = -MAX_REAL
+	    ymin = MAX_REAL
+	    ymax = -MAX_REAL
+	    do i = 1, npts {
+		if (pts[i,ztype] != pts[refpt,ztype])
+		    next
+		xmin = min (xmin, pts[i,x])
+		xmax = max (xmax, pts[i,x])
+		ymin = min (ymin, pts[i,y])
+		ymax = max (ymax, pts[i,y])
+	    }
+	    call gswind (gp, xmin, xmax, ymin, ymax)
+	}
+
+	call gt_swind (gp, gt)
+	call gt_labax (gp, gt)
+
+	if (IS_INDEFI (ztype)) {
+	    do i = 1, npts {
+	        if (pts[i,W] == 0.)
+		    call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+	        else
+		    call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+	    }
+	} else {
+	    do i = 1, npts {
+		if (pts[i,ztype] != pts[refpt,ztype])
+		    next
+	        if (pts[i,W] == 0.)
+		    call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+	        else
+		    call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+	    }
+	}
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsinit.x b/noao/twodspec/longslit/transform/igsfit/igsinit.x
new file mode 100644
index 00000000..f084e7ff
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsinit.x
@@ -0,0 +1,21 @@
+include	<pkg/igsfit.h>
+
+# IGS_INIT -- Initialize the surface fitting parameters.
+
+procedure igs_init (function, xorder, yorder, xmin, xmax, ymin, ymax)
+
+char	function[ARB]		# Function
+int	xorder			# X order
+int	yorder			# Y order
+real	xmin, xmax		# X range
+real	ymin, ymax		# Y range
+
+begin
+	call igs_sets (IGS_FUNCTION, function)
+	call igs_seti (IGS_XORDER, xorder)
+	call igs_seti (IGS_YORDER, yorder)
+	call igs_setr (IGS_XMIN, xmin)
+	call igs_setr (IGS_XMAX, xmax)
+	call igs_setr (IGS_YMIN, ymin)
+	call igs_setr (IGS_YMAX, ymax)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsnearest.x b/noao/twodspec/longslit/transform/igsfit/igsnearest.x
new file mode 100644
index 00000000..69888509
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsnearest.x
@@ -0,0 +1,51 @@
+include	<mach.h>
+include	<gset.h>
+include	<pkg/igsfit.h>
+
+int procedure igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+
+pointer	gp			# GIO pointer
+int	ztype			# Zoom type
+int	refpt			# Reference point
+int	axis[2]			# Axes
+real	pts[npts, ARB]		# Data points
+int	npts			# Number of data points
+real	wx, wy			# Cursor coordinates
+int	wcs			# WCS
+
+int	i, j, x, y
+real	r2, r2min, x0, y0
+
+begin
+	x = axis[1]
+	y = axis[2]
+
+	call gctran (gp, wx, wy, wx, wy, wcs, 0)
+	r2min = MAX_REAL
+	j = 0
+
+	if (IS_INDEFI (ztype)) {
+	    do i = 1, npts {
+	        call gctran (gp, pts[i,x], pts[i,y], x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+	} else {
+	    do i = 1, npts {
+		if (pts[i,ztype] != pts[refpt,ztype])
+		    next
+	        call gctran (gp, pts[i,x], pts[i,y], x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+	}
+
+	call gscur (gp, real (pts[j,x]), real (pts[j,y]))
+	return (j)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsparams.x b/noao/twodspec/longslit/transform/igsfit/igsparams.x
new file mode 100644
index 00000000..9ecdd422
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsparams.x
@@ -0,0 +1,23 @@
+include	<pkg/gtools.h>
+
+# IGS_PARAMS -- Set the GTOOLS parameter string.
+
+procedure igs_params (gt)
+
+pointer	gt		# GTOOLS pointer
+
+pointer	params
+
+include	"igsfit.com"
+
+begin
+	call malloc (params, SZ_LINE, TY_CHAR)
+	call sprintf (Memc[params], SZ_LINE,
+	    "Function = %s, xorder = %d, yorder = %d, rms = %.4g")
+	    call pargstr (function)
+	    call pargi (xorder)
+	    call pargi (yorder)
+	    call pargr (rms)
+	call gt_sets (gt, GTPARAMS, Memc[params])
+	call mfree (params, TY_CHAR)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsset.x b/noao/twodspec/longslit/transform/igsfit/igsset.x
new file mode 100644
index 00000000..ea74e8c9
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsset.x
@@ -0,0 +1,59 @@
+include	<pkg/igsfit.h>
+
+# IGS_SETS -- Set the values of string valued fitting parameters.
+
+procedure igs_sets (param, str)
+
+int	param			# Parameter to be set
+char	str[ARB]		# String value
+
+include	"igsfit.com"
+
+begin
+	switch (param) {
+	case IGS_FUNCTION:
+	    call strcpy (str, function, SZ_LINE)
+	}
+end
+
+
+# IGS_SETI -- Set the values of integer valued fitting parameters.
+
+procedure igs_seti (param, ival)
+
+int	param			# Parameter to be set
+int	ival			# Integer value
+
+include	"igsfit.com"
+
+begin
+	switch (param) {
+	case IGS_XORDER:
+	    xorder = ival
+	case IGS_YORDER:
+	    yorder = ival
+	}
+end
+
+
+# IGS_SETR -- Set the values of real valued fitting parameters.
+
+procedure igs_setr (param, rval)
+
+int	param			# Parameter to be set
+real	rval			# Real value
+
+include	"igsfit.com"
+
+begin
+	switch (param) {
+	case IGS_XMIN:
+	    xmin = rval
+	case IGS_XMAX:
+	    xmax = rval
+	case IGS_YMIN:
+	    ymin = rval
+	case IGS_YMAX:
+	    ymax = rval
+	}
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igssolve.x b/noao/twodspec/longslit/transform/igsfit/igssolve.x
new file mode 100644
index 00000000..a7e39354
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igssolve.x
@@ -0,0 +1,173 @@
+include	<math/gsurfit.h>
+
+
+# IGS_SOLVE1 -- Fit z = f(x, y).
+
+define	SFTYPES		"|chebyshev|legendre|"		# Surface types
+
+procedure igs_solve1 (sf, x, y, z, w, npts)
+
+pointer	sf			# GSURFIT pointer
+real	x[npts]			# X points
+real	y[npts]			# Y points
+real	z[npts]			# Z points
+real	w[npts]			# Weights
+int	npts			# Number of points
+
+int	i, nfunc, ix, iy
+pointer	sf1, sf2, resids
+
+int	strdic()
+
+include	"igsfit.com"
+
+begin
+	# Determine the function type.
+
+	nfunc = strdic (function, function, SZ_LINE, SFTYPES)
+
+	# Fit the first surface.
+
+	ix = min (2, xorder)
+	iy = min (2, yorder)
+	call xgsinit (sf1, nfunc, ix, iy, NO, xmin, xmax, ymin, ymax)
+	call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
+
+	switch (i) {
+	case SINGULAR:
+	    call eprintf ("Singular solution\n")
+	case NO_DEG_FREEDOM:
+	    call error (0, "No degrees of freedom")
+	}
+
+	# Evaluate the first surface and fit the residuals.
+
+	call malloc (resids, npts, TY_REAL)
+	call xgsvector (sf1, x, y, Memr[resids], npts)
+	call asubr (z, Memr[resids], Memr[resids], npts)
+
+	call xgsinit (sf2, nfunc, xorder, yorder, YES, xmin,xmax,ymin,ymax)
+	call xgsfit (sf2, x, y, Memr[resids], w, npts, WTS_USER, i)
+
+	switch (i) {
+	case SINGULAR:
+	    call eprintf ("Singular solution\n")
+	case NO_DEG_FREEDOM:
+	    call error (0, "No degrees of freedom")
+	}
+
+	# Add the two surfaces and free memory.
+
+	call xgsadd (sf1, sf2, sf)
+	call xgsfree (sf1)
+	call xgsfree (sf2)
+	call mfree (resids, TY_REAL)
+end
+
+
+# IGS_SOLVE2 -- Fit z = x + f(y).
+
+
+procedure igs_solve2 (sf, x, y, z, w, npts)
+
+pointer	sf			# GSURFIT pointer
+real	x[npts]			# X points
+real	y[npts]			# Y points
+real	z[npts]			# Z points
+real	w[npts]			# Weights
+int	npts			# Number of points
+
+int	i, nfunc
+real	a
+pointer	sf1
+
+int	strdic()
+real	xgsgcoeff()
+
+include	"igsfit.com"
+
+begin
+	nfunc = strdic (function, function, SZ_LINE, SFTYPES)
+	call xgsinit (sf1, nfunc, 1, yorder, NO, xmin, xmax, ymin, ymax)
+
+	call asubr (z, x, z, npts)
+	call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
+	call aaddr (z, x, z, npts)
+
+	switch (i) {
+	case SINGULAR:
+	    call eprintf ("Singular solution\n")
+	case NO_DEG_FREEDOM:
+	    call error (0, "No degrees of freedom")
+	}
+
+	call xgsfree (sf)
+	call xgsinit (sf, nfunc, 2, yorder, NO, xmin, xmax, ymin, ymax)
+	a = xgsgcoeff (sf1, 1, 1)
+
+	a = a + (xmin + xmax) / 2
+	call xgsscoeff (sf, 1, 1, a)
+
+	a = (xmax - xmin) / 2
+	call xgsscoeff (sf, 2, 1, a)
+
+	do i = 2, yorder {
+	    a = xgsgcoeff (sf1, 1, i)
+	    call xgsscoeff (sf, 1, i, a)
+	}
+
+	call xgsfree (sf1)
+end
+
+# IGS_SOLVE3 -- Fit z = y + f(x).
+
+procedure igs_solve3 (sf, x, y, z, w, npts)
+
+pointer	sf			# GSURFIT pointer
+real	x[npts]			# X points
+real	y[npts]			# Y points
+real	z[npts]			# Z points
+real	w[npts]			# Weights
+int	npts			# Number of points
+
+int	i, nfunc
+real	a
+pointer	sf1
+
+int	strdic()
+real	xgsgcoeff()
+
+include	"igsfit.com"
+
+begin
+	nfunc = strdic (function, function, SZ_LINE, SFTYPES)
+	call xgsinit (sf1, nfunc, xorder, 1, NO, xmin, xmax, ymin, ymax)
+
+	call asubr (z, y, z, npts)
+	call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
+	call aaddr (z, y, z, npts)
+
+	switch (i) {
+	case SINGULAR:
+	    call eprintf ("Singular solution\n")
+	case NO_DEG_FREEDOM:
+	    call error (0, "No degrees of freedom")
+	}
+
+	call xgsfree (sf)
+	call xgsinit (sf, nfunc, xorder, 2, NO, xmin, xmax, ymin, ymax)
+	a = xgsgcoeff (sf1, 1, 1)
+
+	a = a + (ymin + ymax) / 2
+	call xgsscoeff (sf, 1, 1, a)
+
+	a = (ymax - ymin) / 2
+	call xgsscoeff (sf, 1, 2, a)
+
+	do i = 2, xorder {
+	    a = xgsgcoeff (sf1, i, 1)
+	    call xgsscoeff (sf, i, 1, a)
+	}
+
+	call xgsfree (sf1)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsundelete.x b/noao/twodspec/longslit/transform/igsfit/igsundelete.x
new file mode 100644
index 00000000..dc7b802e
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsundelete.x
@@ -0,0 +1,107 @@
+include	<mach.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	<pkg/igsfit.h>
+
+int procedure igs_nearestu (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+
+pointer	gp			# GIO pointer
+int	ztype			# Zoom type
+int	refpt			# Reference point
+int	axis[2]			# Axes
+real	pts[npts, ARB]		# Data points
+int	npts			# Number of data points
+real	wx, wy			# Cursor coordinates
+int	wcs			# WCS
+
+int	i, j, x, y
+real	r2, r2min, x0, y0
+
+begin
+	x = axis[1]
+	y = axis[2]
+
+	call gctran (gp, wx, wy, wx, wy, wcs, 0)
+	r2min = MAX_REAL
+	j = 0
+
+	if (IS_INDEFI (ztype)) {
+	    do i = 1, npts {
+	        if (pts[i,W] != 0.)
+		    next
+	        call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+	} else {
+	    do i = 1, npts {
+	        if ((pts[i,ztype] != pts[refpt,ztype]) || (pts[i,W] != 0.))
+		    next
+	        call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+	}
+
+	return (j)
+end
+
+
+# IGS_UNDELETE - Undelete point or subset.
+
+procedure igs_undelete (gp, gt, ztype, refpt, axis, pts, wts, npts, dtype)
+
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+int	ztype			# Zoom type
+int	refpt			# Reference point for undeletion
+int	axis[2]			# Axes
+real	pts[npts, ARB]		# Data points
+real	wts[npts]		# Original weights
+int	npts			# Number of data points
+int	dtype			# Undeletion type
+
+int	i, x, y
+real	xsize, ysize
+
+real	gt_getr()
+
+begin
+	x = axis[1]
+	y = axis[2]
+
+	xsize = gt_getr (gt, GTXSIZE)
+	ysize = gt_getr (gt, GTYSIZE)
+
+	switch (dtype) {
+	case X, Y, Z:
+	    do i = 1, npts {
+		if (!IS_INDEFI (ztype))
+		    if (pts[refpt,ztype] != pts[i,ztype])
+		        next
+		if (pts[refpt,dtype] != pts[i,dtype])
+		    next
+	        call gseti (gp, G_PMLTYPE, 0)
+	        call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+	        call gseti (gp, G_PMLTYPE, 1)
+	        call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+	        if (wts[i] == 0)
+		    wts[i] = 1
+	        pts[i,W] = wts[i]
+	    }
+	default:
+	    call gseti (gp, G_PMLTYPE, 0)
+	    call gmark (gp, pts[refpt,x], pts[refpt,y], GM_CROSS, xsize, ysize)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, pts[refpt,x], pts[refpt,y], GM_PLUS, xsize, ysize)
+	    if (wts[refpt] == 0)
+		wts[refpt] = 1
+	    pts[refpt,W] = wts[refpt]
+	}
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/mkpkg b/noao/twodspec/longslit/transform/igsfit/mkpkg
new file mode 100644
index 00000000..ac5a6ca9
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/mkpkg
@@ -0,0 +1,21 @@
+# Interactive General Surface Fitting Package
+
+$checkout libpkg.a ../../
+$update	  libpkg.a
+$checkin  libpkg.a ../../
+$exit
+
+libpkg.a:
+	igscolon.x	igsfit.com <gset.h>
+	igsdelete.x	<gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+	igsfit.x	igsfit.com <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+	igsget.x	igsfit.com <pkg/igsfit.h>
+	igsgraph.x	<gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+	igsinit.x	<pkg/igsfit.h>
+	igsnearest.x	<gset.h> <mach.h> <pkg/igsfit.h>
+	igsparams.x	igsfit.com <pkg/gtools.h>
+	igsset.x	igsfit.com <pkg/igsfit.h>
+	igssolve.x	igsfit.com <math/gsurfit.h>
+	igsundelete.x	<gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+	xgs.x		<math/gsurfit.h>
+	;
diff --git a/noao/twodspec/longslit/transform/igsfit/xgs.x b/noao/twodspec/longslit/transform/igsfit/xgs.x
new file mode 100644
index 00000000..7d2ea331
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/xgs.x
@@ -0,0 +1,243 @@
+include <math/gsurfit.h>
+
+# XGS -- These routines provide an interface between real input data and
+# the double precision surface fitting.  Rather than make the input data
+# be double precision we only want the internal surface fitting arithmetic
+# to be double.  But the surface fitting package only provides real
+# arithmetic for real input and double precision arithmetic for double
+# precision input.  Hence these interfaces.  Note that the save and restore
+# functions use double precision.
+
+# XGSINIT --  Procedure to initialize the surface descriptor.
+
+procedure xgsinit (sf, surface_type, xorder, yorder, xterms, xmin, xmax,
+    ymin, ymax)
+
+pointer	sf		# surface descriptor
+int	surface_type	# type of surface to be fitted
+int	xorder		# x order of surface to be fit
+int	yorder		# y order of surface to be fit
+int	xterms		# presence of cross terms
+real    xmin		# minimum value of x
+real	xmax		# maximum value of x
+real	ymin		# minimum value of y
+real	ymax		# maximum value of y
+
+begin
+	call dgsinit (sf, surface_type, xorder, yorder, xterms, double (xmin),
+	    double (xmax), double (ymin), double (ymax))
+end
+
+
+# XGSFIT -- Procedure to solve the normal equations for a surface.
+
+procedure xgsfit (sf, x, y, z, w, npts, wtflag, ier)
+
+pointer	sf		# surface descriptor
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	z[npts]		# data array
+real	w[npts]		# array of weights
+int	npts		# number of data points
+int	wtflag		# type of weighting
+int	ier		# ier = OK, everything OK
+			# ier = SINGULAR, matrix is singular, 1 or more
+			# coefficients are 0.
+			# ier = NO_DEG_FREEDOM, too few points to solve matrix
+
+pointer	sp, xd, yd, zd, wd
+errchk	salloc
+
+begin
+	call smark (sp)
+	call salloc (xd, npts, TY_DOUBLE)
+	call salloc (yd, npts, TY_DOUBLE)
+	call salloc (zd, npts, TY_DOUBLE)
+	call salloc (wd, npts, TY_DOUBLE)
+	call achtrd (x, Memd[xd], npts)
+	call achtrd (y, Memd[yd], npts)
+	call achtrd (z, Memd[zd], npts)
+	call achtrd (w, Memd[wd], npts)
+	call dgsfit (sf, Memd[xd], Memd[yd], Memd[zd], Memd[wd], npts,
+	    wtflag, ier)
+	call sfree (sp)
+end
+
+
+# XGSVECTOR -- Procedure to evaluate the fitted surface at an array of points.
+
+procedure xgsvector (sf, x, y, zfit, npts)
+
+pointer	sf		# pointer to surface descriptor structure
+real	x[ARB]		# x value
+real	y[ARB]		# y value
+real	zfit[ARB]	# fits surface values
+int	npts		# number of data points
+
+pointer	sp, xd, yd, zd
+errchk	salloc
+
+begin
+	call smark (sp)
+	call salloc (xd, npts, TY_DOUBLE)
+	call salloc (yd, npts, TY_DOUBLE)
+	call salloc (zd, npts, TY_DOUBLE)
+	call achtrd (x, Memd[xd], npts)
+	call achtrd (y, Memd[yd], npts)
+	call dgsvector (sf, Memd[xd], Memd[yd], Memd[zd], npts)
+	call achtdr (Memd[zd], zfit, npts)
+	call sfree (sp)
+end
+
+
+# XGSEVAL -- Procedure to evaluate the fitted surface at a single point.
+
+real procedure xgseval (sf, x, y)
+
+pointer	sf		# pointer to surface descriptor structure
+real	x		# x value
+real	y		# y value
+
+double	dgseval()
+
+begin
+	return (real (dgseval (sf, double (x), double (y))))
+end
+
+
+# XGSADD -- Procedure to add the fits from two surfaces together.
+
+procedure xgsadd (sf1, sf2, sf3)
+
+pointer	sf1		# pointer to the first surface
+pointer	sf2		# pointer to the second surface
+pointer	sf3		# pointer to the output surface
+
+begin
+	call dgsadd (sf1, sf2, sf3)
+end
+
+
+# XGSFREE -- Procedure to free the surface descriptor
+
+procedure xgsfree (sf)
+
+pointer	sf	# the surface descriptor
+
+begin
+	call dgsfree (sf)
+end
+
+
+# XGSGCOEFF -- Procedure to fetch a particular coefficient.
+
+real procedure xgsgcoeff (sf, xorder, yorder)
+
+pointer	sf		# pointer to the surface fitting descriptor
+int	xorder		# X order of desired coefficent
+int	yorder		# Y order of desired coefficent
+
+double	dgsgcoeff()
+
+begin
+	return (real (dgsgcoeff (sf, xorder, yorder)))
+end
+
+
+# XGSSCOEFF -- Procedure to set a particular coefficient.
+
+procedure xgsscoeff (sf, xorder, yorder, coeff)
+
+pointer	sf		# pointer to the surface fitting descriptor
+int	xorder		# X order of desired coefficent
+int	yorder		# Y order of desired coefficent
+real	coeff		# Coefficient value
+
+begin
+	call dgsscoeff (sf, xorder, yorder, double (coeff))
+end
+
+
+# XGSGETR -- Procedure to fetch a real gsurfit parameter
+
+real procedure xgsgetr (sf, parameter)
+
+pointer	sf		# pointer to the surface fit
+int	parameter	# parameter to be fetched
+
+double	dgsgetd()
+
+begin
+	return (real (dgsgetd (sf, parameter)))
+end
+
+
+# XGSGETI -- Procedure to fetch an integer parameter
+
+int procedure xgsgeti (sf, parameter)
+
+pointer sf		# pointer to the surface fit
+int	parameter	# integer parameter
+
+int	dgsgeti()
+
+begin
+	return (dgsgeti (sf, parameter))
+end
+
+
+# XGSSAVE -- Procedure to save the surface fit for later use by the
+# evaluate routines.
+#
+# NOTE THAT THIS USES DOUBLE PRECISION FOR THE COEFFICIENTS.
+
+procedure xgssave (sf, fit)
+
+pointer	sf		# pointer to the surface descriptor
+double	fit[ARB]	# array for storing fit
+
+begin
+	call dgssave (sf, fit)
+end
+
+
+# XGSRESTORE -- Procedure to restore the surface fit stored by GSSAVE
+# to the surface descriptor for use by the evaluating routines.
+#
+# NOTE THAT THIS USES DOUBLE PRECISION FOR THE COEFFICIENTS.
+
+procedure xgsrestore (sf, fit)
+
+pointer	sf		# surface descriptor
+double	fit[ARB]	# array containing the surface parameters and
+
+begin
+	call dgsrestore (sf, fit)
+end
+
+
+# XGSDER -- Procedure to calculate a new surface which is a derivative of
+# the previous surface
+
+procedure xgsder (sf1, x, y, zfit, npts, nxd, nyd)
+
+pointer	sf1		# pointer to the previous surface
+real	x[npts]		# x values
+real	y[npts]		# y values
+real	zfit[npts]	# fitted values
+int	npts		# number of points
+int	nxd, nyd	# order of the derivatives in x and y
+
+pointer	sp, xd, yd, zd
+
+begin
+	call smark (sp)
+	call salloc (xd, npts, TY_DOUBLE)
+	call salloc (yd, npts, TY_DOUBLE)
+	call salloc (zd, npts, TY_DOUBLE)
+	call achtrd (x, Memd[xd], npts)
+	call achtrd (y, Memd[yd], npts)
+	call dgsder (sf1, Memd[xd], Memd[yd], Memd[zd], npts, nxd, nyd)
+	call achtdr (Memd[zd], zfit, npts)
+	call sfree (sp)
+end
diff --git a/noao/twodspec/longslit/transform/mkpkg b/noao/twodspec/longslit/transform/mkpkg
new file mode 100644
index 00000000..8ea1b584
--- /dev/null
+++ b/noao/twodspec/longslit/transform/mkpkg
@@ -0,0 +1,20 @@
+# Coordinate Transformation Tasks
+
+$checkout libpkg.a ../
+$update	  libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	@igsfit
+
+	fcdbio.x	<error.h> <math/gsurfit.h> <pkg/dttext.h> <units.h>
+	fcdlist.x	<error.h> <mach.h>
+	fcfitcoords.x	<pkg/gtools.h> <pkg/igsfit.h> <pkg/xtanswer.h>
+	fcgetcoords.x	<imio.h> <mach.h> <pkg/dttext.h> <pkg/igsfit.h>
+	fcgetim.x	
+	fitcoords.x	<error.h> <pkg/igsfit.h> <pkg/xtanswer.h>
+	trsetup.x	<math.h> <math/gsurfit.h> <math/iminterp.h>
+	t_fceval.x
+	t_transform.x	transform.com <imhdr.h> <math/iminterp.h> <units.h>
+	;
diff --git a/noao/twodspec/longslit/transform/t_fceval.x b/noao/twodspec/longslit/transform/t_fceval.x
new file mode 100644
index 00000000..a9c5cc75
--- /dev/null
+++ b/noao/twodspec/longslit/transform/t_fceval.x
@@ -0,0 +1,107 @@
+# T_FCEVAL -- Evaluate FITCOORDS solutions.
+# Input consists of a text file of pixel coordinates to be evaluated and the
+# user coordinate surfaces from FITCOORDS.  The output is a text file of the
+# input coordinates followed by the output coordinates.  When there is no fit
+# for an axis the unit transformation is used and when there is more than one
+# fit for an axis the average is used.
+
+procedure t_fceval ()
+
+pointer	input			# File of input coordinates
+pointer	output			# File of output coordinates
+int	fitnames		# List of user coordinate fits
+pointer	database		# Database
+
+int	i, j, in, out, nsf[2]
+double	x[2], y[2]
+pointer	sp, fitname, sf[2], un[2], sf1, un1
+
+bool	un_compare()
+int	open(), fscan(), nscan()
+int	clpopnu(), clplen(), clgfil()
+double	dgseval()
+errchk	open, lm_dbread
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (database, SZ_FNAME, TY_CHAR)
+	call salloc (fitname, SZ_FNAME, TY_CHAR)
+
+	# Get parameters.
+	call clgstr ("input", Memc[input], SZ_FNAME)
+	call clgstr ("output", Memc[output], SZ_FNAME)
+	fitnames = clpopnu ("fitnames")
+	call clgstr ("database", Memc[database], SZ_FNAME)
+
+	# Open the input and output files.
+	in = open (Memc[input], READ_ONLY, TEXT_FILE)
+	out = open (Memc[output], NEW_FILE, TEXT_FILE)
+
+	# Read the solutions.
+	i = max (1, clplen (fitnames))
+	call salloc (sf[1], i, TY_INT)
+	call salloc (sf[2], i, TY_INT)
+
+	nsf[1] = 0; nsf[2] = 0; un[1] = NULL; un[2] = NULL
+	while (clgfil (fitnames, Memc[fitname], SZ_FNAME) != EOF) {
+	    call lm_dbread (Memc[database], Memc[fitname], j, un1, sf1)
+	    if (un1 != NULL) {
+		if (un[j] == NULL)
+		    un[j] = un1
+		else if (un_compare (un1, un[j]))
+		    call un_close (un1)
+		else
+		    call error (1, "Input units disagree")
+	    }
+
+	    if (sf1 != NULL) {
+		Memi[sf[j]+nsf[j]] = sf1
+		nsf[j] = nsf[j] + 1
+	    }
+	}
+
+	if (nsf[1] + nsf[2] == 0)
+	    call error (0, "No user coordinates")
+
+	# Evaluate the fits at each input coordinate.
+	while (fscan (in) != EOF) {
+	    call gargd (x[1])
+	    call gargd (x[2])
+	    if (nscan() != 2)
+		next
+
+	    do j = 1, 2 {
+		if (nsf[j] == 0)
+		    y[j] = x[j]
+		else {
+		    y[j] = dgseval (Memi[sf[j]], x[1], x[2])
+		    do i = 2, nsf[1]
+			y[j]  = y[j]  + dgseval (Memi[sf[j]+i-1], x[1], y[2])
+		    y[j]  = y[j]  / nsf[j]
+		}
+	    }
+
+	    call fprintf (out, "%g %g %g %g\n")
+		call pargd (x[1])
+		call pargd (x[2])
+		call pargd (y[1])
+		call pargd (y[2])
+	    call flush (out)
+	}
+
+	# Free the surfaces and units structures.
+	do j = 1, 2 {
+	    for (i=1; i<=nsf[j]; i=i+1)
+		call dgsfree (Memi[sf[j]+i-1])
+	    if (un[j] != NULL)
+		call un_close (un[j])
+	}
+
+	# Finish up.
+	call clpcls (fitnames)
+	call close (out)
+	call close (in)
+	call sfree (sp)
+end
diff --git a/noao/twodspec/longslit/transform/t_transform.x b/noao/twodspec/longslit/transform/t_transform.x
new file mode 100644
index 00000000..5610858e
--- /dev/null
+++ b/noao/twodspec/longslit/transform/t_transform.x
@@ -0,0 +1,741 @@
+include	<imhdr.h>
+include	<math/iminterp.h>
+include	<units.h>
+
+define	ITYPES	"|nearest|linear|poly3|poly5|spline3|"
+
+# T_TRANSFORM -- Transform longslit images.
+# Input consists of images to be transformed, the user coordinate surfaces
+# describing the output coordinates in terms of the input coordinates,
+# and the desired coordinates for the output images.  The type of image
+# interpolation is also input.  There is a log output as well as the
+# transformed images.  The output image may replace the input image.
+
+procedure t_transform ()
+
+int	input			# List of input images
+int	output			# List of output images
+int	minput			# List of input masks
+int	moutput			# List of output masks
+int	fitnames		# List of user coordinate fits
+pointer	database		# Database
+char	interp[10]		# Interpolation type
+int	logfiles		# List of log files
+
+int	itypes[II_NTYPES2D], logfd, nusf, nvsf
+pointer	in, out, pmin, pmout
+pointer	un[2], mw, ct, usf, vsf, xmsi, ymsi, jmsi, xout, yout, dxout, dyout
+pointer	sp, image1, image2, image3, minname, moutname, mname, str
+
+int	clpopnu(), clgfil(), clplen(), clgeti(), clgwrd(), open()
+int	imtopenp(), imtlen(), imtgetim()
+bool	clgetb()
+real	clgetr()
+pointer	immap(), mw_openim(), yt_mappm()
+errchk	tr_gsf, tr_setup, open, mw_openim, yt_mappm
+
+data	itypes /II_BINEAREST, II_BILINEAR, II_BIPOLY3, II_BIPOLY5,
+	II_BISPLINE3, II_SINC, II_LSINC, II_DRIZZLE/
+
+include	"transform.com"
+
+
+begin
+	call smark (sp)
+	call salloc (database, SZ_FNAME, TY_CHAR)
+	call salloc (image1, SZ_FNAME, TY_CHAR)
+	call salloc (image2, SZ_FNAME, TY_CHAR)
+	call salloc (image3, SZ_FNAME, TY_CHAR)
+	call salloc (minname, SZ_FNAME, TY_CHAR)
+	call salloc (moutname, SZ_FNAME, TY_CHAR)
+	call salloc (mname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get and error check the input and output image lists and the other
+	# task parameters.
+
+	input = imtopenp ("input")
+	output = imtopenp ("output")
+	if (imtlen (input) != imtlen (output)) {
+	    call imtclose (input)
+	    call imtclose (output)
+	    call error (1, "Number of input and output images differ")
+	}
+	minput = imtopenp ("minput")
+	moutput = imtopenp ("moutput")
+	if (imtlen (minput) > 1 && imtlen (minput) != imtlen (input)) {
+	    call imtclose (input)
+	    call imtclose (output)
+	    call imtclose (minput)
+	    call imtclose (moutput)
+	    call error (1, "Can't associate input masks with input images")
+	}
+	if (imtlen (moutput) > 0 && imtlen (input) != imtlen (moutput)) {
+	    call imtclose (input)
+	    call imtclose (output)
+	    call imtclose (minput)
+	    call imtclose (moutput)
+	    call error (1, "Number output masks differ from input")
+	}
+
+	fitnames = clpopnu ("fitnames")
+	call clgstr ("database", Memc[database], SZ_FNAME)
+	itype = itypes[clgwrd ("interptype", interp, 10, II_FUNCTIONS)]
+	logfiles = clpopnu ("logfiles")
+
+	u1 = clgetr ("x1")
+	u2 = clgetr ("x2")
+	du = clgetr ("dx")
+	nu = clgeti ("nx")
+	v1 = clgetr ("y1")
+	v2 = clgetr ("y2")
+	dv = clgetr ("dy")
+	nv = clgeti ("ny")
+
+	ulog = clgetb ("xlog")
+	vlog = clgetb ("ylog")
+	flux = clgetb ("flux")
+	blank = clgetr ("blank")
+
+	usewcs = (clplen (fitnames) == 0)
+
+	# Transform each input image to the output image.
+	Memc[minname] = EOS
+	Memc[moutname] = EOS
+	Memc[mname] = EOS
+	xmsi = NULL
+	while ((imtgetim (input, Memc[image1], SZ_FNAME) != EOF)  &&
+	    (imtgetim (output, Memc[image2], SZ_FNAME) != EOF)) {
+
+	    # Get mask names.
+	    if (imtgetim (minput, Memc[image3], SZ_FNAME) != EOF)
+	        call strcpy (Memc[image3], Memc[minname], SZ_FNAME)
+	    if (imtgetim (moutput, Memc[image3], SZ_FNAME) != EOF)
+	        call strcpy (Memc[image3], Memc[moutname], SZ_FNAME)
+
+	    # Map the input and output images.
+	    call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3],SZ_FNAME)
+	    in = immap (Memc[image1], READ_ONLY, 0)
+	    out = immap (Memc[image2], NEW_COPY, in)
+
+	    # Map masks.
+	    pmin = NULL; pmout = NULL
+	    if (Memc[minname] != EOS)
+	        pmin = yt_mappm (Memc[minname], in, "logical", Memc[mname],
+		    SZ_FNAME)
+	    if (Memc[moutname] != EOS) {
+	        call xt_maskname (Memc[moutname], "", NEW_IMAGE,
+		    Memc[moutname], SZ_FNAME)
+		pmout = immap (Memc[moutname], NEW_COPY, in)
+		call imastr (out, "BPM", Memc[moutname])
+	    }
+
+	    # Get the coordinate transformation surfaces from the database
+	    # and setup the transformations.
+	    # Do this only on the first pass.
+
+	    if (xmsi == NULL) {
+		if (usewcs) {
+		    mw = mw_openim (in)
+		    call tr_gwcs (mw, un, IM_LEN(in,1), IM_LEN(in,2), ct,
+			usf, nusf, vsf, nvsf)
+		} else {
+		    mw = NULL
+		    ct = NULL
+		    call tr_gsf (Memc[database], fitnames, un, usf, nusf,
+		        vsf, nvsf)
+		}
+		call tr_setup (ct, usf, nusf, vsf, nvsf, un, xmsi, ymsi, jmsi,
+		    xout, yout, dxout, dyout)
+		if (mw != NULL)
+		    call mw_close (mw)
+	    }
+
+	    # Write log information.
+	    while (clgfil (logfiles, Memc[str], SZ_LINE) != EOF) {
+		logfd = open (Memc[str], APPEND, TEXT_FILE)
+		call sysid (Memc[str], SZ_LINE)
+		call fprintf (logfd, "\n%s\n")
+		    call pargstr (Memc[str])
+		call fprintf (logfd, "  Transform %s to %s.\n")
+		    call pargstr (Memc[image1])
+		    call pargstr (Memc[image3])
+		if (pmout != EOS) {
+		    if (pmin != EOS) {
+		        call fprintf (logfd, "  Transform mask %s to %s.\n")
+			    call pargstr (Memc[mname])
+			    call pargstr (Memc[moutname])
+		    } else {
+		        call fprintf (logfd, "  Output mask is %s.\n")
+			    call pargstr (Memc[moutname])
+		    }
+		}
+		if (flux)
+		    call fprintf (logfd, "  Conserve flux per pixel.\n")
+		if (usewcs)
+		    call fprintf (logfd, "  Transforming using image WCS.\n")
+		else {
+		    call fprintf (logfd, "  User coordinate transformations:\n")
+		    while (clgfil (fitnames, Memc[str], SZ_LINE) != EOF) {
+			call fprintf (logfd, "    %s\n")
+			    call pargstr (Memc[str])
+		    }
+		}
+		call fprintf (logfd, "  Interpolation is %s.\n")
+		    call pargstr (interp)
+		if (!IS_INDEFR(blank)) {
+		    call fprintf (logfd, "  Out of bounds pixel value is %g.\n")
+			call pargr (blank)
+		} else
+		    call fprintf (logfd,
+		    "  Using edge extension for out of bounds pixel values.\n")
+		call fprintf (logfd, "  Output coordinate parameters are:\n")
+		call fprintf (logfd,
+	    "    x1 = %10.4g, x2 = %10.4g, dx = %10.4g, nx = %4d, xlog = %b\n")
+		    call pargr (u1)
+		    call pargr (u2)
+		    call pargr (du)
+		    call pargi (nu)
+		    call pargb (ulog)
+		call fprintf (logfd,
+	    "    y1 = %10.4g, y2 = %10.4g, dy = %10.4g, ny = %4d, ylog = %b\n")
+		    call pargr (v1)
+		    call pargr (v2)
+		    call pargr (dv)
+		    call pargi (nv)
+		    call pargb (vlog)
+		call close (logfd)
+	    }
+	    call clprew (logfiles)
+
+	    call tr_transform (in, out, pmin, pmout, un, xmsi, ymsi, jmsi,
+	        Memr[xout], Memr[yout], Memr[dxout], Memr[dyout])
+
+	    if (pmout != NULL)
+	        call imunmap (pmout)
+	    if (pmin != NULL)
+	        call xt_pmunmap (pmin)
+	    call imunmap (in)
+	    call imunmap (out)
+	    call xt_delimtemp (Memc[image2], Memc[image3])
+
+	    if (usewcs) {
+		call mfree (xout, TY_REAL)
+		call mfree (yout, TY_REAL)
+		call mfree (dxout, TY_REAL)
+		call mfree (dyout, TY_REAL)
+		if (xmsi != NULL)
+		    call msifree (xmsi)
+		if (ymsi != NULL)
+		    call msifree (ymsi)
+		if (jmsi != NULL)
+		    call msifree (jmsi)
+		if (un[1] != NULL)
+		    call un_close (un[1])
+		if (un[2] != NULL)
+		    call un_close (un[2])
+		xmsi = NULL
+	    }
+
+	}
+
+	call mfree (xout, TY_REAL)
+	call mfree (yout, TY_REAL)
+	call mfree (dxout, TY_REAL)
+	call mfree (dyout, TY_REAL)
+	if (xmsi != NULL)
+	    call msifree (xmsi)
+	if (ymsi != NULL)
+	    call msifree (ymsi)
+	if (jmsi != NULL)
+	    call msifree (jmsi)
+	if (un[1] != NULL)
+	    call un_close (un[1])
+	if (un[2] != NULL)
+	    call un_close (un[2])
+	call imtclose (minput)
+	call imtclose (moutput)
+	call imtclose (input)
+	call imtclose (output)
+	call clpcls (fitnames)
+	call clpcls (logfiles)
+	call sfree (sp)
+end
+
+
+# TR_SETOUTPUT -- Set the output coordinates in the common block.
+# This procedure allows the user to specifying a part of the output
+# coordinates and let the rest default based on the full limits of
+# the user coordinate surfaces.
+
+procedure tr_setoutput (xmin, xmax, ymin, ymax, umin, umax, vmin, vmax)
+
+real	xmin, xmax, ymin, ymax
+real	umin, umax, vmin, vmax
+
+int	nua, nva
+real	u1a, u2a, dua, v1a, v2a, dva
+
+include	"transform.com"
+
+begin
+	# Save the original values of the user parameters.
+	u1a = u1
+	u2a = u2
+	dua = du
+	nua = nu
+	v1a = v1
+	v2a = v2
+	dva = dv
+	nva = nv
+
+	# If the output coordinate limits are not defined then use the
+	# transformation surface limits.
+
+	if (IS_INDEF (u1))
+	    u1 = umin
+	if (IS_INDEF (u2))
+	    u2 = umax
+	if (IS_INDEF (v1))
+	    v1 = vmin
+	if (IS_INDEF (v2))
+	    v2 = vmax
+
+	# If the number of output pixels are not defined then use the number
+	# of pixels in the input image.
+
+	if (IS_INDEFI (nu))
+	    nu = xmax - xmin + 1
+	if (IS_INDEFI (nv))
+	    nv = ymax - ymin + 1
+
+	# If the coordinate interval is not defined determine it from the
+	# number of pixels and the coordinate limits.  If the interval is
+	# defined then override the number of pixels.
+
+	if (ulog) {
+	    if (IS_INDEF (du))
+		du = (log10 (u2) - log10 (u1)) / (nu - 1)
+	    else if (IS_INDEFI (nua))
+	        nu = nint ((log10 (u2) - log10 (u1)) / du + 1)
+	    else if (IS_INDEF (u1a))
+		u1 = 10.0 ** (log10 (u2) - du * (nu - 1))
+	    else
+		u2 = 10.0 ** (log10 (u1) + du * (nu - 1))
+	} else {
+	    if (IS_INDEF (du))
+	        du = (u2 - u1) / (nu - 1)
+	    else if (IS_INDEFI (nua))
+	        nu = nint ((u2 - u1) / du + 1)
+	    else if (IS_INDEF (u1a))
+		u1 = u2 - du * (nu - 1)
+	    else
+		u2 = u1 + du * (nu - 1)
+	}
+
+	if (vlog) {
+	    if (IS_INDEF (dv))
+		dv = (log10 (v2) - log10 (v1)) / (nv - 1)
+	    else if (IS_INDEFI (nva))
+	        nv = nint ((log10 (v2) - log10 (v1)) / dv + 1)
+	    else if (IS_INDEF (v1a))
+		v1 = 10.0 ** (log10 (v2) - dv * (nv - 1))
+	    else
+		v2 = 10.0 ** (log10 (v1) + dv * (nv - 1))
+	} else {
+	    if (IS_INDEF (dv))
+	        dv = (v2 - v1) / (nv - 1)
+	    else if (IS_INDEFI (nva))
+	        nv = nint ((v2 - v1) / dv + 1)
+	    else if (IS_INDEF (v1a))
+		v1 = v2 - dv * (nv - 1)
+	    else
+		v2 = v1 + dv * (nv - 1)
+	}
+end
+
+
+define	NBUF		16	# Additional buffer for interpolation
+define	NEDGE		2	# Number of edge lines to add for interpolation
+define	MINTERP		100	# Mask value for input mask interpolation
+define	MTHRESH		10	# Interpolated mask value for bad pixels
+define	MBAD		1	# Mask value for output bad pixels
+define	MBLANK		1	# Mask value for out of bounds pixels
+
+# TR_TRANSFORM -- Perform the image transformation using a user specified
+# image interpolator.  If an input and output mask are included the input
+# mask values are set to MINTERP, interpolated in the same way, and any values
+# greater than MTHRESH are set to MBAD.  Note that currently the input mask
+# values are not used in computing the input data interpolation value.
+# The masks MUST be the same size as the input data and are assumed to
+# be registered in logical pixel coordinates.
+
+procedure tr_transform (in, out, pmin, pmout, un, xmsi, ymsi, jmsi, xout, yout,
+	dxout, dyout)
+
+pointer	in, out			#I IMIO data pointers
+pointer	pmin, pmout		#I IMIO mask pointers (NULL if not used)
+pointer	un[2]			#I Units
+pointer xmsi, ymsi		#I Coordinate interpolation pointers
+pointer	jmsi			#I Jacobian interpolation pointer
+real	xout[ARB], yout[ARB]	#I Output grid relative to interpolation surface
+real	dxout[ARB], dyout[ARB]	#I Output coordinate intervals
+
+int	i, j, nxin, nyin, line1, line2, line3, line4, nlines, laxis, paxis
+bool	xofb, yofb
+real	a, b, c, r[2], w[2], cd[2,2]
+pointer	zmsi, mzmsi, buf, mbuf, bufout
+pointer	sp, xin, yin, jbuf, xin1, yin1, y, mw
+
+pointer	mw_open(), impl2r()
+errchk	get_daxis
+
+include	"transform.com"
+
+begin
+	# Initialize the output image header.
+
+	IM_LEN(out, 1) = nu
+	IM_LEN(out, 2) = nv
+	if (pmout != NULL) {
+	    IM_LEN(pmout, 1) = nu
+	    IM_LEN(pmout, 2) = nv
+	}
+
+	mw = mw_open (NULL, 2)
+	call mw_newsystem (mw, "world", 2)
+	do i = 1, 2 {
+	    call mw_swtype (mw, i, 1, "linear", "")
+	    if (un[i] != NULL) {
+		call mw_swattrs (mw, i, "label", UN_LABEL(un[i]))
+		call mw_swattrs (mw, i, "units", UN_UNITS(un[i]))
+	    }
+	}
+
+	r[1] = 1.
+	if (ulog)
+	    w[1] = log10 (u1)
+	else
+	    w[1] = u1
+	cd[1,1] = du
+	cd[1,2] = 0.
+	r[2] = 1.
+	if (vlog)
+	    w[2] = log10 (v1)
+	else
+	    w[2] = v1
+	cd[2,2] = dv
+	cd[2,1] = 0.
+	call mw_swtermr (mw, r, w, cd, 2)
+
+	# The following image parameters are for compatibility with the
+	# ONEDSPEC package if using database solutions.
+
+	if (!usewcs) {
+	    call imastr (out, "DCLOG1", "Transform")
+	    iferr (call imdelf (out, "REFSPEC1"))
+		;
+	    iferr (call imdelf (out, "REFSPEC2"))
+		;
+	    call get_daxis (in, laxis, paxis)
+	    call imaddi (out, "dispaxis", laxis)
+	    switch (laxis) {
+	    case 1:
+		if (ulog)
+		    call imaddi (out, "dc-flag", 1)
+		else
+		    call imaddi (out, "dc-flag", 0)
+		if (un[laxis] == NULL) {
+		    call mw_swattrs (mw, laxis, "label", "Wavelength")
+		    call mw_swattrs (mw, laxis, "units", "Angstroms")
+		}
+	    case 2:
+		if (vlog)
+		    call imaddi (out, "dc-flag", 1)
+		else
+		    call imaddi (out, "dc-flag", 0)
+		if (un[laxis] == NULL) {
+		    call mw_swattrs (mw, laxis, "label", "Wavelength")
+		    call mw_swattrs (mw, laxis, "units", "Angstroms")
+		}
+	    }
+	}
+	call mw_saveim (mw, out)
+	if (pmout != NULL)
+	    call mw_saveim (mw, pmout)
+	call mw_close (mw)
+
+	# Allocate memory for the input coordinates and a vector for the
+	# output y coordinates.  Also initialize the image data buffer.
+
+	call smark (sp)
+	call salloc (xin, nu, TY_REAL)
+	call salloc (yin, nu, TY_REAL)
+	call salloc (y, nu, TY_REAL)
+	if (flux)
+	    call salloc (jbuf, nu, TY_REAL)
+	if (!IS_INDEFR(blank) || pmout != NULL) {
+	    call salloc (xin1, nu, TY_REAL)
+	    call salloc (yin1, nu, TY_REAL)
+	}
+
+	buf = NULL
+	mbuf = NULL
+	nlines = 0
+
+	# Initialize the interpolator.
+
+	call msiinit (zmsi, itype)
+	if (pmin != NULL)
+	    call msiinit (mzmsi, itype)
+
+	# Do each line of the output image.
+
+	nxin = IM_LEN(in, 1)
+	nyin = IM_LEN(in, 2)
+
+	do i = 1, nv {
+
+	    # Evaluate the input coordinates at the output grid for a line
+	    # of the output image using the interpolation surfaces.
+
+	    call amovkr (yout[i], Memr[y], nu)
+	    if (!IS_INDEFR(blank) || pmout != NULL) {
+		call msivector (xmsi, xout, Memr[y], Memr[xin1], nu)
+		call msivector (ymsi, xout, Memr[y], Memr[yin1], nu)
+		call amovr (Memr[xin1], Memr[xin], nu)
+		call amovr (Memr[yin1], Memr[yin], nu)
+	    } else {
+		call msivector (xmsi, xout, Memr[y], Memr[xin], nu)
+		call msivector (ymsi, xout, Memr[y], Memr[yin], nu)
+	    }
+
+	    # Determine the coordinate ranges and check for out of bounds.
+
+	    call alimr (Memr[xin], nu, a, b)
+	    xofb = (a < 1 || b > nxin)
+	    if (xofb) {
+		if (a < 1)
+		    call arltr (Memr[xin], nu, 1., 1.)
+		if (b > nxin)
+		    call argtr (Memr[xin], nu, real (nxin), real (nxin))
+	    }
+
+	    call alimr (Memr[yin], nu, a, b)
+	    yofb = (a < 1 || b > nyin)
+	    if (yofb) {
+		if (a < 1) {
+		    call arltr (Memr[yin], nu, 1., 1.)
+		    a = 1.
+		    b = max (a, b)
+		}
+		if (b > nyin) {
+		    call argtr (Memr[yin], nu, real (nyin), real (nyin))
+		    b = nyin
+		    a = min (a, b)
+		}
+	    }
+
+	    # Get the input image data and fit an interpolator to the data.
+
+	    if ((buf == NULL) || (b > line2) || (a < line1)) {
+		nlines = max (nlines, int (b - a + 2 + NBUF))
+		if (buf == NULL) {
+		    if (a < nyin / 2) {
+		        line1 = max (1, int (a))
+		        line2 = min (nyin, line1 + nlines - 1)
+		    } else {
+		        line2 = min (nyin, int (b+1.))
+		        line1 = max (1, line2 - nlines + 1)
+		    }
+		} else if (b > line2) {
+		    line1 = max (1, int (a))
+		    line2 = min (nyin, line1 + nlines - 1)
+		    line1 = max (1, line2 - nlines + 1)
+		} else {
+		    line2 = min (nyin, int (b+1.))
+		    line1 = max (1, line2 - nlines + 1)
+		    line2 = min (nyin, line1 + nlines - 1)
+		}
+		line3 = max (1, line1 - NEDGE)
+		line4 = min (nyin, line2 + NEDGE)
+	        call tr_bufl2r (in, pmin, line3, line4, buf, mbuf)
+	        call msifit (zmsi, Memr[buf], nxin, line4 - line3 + 1, nxin)
+		if (pmin != NULL)
+		    call msifit (mzmsi, Memr[mbuf], nxin, line4 - line3 + 1,
+		        nxin)
+	    }
+
+	    # The input coordinates must be offset to interpolation data grid.
+	    call asubkr (Memr[yin], real (line3 - 1), Memr[yin], nu)
+
+	    # Evaluate output image pixels, conserve flux (if requested) using
+	    # the Jacobian, and set the out of bounds values.
+
+	    bufout = impl2r (out, i)
+	    call msivector (zmsi, Memr[xin], Memr[yin], Memr[bufout], nu)
+	    if (flux) {
+	        call msivector (jmsi, xout, Memr[y], Memr[jbuf], nu)
+		call amulr (dxout, Memr[jbuf], Memr[jbuf], nu)
+		call amulkr (Memr[jbuf], dyout[i], Memr[jbuf], nu)
+	        call amulr (Memr[bufout], Memr[jbuf], Memr[bufout], nu)
+	    }
+	    if (!IS_INDEFR(blank)) {
+		if (xofb) {
+		    do j = 0, nu-1 {
+			if (Memr[xin1+j] < 1 || Memr[xin1+j] > nxin)
+			    Memr[bufout+j] = blank
+		    }
+		}
+		if (yofb) {
+		    do j = 0, nu-1 {
+			if (Memr[yin1+j] < 1 || Memr[yin1+j] > nyin)
+			    Memr[bufout+j] = blank
+		    }
+		}
+	    }
+
+	    # Evaluate output mask pixels and set output bad values.
+
+	    if (pmout != NULL) {
+		bufout = impl2r (pmout, i)
+		if (pmin != NULL) {
+		    call msivector (mzmsi, Memr[xin], Memr[yin], Memr[bufout],
+		        nu)
+		    do j = 0, nu-1 {
+			c = Memr[bufout+j]
+			if (Memr[xin1+j] < 1 || Memr[xin1+j] > nxin ||
+			    Memr[yin1+j] < 1 || Memr[yin1+j] > nyin)
+			    Memr[bufout+j] = MBLANK
+			else if (c > 0.) {
+			    if (c > MTHRESH)
+				Memr[bufout+j] = MBAD
+			    else
+				Memr[bufout+j] = 0
+			}
+		    }
+		} else {
+		    call aclrr (Memr[bufout], nu)
+		    if (xofb) {
+			do j = 0, nu-1 {
+			    if (Memr[xin1+j] < 1 || Memr[xin1+j] > nxin)
+				Memr[bufout+j] = MBLANK
+			}
+		    }
+		    if (yofb) {
+			do j = 0, nu-1 {
+			    if (Memr[yin1+j] < 1 || Memr[yin1+j] > nyin)
+				Memr[bufout+j] = MBLANK
+			}
+		    }
+		}
+	    }
+	}
+
+	# Free memory.
+
+	call mfree (buf, TY_REAL)
+	call mfree (mbuf, TY_REAL)
+	call msifree (zmsi)
+	if (pmin != NULL)
+	    call msifree (mzmsi)
+	call sfree (sp)
+end
+
+
+# TR_BUFL2R -- Maintain buffer of image lines.  A new buffer is created when
+# the buffer pointer is null or if the number of lines requested is changed.
+# The minimum number of image reads is used.
+
+procedure tr_bufl2r (im, pmin, line1, line2, buf, mbuf)
+
+pointer	im		#I Image pointer
+pointer	pmin		#I Mask pointer
+int	line1		#I First image line of buffer
+int	line2		#I Last image line of buffer
+pointer	buf		#U Output data buffer
+pointer	mbuf		#U Output mask buffer
+
+int	i, nlines, nx, last1, last2, nlast
+pointer	buf1, buf2
+
+pointer	imgl2r()
+
+begin
+	nlines = line2 - line1 + 1
+
+	# If the buffer pointer is undefined then allocate memory for the
+	# buffer.  If the number of lines requested changes reallocate
+	# the buffer.  Initialize the last line values to force a full
+	# buffer image read.
+
+	if (buf == NULL) {
+	    nx = IM_LEN(im, 1)
+	    call malloc (buf, nx * nlines, TY_REAL)
+	    if (pmin != NULL)
+		call malloc (mbuf, nx * nlines, TY_REAL)
+	    last1 = line1 - nlines
+	    last2 = line2 - nlines
+	} else if (nlines != nlast) {
+	    call realloc (buf, nx * nlines, TY_REAL)
+	    if (pmin != NULL)
+		call realloc (mbuf, nx * nlines, TY_REAL)
+	    last1 = line1 - nlines
+	    last2 = line2 - nlines
+	}
+
+	# Read only the image lines with are different from the last buffer.
+
+	if (line1 < last1) {
+	    do i = line2, line1, -1 {
+		if (i > last1)
+		    buf1 = buf + (i - last1) * nx
+		else
+		    buf1 = imgl2r (im, i)
+		    
+		buf2 = buf + (i - line1) * nx
+		call amovr (Memr[buf1], Memr[buf2], nx)
+	    }
+	} else if (line2 > last2) {
+	    do i = line1, line2 {
+		if (i < last2)
+		    buf1 = buf + (i - last1) * nx
+		else
+		    buf1 = imgl2r (im, i)
+		    
+		buf2 = buf + (i - line1) * nx
+		call amovr (Memr[buf1], Memr[buf2], nx)
+	    }
+	}
+	if (pmin != NULL) {
+	    if (line1 < last1) {
+		do i = line2, line1, -1 {
+		    if (i > last1)
+			buf1 = mbuf + (i - last1) * nx
+		    else
+			buf1 = imgl2r (pmin, i)
+			
+		    buf2 = mbuf + (i - line1) * nx
+		    call amovr (Memr[buf1], Memr[buf2], nx)
+		    call argtr (Memr[buf2], nx, 0.1, real(MINTERP))
+		}
+	    } else if (line2 > last2) {
+		do i = line1, line2 {
+		    if (i < last2)
+			buf1 = mbuf + (i - last1) * nx
+		    else
+			buf1 = imgl2r (pmin, i)
+			
+		    buf2 = mbuf + (i - line1) * nx
+		    call amovr (Memr[buf1], Memr[buf2], nx)
+		    call argtr (Memr[buf2], nx, 0.1, real(MINTERP))
+		}
+	    }
+	}
+
+	# Save the buffer parameters.
+
+	last1 = line1
+	last2 = line2
+	nlast = nlines
+end
diff --git a/noao/twodspec/longslit/transform/transform.com b/noao/twodspec/longslit/transform/transform.com
new file mode 100644
index 00000000..baaae3ab
--- /dev/null
+++ b/noao/twodspec/longslit/transform/transform.com
@@ -0,0 +1,14 @@
+# TRANSFORM -- Common task parameters.
+
+int	itype			# Interpolation type
+real	u1, v1			# Starting coordinates
+real	u2, v2			# Ending coordinates
+real	du, dv			# Coordinate intervals
+int	nu, nv			# Number of pixels
+bool	ulog, vlog		# Logrithmic coordinates?
+bool	flux			# Conserve flux per pixel?
+bool	usewcs			# Use WCS?
+real	blank			# Blank value
+
+common	/trcom/ u1, v1, u2, v2, du, dv, nu, nv, itype, ulog, vlog,
+		flux, usewcs, blank
diff --git a/noao/twodspec/longslit/transform/trsetup.x b/noao/twodspec/longslit/transform/trsetup.x
new file mode 100644
index 00000000..72db570d
--- /dev/null
+++ b/noao/twodspec/longslit/transform/trsetup.x
@@ -0,0 +1,663 @@
+include	<math.h>
+include	<math/gsurfit.h>
+include	<math/iminterp.h>
+
+# Wrapper for MWCS CT pointer to include the image pixel range.
+
+define	CT_LW		Memi[$1]		# MWCS CT (logical -> world)
+define	CT_WL		Memi[$1+1]		# MWCS CT (world -> logical)
+define	CT_NX		Memi[$1+2]		# Number of pixels in X
+define	CT_NY		Memi[$1+3]		# Number of pixels Y
+
+
+# TR_GSF -- Get coordinate surface fits from the database.
+
+procedure tr_gsf (database, sflist, un, usf, nusf, vsf, nvsf)
+
+char	database		#I Database containing coordinate surfaces
+int	sflist			#I List of user coordinate surfaces
+pointer	un[2]			#O Units pointers
+pointer	usf			#O Pointer to array of U surface fits
+int	nusf			#O Number of U surface fits
+pointer	vsf			#O Pointer to array of V surface fits
+int	nvsf			#O Number of U surface fits
+
+int	i, nsf
+pointer	sp, sfname, un1, sf
+
+bool	un_compare()
+int	clgfil(), clplen()
+
+begin
+	# Get the user coordinate surfaces and separate them into U and V.
+	# Check that all surfaces have the same range of X and Y and determine
+	# the range of U and V.
+
+	call smark (sp)
+	call salloc (sfname, SZ_FNAME, TY_CHAR)
+
+	nsf = max (1, clplen (sflist))
+	call malloc (usf, nsf, TY_INT)
+	call malloc (vsf, nsf, TY_INT)
+
+	un[1] = NULL
+	un[2] = NULL
+	Memi[usf] = NULL
+	Memi[vsf] = NULL
+	nusf = 0
+	nvsf = 0
+	while (clgfil (sflist, Memc[sfname], SZ_FNAME) != EOF) {
+	    call lm_dbread (database, Memc[sfname], i, un1, sf)
+	    if (un1 != NULL) {
+		if (un[i] == NULL)
+		    un[i] = un1
+		else if (un_compare (un1, un[i]))
+		    call un_close (un1)
+		else {
+		    call un_close (un1)
+		    call un_close (un[i])
+		    call sfree (sp)
+		    call error (1, "Input units disagree")
+		}
+	    }
+
+	    if (sf != NULL) {
+		if (i == 1) {
+		    nusf = nusf+1
+		    Memi[usf+nusf-1] = sf
+		} else if (i == 2) {
+		    nvsf = nvsf+1
+		    Memi[vsf+nvsf-1] = sf
+		}
+	    }
+	}
+	call clprew (sflist)
+
+	if (nusf + nvsf == 0)
+	    call error (0, "No user coordinates")
+
+	call sfree (sp)
+end
+
+
+# TR_GWCS -- Get WCS.
+
+procedure tr_gwcs (mw, un, nx, ny, ct, usf, nusf, vsf, nvsf)
+
+pointer	mw			#I MWCS pointer
+pointer	un[2]			#O Units pointers
+int	nx, ny			#I Image size
+
+pointer	ct			#O CT pointer
+pointer	usf			#O Pointer to array of U surface fits
+int	nusf			#O Number of U surface fits
+pointer	vsf			#O Pointer to array of V surface fits
+int	nvsf			#O Number of U surface fits
+
+int	i
+pointer	sp, units, un_open(), mw_sctran()
+errchk	un_open
+
+begin
+	call smark (sp)
+	call salloc (units, SZ_FNAME, TY_CHAR)
+
+	call malloc (ct, 4, TY_STRUCT)
+	nusf = 1
+	call calloc (usf, nusf, TY_INT) 
+	nvsf = 1
+	call calloc (vsf, nvsf, TY_INT) 
+
+	CT_LW(ct) =  mw_sctran (mw, "logical", "world", 3)
+	CT_WL(ct) =  mw_sctran (mw, "world", "logical", 3)
+	CT_NX(ct) = nx
+	CT_NY(ct) = ny
+
+	do i = 1, 2 {
+	    ifnoerr (call mw_gwattrs (mw, i, "units", Memc[units], SZ_FNAME))
+	        un[i] = un_open (Memc[units])
+	    else
+		un[i] = NULL
+	}
+end
+
+
+# TR_SETUP -- Setup the transformation interpolation.
+# 
+# At each point (U,V) in the output image we need to know the coordinate
+# (X,Y) of the input images to be interpolated.  This means we need
+# to determine X(U,V) and Y(U,V).  The input user coordinate surfaces,
+# however, are U(X,Y) and V(X,Y) (a missing surface implies a one to one
+# mapping of U=X or V=Y).  This requires simultaneously inverting the user
+# coordinate surfaces.  This is a slow process using a gradient following
+# iterative technique.
+# 
+# Note that when an WCS is used, the MWCS routines already provide the
+# inverse mapping.  But even in this case it may be slow and so we use the
+# same sampling and surface fitting technique for setting up the inversion
+# mapping.
+# 
+# The inverted coordinates are determined on a evenly subsampled grid of
+# linear output coordinates.  A linear interpolation surface can then be fit
+# to this grid which is much faster to evaluate at each output coordinate.
+# These interpolation surfaces are returned.  If flux is to be conserved a
+# similar interpolation surface for the Jacobian, J(U,V) is also returned.
+# There may also be a mapping of the output image into logrithmic intervals
+# which maps to the linearly sampled interpolation surfaces.  The mappings
+# of the output U and V intervals to the subsampled interpolation coordinates
+# are also returned.
+# 
+# 1. Set the output coordinate system based on the ranges of X, Y, U, and V.
+# 2. Determine X(U,V), Y(U,V), and J(U,V) on a evenly subsampled grid of
+#    U and V.
+# 3. Fit linear interpolation surfaces to these data.
+# 4. Compute the mapping between output coordinates along each axis, which
+#    may be logrithmic, into the subsampling interpolation coordinates.
+
+procedure tr_setup (ct, usf, nusf, vsf, nvsf, un, xmsi, ymsi, jmsi,
+	uout, vout, duout, dvout)
+
+pointer	ct			#I CT pointer
+pointer	usf			#U Pointers to U surface fits: freed upon return
+int	nusf			#I Number of U surface fits
+pointer	vsf			#U Pointers to V surface fits: freed upon return
+int	nvsf			#I Number of V surface fits
+pointer	un[2]			#O Units pointers
+pointer	xmsi, ymsi, jmsi	#O Surface interpolators for X, Y and Jacobian
+pointer	uout, vout		#O Output coordinates relative to interpolator 
+pointer	duout, dvout		#O Output coordinate intervals
+
+int	i, j, step, nu1, nv1
+real	xmin, xmax, ymin, ymax, umin, umax, vmin, vmax
+real	u, v, x, y, du1, dv1, der[8]
+double	dval
+pointer	xgrid, ygrid, zgrid, ptr1, ptr2, ptr3
+
+real	tr_getr(), tr_eval()
+
+include	"transform.com"
+
+begin
+	#step = clgeti ("step")
+	step = 10
+
+	xmin = INDEF
+	xmax = INDEF
+	ymin = INDEF
+	ymax = INDEF
+	umin = INDEF
+	umax = INDEF
+	vmin = INDEF
+	vmax = INDEF
+	do i = 1, nusf {
+	    if (IS_INDEF (xmin)) {
+		xmin = tr_getr (ct, Memi[usf+i-1], GSXMIN)
+		xmax = tr_getr (ct, Memi[usf+i-1], GSXMAX)
+		ymin = tr_getr (ct, Memi[usf+i-1], GSYMIN)
+		ymax = tr_getr (ct, Memi[usf+i-1], GSYMAX)
+	    } else {
+		if ((xmin != tr_getr (ct, Memi[usf+i-1], GSXMIN)) ||
+		    (xmax != tr_getr (ct, Memi[usf+i-1], GSXMAX)) ||
+		    (ymin != tr_getr (ct, Memi[usf+i-1], GSYMIN)) ||
+		    (ymax != tr_getr (ct, Memi[usf+i-1], GSYMAX)))
+		    call error (0, "tr_setup: Inconsistent coordinate fits")
+	    }
+
+	    if (IS_INDEF (umin)) {
+	        umin = tr_eval (ct, Memi[usf+i-1], 1, xmin, ymin)
+	        umax = umin
+	    }
+	    u = tr_eval (ct, Memi[usf+i-1], 1, xmin, ymin)
+	    umin = min (u, umin)
+	    umax = max (u, umax)
+	    u = tr_eval (ct, Memi[usf+i-1], 1, xmax, ymin)
+	    umin = min (u, umin)
+	    umax = max (u, umax)
+	    u = tr_eval (ct, Memi[usf+i-1], 1, xmin, ymax)
+	    umin = min (u, umin)
+	    umax = max (u, umax)
+	    u = tr_eval (ct, Memi[usf+i-1], 1, xmax, ymax)
+	    umin = min (u, umin)
+	    umax = max (u, umax)
+	}
+	do i = 1, nvsf {
+	    if (IS_INDEF (xmin)) {
+		xmin = tr_getr (ct, Memi[vsf+i-1], GSXMIN)
+		xmax = tr_getr (ct, Memi[vsf+i-1], GSXMAX)
+		ymin = tr_getr (ct, Memi[vsf+i-1], GSYMIN)
+		ymax = tr_getr (ct, Memi[vsf+i-1], GSYMAX)
+	    } else {
+		if ((xmin != tr_getr (ct, Memi[vsf+i-1], GSXMIN)) ||
+		    (xmax != tr_getr (ct, Memi[vsf+i-1], GSXMAX)) ||
+		    (ymin != tr_getr (ct, Memi[vsf+i-1], GSYMIN)) ||
+		    (ymax != tr_getr (ct, Memi[vsf+i-1], GSYMAX)))
+		    call error (0, "tr_setup: Inconsistent coordinate fits")
+	    }
+
+	    if (IS_INDEF (vmin)) {
+	        vmin = tr_eval (ct, Memi[vsf+i-1], 2, xmin, ymin)
+	        vmax = vmin
+	    }
+	    v = tr_eval (ct, Memi[vsf+i-1], 2, xmin, ymin)
+	    vmin = min (v, vmin)
+	    vmax = max (v, vmax)
+	    v = tr_eval (ct, Memi[vsf+i-1], 2, xmax, ymin)
+	    vmin = min (v, vmin)
+	    vmax = max (v, vmax)
+	    v = tr_eval (ct, Memi[vsf+i-1], 2, xmin, ymax)
+	    vmin = min (v, vmin)
+	    vmax = max (v, vmax)
+	    v = tr_eval (ct, Memi[vsf+i-1], 2, xmax, ymax)
+	    vmin = min (v, vmin)
+	    vmax = max (v, vmax)
+	}
+	if (IS_INDEF (umin)) {
+	    umin = xmin
+	    umax = xmax
+	}
+	if (IS_INDEF (vmin)) {
+	    vmin = ymin
+	    vmax = ymax
+	}
+
+	# Set the output coordinate system which is in a common block.
+	call tr_setoutput (xmin, xmax, ymin, ymax, umin, umax, vmin, vmax)
+
+	# Subsample the inverted coordinates and fit an interpolation
+	# surface.  The grid is evaluated in a back and forth pattern to
+	# use the last point evaluated and the starting point for the next
+	# point.  This allows the interative inversion routine to work most
+	# efficiently with typically only two evaluations per step.
+
+	nu1 = max (2, nu / step)
+	nv1 = max (2, nv / step)
+	du1 = (u2 - u1) / (nu1 - 1)
+	dv1 = (v2 - v1) / (nv1 - 1)
+
+	call malloc (xgrid, nu1 * nv1, TY_REAL)
+	call malloc (ygrid, nu1 * nv1, TY_REAL)
+	call malloc (zgrid, nu1 * nv1, TY_REAL)
+
+	call tr_init (ct, Memi[usf], nusf, Memi[vsf], nvsf, xmin, ymin, der)
+	do i = 1, nv1, 2 {
+	    # Do this line from left to right.
+	    ptr1 = xgrid + (i - 1) * nu1 - 1
+	    ptr2 = ygrid + (i - 1) * nu1 - 1
+	    ptr3 = zgrid + (i - 1) * nu1 - 1
+	    v = v1 + (i - 1) * dv1
+	    do j = 1, nu1 {
+		u = u1 + (j - 1) * du1
+		call tr_invert (ct, Memi[usf], nusf, Memi[vsf], nvsf, u, v,
+		    x, y, der, xmin, xmax, ymin, ymax)
+		# V2.10.2
+		#Memr[ptr1+j] = der[1]
+		#Memr[ptr2+j] = der[2]
+		# After V2.10.3
+		Memr[ptr1+j] = x
+		Memr[ptr2+j] = y
+
+		Memr[ptr3+j] = 1. / abs (der[4] * der[8] - der[5] * der[7])
+	    }
+	    if (i == nv1)
+		break
+
+	    # Do the next line from right to left.
+	    ptr1 = xgrid + i * nu1 - 1
+	    ptr2 = ygrid + i * nu1 - 1
+	    ptr3 = zgrid + i * nu1 - 1
+	    v = v1 + i * dv1
+	    do j = nu1, 1, -1 {
+		u = u1 + (j - 1) * du1
+		call tr_invert (ct, Memi[usf], nusf, Memi[vsf], nvsf, u, v,
+		    x, y, der, xmin, xmax, ymin, ymax)
+		# V2.10.2
+		#Memr[ptr1+j] = der[1]
+		#Memr[ptr2+j] = der[2]
+		# V2.10.3
+		Memr[ptr1+j] = x
+		Memr[ptr2+j] = y
+		Memr[ptr3+j] = 1. / abs (der[4] * der[8] - der[5] * der[7])
+	    }
+	}
+
+	# Free the surfaces since we are now done with them.
+	if (ct != NULL)
+	    call mfree (ct, TY_STRUCT)
+	for (i=1; i<=nusf; i=i+1)
+	    if (Memi[usf+i-1] != NULL)
+		call xgsfree (Memi[usf+i-1])
+	call mfree (usf, TY_POINTER)
+	for (i=1; i<=nvsf; i=i+1)
+	    if (Memi[vsf+i-1] != NULL)
+		call xgsfree (Memi[vsf+i-1])
+	call mfree (vsf, TY_POINTER)
+
+	# Fit a linear interpolator to the subsampled grids of X(U,V), Y(U,V),
+	# and J(U,V) to avoid having to evaluate the inverse at each point in
+	# the output image.  The inversion is slow because of the many
+	# evaluations of the surfaces coordinates.  Also compute an return
+	# arrays mapping the output coordinates to the subsampled coordinates.
+	# This may include a transformation to logrithmic intervals.
+
+	call msiinit (xmsi, II_BILINEAR)
+	call msifit (xmsi, Memr[xgrid], nu1, nv1, nu1)
+	call mfree (xgrid, TY_REAL)
+
+	call msiinit (ymsi, II_BILINEAR)
+	call msifit (ymsi, Memr[ygrid], nu1, nv1, nu1)
+	call mfree (ygrid, TY_REAL)
+
+	if (flux) {
+	    call msiinit (jmsi, II_BILINEAR)
+	    call msifit (jmsi, Memr[zgrid], nu1, nv1, nu1)
+	}
+	call mfree (zgrid, TY_REAL)
+
+	# Compute the mapping between output coordinates and the subsampled
+	# interpolation surface.  Also compute the intervals used to define
+	# the pixel areas for conserving flux.
+
+	call malloc (uout, nu, TY_REAL)
+	call malloc (duout, nu, TY_REAL)
+	if (ulog) {
+	    dval = log10 (double(u1))
+	    do i = 0, nu - 1
+		Memr[uout+i] = 10.**(dval+i*du)
+	    call amulkr (Memr[uout], du * LN_10, Memr[duout], nu)
+	} else {
+	    do i = 0, nu - 1
+		Memr[uout+i] = u1 + i * du
+	    call amovkr (du, Memr[duout], nu)
+	}
+	u2 = Memr[uout+nu-1]
+
+	call malloc (vout, nv, TY_REAL)
+	call malloc (dvout, nv, TY_REAL)
+	if (vlog) {
+	    dval = log10 (double(v1))
+	    do i = 0, nv - 1
+		Memr[vout+i] = 10.**(dval+i*dv)
+	    call amulkr (Memr[vout], dv * LN_10, Memr[dvout], nv)
+	} else {
+	    do i = 0, nv - 1
+		Memr[vout+i] = v1 + i * dv
+	    call amovkr (dv, Memr[dvout], nv)
+	}
+	v2 = Memr[vout+nv-1]
+
+	# Convert to interpolation coordinates.
+	umin = 1.; umax = nu
+	do i = 0, nu - 1
+	    Memr[uout+i] = max (umin, min (umax, (Memr[uout+i]-u1)/du1+1))
+	vmin = 1.; vmax = nv
+	do i = 0, nv - 1
+	    Memr[vout+i] = max (vmin, min (vmax, (Memr[vout+i]-v1)/dv1+1))
+end
+
+
+define	MAX_ITERATE	10
+define	ERROR		0.05
+define	FUDGE		0.5
+
+# TR_INVERT -- Given user coordinate surfaces U(X,Y) and V(X,Y)
+# (if none use one-to-one mapping and if more than one average)
+# corresponding to a given U and V and also the various partial
+# derivatives.  This is done using a gradient following interative
+# method based on evaluating the partial derivative at each point
+# and solving the linear Taylor expansions simultaneously.  The last
+# point sampled is used as the starting point.  Thus, if the
+# input U and V progress smoothly then the number of iterations
+# can be small.  The output is returned in x and y and in the derivative array
+# DER.  A point outside of the surfaces is returned as the nearest
+# point at the edge of the surfaces in the DER array.
+#
+# If a WCS is used then we let MWCS do the inversion and compute the
+# derivatives numerically.
+
+procedure tr_invert (ct, usf, nusf, vsf, nvsf, u, v, x, y, der,
+	xmin, xmax, ymin, ymax)
+
+pointer	ct			#I CT pointer
+pointer	usf[ARB], vsf[ARB]	#I User coordinate surfaces U(X,Y) and V(X,Y)
+int	nusf, nvsf		#I Number of surfaces for each coordinate
+real	u, v			#I Input U and V to determine X and Y
+real	x, y			#O Output X and Y
+real	der[8]			#U Last result as input, new result as output 
+				#  1=X, 2=Y, 3=U, 4=DUDX, 5=DUDY, 6=V,
+				#  7=DVDX, 8=DVDY
+real	xmin, xmax, ymin, ymax	#I Limits of coordinate surfaces.
+
+int	i, j, nedge
+real	fudge, du, dv, dx, dy, a, b, tmp[4]
+
+begin
+	# If using a WCS we let MWCS do the inversion.
+	if (ct != NULL) {
+	    call mw_c2tranr (CT_WL(ct), u, v, x, y)
+	    call mw_c2tranr (CT_LW(ct), x-0.5, y, tmp[1], tmp[3])
+	    call mw_c2tranr (CT_LW(ct), x+0.5, y, tmp[2], tmp[4])
+	    der[4] = tmp[2] - tmp[1]
+	    der[7] = tmp[4] - tmp[3]
+	    call mw_c2tranr (CT_LW(ct), x, y-0.5, tmp[1], tmp[3])
+	    call mw_c2tranr (CT_LW(ct), x, y+0.5, tmp[2], tmp[4])
+	    der[5] = tmp[2] - tmp[1]
+	    der[8] = tmp[4] - tmp[3]
+	    return
+	}
+
+	# Use the last result as the starting point for the next position.
+	# If this is near the desired value then the interation will converge
+	# quickly.  Allow a iteration to go off the surface twice.
+	# Quit when DX and DY are within ERROR.
+
+	nedge = 0
+	do i = 1, MAX_ITERATE {
+	    du = u - der[3]
+	    dv = v - der[6]
+	    a = der[8] * du - der[5] * dv
+	    b = der[8] * der[4] - der[5] * der[7]
+	    if (b == 0.) {
+		if (a < 0.)
+		    dx = -2.
+		else
+		    dx = 2.
+	    } else
+	        dx = a / b
+	    a = dv - der[7] * dx
+	    b = der[8]
+	    if (b == 0.) {
+		if (a < 0.)
+		    dy = -2.
+		else
+		    dy = 2.
+	    } else
+	        dy = a / b
+	    fudge = 1 - FUDGE / i
+	    x = der[1] + fudge * dx
+	    y = der[2] + fudge * dy
+	    der[1] = max (xmin, min (xmax, x))
+	    der[2] = max (ymin, min (ymax, y))
+#	    if (x < xmin || x > xmax)
+#	        nedge = nedge + 1
+#	    if (y < ymin || y > ymax)
+#	        nedge = nedge + 1
+#	    if (nedge > 2)
+#	        break
+	    if ((abs (dx) < ERROR) && (abs (dy) < ERROR))
+	        break
+
+	    if (nusf == 0)
+		der[3] = der[1]
+	    else if (nusf == 1) {
+	        call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
+	        call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
+	        call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
+	    } else {
+	        call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
+	        call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
+	        call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
+		do j = 2, nusf {
+	            call xgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0)
+	            call xgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0)
+	            call xgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1)
+		    der[3] = der[3] + tmp[1]
+		    der[4] = der[4] + tmp[2]
+		    der[5] = der[5] + tmp[3]
+		}
+		der[3] = der[3] / nusf
+		der[4] = der[4] / nusf
+		der[5] = der[5] / nusf
+	    }
+
+	    if (nvsf == 0)
+		der[6] = der[2]
+	    else if (nvsf == 1) {
+	        call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
+	        call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
+	        call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
+	    } else {
+	        call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
+	        call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
+	        call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
+		do j = 2, nvsf {
+	            call xgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0)
+	            call xgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0)
+	            call xgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1)
+		    der[6] = der[6] + tmp[1]
+		    der[7] = der[7] + tmp[2]
+		    der[8] = der[8] + tmp[3]
+		}
+		der[6] = der[6] / nvsf
+		der[7] = der[7] / nvsf
+		der[8] = der[8] / nvsf
+	    }
+	}
+end
+
+
+# TR_INIT -- Since the inversion iteration always begins from the last
+# point we need to initialize before the first call to TR_INVERT.
+# When using a WCS this simply returns.
+
+procedure tr_init (ct, usf, nusf, vsf, nvsf, x, y, der)
+
+pointer	ct			#I CT pointer
+pointer	usf[ARB], vsf[ARB]	#I User coordinate surfaces
+int	nusf, nvsf		#I Number of surfaces for each coordinate
+real	x, y			#I Starting X and Y
+real	der[8]			#O Inversion data
+
+int	j
+real	tmp[3]
+
+begin
+	if (ct != NULL)
+	    return
+
+	der[1] = x
+	der[2] = y
+	if (nusf == 0) {
+	    der[3] = der[1]
+	    der[4] = 1.
+	    der[5] = 0.
+	} else if (nusf == 1) {
+	    call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
+	    call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
+	    call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
+	} else {
+	    call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
+	    call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
+	    call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
+	    do j = 2, nusf {
+	        call xgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0)
+	        call xgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0)
+	        call xgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1)
+		der[3] = der[3] + tmp[1]
+		der[4] = der[4] + tmp[2]
+		der[5] = der[5] + tmp[3]
+	    }
+	    der[3] = der[3] / nusf
+	    der[4] = der[4] / nusf
+	    der[5] = der[5] / nusf
+	}
+
+	if (nvsf == 0) {
+	    der[6] = der[2]
+	    der[7] = 0.
+	    der[8] = 1.
+	} else if (nvsf == 1) {
+	    call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
+	    call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
+	    call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
+	} else {
+	    call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
+	    call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
+	    call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
+	    do j = 2, nvsf {
+	        call xgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0)
+	        call xgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0)
+	        call xgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1)
+		der[6] = der[6] + tmp[1]
+		der[7] = der[7] + tmp[2]
+		der[8] = der[8] + tmp[3]
+	    }
+	    der[6] = der[6] / nvsf
+	    der[7] = der[7] / nvsf
+	    der[8] = der[8] / nvsf
+	}
+end
+
+
+# TR_EVAL -- Evalute coordinate function.
+#
+# This is an interface routine to allow using either an MWCS CT (coordinate
+# transform) pointer or a GSURFIT SF (2D surface function) pointer.  The
+# surface method is used with a FITCOORDS database.  The MWCS method is
+# used to retransform an image with a WCS.
+
+real procedure tr_eval (ct, sf, axis, x, y)
+
+pointer	ct			#I CT pointer
+pointer	sf			#I SF pointer
+int	axis			#I World coordinate axis to return
+real	x, y			#I Pixel coordinate to transform
+
+real	w[2], xgseval()
+
+begin
+	if (sf != NULL)
+	    return (xgseval (sf, x, y))
+	
+	call mw_c2tranr (CT_LW(ct), x, y, w[1], w[2])
+	return (w[axis])
+end
+
+
+# TR_GETR -- Get real valued parameter.
+#
+# This is an interface routine to allow using either an MWCS CT (coordinate
+# transform) pointer or a GSURFIT SF (2D surface function) pointer.  The
+# surface method is used with a FITCOORDS database.  The MWCS method is
+# used to retransform an image with a WCS.
+
+real procedure tr_getr (ct, sf, param)
+
+pointer	ct			#I CT pointer
+pointer	sf			#I SF pointer
+int	param			#I Parameter code
+
+real	xgsgetr()
+
+begin
+	if (sf != NULL)
+	    return (xgsgetr (sf, param))
+
+	switch (param) {
+	case GSXMIN, GSYMIN: 
+	    return (real (1))
+	case GSXMAX:
+	    return (real (CT_NX(ct)))
+	case GSYMAX:
+	    return (real (CT_NY(ct)))
+	}
+end