Initial commit

33 files changed, 4311 insertions, 0 deletions

diff --git a/pkg/xtools/skywcs/doc/README b/pkg/xtools/skywcs/doc/README
new file mode 100644
index 00000000..b0998629
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/README
@@ -0,0 +1,301 @@
+                SKYWCS: The Sky Coordinates Package
+
+1. Introduction
+
+    The skywcs package contains a simple set of routines for managing sky
+coordinate information and for transforming from one sky coordinate system to
+another. The sky coordinate system is defined either by a system name, e.g.
+"J2000", "galactic", etc., or by an image system name, e.g. "dev$ypix" or
+"dev$ypix world".
+
+    The skywcs routine are layered on the Starlink Positional Astronomy library
+SLALIB which is installed in the IRAF MATH package. Type "help slalib option=
+sys" for more information about SLALIB.
+
+
+2. The Interface Routines
+
+The package prefix is sk. The interface routines are listed below.
+
+   stat = sk_decwcs (ccsystem, mw, coo, imcoo)
+   stat = sk_decwstr (ccsystem, coo, imcoo)
+   stat = sk_decim (im, wcs, mw, coo)
+	  sk_enwcs (coo, ccsystem, maxch)
+ newcoo = sk_copy (coo)
+          sk_iiprint (label, imagesys, mw, coo) 
+          sk_iiwrite (fd, label, imagesys, mw, coo) 
+[id]val = sk_stat[id] (coo, param)
+          sk_stats (coo, param, str, maxch)
+	  sk_set[id] (coo, param, [id]val)
+          sk_sets (coo, param, str)
+	  sk_ultran (incoo, outcoo, ilng, ilat, olng, olat, npts)
+	  sk_lltran (incoo, outoo, ilng, ilat, ipmlng, ipmlat, px, rv,
+             olng, olat)
+	  sk_equatorial (incoo, outcoo, ilng, ilat, ipmlng, ipmlat, px,
+	     rv, olng, olat)
+	  sk_saveim (coo, mw, im)
+	  sk_close (coo)
+
+
+3. Notes
+
+    An "include <pkg/skywcs.h>" statement must appear in the calling program
+to make the skywcs package parameter definitions visible to the calling
+program.
+
+    A "-lxtools -lslalib" must be included in the calling program link line
+to link in the skywcs and the slalib routines.
+
+    The sky coordinate descriptor is created with a call to one of the
+sk_decwcs, sk_decwstr, or sk_imwcs routines. If the source of the sky
+coordinate descriptor is an image then an IRAF MWCS descriptor will be returned
+with the sky oordinate descriptor. The sky coordinate descriptor is freed with a
+call to sk_close. A separate call to mw_close must be made to free the MWCS
+descriptor if one was allocated.
+
+    By default the main skywcs coordinate transformation routine sk_ultran
+assumes that the input and output sky coordinates are in hours and degrees
+if the input and output coordinate systems are equatorial, otherwise the
+coordinates are assumed to be in degrees and degrees. The default input and
+output sky coordinate units can be reset with calls to sk_seti. Two lower level
+coordinate transformations for handling proper motions sk_lltran and
+sk_equatorial are also available. These routines assume that the input and
+output coordinates and proper motions are in radians.
+
+    Calling programs working with both sky coordinate and MWCS descriptors
+need to be aware that the MWCS routines assume that all sky coordinates
+will be input and output in degrees and adjust their code accordingly.
+
+    The skywcs routine sk_saveim can be used to update an image header.
+
+
+3. Examples
+
+Example 1: Convert from B1950 coordinates to J2000 coordinates.
+
+    include <skywcs.h>
+
+    ....
+
+    # Open input coordinate system.
+    instat = sk_decwstr ("B1950", incoo, NULL)
+    if (instat == ERR) {
+	call sk_close (incoo)
+	return
+    }
+
+    # Open output coordinate system.
+    outstat = sk_decwstr ("J2000", outcoo, NULL)
+    if (outstat == ERR) {
+	call sk_close (outcoo)
+	return
+    }
+
+    # Do the transformation assuming the input coordinates are in hours
+    # and degrees. The output coordinates will be in hours and degrees
+    # as well.
+    call sk_ultran (incoo, outcoo, rain, decin, raout, decout, npts) 
+
+    # Close the coordinate descriptors.
+    call sk_close (incoo)
+    call sk_close (outcoo)
+
+    ...
+
+
+Example 2: Repeat example 1 but convert to galactic coordinates.
+
+    include <skywcs.h>
+
+    ....
+
+    # Open the input coordinate system.
+    instat = sk_decwstr ("B1950", incoo, NULL)
+    if (instat == ERR) {
+	call sk_close (incoo)
+	return
+    }
+
+    # Open the output coordinate system.
+    outstat = sk_decwstr ("galactic", outcoo, NULL)
+    if (outstat == ERR) {
+	call sk_close (outcoo)
+	return
+    }
+
+    # Dd the transformation assuming the input coordinates are in hours and
+    # degrees. The output coordinates will be in degrees and degrees.
+    call sk_ultran (incoo, outcoo, rain, decin, raout, decout, npts) 
+
+    # Close the coordinate descriptors.
+    call sk_close (incoo)
+    call sk_close (outcoo)
+
+    ...
+
+Example 3: Convert a grid of pixel coordinates in the input image to the
+equivalent pixel coordinate in the output image using the image world
+coordinate systems to connect the two.
+
+    include <skywcs.h>
+
+    ....
+
+    # Mwref will be defined because the input system is an image.
+    refstat = sk_decwcs ("refimage logical", mwref, refcoo, NULL)
+    if (refstat == ERR || mwref == NULL) {
+	if (mwref != NULL)
+	    call mw_close (mwref)
+	call sk_close (refcoo)
+	return
+    }
+
+    # Set the reference coordinate descriptor so it expects input in degrees
+    # and degrees.
+    call sk_seti (refcoo, S_NLNGUNUTS, SKY_DEGREES)
+    call sk_seti (refcoo, S_NLATUNUTS, SKY_DEGREES)
+
+    # Mwout will be defined because the output system is an image.
+    outstat = sk_decwcs ("image logical", mwout, outcoo, NULL)
+    if (outstat == ERR || mwout == NULL) {
+	if (mwout != NULL)
+	    call mw_close (mwout)
+	call sk_close (outcoo)
+	call mw_close (mwref)
+	call sk_close (refcoo)
+	return
+    }
+
+    # Set the output coordinate descriptor so it will output coordinates
+    # in degrees and degrees.
+    call sk_seti (outcoo, S_NLNGUNUTS, SKY_DEGREES)
+    call sk_seti (outcoo, S_NLATUNUTS, SKY_DEGREES)
+
+    # Compute pixel grid in refimage and store coordinate in the arrays
+    # xref and yref.
+    npts = 0
+    do j = 1, IM_LEN(im,2), 100 {
+        do i = 1, IM_LEN(im,1), 100 {
+	    npts = npts + 1
+	    xref[npts] = i
+	    yref[npts] = j
+	}
+    }
+
+    # Convert xref and yref to celestial coordinates raref and decref using
+    # mwref. The output coordinates will be in degrees and degrees.
+    ctref = mw_sctran (mwref, "logical", "world", 03B)
+    do i = 1, npts
+	call mw_c2trand (ctref, xref[i], yref[i], raref[i], decref[i])
+    call ct_free (ctref)
+
+    # Convert the reference celestial coordinates to the output celestial
+    # coordinate system using the coordinate descriptors.
+    call sk_ultran (refcoo, outcoo, raref, decref, raout, decout, npts)
+
+    # Convert the output celestial coordinates to  pixel coordinates in
+    # the other image using mwout.
+    ctout = mw_sctran (mwout, "world", "logical", 03B)
+    do i = 1, npts
+	call mw_c2trand (ctout, raout[i], decout[i], xout[i], yout[i])
+    call ct_free (ctout)
+
+    # Print the input and output pixel coordinates.
+    do i = 1, npts {
+	call printf ("%10.3f %10.3f  %10.3f %10.3f\n")
+	    call pargd (xref[i])
+	    call pargd (yref[i])
+	    call pargd (xout[i])
+	    call pargd (yout[i])
+    }
+
+    # Tidy up.
+    call mw_close (mwref)
+    call mw_close (mwout)
+    call sk_close (refcoo)
+    call sk_close (outcoo)
+
+
+Example 4: Convert a 2D image with an J2000 tangent plane projection wcs to the
+equivalent galactic wcs. The transformation requires a shift in origin and a
+rotation. Assume that the ra axis is 1 and the dec axis is 2. The details of
+how to compute the rotation are not shown here. See the imcctran task for
+details.
+
+    include <mwset.h> 
+    include <skywcs.h>
+
+    ...
+
+    # Open image.
+    im = immap (image, READ_WRITE, 0)
+
+    # Open the image coordinate system.
+    instat = sk_decim (im, "logical", mwin, cooin)
+    if (instat == ERR || mwin == NULL) {
+	...
+	call sk_close (cooin)
+	...
+    }
+
+    # Get the dimensions of the mwcs descriptor. This should be 2.
+    ndim = mw_ndim (mwin, MW_NPHYSDIM)
+
+    # Get the default coordinates to degrees and degreees.
+    call sk_seti (cooin, S_NLNGUNITS, SKY_DEGREES)
+    call sk_seti (cooin, S_NATGUNITS, SKY_DEGREES)
+
+    # Open the output coordinate system. Mwout is NULL because this system
+    # is not an image.
+    outstat = sk_decwstr ("galactic", mwout, cooout, cooin)
+    if (outstat == ERR) {
+        ...
+	call sk_close (outstat)
+	...
+    }
+
+    # Make a copy of the mwcs descriptor.
+    mwout = mw_newcopy (mwin)
+
+    # Allocate space for the r and w vectors and cd matrix.
+    call malloc (r, ndim, TY_DOUBLE)
+    call malloc (w, ndim, TY_DOUBLE)
+    call malloc (cd, ndim * ndim, TY_DOUBLE)
+    call malloc (newcd, ndim * ndim, TY_DOUBLE)
+
+    # Assume for simplicty that the MWCS LTERM is the identify transform.
+    # so we don't have to worry about it. Get the WTERM which consists
+    # of r the reference point in pixels, w the reference point in degrees,
+    # and the cd matrix in degrees per pixel.
+    call mw_gwtermd (mwin, Memd[r], Memd[w], Memd[cd], ndim)
+
+    # Convert the world coordinates zero point. The pixel zero point
+    # remains the same.
+    tilng = Memd[w]
+    tilat = Memd[w+1]
+    call sk_ultran (incoo, outcoo, tilng, tilat, tolng, tolat, 1)
+    Memd[w] = tolng
+    Memd[w+1] = tolat
+
+    # Figure out how much to rotate the coordinate system and edit the
+    # compute a new CD matrix. Call it newcd.
+    ...
+
+    # Enter the new CD matrix and zero point.
+    call mw_swterm (mwout, Memd[r], Memd[w], Memd[newcd], ndim)
+
+    # Update the header. 
+    call sk_saveim (cooout, mwout, im)
+    call mw_saveim (mwout, im)
+    ...
+
+    # Tidy up.
+    call mfree (r, TY_DOUBLE)
+    call mfree (w, TY_DOUBLE)
+    call mfree (cd, TY_DOUBLE)
+    call mfree (newcd, TY_DOUBLE)
+    call mw_close (mwin)
+    call mw_close (mwout)
+    call sk_close (cooin)
+    call sk_close (cooout)
+    call imunmap (im)
diff --git a/pkg/xtools/skywcs/doc/ccsystems.hlp b/pkg/xtools/skywcs/doc/ccsystems.hlp
new file mode 100644
index 00000000..63a2fde6
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/ccsystems.hlp
@@ -0,0 +1,134 @@
+.help ccsystems Mar00 Skywcs
+.ih
+NAME
+ccsystems -- list and describe the supported sky coordinate systems
+.ih
+USAGE
+help ccsystems
+
+.ih
+SKY COORDINATE SYSTEMS
+
+The sky package supports the equatorial ("fk4", "fk4-noe", "fk5", "icrs"),
+ecliptic, galactic, and supergalactic celestial coordinate systems. In most
+cases and unless otherwise noted users can input their coordinates in
+any one of these systems as long as they specify the coordinate system
+correctly.
+
+Considerable flexibility is permitted in how the coordinate systems are
+specified, e.g. J2000.0, j2000.0, 2000.0, fk5, fk5 J2000, and fk5 2000.0
+all specify the mean place post-IAU 1976 or FK5 system. Missing equinox and
+epoch fields assume reasonable defaults. In most cases the
+systems of most interest to users are "icrs", "j2000", and "b1950"
+which stand for the ICRS J2000.0, FK5 J2000.0 and FK4 B1950.0 celestial
+coordinate systems respectively. The full set of options are listed below:
+
+.ls equinox [epoch]
+The equatorial mean place post-IAU 1976 (FK5) system if equinox is a
+Julian epoch, e.g. J2000.0 or 2000.0, or the equatorial mean place
+pre-IAU 1976 system (FK4) if equinox is a Besselian epoch, e.g. B1950.0
+or 1950.0. Julian equinoxes are prefixed by a J or j, Besselian equinoxes
+by a B or b. Equinoxes without the J / j or B / b prefix are treated as
+Besselian epochs if they are < 1984.0, Julian epochs if they are >= 1984.0.
+Epoch is the epoch of the observation and may be a Julian
+epoch, a Besselian epoch, or a Julian date. Julian epochs
+are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to the epoch type of
+equinox if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date.  If undefined epoch defaults to equinox.
+.le
+.ls icrs [equinox] [epoch]
+The International Celestial Reference System where equinox is
+a Julian or Besselian epoch e.g. J2000.0  or B1980.0.
+Equinoxes without the J / j or B / b prefix are treated as Julian epochs.
+The default value of equinox is J2000.0.
+Epoch is a Besselian epoch, a Julian epoch, or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Julian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date.  If undefined epoch defaults to equinox.
+.le
+.ls fk5 [equinox] [epoch]
+The equatorial mean place post-IAU 1976 (FK5) system where equinox is
+a Julian or Besselian epoch e.g. J2000.0  or B1980.0.
+Equinoxes without the J / j or B / b prefix are treated as Julian epochs.
+The default value of equinox is J2000.0.
+Epoch is a Besselian epoch, a Julian epoch, or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Julian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date.  If undefined epoch defaults to equinox.
+.le
+.ls fk4 [equinox] [epoch]
+The equatorial mean place pre-IAU 1976 (FK4) system where equinox is a
+Besselian or Julian epoch e.g. B1950.0  or J2000.0,
+and epoch is the Besselian epoch, the Julian epoch, or the Julian date of the
+observation.
+Equinoxes without the J / j or B / b prefix are treated
+as Besselian epochs. The default value of equinox is B1950.0. Epoch
+is a Besselian epoch, a Julian epoch, or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Besselian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date.  If undefined epoch defaults to equinox.
+.le
+.ls noefk4 [equinox] [epoch]
+The equatorial mean place pre-IAU 1976 (FK4) system but without the E-terms
+where equinox is a Besselian or Julian epoch e.g. B1950.0 or J2000.0,
+and epoch is the Besselian epoch, the Julian epoch, or the Julian date of the
+observation.
+Equinoxes without the J / j or B / b prefix are treated
+as Besselian epochs. The default value of equinox is B1950.0.
+Epoch is a Besselian epoch, a Julian epoch, or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Besselian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian day.  If undefined epoch defaults to equinox.
+.le
+.ls apparent epoch
+The equatorial geocentric apparent place post-IAU 1976 system where
+epoch is the epoch of observation.
+Epoch is a Besselian epoch, a Julian epoch or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Besselian
+epochs if the epoch value < 1984.0, Julian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date.
+.le
+.ls ecliptic epoch
+The ecliptic coordinate system where epoch is the epoch of observation.
+Epoch is a Besselian epoch, a Julian epoch, or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Besselian epochs
+if the epoch values < 1984.0, Julian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian day.
+.le
+.ls galactic [epoch]
+The IAU 1958 galactic coordinate system.
+Epoch is a Besselian epoch, a Julian epoch or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Besselian
+epochs if the epoch value < 1984.0, Julian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date. The default value of epoch is B1950.0.
+.le
+.ls supergalactic [epoch]
+The deVaucouleurs supergalactic coordinate system.
+Epoch is a Besselian epoch, a Julian epoch or a Julian date.
+Julian epochs are prefixed by a J or j, Besselian epochs by a B or b.
+Epochs without the J / j or B / b prefix default to Besselian
+epochs if the epoch value < 1984.0, Julian epochs
+if the epoch value <= 3000.0, otherwise epoch is interpreted as
+a Julian date. The default value of epoch is B1950.0.
+.le
+
+Fields enclosed in [] are optional with the defaults as described. The epoch
+field for the "icrs" , "fk5", "galactic", and "supergalactic" coordinate
+systems is only used if the input coordinates are in the equatorial fk4,
+noefk4, fk5, or icrs systems and proper motions are used to transform from
+coordinate system to another.
+
+.ih
+SEE ALSO
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skclose.hlp b/pkg/xtools/skywcs/doc/skclose.hlp
new file mode 100644
index 00000000..191b08b5
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skclose.hlp
@@ -0,0 +1,23 @@
+.help skclose Mar00 Skywcs
+.ih
+NAME
+skclose -- free the sky coordinate descriptor
+.ih
+SYNOPSIS
+call sk_close (coo)
+
+.nf
+pointer	coo		# the sky coordinate descriptor
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor to be freed.
+.le
+.ih
+DESCRIPTION
+Sk_close frees a previously allocated sky coordinate descriptor.
+.ih
+SEE ALSO
+skdecwcs, skdecwstr, skdecim, skcopy
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skcopy.hlp b/pkg/xtools/skywcs/doc/skcopy.hlp
new file mode 100644
index 00000000..68219c0d
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skcopy.hlp
@@ -0,0 +1,24 @@
+.help skcopy Mar00 Skywcs
+.ih
+NAME
+skcopy -- copy a sky coordinate descriptor
+.ih
+SYNOPSIS
+newcoo = sk_copy (coo)
+
+.nf
+pointer	coo		# the sky coordinate descriptor
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor to be copied.
+.le
+.ih
+DESCRIPTION
+Sk_copy is a pointer function which returns a copy of the input sky coordinate
+descriptor as its function value.
+.ih
+SEE ALSO
+skdecwcs, skdecwstr, skdecim, skclose
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skdecim.hlp b/pkg/xtools/skywcs/doc/skdecim.hlp
new file mode 100644
index 00000000..6e570e47
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skdecim.hlp
@@ -0,0 +1,56 @@
+.help skdecim Mar00 Skywcs
+.ih
+NAME
+skdecim -- open a sky coordinate descriptor using an image descriptor
+.ih
+SYNOPSIS
+stat = sk_decim (im, mw, coo, imcoo)
+
+.nf
+pointer	im		# the input image descriptor
+pointer	mw		# the output mwcs descriptor
+pointer	coo		# the output sky coordinate descriptor
+pointer	imcoo		# the input image sky coordinate descriptor 
+.fi
+.ih
+ARGUMENTS
+.ls  im
+The input image descriptor.
+.le
+.ls mw
+The output mwcs descriptor.  A NULL value for mw is returned if the image
+world coordinate system cannot be read.
+.le
+.ls coo
+The output sky coordinate descriptor.
+.le
+.ls imcoo
+The parent image sky coordinate descriptor. Imcoo is set to NULL
+except in cases where the sky coordinate descriptor for an image is
+transformed and written back to the same image.
+.le
+.ih
+DESCRIPTION
+Sk_decim is an integer function which returns OK or ERR as its function
+value. ERR is returned if a valid sky coordinate system cannot be opened,
+OK otherwise.
+
+Sk_decim returns the image MWCS descriptor mw. The MWCS descriptor is used
+to convert from pixel coordinates to world coordinates and vice versa.
+The MWCS descriptor must be freed with a call to the MWCS routine
+mw_close before task termination.
+
+Sk_decim returns the sky descriptor coo. The sky coordinate descriptor 
+is defined even if an error is detected in reading the image celestial 
+coordinate system, and must be freed with a call to sk_close before
+task termination.
+
+.ih
+NOTES
+Type "help ccsystems" to see the list of the supported sky coordinate systems.
+
+Type "help mwcs$MWCS.hlp fi+" to find out more about the IRAF image world
+coordinate system library MWCS.
+SEE ALSO
+skdecwcs, skdecwstr, skcopy, skclose
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skdecwcs.hlp b/pkg/xtools/skywcs/doc/skdecwcs.hlp
new file mode 100644
index 00000000..2081fd50
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skdecwcs.hlp
@@ -0,0 +1,62 @@
+.help skdecwcs Mar00 Skywcs 
+.ih
+NAME
+skdecwcs -- open a sky coordinate descriptor using an image or system name
+.ih
+SYNOPSIS
+stat = sk_decwcs (ccsystem, mw, coo, imcoo)
+
+.nf
+char	ccsystem	# the input celestial coordinate system name
+pointer	mw		# the output mwcs descriptor
+pointer	coo		# the output sky coordinate descriptor
+pointer	imcoo		# the input image sky coordinate descriptor 
+.fi
+.ih
+ARGUMENTS
+.ls  ccsystem.   
+The celestial coordinate system name. Ccsystem is a either an image system
+name, e.g.  "dev$ypix logical" or "dev$ypix world" or a system name, e.g.
+"J2000" or "galactic".
+.le
+.ls mw
+The output mwcs descriptor. A NULL value for mw is returned if the
+image world coordinate system cannot be read or ccsystem is not an image
+system name. 
+.le
+.ls coo
+The output sky coordinate descriptor.
+.le
+.ls imcoo
+The parent image coordinate descriptor. Imcoo is set to NULL
+except in cases where the sky coordinate descriptor for an image is
+transformed and written back to the same image.
+.le
+.ih
+DESCRIPTION
+Sk_decwcs is an integer function which returns OK or ERR as its function
+value. ERR is returned if a valid sky coordinate system cannot be opened,
+OK otherwise.
+
+Sk_decwcs returns the image MWCS descriptor mw if ccsystem is an image
+system, otherwise it returns NULL. The MWCS descriptor is used
+to convert from pixel coordinates to world coordinates and vice versa.
+The MWCS descriptor must be freed with a call to the MWCS routine
+mw_close before task termination.
+
+Sk_decwcs returns the sky descriptor coo. The sky coordinate descriptor
+is defined even if an error is detected in reading the image celestial
+coordinate system, and must be freed with a call to sk_close before
+task termination.
+
+.ih
+NOTES
+Type "help ccsystems" to see the list of the supported sky coordinate systems.
+
+Type "help mwcs$MWCS.hlp fi+" to find out more about the IRAF image world
+coordinate system library MWCS.
+
+
+SEE ALSO
+skdecwstr, skdecim
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skdecwstr.hlp b/pkg/xtools/skywcs/doc/skdecwstr.hlp
new file mode 100644
index 00000000..0edf7fa0
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skdecwstr.hlp
@@ -0,0 +1,46 @@
+.help skdecwstr Mar00 Skywcs
+.ih
+NAME
+skdecwstr -- open a sky coordinate descriptor using a system name
+.ih
+SYNOPSIS
+stat = sk_decwstr (csystem, coo, imcoo)
+
+.nf
+char	csystem		# the input celestial coordinate system name
+pointer	coo		# the output sky coordinate descriptor
+pointer	imcoo		# the input image sky coordinate descriptor 
+.fi
+.ih
+ARGUMENTS
+.ls  csystem 
+The sky coordinates definition. Ccsystem is a system name,  e.g. "J2000"
+or "galactic".
+.le
+.ls coo
+The output sky coordinate descriptor.
+.le
+.ls imcoo
+The parent image coordinate descriptor. Imcoo is set to NULL
+except in cases where the sky coordinate descriptor for an image is
+transformed and written back to the same image.
+.le
+.ih
+DESCRIPTION
+Sk_decwstr is an integer function which returns OK or ERR as its function
+value. ERR is returned if a valid sky coordinate system cannot be opened,
+OK otherwise.
+
+Sk_decwstr returns the sky descriptor coo. The sky coordinate descriptor
+is defined even if an error is detected in reading the image celestial
+coordinate system, and must be freed with a call to sk_close before
+task termination.
+
+.ih
+NOTES
+
+Type "help ccsystems" to get a list of the supported sky coordinate systems.
+
+SEE ALSO
+skdecwcs, skdecim, skcopy, skclose
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skenwcs.hlp b/pkg/xtools/skywcs/doc/skenwcs.hlp
new file mode 100644
index 00000000..cc388108
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skenwcs.hlp
@@ -0,0 +1,32 @@
+.help skenwcs Mar00 Skywcs
+.ih
+NAME
+skenwcs -- encode a system name using a sky coordinate descriptor
+.ih
+SYNOPSIS
+
+call sk_enwcs (coo, csystem, maxch)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+char	csystem		# the output system name
+int	maxch		# the maximum size of the output system name
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The input sky coordinate descriptor
+.le
+.ls csystem
+The output system name, e.g. "galactic".
+.le
+.ls maxch
+The maximum size of the output system name.
+.le
+.ih
+DESCRIPTION
+Sk_enwcs returns the sky coordinate system name.
+.ih
+SEE ALSO
+skdecwcs, skdecwstr
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skequatorial.hlp b/pkg/xtools/skywcs/doc/skequatorial.hlp
new file mode 100644
index 00000000..4500b881
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skequatorial.hlp
@@ -0,0 +1,58 @@
+.help skequatorial Mar00 Skywcs
+.ih
+NAME
+skequatorial -- apply pm and transform between equatorial coordinate systems
+.ih
+SYNOPSIS
+call sk_equatorial (incoo, outcoo, ilng, ilat, ipmlng, ipmlat, px, rv,
+	            olng, olat)
+
+.nf
+pointer	incoo		# the input sky coordinate descriptor
+pointer	outcoo		# the output sky coordinate descriptor
+double ilng, ilat	# the input sky coordinates in radians
+double ipmlng, ipmlat	# the input proper motions in radians  / year
+double	px		# the input parallax in arcsec
+double	rv		# the input radial velocity in km / sec (+ve receding)
+double olng, olat	# the output sky coordinates in radians
+.fi
+.ih
+ARGUMENTS
+.ls  incoo    
+The input sky coordinate descriptor.
+.le
+.ls  outcoo    
+The output sky coordinate descriptor.
+.le
+.ls ilng, ilat
+The input sky coordinates in radians.
+.le
+.ls ipmlng, ipmlat
+The input proper motions. If proper motions are unknown do not set ipmlng
+and ipmlat to 0.0, use sk_ultran instead. Note that the ra proper motion
+is in dra  not cos (dec) * dra units.
+.le
+.ls px
+The parallax in arcseconds. Use 0.0 if the proper motion is unknown unknown. 
+The parallax value is used only if proper motions are defined.
+.le
+.ls rv
+The radial velocity in km / sec. Use 0.0 if the radial velocity is unknown.
+The radial velocity value is used only if proper motions are defined.
+.le
+.ls olng, olat
+The output sky coordinates in radians.
+.le
+.ih
+DESCRIPTION
+The coordinates in the input sky coordinate system are converted to 
+coordinates in the output sky coordinate system.
+.ih
+NOTES
+If the proper motions are undefined use the routine sk_ultran. Zero valued
+proper motions are not the same as undefined proper motions.
+
+.ih
+SEE ALSO
+sk_lltran, sk_ultran
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skiiprint.hlp b/pkg/xtools/skywcs/doc/skiiprint.hlp
new file mode 100644
index 00000000..217819c2
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skiiprint.hlp
@@ -0,0 +1,39 @@
+.help skiiprint Mar00 Skywcs
+.ih
+NAME
+skiiprint -- print the sky coordinate system summary
+.ih
+SYNOPSIS
+
+call sk_iprint (label, imagesys, mw, coo)
+
+.nf
+char	label		# the input user label
+char	imagesys	# the input image system
+pointer	mw		# the input mwcs descriptor 
+pointer	coo		# the sky coordinate descriptor
+.fi
+.ih
+ARGUMENTS
+.ls label
+The input user supplied label, e.g. "Input System", "Ref System",
+"Output System" etc.
+.le
+.ls imagesys
+The input image system, e.g. "dev$ypix logical", "dev$ypix world", etc.
+.le
+.ls mwcs
+The input image mwcs descriptor if defined. If mwcs is defined then 
+information about which sky coordinate  corresponds to which image
+axis etc is read from the mwcs descriptor.
+.le
+.ls  coo    
+The input sky coordinate descriptor.
+.le
+.ih
+DESCRIPTION
+A summary of the sky coordinate system is printed on the standard output.
+.ih
+SEE ALSO
+skiiwrite
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skiiwrite.hlp b/pkg/xtools/skywcs/doc/skiiwrite.hlp
new file mode 100644
index 00000000..c82472f4
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skiiwrite.hlp
@@ -0,0 +1,43 @@
+.help skiiwrite Mar00 Skywcs
+.ih
+NAME
+skiiwrite -- write the sky coordinate system summary to a file
+.ih
+SYNOPSIS
+
+call sk_iiwrite (outfd, label, imagesys, mw, coo)
+
+.nf
+int	outfd		# the input file descriptor
+char	label		# the input user label
+char	imagesys	# the input image system
+pointer	mw		# the input mwcs descriptor
+pointer	coo		# the sky coordinate descriptor
+.fi
+.ih
+ARGUMENTS
+.ls outfd
+The input file descriptor.
+.le
+.ls label
+The input user supplied label, e.g. "Input System", "Ref System",
+"Output System" etc.
+.le
+.ls imagesys
+The input image system, e.g. "dev$ypix logical", "dev$ypix world", etc.
+.le
+.ls mwcs
+The input image mwcs descriptor if defined. If mwcs is defined then
+information about which sky coordinate  corresponds to which image
+axis etc is read from the mwcs descriptor.
+.le
+.ls  coo    
+The input sky coordinate descriptor.
+.le
+.ih
+DESCRIPTION
+A summary of the sky coordinate system is written to a file.
+.ih
+SEE ALSO
+skiiprint
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/sklltran.hlp b/pkg/xtools/skywcs/doc/sklltran.hlp
new file mode 100644
index 00000000..b45f3ea4
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/sklltran.hlp
@@ -0,0 +1,59 @@
+.help sklltran Mar00 Skywcs
+.ih
+NAME
+sklltran -- apply pm and transform between coordinate systems
+.ih
+SYNOPSIS
+call sk_lltran (incoo, outcoo, ilng, ilat, ipmlng, ipmlat, px, rv, olng, olat)
+
+.nf
+pointer	incoo		# the input sky coordinate descriptor
+pointer	outcoo		# the output sky coordinate descriptor
+double ilng, ilat	# the input sky coordinates in radians
+double ipmlng, ipmlat	# the input proper motions in radians  / year
+double	px		# the input parallax in arcsec
+double	rv		# the input radial velocity in km / sec (+ve receding)
+double olng, olat	# the output sky coordinates in radians
+.fi
+.ih
+ARGUMENTS
+.ls  incoo    
+The input sky coordinate descriptor.
+.le
+.ls  outcoo    
+The output sky coordinate descriptor.
+.le
+.ls ilng, ilat
+The input sky coordinates in radians.
+.le
+.ls ipmlng, ipmlat
+The input proper motions. For these to be applied the input coordinate
+system must be an equatorial coordinate system. If proper motions are
+unknown do not set ipmlng and ipmlat to 0.0, use sk_ultran instead. Note that
+the ra proper motion is in dra  not cos (dec) * dra units.
+.le
+.ls px
+The parallax in arcseconds. Use 0.0 if the proper motion is unknown unknown.
+The parallax value is used only if proper motions are defined.
+.le
+.ls rv
+The radial velocity in km / sec. Use 0.0 if the radial velocity is unknown.
+The radial velocity value is used only if proper motions are defined.
+.le
+.ls olng, olat
+The onput sky coordinates in radians.
+.le
+
+.ih
+DESCRIPTION
+The coordinates in the input sky coordinate system are converted to
+coordinates in the output sky coordinate system.
+.ih
+NOTES
+If the proper motions are undefined use the routine sk_ultran. Zero valued
+proper motions are not the same as undefined proper motions.
+
+.ih
+SEE ALSO
+sk_ultran, sk_equatorial
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/sksaveim.hlp b/pkg/xtools/skywcs/doc/sksaveim.hlp
new file mode 100644
index 00000000..82c16f3f
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/sksaveim.hlp
@@ -0,0 +1,39 @@
+.help sksaveim Mar00 Skywcs
+.ih
+NAME
+sksaveim -- update the image header using a sky coordinate descriptor
+.ih
+SYNOPSIS
+call sk_saveim (coo, mw, im)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+pointer	mw		# the input mwcs descriptor
+pointer	im		# the input image descriptor
+.fi
+.ih
+ARGUMENTS
+.ls coo
+The input sky coordinate descriptor.
+.le
+.ls mw
+The IRAF mwcs descriptor.
+.le
+.ls  im
+The input image descriptor.
+.le
+.ih
+DESCRIPTION
+The image world coordinate system is updated using information in
+the sky coordinate descriptor and the mwcs descriptor.
+
+.ih
+NOTES
+Note that the sk_saveim call does not include a call to the MWCS mw_saveim
+routine. This call must be made separately.
+
+Type "help mwcs$MWCS.hlp fi+" to find out more about the IRAF image world
+coordinate system code.
+SEE ALSO
+skdecwcs, skdecim
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/sksetd.hlp b/pkg/xtools/skywcs/doc/sksetd.hlp
new file mode 100644
index 00000000..f518d71c
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/sksetd.hlp
@@ -0,0 +1,53 @@
+.help sksetd Mar00 Skywcs
+.ih
+NAME
+sksetd -- set a double sky coordinate descriptor parameter
+.ih
+SYNOPSIS
+include <skywcs.h>
+
+call sk_setd (coo, parameter, dval)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+int	parameter	# the double parameter to be set
+double	dval		# the value of the parameter to be set
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor.
+.le
+.ls parameter
+The parameter to be set. The double parameter definitions in skywcs.h are:
+.nf
+	S_VXOFF		# the logical ra / longitude offset in pixels
+	S_VYOFF		# the logical dec / latitude offset in pixels
+	S_VXSTEP	# the logical ra / longitude step size in pixels
+	S_VYSTEP	# the logical dec / latitude step size in pixels
+	S_EQUINOX	# the equinox in years
+	S_EPOCH		# the MJD of the observation 
+.fi
+.le
+.ls dval
+The value of the parameter to be set.
+.le
+.ih
+DESCRIPTION
+Sk_setd sets the values of double sky coordinate descriptor parameters.
+.ih
+NOTES
+The offsets and step sizes default to 0 and 1 for both axes. However
+if the sky coordinate descriptor was derived from an input image section, e.g.
+"dev$ypix[100:300,100:300]" these numbers may assume other values in some
+circumstances.
+
+The equinox and epoch of observation are normally set by the calling program
+when the sky coordinate descriptor is initialized, e.g. they default
+to 2000.0 and  51544.50000 if the input coordinate system was "fk5".
+
+In most cases these parameters should not be set by the user.
+.ih
+SEE ALSO
+skseti, sksets
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skseti.hlp b/pkg/xtools/skywcs/doc/skseti.hlp
new file mode 100644
index 00000000..b08be476
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skseti.hlp
@@ -0,0 +1,93 @@
+.help skseti Mar00 Skywcs
+.ih
+NAME
+skseti -- set an integer sky coordinate descriptor parameter
+.ih
+SYNOPSIS
+include <skywcs.h>
+
+call sk_seti (coo, parameter, ival)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+int	parameter	# the integer parameter to be set
+int	ival		# the value of the parameter to be set
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor.
+.le
+.ls parameter
+The parameter to be set. The double parameter definitions in skywcs.h are:
+.nf
+	S_CTYPE		# the celestial coordinate system type
+	S_RADECSYS	# the equatorial system type
+	S_NLNGUNITS	# the ra / longitude units
+	S_NLATUNITS	# the dec/ latitude units
+	S_WTYPE		# the projection type
+	S_PLNGAX	# the physical ra / longitude axis
+	S_PLATAX	# the physical dec / latitude axis
+	S_XLAX		# the logical ra / longitude axis
+	S_YLAX		# the logical dec / latitude axis
+	S_PIXTYPE	# the IRAF pixel coordinate system type
+	S_NLNGAX	# the length of ra / longitude axis
+	S_NLATAX	# the length of dec / latitude axis
+	S_STATUS	# the coordinate system status 
+.fi
+.le
+.ls ival
+The value of the parameter to be set.
+.le
+.ih
+DESCRIPTION
+Sk_seti sets the values of integer sky coordinate descriptor parameters.
+.ih
+NOTES
+Permitted values of S_CTYPE are CTYPE_EQUATORIAL, CTYPE_ECLIPTIC, 
+CTYPE_GALACTIC, and CTYPE_SUPERGALACTIC. The corresponding string dictionary
+is CTYPE_LIST.
+
+Permitted types of S_RADECSYS are EQTYPE_FK4, EQTYPE_FK4NOE,
+EQTYPE_FK5, EQTYPE, ICRS, and EQTYPE_GAPPT. The corresponding string
+dictionary is EQTYPE_LIST.
+
+Permitted values of S_WTYPE are WTYPE_LIN, WTYPE_AZP, WTYPE_TAN,  WTYPE_SIN,
+WTYPE_STG, WTYPE_ARC, WTYPE_ZPN, WTYPE_ZEA, WTYPE_AIR, WTYPE_CYP, WTYPE_CAR,
+WTYPE_MER, WTYPE_CEA, WTYPE_COP, WTYPE_COD, WTYPE_COE, WTYPE_COO, WTYPE_BON,
+WTYPE_PCO, WTYPE_GLS, WTYPE_PAR, WTYPE_AIT, WTYPE_MOL, WTYPE_CSC, WTYPE_QSC,
+WTYPE_TSC, WTYPE_TNX, WTYPE_ZPX. The corresponding string dictionary is
+WTYPE_LIST.
+
+Permitted values of S_PIXTYPE are  PIXTYPE_LOGICAL, PIXTYPE_TV,
+PIXTYPE_PHYSICAL.  and PIXTPE_WORLD. The corresponding string dictionary
+is PIXTYPE_LIST.
+
+Permitted values of S_NLNGUNITS are SKY_HOURS, SKY_DEGREES, and SKY_RADIANS.
+The corresponding string dictionary is SKY_LNG_UNITLIST.
+Permitted values of S_NLATUNITS are SKY_DEGREES, and SKY_RADIANS.
+The corresponding string dictionary is SKY_LAT_UNITLIST.
+
+The parameters S_CTYPE, S_RADECSYS, S_NLNGUNITS, and S_NLATUNITS are
+important for all sky coordinate descriptors regardless of the source.
+The parameters S_WTYPE, S_PLNGAX, S_PLATAX, S_XLAX, S_YLAX, S_PIXTYPE,
+S_NLNGAX, and S_NLATAX are only important for sky coordinate descriptors
+derived from an image sky coordinate systems. S_STATUS is OK if the sky
+coordinate descriptor describes a valid celestial coordinate system, ERR
+otherwise.
+
+In most cases these parameters should not be modified by the user. The
+major exceptions are the units parameters S_NLNGUNITS and N_LATUNITS
+which assumes default values fo hours and degrees for equatorial sky
+coordinate systems and degrees and degrees for other sky coordinate systems.
+If the user input and output units are different from the normal defaults
+then the units parameters should be set appropriately.
+
+Parameters that occasionally need to be reset when a coordinate system
+is created, edited, or saved to an image are S_WTYPE, S_PIXTYPE, S_PLNGAX,
+and S_PLATAX.
+
+.ih
+SEE ALSO
+sksetd, sksets
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/sksets.hlp b/pkg/xtools/skywcs/doc/sksets.hlp
new file mode 100644
index 00000000..8e4179b4
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/sksets.hlp
@@ -0,0 +1,36 @@
+.help sksets Mar00 Skywcs
+.ih
+NAME
+sksets -- set a string sky coordinate descriptor parameter
+.ih
+SYNOPSIS
+include <skywcs.h>
+
+call sk_sets (coo, parameter, str)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+int	parameter	# the string parameter to be set
+char	str		# the value of the string parameter to be set
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor.
+.le
+.ls parameter
+The parameter to be set. The string parameter definitions in skywcs.h are:
+.nf
+	S_COOSYSTEM	# the celestial coordinate system name
+.fi
+.le
+.ls str
+The value of the parameter to be set.
+.le
+.ih
+DESCRIPTION
+Sk_sets sets the values of string sky coordinate descriptor parameters. 
+.ih
+SEE ALSO
+sksetd, skseti
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skstatd.hlp b/pkg/xtools/skywcs/doc/skstatd.hlp
new file mode 100644
index 00000000..52dc0c70
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skstatd.hlp
@@ -0,0 +1,49 @@
+.help skstatd Mar00 Skywcs
+.ih
+NAME
+skstatd -- get a double sky coordinate descriptor parameter
+.ih
+SYNOPSIS
+include <skywcs.h>
+
+dval = sk_statd (coo, parameter)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+int	parameter	# the double parameter to be returned
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor.
+.le
+.ls parameter
+The oarameter to be returned. The double parameter definitions in skywcs.h are:
+.nf
+	S_VXOFF		# the logical ra / longitude offset in pixels
+	S_VYOFF		# the logical dec / latitude offset in pixels
+	S_VXSTEP	# the logical ra / longitude step size in pixels
+	S_VYSTEP	# the logical dec / latitude step size in pixels
+	S_EQUINOX	# the equinox in years
+	S_EPOCH		# the MJD of the observation
+.fi
+.le
+.ih
+DESCRIPTION
+Sk_statd returns the values of double sky coordinate descriptor parameters.
+
+.ih
+NOTES
+The offsets and step sizes default to 0 and 1 for both axes. However
+if the sky coordinate descriptor was derived from an input image section, e.g.
+"dev$ypix[100:300,100:300]" these numbers may assume other values in some
+circumstances.
+
+The equinox and epoch of observation are normally set by the calling program
+when the sky coordinate descriptor is initialized, e.g. they default
+to 2000.0 and  51544.50000 if the input coordinate system was "fk5".
+
+.ih
+SEE ALSO
+skstati, skstats
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skstati.hlp b/pkg/xtools/skywcs/doc/skstati.hlp
new file mode 100644
index 00000000..90d33eb1
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skstati.hlp
@@ -0,0 +1,79 @@
+.help skstati Mar00 Skywcs
+.ih
+NAME
+skstati -- get an integer sky coordinate descriptor parameter
+.ih
+SYNOPSIS
+include <skywcs.h>
+
+ival = sk_stati (coo, parameter)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+int	parameter	# the integer parameter to be returned
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor.
+.le
+.ls parameter
+Parameter to be returned. The integer parameter definitions in skywcs.h are:
+.nf
+	S_CTYPE		# the celestial coordinate system type
+	S_RADECSYS	# the equatorial system type
+	S_NLNGUNITS	# the ra / longitude units
+	S_NLATUNITS	# the dec/ latitude units
+	S_WTYPE		# the projection type
+	S_PLNGAX	# the physical ra / longitude axis
+	S_PLATAX	# the physical dec / latitude axis
+	S_XLAX		# the logical ra / longitude axis
+	S_YLAX		# the logical dec / latitude axis
+	S_PIXTYPE	# the IRAF pixel coordinate system type
+	S_NLNGAX	# the length of the ra / longitude axis
+	S_NLATAX	# the length of the dec / latitude axis
+	S_STATUS	# the coordinate system status 
+.fi
+.le
+.ih
+DESCRIPTION
+Sk_stati returns the values of integer sky coordinate descriptor parameters.
+
+.ih
+NOTES
+Permitted values of S_CTYPE are CTYPE_EQUATORIAL, CTYPE_ECLIPTIC, 
+CTYPE_GALACTIC, and CTYPE_SUPERGALACTIC. The corresponding string dictionary
+is CTYPE_LIST.
+
+Permitted types of S_RADECSYS are EQTYPE_FK4, EQTYPE_FK4NOE,
+EQTYPE_FK5, EQTYPE, ICRS, and EQTYPE_GAPPT. The corresponding string
+dictionary is EQTYPE_LIST.
+
+Permitted values of S_WTYPE are WTYPE_LIN, WTYPE_AZP, WTYPE_TAN,  WTYPE_SIN,
+WTYPE_STG, WTYPE_ARC, WTYPE_ZPN, WTYPE_ZEA, WTYPE_AIR, WTYPE_CYP, WTYPE_CAR,
+WTYPE_MER, WTYPE_CEA, WTYPE_COP, WTYPE_COD, WTYPE_COE, WTYPE_COO, WTYPE_BON,
+WTYPE_PCO, WTYPE_GLS, WTYPE_PAR, WTYPE_AIT, WTYPE_MOL, WTYPE_CSC, WTYPE_QSC,
+WTYPE_TSC, WTYPE_TNX, WTYPE_ZPX. The corresponding string dictionary is
+WTYPE_LIST.
+
+Permitted values of S_PIXTYPE are  PIXTYPE_LOGICAL, PIXTYPE_TV,
+PIXTYPE_PHYSICAL.  and PIXTPE_WORLD. The corresponding string dictionary
+is PIXTYPE_LIST.
+
+Permitted values of S_NLNGUNITS are SKY_HOURS, SKY_DEGREES, and SKY_RADIANS.
+The corresponding string dictionary is SKY_LNG_UNITLIST.
+Permitted values of S_NLATUNITS are SKY_DEGREES, and SKY_RADIANS.
+The corresponding string dictionary is SKY_LAT_UNITLIST.
+
+The parameters S_CTYPE, S_RADECSYS, S_NLNGUNITS, and S_NLATUNITS are
+important for all sky coordinate descriptors regardless of the source.
+The parameters S_WTYPE, S_PLNGAX, S_PLATAX, S_XLAX, S_YLAX, S_PIXTYPE,
+S_NLNGAX, and S_NLATAX are only important for sky coordinate descriptors
+derived from an image sky coordinate systems. S_STATUS is OK if the sky
+coordinate descriptor describes a valid celestial coordinate system, ERR
+otherwise.
+
+.ih
+SEE ALSO
+skstatd, skstats
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skstats.hlp b/pkg/xtools/skywcs/doc/skstats.hlp
new file mode 100644
index 00000000..483ed3e5
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skstats.hlp
@@ -0,0 +1,40 @@
+.help skstats Mar00 Skywcs
+.ih
+NAME
+skstats -- get a string sky coordinate descriptor parameter
+.ih
+SYNOPSIS
+include <skywcs.h>
+
+call sk_stats (coo, parameter, str, maxch)
+
+.nf
+pointer	coo		# the input sky coordinate descriptor
+int	parameter	# the string parameter to be returned
+char	str		# the returned string parameter value
+int	maxch		# the maximum size of the returned string parameter
+.fi
+.ih
+ARGUMENTS
+.ls  coo    
+The sky coordinate descriptor.
+.le
+.ls parameter
+The parameter to be returned. The string parameter definitions in skywcs.h are:
+.nf
+	S_COOSYSTEM	# the celestial coordinate system name
+.fi
+.le
+.ls str
+The value of the returned string.
+.le
+.ls maxch
+The maximum size of the returned string.
+.le
+.ih
+DESCRIPTION
+Sk_stats returns the values of string sky coordinate descriptor parameters.
+.ih
+SEE ALSO
+skstati, skstatd
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skultran.hlp b/pkg/xtools/skywcs/doc/skultran.hlp
new file mode 100644
index 00000000..ca02385e
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skultran.hlp
@@ -0,0 +1,50 @@
+.help skultran Mar00 Skywcs
+.ih
+NAME
+skultran -- transform between coordinate systems
+.ih
+SYNOPSIS
+call sk_ultran (incoo, outcoo, ilng, ilat, olng, olat, npts)
+
+.nf
+pointer	incoo		# the input sky coordinate descriptor
+pointer	outcoo		# the output sky coordinate descriptor
+double ilng, ilat	# the input celestial coordinates in expected units
+double olng, olat	# the output celestial coordinates in expected units
+int	npts		# the number of input and output coordinate pairs
+.fi
+.ih
+ARGUMENTS
+.ls  incoo    
+The input sky coordinate descriptor.
+.le
+.ls  outcoo    
+The output sky coordinate descriptor.
+.le
+.ls ilng, ilat
+The input sky coordinates in the units defined by the integer parameters
+S_NLNGUNITS and S_NLATUNITS.
+.le
+.ls olng, olat
+The output sky coordinates in the units defined by the integer parameters
+S_NLNGUNITS and S_NLATUNITS.
+.le
+.ls npts
+The number of input and output coordinate pairs.
+.le
+.ih
+DESCRIPTION
+The coordinates in the input coordinate system are converted to 
+coordinates in the output coordinates system.
+
+If the calling program has not set the S_NLNGUNITS and S_NLATUNITS parameters
+in either system the expected coordinates are hours and degrees for
+equatorial sky coordinate systems and degrees and degrees for other sky
+coordinate systems. The calling program must either perform the necessary
+coordinate conversions or set the units parameters in the input and output
+sky coordinate descriptors appropriately.
+
+.ih
+SEE ALSO
+sk_lltran, sk_equatorial
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skywcs.hd b/pkg/xtools/skywcs/doc/skywcs.hd
new file mode 100644
index 00000000..74bac140
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skywcs.hd
@@ -0,0 +1,25 @@
+# Help directory for the SKYWCS library
+
+$doc		= "./"
+$source		= "../"
+
+skdecwcs	hlp=doc$skdecwcs.hlp,		src=source$skdecode.x
+skdecwstr	hlp=doc$skdecwstr.hlp,		src=source$skdecode.x
+skdecim		hlp=doc$skdecim.hlp,		src=source$skdecode.x
+skenwcs		hlp=doc$skenwcs.hlp,		src=source$skdecode.x
+skcopy		hlp=doc$skcopy.hlp,		src=source$skdecode.x
+skiiprint	hlp=doc$skiiprint.hlp,		src=source$skwrite.x
+skiiwrite	hlp=doc$skiiwrite.hlp,		src=source$skwrite.x
+skstati		hlp=doc$skstati.hlp,		src=source$skstat.x
+skstatd		hlp=doc$skstatd.hlp,		src=source$skstat.x
+skstats		hlp=doc$skstats.hlp,		src=source$skstat.x
+skseti		hlp=doc$skseti.hlp,		src=source$skset.x
+sksetd		hlp=doc$sksetd.hlp,		src=source$skset.x
+sksets		hlp=doc$sksets.hlp,		src=source$skset.x
+skultran	hlp=doc$skultran.hlp,		src=source$skytransform.x
+sklltran	hlp=doc$sklltran.hlp,		src=source$skytransform.x
+skequatorial	hlp=doc$skequatorial.hlp,	src=source$skytransform.x
+sksaveim	hlp=doc$sksaveim.hlp,		src=source$sksaveim.x
+skclose		hlp=doc$skclose.hlp,		src=source$skdecode.x
+
+ccsystems	hlp=doc$ccsystems.hlp
diff --git a/pkg/xtools/skywcs/doc/skywcs.hlp b/pkg/xtools/skywcs/doc/skywcs.hlp
new file mode 100644
index 00000000..d02f4d2f
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skywcs.hlp
@@ -0,0 +1,306 @@
+.help skywcs Oct00 xtools
+.ih
+NAME
+skywcs -- sky coordinates package
+.ih
+SYNOPSIS
+
+.nf
+   stat = sk_decwcs (ccsystem, mw, coo, imcoo)
+   stat = sk_decwstr (ccsystem, coo, imcoo)
+   stat = sk_decim (im, wcs, mw, coo)
+	  sk_enwcs (coo, ccsystem, maxch)
+ newcoo = sk_copy (coo)
+          sk_iiprint (label, imagesys, mw, coo) 
+          sk_iiwrite (fd, label, imagesys, mw, coo) 
+[id]val = sk_stat[id] (coo, param)
+          sk_stats (coo, param, str, maxch)
+	  sk_set[id] (coo, param, [id]val)
+          sk_sets (coo, param, str)
+	  sk_ultran (incoo, outcoo, ilng, ilat, olng, olat, npts)
+	  sk_lltran (incoo, outoo, ilng, ilat, ipmlng, ipmlat, px, rv,
+             olng, olat)
+	  sk_equatorial (incoo, outcoo, ilng, ilat, ipmlng, ipmlat, px,
+	     rv, olng, olat)
+	  sk_saveim (coo, mw, im)
+	  sk_close (coo)
+
+.fi
+.ih
+DESCRIPTION
+
+The skywcs package contains a simple set of routines for managing
+sky coordinate information and for transforming from one sky coordinate
+system to another. The sky coordinate system is defined either by a system
+name, e.g. "J2000", "galactic", etc. or by an image system name, e.g.
+"dev$ypix" or "dev$ypix world".
+
+The skywcs routine are layered on the Starlink Positional Astronomy library
+SLALIB which is installed in the IRAF MATH package.  Type "help slalib
+option=sys" for more information about SLALIB.
+
+
+.ih
+NOTES
+
+An "include <skywcs.h>" statement must be included in the calling program
+to make the skywcs package parameter definitions visible to the calling
+program.
+
+The sky coordinate descriptor is created with a call to one of the sk_decwcs
+sk_decwstr or sk_imwcs routines. If the source of sky coordinate descriptor
+is an image then an IRAF MWCS descriptor will be returned with the sky 
+oordinate descriptor. The sky coordinate descriptor is freed with a
+call to sk_close. A separate call to mw_close must be made to free the
+MWCS descriptor if one was allocated.
+
+By default the main skywcs coordinate transformation routine sk_ultran
+assumes that the input and output sky coordinates are in hours and degrees
+if the input and output coordinate systems are equatorial, otherwise the
+coordinates are assumed to be in degrees and degrees. The default input and
+output sky coordinate units can be reset with calls to sk_seti. Two lower level
+coordinate transformations for handling proper motions sk_lltran and
+sk_equatorial are also available. These routines expect the input and output
+coordinates and proper motions to be in radians.
+
+Calling programs working with both sky coordinate and MWCS descriptors
+need to be aware that the MWCS routines assume that all sky coordinates
+must be input in degrees and will be output in degrees and adjust their
+code accordingly.
+
+The skywcs routine sk_saveim can be used to update an image header.
+
+
+.ih
+EXAMPLES
+.nf
+Example 1: Convert from B1950 coordinates to J2000 coordinates.
+
+    include <skywcs.h>
+
+    ....
+
+    # Open input coordinate system.
+    instat = sk_decwstr ("B1950", incoo, NULL)
+    if (instat == ERR) {
+	call sk_close (incoo)
+	return
+    }
+
+    # Open output coordinate system.
+    outstat = sk_decwstr ("J2000", outcoo, NULL)
+    if (outstat == ERR) {
+	call sk_close (outcoo)
+	return
+    }
+
+    # Do the transformation assuming the input coordinates are in hours
+    # and degrees. The output coordinates will be in hours and degrees
+    # as well.
+    call sk_ultran (incoo, outcoo, rain, decin, raout, decout, npts) 
+
+    # Close the coordinate descriptors.
+    call sk_close (incoo)
+    call sk_close (outcoo)
+
+    ...
+
+
+Example 2: Repeat example 1 but convert to galactic coordinates.
+
+    include <skywcs.h>
+
+    ....
+
+    # Open the input coordinate system.
+    instat = sk_decwstr ("B1950", incoo, NULL)
+    if (instat == ERR) {
+	call sk_close (incoo)
+	return
+    }
+
+    # Open the output coordinate system.
+    outstat = sk_decwstr ("galactic", outcoo, NULL)
+    if (outstat == ERR) {
+	call sk_close (outcoo)
+	return
+    }
+
+    # Do the transformation assuming the input coordinates are in hours and
+    # degrees. The output coordinates will be in degrees and degrees.
+    call sk_ultran (incoo, outcoo, rain, decin, raout, decout, npts) 
+
+    # Close the coordinate descriptors.
+    call sk_close (incoo)
+    call sk_close (outcoo)
+
+    ...
+
+Example 3: Convert a grid of pixel coordinates in the input image to the
+	   equivalent pixel coordinate in the output image using the
+	   image world coordinate systems to connect the two.
+
+    include <skywcs.h>
+
+    ....
+
+    # Mwref will be defined because the input system is an image.
+    refstat = sk_decwcs ("refimage logical", mwref, refcoo, NULL)
+    if (refstat == ERR || mwref == NULL) {
+	if (mwref != NULL)
+	    call mw_close (mwref)
+	call sk_close (refcoo)
+	return
+    }
+
+    # Set the reference coordinate descriptor so it expects input in degrees
+    # and degrees.
+    call sk_seti (refcoo, S_NLNGUNUTS, SKY_DEGREES)
+    call sk_seti (refcoo, S_NLATUNUTS, SKY_DEGREES)
+
+    # Mwout will be defined because the output system is an image.
+    outstat = sk_decwcs ("image logical", mwout, outcoo, NULL)
+    if (outstat == ERR || mwout == NULL) {
+	if (mwout != NULL)
+	    call mw_close (mwout)
+	call sk_close (outcoo)
+	call mw_close (mwref)
+	call sk_close (refcoo)
+	return
+    }
+
+    # Set the output coordinate descriptor so it will output coordinates
+    # in degrees and degrees.
+    call sk_seti (outcoo, S_NLNGUNUTS, SKY_DEGREES)
+    call sk_seti (outcoo, S_NLATUNUTS, SKY_DEGREES)
+
+    # Compute pixel grid in refimage and store coordinate in the arrays
+    # xref and yref.
+    npts = 0
+    do j = 1, IM_LEN(im,2), 100 {
+        do i = 1, IM_LEN(im,1), 100 {
+	    npts = npts + 1
+	    xref[npts] = i
+	    yref[npts] = j
+	}
+    }
+
+    # Convert xref and yref to celestial coordinates raref and decref using
+    # mwref. The output coordinates will be in degrees and degrees.
+    ctref = mw_sctran (mwref, "logical", "world", 03B)
+    do i = 1, npts
+	call mw_c2trand (ctref, xref[i], yref[i], raref[i], decref[i])
+    call ct_free (ctref)
+
+    # Convert the reference celestial coordinates to the output celestial
+    # coordinate system using the coordinate descriptors.
+    call sk_ultran (refcoo, outcoo, raref, decref, raout, decout, npts)
+
+    # Convert the output celestial coordinates to  pixel coordinates in
+    # the other image using mwout.
+    ctout = mw_sctran (mwout, "world", "logical", 03B)
+    do i = 1, npts
+	call mw_c2trand (ctout, raout[i], decout[i], xout[i], yout[i])
+    call ct_free (ctout)
+
+    # Print the input and output pixel coordinates.
+    do i = 1, npts {
+	call printf ("%10.3f %10.3f  %10.3f %10.3f\n")
+	    call pargd (xref[i])
+	    call pargd (yref[i])
+	    call pargd (xout[i])
+	    call pargd (yout[i])
+    }
+
+    # Tidy up.
+    call mw_close (mwref)
+    call mw_close (mwout)
+    call sk_close (refcoo)
+    call sk_close (outcoo)
+
+
+Example 4: Convert a 2D image with an J2000 tangent plane projection
+	   wcs to the equivalent galactic wcs. The transformation
+	   requires a shift in origin and a rotation. Assume that the ra
+	   axis is 1 and the dec axis is 2. The details of how to compute
+	   the rotation are not shown here. See the imcctran task for details.
+
+    include <mwset.h> 
+    include <skywcs.h>
+
+    ...
+
+    # Open image.
+    im = immap (image, READ_WRITE, 0)
+
+    # Open the image coordinate system.
+    instat = sk_decim (im, "logical", mwin, cooin)
+    if (instat == ERR || mwin == NULL) {
+	...
+	call sk_close (cooin)
+	...
+    }
+
+    # Get the dimensions of the mwcs descriptor. This should be 2.
+    ndim = mw_ndim (mwin, MW_NPHYSDIM)
+
+    # Get the default coordinates to degrees and degreees.
+    call sk_seti (cooin, S_NLNGUNITS, SKY_DEGREES)
+    call sk_seti (cooin, S_NATGUNITS, SKY_DEGREES)
+
+    # Open the output coordinate system. Mwout is NULL because this system
+    # is not an image.
+    outstat = sk_decwstr ("galactic", mwout, cooout, cooin)
+    if (outstat == ERR) {
+        ...
+	call sk_close (outstat)
+	...
+    }
+
+    # Make a copy of the mwcs descriptor.
+    mwout = mw_newcopy (mwin)
+
+    # Allocate space for the r and w vectors and cd matrix.
+    call malloc (r, ndim, TY_DOUBLE)
+    call malloc (w, ndim, TY_DOUBLE)
+    call malloc (cd, ndim * ndim, TY_DOUBLE)
+    call malloc (newcd, ndim * ndim, TY_DOUBLE)
+
+    # Assume for simplicty that the MWCS LTERM is the identify transform.
+    # so we don't have to worry about it. Get the WTERM which consists
+    # of r the reference point in pixels, w the reference point in degrees,
+    # and the cd matrix in degrees per pixel.
+    call mw_gwtermd (mwin, Memd[r], Memd[w], Memd[cd], ndim)
+
+    # Convert the world coordinates zero point. The pixel zero point
+    # remains the same.
+    tilng = Memd[w]
+    tilat = Memd[w+1]
+    call sk_ultran (incoo, outcoo, tilng, tilat, tolng, tolat, 1)
+    Memd[w] = tolng
+    Memd[w+1] = tolat
+
+    # Figure out how much to rotate the coordinate system and edit the
+    # compute a new CD matrix. Call it newcd.
+    ...
+
+    # Enter the new CD matrix and zero point.
+    call mw_swterm (mwout, Memd[r], Memd[w], Memd[newcd], ndim)
+
+    # Update the header. 
+    call sk_saveim (cooout, mwout, im)
+    call mw_saveim (mwout, im)
+    ...
+
+    # Tidy up.
+    call mfree (r, TY_DOUBLE)
+    call mfree (w, TY_DOUBLE)
+    call mfree (cd, TY_DOUBLE)
+    call mfree (newcd, TY_DOUBLE)
+    call mw_close (mwin)
+    call mw_close (mwout)
+    call sk_close (cooin)
+    call sk_close (cooout)
+    call imunmap (im)
+.fi
+.endhelp
diff --git a/pkg/xtools/skywcs/doc/skywcs.men b/pkg/xtools/skywcs/doc/skywcs.men
new file mode 100644
index 00000000..7502bcd0
--- /dev/null
+++ b/pkg/xtools/skywcs/doc/skywcs.men
@@ -0,0 +1,15 @@
+    skdecwcs - Open a sky coordinate descriptor using an image or system name
+   skdecwstr - Open a sky coordinate descriptor using a system name
+     skdecim - Open a sky coordinate descriptor using an image descriptor
+     skenwcs - Encode a system name using a sky coordinate descriptor
+      skcopy - Copy a sky coordinate descriptor
+ skstat[ids] - Get a sky coordinate descriptor parameter value
+  skset[ids] - Set a sky coordinate descriptor parameter value
+   skiiprint - Print a sky coordinate descriptor summary
+   skiiwrite - Write a sky coordinate descriptor summary
+    skultran - Transform between coordinate systems
+    sklltran - Apply pm and transform between coordinates systems 
+skequatorial - Apply pm and transform between equatorial coordinate systems
+    sksaveim - Update image header using sky coordinate descriptor
+     skclose - Close the sky coordinate descriptor
+   ccsystems - Describe the supported celestial coordinate systems
diff --git a/pkg/xtools/skywcs/mkpkg b/pkg/xtools/skywcs/mkpkg
new file mode 100644
index 00000000..9a46ce5a
--- /dev/null
+++ b/pkg/xtools/skywcs/mkpkg
@@ -0,0 +1,16 @@
+# Libary for the celestial coordinate sytem pacakge
+
+$checkout libxtools.a lib$
+$update   libxtools.a
+$checkin  libxtools.a lib$
+$exit
+
+libxtools.a:
+	skdecode.x	<imio.h> <imhdr.h> <mwset.h> "skywcsdef.h" "skywcs.h"
+	skwrite.x	"skywcsdef.h" "skywcs.h"
+	skstat.x	"skywcsdef.h" "skywcs.h"
+	skset.x		"skywcsdef.h" "skywcs.h"
+	sktransform.x	<math.h> "skywcsdef.h" "skywcs.h"
+	sksaveim.x	"skywcsdef.h" "skywcs.h"
+	skwrdstr.x
+	;
diff --git a/pkg/xtools/skywcs/skdecode.x b/pkg/xtools/skywcs/skdecode.x
new file mode 100644
index 00000000..5fa88f3b
--- /dev/null
+++ b/pkg/xtools/skywcs/skdecode.x
@@ -0,0 +1,999 @@
+include <imio.h>
+include <imhdr.h>
+include <mwset.h>
+include "skywcs.h"
+include "skywcsdef.h"
+
+# SK_DECWCS -- Decode the wcs string which may be either an image name
+# plus wcs, e.g. "dev$pix logical" or a string describing the celestial
+# coordinate system, e.g. "J2000" or "galactic" into a celestial coordinate
+# structure. If the input wcs is an image wcs then a non-NULL pointer to
+# the image wcs structure is also returned. ERR is returned if a valid
+# celestial coordinate structure cannot be created.
+
+int procedure sk_decwcs (instr, mw, coo, imcoo)
+
+char	instr[ARB]		#I the input wcs string
+pointer	mw			#O the pointer to the image wcs structure
+pointer	coo			#O the pointer to the coordinate structure
+pointer	imcoo			#I pointer to an existing coordinate structure 
+
+int	stat
+pointer	sp, str1, str2, laxno, paxval, im
+int	sk_strwcs(), sk_decim()
+pointer	immap()
+errchk	immap()
+
+begin
+	call calloc (coo, LEN_SKYCOOSTRUCT, TY_STRUCT)
+	call strcpy (instr, SKY_COOSYSTEM(coo), SZ_FNAME)
+
+	# Allocate some working space.
+	call smark (sp)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (str2, SZ_LINE, TY_CHAR)
+	call salloc (laxno, IM_MAXDIM, TY_INT)
+	call salloc (paxval, IM_MAXDIM, TY_INT)
+
+	# Decode the wcs.
+	call sscan (instr)
+	    call gargwrd (Memc[str1], SZ_LINE)
+	    call gargwrd (Memc[str2], SZ_LINE)
+
+	# First try to open an image wcs.
+	iferr {
+	    im = immap (Memc[str1], READ_ONLY, 0)
+
+	# Decode the user wcs.
+	} then {
+
+	    # Initialize.
+	    mw = NULL
+	    if (imcoo == NULL) {
+	        SKY_NLNGAX(coo) = 2048
+	        SKY_NLATAX(coo) = 2048
+	        SKY_PLNGAX(coo) = 1
+	        SKY_PLATAX(coo) = 2
+	        SKY_XLAX(coo) = 1
+	        SKY_YLAX(coo) = 2
+	        SKY_VXOFF(coo) = 0.0d0
+	        SKY_VYOFF(coo) = 0.0d0
+	        SKY_VXSTEP(coo) = 1.0d0
+	        SKY_VYSTEP(coo) = 1.0d0
+	        SKY_WTYPE(coo) = 0
+	    } else {
+	        SKY_NLNGAX(coo) = SKY_NLNGAX(imcoo)
+	        SKY_NLATAX(coo) = SKY_NLATAX(imcoo)
+	        SKY_PLNGAX(coo) = SKY_PLNGAX(imcoo)
+	        SKY_PLATAX(coo) = SKY_PLATAX(imcoo)
+	        SKY_XLAX(coo) = SKY_XLAX(imcoo)
+	        SKY_YLAX(coo) = SKY_YLAX(imcoo)
+	        SKY_VXOFF(coo) = SKY_VXOFF(imcoo)
+	        SKY_VYOFF(coo) = SKY_VYOFF(imcoo)
+	        SKY_VXSTEP(coo) = SKY_VXSTEP(imcoo)
+	        SKY_VYSTEP(coo) = SKY_VYSTEP(imcoo)
+	        SKY_WTYPE(coo) = SKY_WTYPE(imcoo)
+	    }
+	    SKY_PIXTYPE(coo) = PIXTYPE_WORLD
+
+	    # Decode the actual wcs.
+	    stat = sk_strwcs (instr, SKY_CTYPE(coo), SKY_RADECSYS(coo),
+	        SKY_EQUINOX(coo), SKY_EPOCH(coo))
+	    switch (SKY_CTYPE(coo)) {
+	    case CTYPE_EQUATORIAL:
+		SKY_NLNGUNITS(coo) = SKY_HOURS
+		SKY_NLATUNITS(coo) = SKY_DEGREES
+	    default:
+		SKY_NLNGUNITS(coo) = SKY_DEGREES
+		SKY_NLATUNITS(coo) = SKY_DEGREES
+	    }
+
+	# Decode the image wcs.
+	} else {
+	    stat = sk_decim (im, Memc[str2], mw, coo)
+	    call imunmap (im)
+	}
+
+	call sfree (sp)
+
+	SKY_STATUS(coo) = stat
+	return (stat)
+end
+
+
+# SK_DECWSTR -- Decode the wcs string coordinate system, e.g. "J2000" or
+# "galactic" into a celestial coordinate structure. ERR is returned if a
+# valid celestial coordinate structure cannot be created.
+
+int procedure sk_decwstr (instr, coo, imcoo)
+
+char	instr[ARB]		#I the input wcs string
+pointer	coo			#O the pointer to the coordinate structure
+pointer	imcoo			#I pointer to an existing coordinate structure 
+
+int	stat
+int	sk_strwcs()
+
+begin
+	call calloc (coo, LEN_SKYCOOSTRUCT, TY_STRUCT)
+	call strcpy (instr, SKY_COOSYSTEM(coo), SZ_FNAME)
+
+	# Initialize.
+	if (imcoo == NULL) {
+	    SKY_NLNGAX(coo) = 2048
+	    SKY_NLATAX(coo) = 2048
+	    SKY_PLNGAX(coo) = 1
+	    SKY_PLATAX(coo) = 2
+	    SKY_XLAX(coo) = 1
+	    SKY_YLAX(coo) = 2
+	    SKY_VXOFF(coo) = 0.0d0
+	    SKY_VYOFF(coo) = 0.0d0
+	    SKY_VXSTEP(coo) = 1.0d0
+	    SKY_VYSTEP(coo) = 1.0d0
+	    SKY_WTYPE(coo) = 0
+	} else {
+	    SKY_NLNGAX(coo) = SKY_NLNGAX(imcoo)
+	    SKY_NLATAX(coo) = SKY_NLATAX(imcoo)
+	    SKY_PLNGAX(coo) = SKY_PLNGAX(imcoo)
+	    SKY_PLATAX(coo) = SKY_PLATAX(imcoo)
+	    SKY_XLAX(coo) = SKY_XLAX(imcoo)
+	    SKY_YLAX(coo) = SKY_YLAX(imcoo)
+	    SKY_VXOFF(coo) = SKY_VXOFF(imcoo)
+	    SKY_VYOFF(coo) = SKY_VYOFF(imcoo)
+	    SKY_VXSTEP(coo) = SKY_VXSTEP(imcoo)
+	    SKY_VYSTEP(coo) = SKY_VYSTEP(imcoo)
+	    SKY_WTYPE(coo) = SKY_WTYPE(imcoo)
+	}
+	SKY_PIXTYPE(coo) = PIXTYPE_WORLD
+
+	# Decode the actual wcs.
+	stat = sk_strwcs (instr, SKY_CTYPE(coo), SKY_RADECSYS(coo),
+            SKY_EQUINOX(coo), SKY_EPOCH(coo))
+	switch (SKY_CTYPE(coo)) {
+	case CTYPE_EQUATORIAL:
+	    SKY_NLNGUNITS(coo) = SKY_HOURS
+	    SKY_NLATUNITS(coo) = SKY_DEGREES
+	default:
+	    SKY_NLNGUNITS(coo) = SKY_DEGREES
+	    SKY_NLATUNITS(coo) = SKY_DEGREES
+	}
+
+	SKY_STATUS(coo) = stat
+
+	return (stat)
+end
+
+
+# SK_DECIM -- Given an image descriptor and an image wcs string create a
+# celstial coordinate structure. A non-NULL pointer to the image wcs structure
+# is also returned. ERR is returned if a valid  celestial coordinate descriptor
+# cannot be created.
+
+
+int procedure sk_decim (im, wcs, mw, coo)
+
+pointer	im			#I the pointer to the input image
+char	wcs[ARB]		#I the wcs string [logical|tv|physical|world]
+pointer	mw			#O the pointer to the image wcs structure
+pointer	coo			#O the pointer to the coordinate structure
+
+int	stat
+pointer	sp, str1, laxno, paxval
+int	sk_imwcs(), strdic(), mw_stati()
+pointer	mw_openim()
+errchk	mw_openim()
+
+begin
+	call malloc (coo, LEN_SKYCOOSTRUCT, TY_STRUCT)
+	call sprintf (SKY_COOSYSTEM(coo), SZ_FNAME, "%s %s")
+	    call pargstr (IM_HDRFILE(im))
+	    call pargstr (wcs)
+
+	call smark (sp)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (laxno, IM_MAXDIM, TY_INT)
+	call salloc (paxval, IM_MAXDIM, TY_INT)
+
+	# Try to open the image wcs.
+	iferr {
+	    mw = mw_openim (im)
+
+	# Set up a dummy wcs.
+	} then {
+
+	    #Initialize.
+	    SKY_CTYPE(coo) = 0
+	    SKY_RADECSYS(coo) = 0
+	    SKY_EQUINOX(coo) = INDEFD
+	    SKY_EPOCH(coo) = INDEFD
+	    mw = NULL
+	    SKY_PLNGAX(coo) = 1
+	    SKY_PLATAX(coo) = 2
+	    SKY_XLAX(coo) = 1
+	    SKY_YLAX(coo) = 2
+	    SKY_NLNGAX(coo) = 2048
+	    SKY_NLATAX(coo) = 2048
+	    SKY_VXOFF(coo) = 0.0d0
+	    SKY_VYOFF(coo) = 0.0d0
+	    SKY_VXSTEP(coo) = 1.0d0
+	    SKY_VYSTEP(coo) = 1.0d0
+	    SKY_WTYPE(coo) = 0
+	    SKY_PIXTYPE(coo) = PIXTYPE_LOGICAL
+	    SKY_NLNGUNITS(coo) = SKY_DEGREES
+	    SKY_NLATUNITS(coo) = SKY_DEGREES
+	    stat = ERR
+
+	# Decode the wcs.
+	} else {
+	    SKY_PIXTYPE(coo) = strdic (wcs, Memc[str1], SZ_LINE, PIXTYPE_LIST)
+	    if (SKY_PIXTYPE(coo) <= 0)
+	        SKY_PIXTYPE(coo) = PIXTYPE_LOGICAL
+	    if (sk_imwcs (im, mw, SKY_CTYPE(coo), SKY_PLNGAX(coo),
+		    SKY_PLATAX(coo), SKY_WTYPE(coo), SKY_RADECSYS(coo),
+		    SKY_EQUINOX(coo), SKY_EPOCH(coo)) == OK) {
+	    	switch (SKY_CTYPE(coo)) {
+	    	case CTYPE_EQUATORIAL:
+	    	    SKY_NLNGUNITS(coo) = SKY_HOURS
+	    	    SKY_NLATUNITS(coo) = SKY_DEGREES
+	    	default:
+	    	    SKY_NLNGUNITS(coo) = SKY_DEGREES
+	    	    SKY_NLATUNITS(coo) = SKY_DEGREES
+	    	}
+		call mw_gaxmap (mw, Memi[laxno], Memi[paxval], mw_stati(mw,
+		    MW_NPHYSDIM))
+		if (Memi[laxno+SKY_PLNGAX(coo)-1] <
+		    Memi[laxno+SKY_PLATAX(coo)-1]) {
+		    SKY_XLAX(coo) = Memi[laxno+SKY_PLNGAX(coo)-1] 
+		    SKY_YLAX(coo) = Memi[laxno+SKY_PLATAX(coo)-1] 
+		} else {
+		    SKY_XLAX(coo) = Memi[laxno+SKY_PLATAX(coo)-1] 
+		    SKY_YLAX(coo) = Memi[laxno+SKY_PLNGAX(coo)-1] 
+		}
+		if (SKY_XLAX(coo) <= 0 || SKY_YLAX(coo) <= 0) {
+		    SKY_VXOFF(coo) = 0.0d0
+		    SKY_VYOFF(coo) = 0.0d0
+		    SKY_VXSTEP(coo) = 1.0d0
+		    SKY_VYSTEP(coo) = 1.0d0
+	            SKY_NLNGAX(coo) = 2048
+	            SKY_NLATAX(coo) = 2048
+		    stat = ERR
+		} else {
+		    SKY_VXOFF(coo) = IM_VOFF(im,IM_VMAP(im,SKY_XLAX(coo)))
+		    SKY_VYOFF(coo) = IM_VOFF(im,IM_VMAP(im,SKY_YLAX(coo)))
+		    SKY_VXSTEP(coo) = IM_VSTEP(im,SKY_XLAX(coo))
+		    SKY_VYSTEP(coo) = IM_VSTEP(im,SKY_YLAX(coo))
+	            SKY_NLNGAX(coo) = IM_LEN(im,SKY_XLAX(coo))
+	            SKY_NLATAX(coo) = IM_LEN(im,SKY_YLAX(coo))
+		    stat = OK
+		}
+	    } else {
+		call mw_close (mw)
+		mw = NULL
+	        SKY_XLAX(coo) = 1
+	        SKY_YLAX(coo) = 2
+	        SKY_NLNGAX(coo) = 2048
+	        SKY_NLATAX(coo) = 2048
+	        SKY_VXOFF(coo) = 0.0d0
+	        SKY_VYOFF(coo) = 0.0d0
+	        SKY_VXSTEP(coo) = 1.0d0
+	        SKY_VYSTEP(coo) = 1.0d0
+	    	SKY_NLNGUNITS(coo) = SKY_DEGREES
+	    	SKY_NLATUNITS(coo) = SKY_DEGREES
+		stat = ERR
+	    }
+	}
+
+	call sfree (sp)
+
+	SKY_STATUS(coo) = stat
+	return (stat)
+end
+
+
+# SK_STRWCS -- Decode the sky coordinate system from an input string.
+# The string syntax is [ctype] equinox [epoch]. The various options
+# have been placed case statements. Although there is considerable
+# duplication of code in the case statements, there are minor differences
+# and I found it clearer to write it out rather than trying to be
+# concise. I might want to clean this up a bit later.
+
+int procedure sk_strwcs (instr, ctype, radecsys, equinox, epoch)
+
+char	instr[ARB]		#I the input wcs string
+int	ctype			#O the output coordinate type
+int	radecsys		#O the output equatorial reference system
+double	equinox			#O the output equinox
+double	epoch			#O the output epoch of the observation
+
+int	ip, nitems, sctype, sradecsys, stat
+pointer	sp, str1, str2
+int	strdic(), nscan(), ctod()
+double	sl_ej2d(), sl_epb(), sl_eb2d(), sl_epj()
+
+begin
+	# Initialize.
+	ctype = 0
+	radecsys = 0
+	equinox = INDEFD
+	epoch = INDEFD
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str1, SZ_LINE, TY_CHAR)
+	call salloc (str2, SZ_LINE, TY_CHAR)
+
+	# Determine the coordinate string.
+	call sscan (instr)
+	    call gargwrd (Memc[str1], SZ_LINE)
+
+	# Return with an error if the string is blank.
+	if (Memc[str1] == EOS || nscan() < 1) {
+	    call sfree (sp)
+	    return (ERR)
+	} else
+	    nitems = 1
+	    
+	# If the coordinate type is undefined temporarily default it to
+	# equatorial.
+	sctype = strdic (Memc[str1], Memc[str2], SZ_LINE, FTYPE_LIST) 
+	if (sctype <= 0) {
+	    ctype = CTYPE_EQUATORIAL
+	} else {
+	    switch (sctype) {
+	    case FTYPE_FK4:
+	        ctype = CTYPE_EQUATORIAL
+		radecsys = EQTYPE_FK4
+	    case FTYPE_FK4NOE:
+	        ctype = CTYPE_EQUATORIAL
+		radecsys = EQTYPE_FK4NOE
+	    case FTYPE_FK5:
+	        ctype = CTYPE_EQUATORIAL
+		radecsys = EQTYPE_FK5
+	    case FTYPE_ICRS:
+	        ctype = CTYPE_EQUATORIAL
+		radecsys = EQTYPE_ICRS
+	    case FTYPE_GAPPT:
+	        ctype = CTYPE_EQUATORIAL
+		radecsys = EQTYPE_GAPPT
+	    case FTYPE_ECLIPTIC:
+	        ctype = CTYPE_ECLIPTIC
+	    case FTYPE_GALACTIC:
+	        ctype = CTYPE_GALACTIC
+	    case FTYPE_SUPERGALACTIC:
+	        ctype = CTYPE_SUPERGALACTIC
+	    }
+	    call gargwrd (Memc[str1], SZ_LINE)
+	    if (nscan() > nitems)
+		nitems = nitems + 1
+	}
+	sctype = ctype
+	sradecsys = radecsys
+
+	# Decode the coordinate system.
+	switch (sctype) {
+
+	# Decode the equatorial system, equinox, and epoch.
+	case CTYPE_EQUATORIAL:
+
+	    switch (sradecsys) {
+	    case EQTYPE_FK4, EQTYPE_FK4NOE:
+		if (Memc[str1] == 'J' || Memc[str1] == 'j' ||
+		    Memc[str1] == 'B' || Memc[str1] == 'b')
+		    ip = 2
+		else 
+		    ip = 1
+	        if (ctod (Memc[str1], ip, equinox) <= 0)
+		    equinox = 1950.0d0
+		if (Memc[str1] == 'J' || Memc[str1] == 'j')
+		    equinox = sl_epb (sl_ej2d (equinox))
+
+	        call gargwrd (Memc[str2], SZ_LINE)
+	        if (nscan() <= nitems)
+		    epoch = sl_eb2d (equinox)
+		else {
+		    if (Memc[str2] == 'J' || Memc[str2] == 'j' ||
+		        Memc[str2] == 'B' || Memc[str2] == 'b')
+		        ip = 2
+		    else
+		        ip = 1
+		    if (ctod (Memc[str2], ip, epoch) <= 0)
+		        epoch = sl_eb2d (equinox)
+		    else if (epoch <= 3000.0d0 && (Memc[str2] == 'J' ||
+		        Memc[str2] == 'j'))
+		        epoch = sl_ej2d (epoch)
+		    else if (epoch > 3000.0d0)
+		        epoch = epoch - 2400000.5d0
+		    else
+			epoch = sl_eb2d (epoch)
+		}
+
+	    case EQTYPE_FK5, EQTYPE_ICRS:
+		if (Memc[str1] == 'J' || Memc[str1] == 'j' ||
+		    Memc[str1] == 'B' || Memc[str1] == 'b')
+		    ip = 2
+		else 
+		    ip = 1
+	        if (ctod (Memc[str1], ip, equinox) <= 0)
+		    equinox = 2000.0d0
+		if (Memc[str1] == 'B' || Memc[str1] == 'b')
+		    equinox = sl_epj(sl_eb2d (equinox))
+
+	        call gargwrd (Memc[str2], SZ_LINE)
+	        if (nscan() <= nitems)
+		    epoch = sl_ej2d (equinox)
+		else {
+		    if (Memc[str2] == 'J' || Memc[str2] == 'j' ||
+		        Memc[str2] == 'B' || Memc[str2] == 'b')
+		        ip = 2
+		    else
+		        ip = 1
+		    if (ctod (Memc[str2], ip, epoch) <= 0)
+		        epoch = sl_ej2d (equinox)
+		    else if (epoch <= 3000.0d0 && (Memc[str2] == 'B' ||
+		        Memc[str2] == 'b'))
+		        epoch = sl_eb2d (epoch)
+		    else if (epoch > 3000.0d0)
+		        epoch = epoch - 2400000.5d0
+		    else
+			epoch = sl_ej2d (epoch)
+		}
+
+	    case EQTYPE_GAPPT:
+		equinox = 2000.0d0
+		if (Memc[str1] == 'J' || Memc[str1] == 'j' ||
+		    Memc[str1] == 'B' || Memc[str1] == 'b')
+		    ip = 2
+		else
+		    ip = 1
+	        if (ctod (Memc[str1], ip, epoch) <= 0) {
+		    epoch = INDEFD
+		} else if (epoch <= 3000.0d0) {
+		    if (Memc[str1] == 'B' || Memc[str1] == 'b')
+		        epoch = sl_eb2d (epoch)
+		    else if (Memc[str1] == 'J' || Memc[str1] == 'j')
+		        epoch = sl_ej2d (epoch)
+		    else if (epoch < 1984.0d0)
+		        epoch = sl_eb2d (epoch)
+		    else
+		        epoch = sl_ej2d (epoch)
+		} else {
+		    epoch = epoch - 2400000.5d0
+		} 
+
+	    default:
+		ip = 1
+		if (Memc[str1] == 'B' || Memc[str1] == 'b') {
+		    radecsys = EQTYPE_FK4
+		    ip = ip + 1
+	            if (ctod (Memc[str1], ip, equinox) <= 0)
+			equinox = 1950.0d0
+
+	            call gargwrd (Memc[str2], SZ_LINE)
+	            if (nscan() <= nitems)
+		        epoch = sl_eb2d (equinox)
+		    else {
+		        if (Memc[str2] == 'J' || Memc[str2] == 'j')
+		            ip = 2
+		        else if (Memc[str2] == 'B' || Memc[str2] == 'b')
+		            ip = 2
+		        else
+		            ip = 1
+		        if (ctod (Memc[str2], ip, epoch) <= 0)
+		            epoch = sl_eb2d (equinox)
+		        else if (epoch <= 3000.0d0 && (Memc[str2] == 'J' ||
+			    Memc[str2] == 'j'))
+		            epoch = sl_ej2d (epoch)
+		        else if (epoch > 3000.0d0)
+		            epoch = epoch - 2400000.5d0
+			else
+			    epoch = sl_eb2d (epoch)
+		    }
+
+		} else if (Memc[str1] == 'J' || Memc[str1] == 'j') {
+		    radecsys = EQTYPE_FK5
+		    ip = ip + 1
+	            if (ctod (Memc[str1], ip, equinox) <= 0)
+			equinox = 2000.0d0
+
+	            call gargwrd (Memc[str2], SZ_LINE)
+	            if (nscan() <= nitems)
+		        epoch = sl_ej2d (equinox)
+		    else {
+		        if (Memc[str2] == 'J' || Memc[str2] == 'j' ||
+		            Memc[str2] == 'B' || Memc[str2] == 'b')
+		            ip = 2
+		        else
+		            ip = 1
+		        if (ctod (Memc[str2], ip, epoch) <= 0)
+		            epoch = sl_ej2d (equinox)
+		        else if (epoch <= 3000.0d0 && (Memc[str2] == 'B' ||
+			    Memc[str2] == 'b'))
+		            epoch = sl_eb2d (epoch)
+		        else if (epoch > 3000.0d0)
+		            epoch = epoch - 2400000.5d0
+			else
+			    epoch = sl_ej2d (epoch)
+		    }
+
+		} else if (ctod (Memc[str1], ip, equinox) <= 0) {
+		    ctype = 0
+		    radecsys = 0
+		    equinox = INDEFD
+		    epoch = INDEFD
+
+		} else if (equinox < 1984.0d0) {
+		    radecsys = EQTYPE_FK4
+	            call gargwrd (Memc[str2], SZ_LINE)
+	            if (nscan() <= nitems)
+		        epoch = sl_eb2d (equinox)
+		    else {
+		        if (Memc[str2] == 'J' || Memc[str2] == 'j' ||
+		            Memc[str2] == 'B' || Memc[str2] == 'b')
+		            ip = 2
+		        else
+		            ip = 1
+		        if (ctod (Memc[str2], ip, epoch) <= 0)
+		            epoch = sl_eb2d (equinox)
+		        else if (epoch <= 3000.0d0 && (Memc[str2] == 'J' ||
+			    Memc[str2] == 'j'))
+		            epoch = sl_ej2d (epoch)
+		        else if (epoch > 3000.0d0)
+		            epoch = epoch - 2400000.5d0
+			else
+			    epoch = sl_eb2d (epoch)
+		    }
+
+		} else {
+		    radecsys = EQTYPE_FK5
+	            call gargwrd (Memc[str2], SZ_LINE)
+	            if (nscan() <= nitems)
+		        epoch = sl_ej2d (equinox)
+		    else {
+		        if (Memc[str2] == 'J' || Memc[str2] == 'j' ||
+		            Memc[str2] == 'B' || Memc[str2] == 'b')
+		            ip = 2
+		        else
+		            ip = 1
+		        if (ctod (Memc[str2], ip, epoch) <= 0)
+		            epoch = sl_ej2d (equinox)
+		        else if (epoch <= 3000.0d0 && (Memc[str2] == 'B' ||
+			    Memc[str2] == 'b'))
+		            epoch = sl_eb2d (epoch)
+		        else if (epoch > 3000.0d0)
+		            epoch = epoch - 2400000.5d0
+			else
+			    epoch = sl_ej2d (epoch)
+		    }
+		}
+	    }
+
+	# Decode the ecliptic coordinate system.
+	case CTYPE_ECLIPTIC:
+	    if (Memc[str1] == 'J' || Memc[str1] == 'j' ||
+	        Memc[str1] == 'B' || Memc[str1] == 'b')
+		ip = 2
+	    else
+		ip = 1
+	    if (ctod (Memc[str1], ip, epoch) <= 0) {
+		epoch = INDEFD
+	    } else if (epoch <= 3000.0d0) {
+	        if (Memc[str1] == 'B' || Memc[str1] == 'b')
+		    epoch = sl_eb2d (epoch)
+	        else if (Memc[str1] == 'J' || Memc[str1] == 'j')
+		    epoch = sl_ej2d (epoch)
+		else if (epoch < 1984.0d0)
+		    epoch = sl_eb2d (epoch)
+		else
+		    epoch = sl_ej2d (epoch)
+	    } else {
+		epoch = epoch - 2400000.5d0
+	    }
+
+	# Decode the galactic and supergalactic coordinate system.
+	case CTYPE_GALACTIC, CTYPE_SUPERGALACTIC:
+	    if (Memc[str1] == 'J' || Memc[str1] == 'j' ||
+		Memc[str1] == 'B' || Memc[str1] == 'b')
+		ip = 2
+	    else
+		ip = 1
+	    if (ctod (Memc[str1], ip, epoch) <= 0) {
+		epoch = sl_eb2d (1950.0d0)
+	    } else if (epoch <= 3000.0d0) {
+	        if (Memc[str1] == 'J' || Memc[str1] == 'j')
+		    epoch = sl_ej2d (epoch)
+	        else if (Memc[str1] == 'B' || Memc[str1] == 'b')
+		    epoch = sl_eb2d (epoch)
+	        else if (epoch < 1984.0d0) 
+		    epoch = sl_eb2d (epoch)
+		else
+		    epoch = sl_ej2d (epoch)
+	    } else {
+		epoch = epoch - 2400000.5d0
+	    }
+	}
+
+	# Return the appropriate error status.
+	if (ctype == 0)
+	    stat = ERR
+	else if (ctype == CTYPE_EQUATORIAL && (radecsys == 0 ||
+	    IS_INDEFD(equinox) || IS_INDEFD(epoch)))
+	    stat = ERR
+	else if (ctype == CTYPE_ECLIPTIC && IS_INDEFD(epoch))
+	    stat = ERR
+	else
+	    stat = OK
+
+	call sfree (sp)
+
+	return (stat)
+end
+
+
+# SK_IMWCS -- Decode the sky coordinate system of the image. Return
+# an error if the sky coordinate system is not one of the supported types
+# or required information is missing from the image header.
+
+int procedure sk_imwcs (im, mw, ctype, lngax, latax, wtype, radecsys,
+	equinox, epoch)
+
+pointer	im			#I the image pointer
+pointer	mw			#I pointer to the world coordinate system
+int	ctype			#O the output coordinate type
+int	lngax			#O the output ra/glon/elon axis
+int	latax			#O the output dec/glat/elat axis
+int	wtype			#O the output projection type
+int	radecsys		#O the output equatorial reference system
+double	equinox			#O the output equinox
+double	epoch			#O the output epoch of the observation
+
+int	i, ndim, axtype, day, month, year, ier, oldfits
+pointer	sp, atval
+double	hours
+double	imgetd(), sl_eb2d(), sl_ej2d()
+int	mw_stati(), strdic(), dtm_decode()
+errchk	mw_gwattrs(), imgstr(), imgetd()
+
+begin
+	call smark (sp)
+	call salloc (atval, SZ_LINE, TY_CHAR)
+
+	# Initialize
+	ctype = 0
+	lngax = 0
+	latax = 0
+	wtype = 0
+	radecsys = 0
+	equinox = INDEFD
+	epoch = INDEFD
+
+	# Determine the sky coordinate system of the image.
+	ndim = mw_stati (mw, MW_NPHYSDIM)
+	do i = 1, ndim {
+	    iferr (call mw_gwattrs (mw, i, "axtype", Memc[atval], SZ_LINE))
+		call strcpy ("INDEF", Memc[atval], SZ_LINE)
+	    axtype = strdic (Memc[atval], Memc[atval], SZ_LINE, AXTYPE_LIST)
+	    switch (axtype) {
+	    case AXTYPE_RA, AXTYPE_DEC: 
+		ctype = CTYPE_EQUATORIAL
+	    case AXTYPE_ELON, AXTYPE_ELAT: 
+		ctype = CTYPE_ECLIPTIC
+	    case AXTYPE_GLON, AXTYPE_GLAT: 
+		ctype = CTYPE_GALACTIC
+	    case AXTYPE_SLON, AXTYPE_SLAT: 
+		ctype = CTYPE_SUPERGALACTIC
+	    default:
+		;
+	    }
+	    switch (axtype) {
+	    case AXTYPE_RA, AXTYPE_ELON, AXTYPE_GLON, AXTYPE_SLON:
+		lngax = i
+	    case AXTYPE_DEC, AXTYPE_ELAT, AXTYPE_GLAT, AXTYPE_SLAT:
+		latax = i
+	    default:
+		;
+	    }
+	}
+
+	# Return if the sky coordinate system cannot be decoded.
+	if (ctype == 0 || lngax == 0 || latax == 0) {
+	    call sfree (sp)
+	    return (ERR)
+	}
+
+	# Decode the sky projection.
+	iferr {
+	    call mw_gwattrs (mw, lngax, "wtype", Memc[atval], SZ_LINE)
+	} then {
+	    iferr (call mw_gwattrs(mw, latax, "wtype", Memc[atval], SZ_LINE))
+		call strcpy ("linear", Memc[atval], SZ_LINE)
+	}
+	wtype = strdic (Memc[atval], Memc[atval], SZ_LINE, WTYPE_LIST)
+
+	# Return if the sky projection system is not supported.
+	if (wtype == 0) {
+	    call sfree (sp)
+	    return (ERR)
+	}
+
+	# Determine the RA/DEC system and equinox.
+	if (ctype == CTYPE_EQUATORIAL) {
+
+	    # Get the equinox of the coordinate system. The EQUINOX keyword
+	    # takes precedence over EPOCH.
+	    iferr {
+	        equinox = imgetd (im, "EQUINOX")
+	    } then {
+		iferr {
+	            equinox = imgetd (im, "EPOCH")
+		} then {
+		    equinox = INDEFD
+		}
+	    }
+
+	    # Determine which equatorial system will be used. The default
+	    # is FK4 if equinox < 1984.0, FK5 if equinox is >= 1984.
+	    iferr {
+	        call imgstr (im, "RADECSYS", Memc[atval], SZ_LINE)
+	    } then {
+	        radecsys = 0
+	    } else {
+		call strlwr (Memc[atval])
+	        radecsys = strdic (Memc[atval], Memc[atval], SZ_LINE,
+		    EQTYPE_LIST)
+	    }
+	    if (radecsys == 0) {
+		if (IS_INDEFD(equinox))
+		    radecsys = EQTYPE_FK5
+		else if (equinox < 1984.0d0)
+		    radecsys = EQTYPE_FK4
+		else
+		    radecsys = EQTYPE_FK5
+	    }
+
+	    # Get the MJD of the observation. If there is no MJD in the
+	    # header use the DATE_OBS keyword value and transform it to
+	    # an MJD.
+	    iferr {
+	        epoch = imgetd (im, "MJD-WCS")
+	    } then {
+	        iferr {
+	            epoch = imgetd (im, "MJD-OBS")
+	        } then {
+		    iferr {
+	                call imgstr (im, "DATE-OBS", Memc[atval], SZ_LINE)
+		    } then {
+		        epoch = INDEFD
+		    } else if (dtm_decode (Memc[atval], year, month, day,
+			hours, oldfits) == OK) {
+		        call sl_cadj (year, month, day, epoch, ier)
+		        if (ier != 0)
+			    epoch = INDEFD
+			else if (! IS_INDEFD(hours) && hours >= 0.0d0 &&
+			    hours <= 24.0d0)
+			    epoch = epoch + hours / 24.0d0
+		    } else
+		        epoch = INDEFD
+	        }
+	    }
+
+	    # Set the default equinox and epoch appropriate for each
+	    # equatorial system if these are undefined.
+	    switch (radecsys) {
+	    case EQTYPE_FK4, EQTYPE_FK4NOE:
+		if (IS_INDEFD(equinox))
+		    equinox = 1950.0d0
+		if (IS_INDEFD(epoch))
+		    epoch = sl_eb2d (1950.0d0)
+	    case EQTYPE_FK5, EQTYPE_ICRS:
+		if (IS_INDEFD(equinox))
+		    equinox = 2000.0d0
+		if (IS_INDEFD(epoch))
+		    epoch = sl_ej2d (2000.0d0)
+	    case EQTYPE_GAPPT:
+		equinox = 2000.0d0
+		;
+	    }
+
+	    # Return if the epoch is undefined. This can only occur if
+	    # the equatorial coordinate system is GAPPT and there is NO
+	    # epoch of observation in the image header.
+	    if (IS_INDEFD(epoch)) {
+		call sfree (sp)
+		return (ERR)
+	    }
+	} 
+
+	# Get the MJD of the observation. If there is no MJD in the
+	# header use the DATE_OBS keyword value and transform it to
+	# an MJD.
+	if (ctype == CTYPE_ECLIPTIC) {
+
+	    iferr {
+	        epoch = imgetd (im, "MJD-WCS")
+	    } then {
+	        iferr {
+	            epoch = imgetd (im, "MJD-OBS")
+	        } then {
+		    iferr {
+	                call imgstr (im, "DATE-OBS", Memc[atval], SZ_LINE)
+		    } then {
+		        epoch = INDEFD
+		    } else if (dtm_decode (Memc[atval], year, month, day,
+			hours, oldfits) == OK) {
+		        call sl_cadj (year, month, day, epoch, ier)
+		        if (ier != 0)
+			    epoch = INDEFD
+			else if (! IS_INDEFD(hours) && hours >= 0.0d0 &&
+			    hours <= 24.0d0)
+			    epoch = epoch + hours / 24.0d0
+		    } else
+		        epoch = INDEFD
+	        }
+	    }
+
+	    # Return if the epoch is undefined.
+	    if (IS_INDEFD(epoch)) {
+		call sfree (sp)
+		return (ERR)
+	    }
+	}
+
+	if (ctype == CTYPE_GALACTIC || ctype == CTYPE_SUPERGALACTIC) {
+
+	    # Get the MJD of the observation. If there is no MJD in the
+	    # header use the DATE_OBS keyword value and transform it to
+	    # an MJD.
+	    iferr {
+	        epoch = imgetd (im, "MJD-WCS")
+	    } then {
+	        iferr {
+	            epoch = imgetd (im, "MJD-OBS")
+	        } then {
+		    iferr {
+	                call imgstr (im, "DATE-OBS", Memc[atval], SZ_LINE)
+		    } then {
+		        epoch = sl_eb2d (1950.0d0)
+		    } else if (dtm_decode (Memc[atval], year, month, day,
+		        hours, oldfits) == OK) {
+		        call sl_cadj (year, month, day, epoch, ier)
+		        if (ier != 0)
+			    epoch = sl_eb2d (1950.0d0)
+			else {
+			    if (! IS_INDEFD(hours) && hours >= 0.0d0 &&
+				hours <= 24.0d0)
+			        epoch = epoch + hours / 24.0d0
+			    #if (epoch < 1984.0d0)
+		    	        #epoch = sl_eb2d (epoch)
+			    #else
+		    	        #epoch = sl_ej2d (epoch)
+			}
+		    } else
+		        epoch = sl_eb2d (1950.0d0)
+	        }
+	    }
+	}
+
+	call sfree (sp)
+
+	return (OK)
+end
+
+
+# SK_ENWCS -- Encode the celestial wcs system.
+
+procedure sk_enwcs (coo, wcsstr, maxch)
+
+pointer	coo			#I the celestial coordinate system descriptor
+char	wcsstr[ARB]		#O the output wcs string
+int	maxch			#I the size of the output string
+
+double	sk_statd(), sl_epj(), sl_epb()
+int	sk_stati()
+
+begin
+	switch (sk_stati (coo, S_CTYPE)) {
+
+	case CTYPE_EQUATORIAL:
+
+	    switch (sk_stati(coo, S_RADECSYS)) {
+
+	    case EQTYPE_GAPPT:
+		if (IS_INDEFD(sk_statd(coo, S_EPOCH))) {
+	            call sprintf (wcsstr, maxch, "apparent")
+		} else {
+	            call sprintf (wcsstr, maxch, "apparent J%0.8f")
+		        call pargd (sl_epj(sk_statd(coo, S_EPOCH)))
+		}
+
+	    case EQTYPE_FK5:
+	        call sprintf (wcsstr, maxch, "fk5 J%0.3f J%0.8f")
+		    call pargd (sk_statd(coo, S_EQUINOX))
+		    call pargd (sl_epj(sk_statd(coo, S_EPOCH)))
+
+	    case EQTYPE_ICRS:
+	        call sprintf (wcsstr, maxch, "icrs J%0.3f J%0.8f")
+		    call pargd (sk_statd(coo, S_EQUINOX))
+		    call pargd (sl_epj(sk_statd(coo, S_EPOCH)))
+
+	    case EQTYPE_FK4:
+	        call sprintf (wcsstr, maxch, "fk4 B%0.3f B%0.8f")
+		    call pargd (sk_statd(coo, S_EQUINOX))
+		    call pargd (sl_epb(sk_statd(coo, S_EPOCH)))
+
+	    case EQTYPE_FK4NOE:
+	        call sprintf (wcsstr, maxch, "fk4noe B%0.3f B%0.8f")
+		    call pargd (sk_statd(coo, S_EQUINOX))
+		    call pargd (sl_epb(sk_statd(coo, S_EPOCH)))
+
+	    default:
+		wcsstr[1] = EOS
+	    }
+
+	case CTYPE_ECLIPTIC:
+	    if (IS_INDEFD(sk_statd(coo, S_EPOCH))) {
+	        call sprintf (wcsstr, maxch, "ecliptic")
+	    } else {
+	        call sprintf (wcsstr, maxch, "ecliptic J%0.8f")
+		    call pargd (sl_epj(sk_statd(coo, S_EPOCH)))
+	    }
+
+	case CTYPE_GALACTIC:
+	    call sprintf (wcsstr, maxch, "galactic J%0.8f")
+		call pargd (sl_epj(sk_statd(coo, S_EPOCH)))
+
+	case CTYPE_SUPERGALACTIC:
+	    call sprintf (wcsstr, maxch, "supergalactic j%0.8f")
+		call pargd (sl_epj(sk_statd(coo, S_EPOCH)))
+	}
+end
+
+
+# SK_COPY -- Copy the coodinate structure.
+
+pointer procedure sk_copy (cooin)
+
+pointer	cooin			#I the pointer to the input structure
+
+pointer	cooout
+
+begin
+	if (cooin == NULL)
+	    cooout = NULL
+	else {
+	    call calloc (cooout, LEN_SKYCOOSTRUCT, TY_STRUCT)
+            SKY_VXOFF(cooout) = SKY_VXOFF(cooin)
+            SKY_VYOFF(cooout) = SKY_VYOFF(cooin)
+            SKY_VXSTEP(cooout) = SKY_VXSTEP(cooin)
+            SKY_VYSTEP(cooout) = SKY_VYSTEP(cooin)
+	    SKY_EQUINOX(cooout) = SKY_EQUINOX(cooin)
+	    SKY_EPOCH(cooout) = SKY_EPOCH(cooin)
+	    SKY_CTYPE(cooout) = SKY_CTYPE(cooin)
+	    SKY_RADECSYS(cooout) = SKY_RADECSYS(cooin)
+            SKY_WTYPE(cooout) = SKY_WTYPE(cooin)
+            SKY_PLNGAX(cooout) = SKY_PLNGAX(cooin)
+            SKY_PLATAX(cooout) = SKY_PLATAX(cooin)
+            SKY_XLAX(cooout) = SKY_XLAX(cooin)
+            SKY_YLAX(cooout) = SKY_YLAX(cooin)
+            SKY_PIXTYPE(cooout) = SKY_PIXTYPE(cooin)
+    	    SKY_NLNGAX(cooout) = SKY_NLNGAX(cooin)
+            SKY_NLATAX(cooout) = SKY_NLATAX(cooin)
+    	    SKY_NLNGUNITS(cooout) = SKY_NLNGUNITS(cooin)
+            SKY_NLATUNITS(cooout) = SKY_NLATUNITS(cooin)
+	    call strcpy (SKY_COOSYSTEM(cooin), SKY_COOSYSTEM(cooout),
+		SZ_FNAME)
+	}
+
+	return (cooout)
+end
+
+
+# SK_CLOSE -- Free the coordinate structure.
+
+procedure sk_close (coo)
+
+pointer	coo			#U the input coordinate structure
+
+begin
+	if (coo != NULL)
+	    call mfree (coo, TY_STRUCT)
+end
diff --git a/pkg/xtools/skywcs/sksaveim.x b/pkg/xtools/skywcs/sksaveim.x
new file mode 100644
index 00000000..77b5a1d9
--- /dev/null
+++ b/pkg/xtools/skywcs/sksaveim.x
@@ -0,0 +1,157 @@
+include "skywcsdef.h"
+include "skywcs.h"
+
+# SK_SAVEIM -- Update the image header keywords that describe the
+# fundamental coordinate system, CTYPE, RADECSYS, EQUINOX (EPOCH), and
+# MJD-WCS.
+
+procedure sk_saveim (coo, mw, im)
+
+pointer	coo			#I pointer to the coordinate structure
+pointer	mw			#I pointer to the mwcs structure
+pointer	im			#I image descriptor
+
+errchk	imdelf()
+
+begin
+	# Move all this to a separate routine
+	switch (SKY_CTYPE(coo)) {
+
+	case CTYPE_EQUATORIAL:
+	    call mw_swattrs (mw, SKY_PLNGAX(coo), "axtype", "ra")
+	    call mw_swattrs (mw, SKY_PLATAX(coo), "axtype", "dec")
+	    switch (SKY_RADECSYS(coo)) {
+	    case EQTYPE_FK4:
+		call imastr (im, "radecsys", "FK4")
+		call imaddd (im, "equinox", SKY_EQUINOX(coo))
+		call imaddd (im, "mjd-wcs", SKY_EPOCH(coo))
+	    case EQTYPE_FK4NOE:
+		call imastr (im, "radecsys", "FK4NOE")
+		call imaddd (im, "equinox", SKY_EQUINOX(coo))
+		call imaddd (im, "mjd-wcs", SKY_EPOCH(coo))
+	    case EQTYPE_FK5:
+		call imastr (im, "radecsys", "FK5")
+		call imaddd (im, "equinox", SKY_EQUINOX(coo))
+	        iferr (call imdelf (im, "mjd-wcs"))
+		    ;
+	    case EQTYPE_ICRS:
+		call imastr (im, "radecsys", "ICRS")
+		call imaddd (im, "equinox", SKY_EQUINOX(coo))
+	        iferr (call imdelf (im, "mjd-wcs"))
+		    ;
+	    case EQTYPE_GAPPT:
+		call imastr (im, "radecsys", "GAPPT")
+		iferr (call imdelf (im, "equinox"))
+		    ;
+		call imaddd (im, "mjd-wcs", SKY_EPOCH(coo))
+	    }
+
+	case CTYPE_ECLIPTIC:
+	    call mw_swattrs (mw, SKY_PLNGAX(coo), "axtype", "elon")
+	    call mw_swattrs (mw, SKY_PLATAX(coo), "axtype", "elat")
+	    iferr (call imdelf (im, "radecsys"))
+		;
+	    iferr (call imdelf (im, "equinox"))
+		;
+	    call imaddd (im, "mjd-wcs", SKY_EPOCH(coo))
+
+	case CTYPE_GALACTIC:
+	    call mw_swattrs (mw, SKY_PLNGAX(coo), "axtype", "glon")
+	    call mw_swattrs (mw, SKY_PLATAX(coo), "axtype", "glat")
+	    iferr (call imdelf (im, "radecsys"))
+		;
+	    iferr (call imdelf (im, "equinox"))
+		;
+	    iferr (call imdelf (im, "mjd-wcs"))
+		;
+
+	case CTYPE_SUPERGALACTIC:
+	    call mw_swattrs (mw, SKY_PLNGAX(coo), "axtype", "slon")
+	    call mw_swattrs (mw, SKY_PLATAX(coo), "axtype", "slat")
+	    iferr (call imdelf (im, "radecsys"))
+		;
+	    iferr (call imdelf (im, "equinox"))
+		;
+	    iferr (call imdelf (im, "mjd-wcs"))
+		;
+	}
+end
+
+
+# SK_CTYPEIM -- Modify the CTYPE keywords appropriately. This step will
+# become unnecessary when MWCS is updated to deal with non-equatorial celestial
+# coordinate systems.
+
+procedure sk_ctypeim (coo, im)
+
+pointer	coo			#I pointer to the coordinate structure
+pointer	im			#I image descriptor
+
+pointer	sp, wtype, key1, key2, attr
+int	sk_wrdstr()
+
+begin
+	call smark (sp)
+	call salloc (key1, 8, TY_CHAR)
+	call salloc (key2, 8, TY_CHAR)
+	call salloc (wtype, 3, TY_CHAR)
+	call salloc (attr, 8, TY_CHAR)
+
+	call sprintf (Memc[key1], 8, "CTYPE%d")
+	    call pargi (SKY_PLNGAX(coo))
+	call sprintf (Memc[key2], 8, "CTYPE%d")
+	    call pargi (SKY_PLATAX(coo))
+
+	if (SKY_WTYPE(coo) <= 0 || SKY_WTYPE(coo) == WTYPE_LIN) {
+	    call imastr (im, Memc[key1], "LINEAR")
+	    call imastr (im, Memc[key2], "LINEAR")
+	    call sfree (sp)
+	    return
+	}
+
+	if (sk_wrdstr (SKY_WTYPE(coo), Memc[wtype], 3, WTYPE_LIST) <= 0)
+	    call strcpy ("tan", Memc[wtype], 3)
+	call strupr (Memc[wtype])
+
+	# Move all this to a separate routine
+	switch (SKY_CTYPE(coo)) {
+
+	case CTYPE_EQUATORIAL:
+	    call sprintf (Memc[attr], 8, "RA---%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key1], Memc[attr])
+	    call sprintf (Memc[attr], 8, "DEC--%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key2], Memc[attr])
+
+	case CTYPE_ECLIPTIC:
+	    call sprintf (Memc[attr], 8, "ELON-%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key1], Memc[attr])
+	    call sprintf (Memc[attr], 8, "ELAT-%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key2], Memc[attr])
+
+	case CTYPE_GALACTIC:
+	    call sprintf (Memc[attr], 8, "GLON-%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key1], Memc[attr])
+	    call sprintf (Memc[attr], 8, "GLAT-%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key2], Memc[attr])
+
+	case CTYPE_SUPERGALACTIC:
+	    call sprintf (Memc[attr], 8, "SLON-%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key1], Memc[attr])
+	    call sprintf (Memc[attr], 8, "SLAT-%3s") 
+		call pargstr (Memc[wtype])
+	    call imastr (im, Memc[key2], Memc[attr])
+
+	default:
+	    call imastr (im, Memc[key1], "LINEAR")
+	    call imastr (im, Memc[key2], "LINEAR")
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/xtools/skywcs/skset.x b/pkg/xtools/skywcs/skset.x
new file mode 100644
index 00000000..9e7191c3
--- /dev/null
+++ b/pkg/xtools/skywcs/skset.x
@@ -0,0 +1,90 @@
+include "skywcsdef.h"
+include "skywcs.h"
+
+
+# SK_SETD -- Set a double precision coordinate parameter.
+
+procedure sk_setd (coo, param, value)
+
+pointer	coo			#I pointer to the coordinate structure
+int	param			#I the input parameter
+double	value			#I the parameter value
+
+begin
+	switch (param) {
+	case S_VXOFF:
+	    SKY_VXOFF(coo) = value
+	case S_VYOFF:
+	    SKY_VYOFF(coo) = value
+	case S_VXSTEP:
+	    SKY_VXSTEP(coo) = value
+	case S_VYSTEP:
+	    SKY_VYSTEP(coo) = value
+	case S_EQUINOX:
+	    SKY_EQUINOX(coo) = value
+	case S_EPOCH:
+	    SKY_EPOCH(coo) = value
+	default:
+	    call error (0, "SKY_SETD: Unknown coordinate system parameter")
+	}
+end
+
+
+# SK_SETI -- Set an integer coordinate parameter.
+
+procedure sk_seti (coo, param, value)
+
+pointer	coo			#I pointer to the coordinate structure
+int	param			#I the input parameter
+int	value			#I the parameter value
+
+begin
+	switch (param) {
+	case S_CTYPE:
+	    SKY_CTYPE(coo) = value
+	case S_RADECSYS:
+	    SKY_RADECSYS(coo) = value
+	case S_WTYPE:
+	    SKY_WTYPE(coo) = value
+	case S_PLNGAX:
+	    SKY_PLNGAX(coo) = value
+	case S_PLATAX:
+	    SKY_PLATAX(coo) = value
+	case S_XLAX:
+	    SKY_XLAX(coo) = value
+	case S_YLAX:
+	    SKY_YLAX(coo) = value
+	case S_PIXTYPE:
+	    SKY_PIXTYPE(coo) = value
+	case S_NLNGAX:
+	    SKY_NLNGAX(coo) = value
+	case S_NLATAX:
+	    SKY_NLATAX(coo) = value
+	case S_NLNGUNITS:
+	    SKY_NLNGUNITS(coo) = value
+	case S_NLATUNITS:
+	    SKY_NLATUNITS(coo) = value
+	case S_STATUS:
+	    SKY_STATUS(coo) = value
+	default:
+	    call error (0, "SKY_SETI: Unknown coordinate system parameter")
+	}
+end
+
+
+# SK_SETS -- Set a character string coordinate parameter.
+
+procedure sk_sets (coo, param, value)
+
+pointer	coo			#I pointer to the coordinate structure
+int	param			#I the input parameter
+char	value[ARB]		#I the parameter value
+
+begin
+	switch (param) {
+	case S_COOSYSTEM:
+	    call strcpy (value, SKY_COOSYSTEM(coo), SZ_FNAME)
+	default:
+	    call error (0, "SKY_SETSTR: Unknown coordinate system parameter")
+	}
+end
diff --git a/pkg/xtools/skywcs/skstat.x b/pkg/xtools/skywcs/skstat.x
new file mode 100644
index 00000000..82d2f1c2
--- /dev/null
+++ b/pkg/xtools/skywcs/skstat.x
@@ -0,0 +1,90 @@
+include "skywcsdef.h"
+include "skywcs.h"
+
+
+# SK_STATD -- Get a double precision coordinate parameter.
+
+double procedure sk_statd (coo, param)
+
+pointer	coo			#I pointer to the coordinate structure
+int	param			#I the input parameter
+
+begin
+	switch (param) {
+	case S_VXOFF:
+	    return (SKY_VXOFF(coo))
+	case S_VYOFF:
+	    return (SKY_VYOFF(coo))
+	case S_VXSTEP:
+	    return (SKY_VXSTEP(coo))
+	case S_VYSTEP:
+	    return (SKY_VYSTEP(coo))
+	case S_EQUINOX:
+	    return (SKY_EQUINOX(coo))
+	case S_EPOCH:
+	    return (SKY_EPOCH(coo))
+	default:
+	    call error (0, "SKY_STATD: Unknown coordinate system parameter")
+	}
+end
+
+
+# SK_STATI -- Get an integer coordinate parameter.
+
+int procedure sk_stati (coo, param)
+
+pointer	coo			#I pointer to the coordinate structure
+int	param			#I the input parameter
+
+begin
+	switch (param) {
+	case S_CTYPE:
+	    return (SKY_CTYPE(coo))
+	case S_RADECSYS:
+	    return (SKY_RADECSYS(coo))
+	case S_WTYPE:
+	    return (SKY_WTYPE(coo))
+	case S_PLNGAX:
+	    return (SKY_PLNGAX(coo))
+	case S_PLATAX:
+	    return (SKY_PLATAX(coo))
+	case S_XLAX:
+	    return (SKY_XLAX(coo))
+	case S_YLAX:
+	    return (SKY_YLAX(coo))
+	case S_PIXTYPE:
+	    return (SKY_PIXTYPE(coo))
+	case S_NLNGAX:
+	    return (SKY_NLNGAX(coo))
+	case S_NLATAX:
+	    return (SKY_NLATAX(coo))
+	case S_NLNGUNITS:
+	    return (SKY_NLNGUNITS(coo))
+	case S_NLATUNITS:
+	    return (SKY_NLATUNITS(coo))
+	case S_STATUS:
+	    return (SKY_STATUS(coo))
+	default:
+	    call error (0, "SKY_STATI: Unknown coordinate system parameter")
+	}
+end
+
+
+
+# SK_STATS -- Get a character string coordinate parameter.
+
+procedure sk_stats (coo, param, value, maxch)
+
+pointer	coo			#I pointer to the coordinate structure
+int	param			#I the input parameter
+char	value			#O the output string
+int	maxch			#I the maximum size of the string
+
+begin
+	switch (param) {
+	case S_COOSYSTEM:
+	    call strcpy (SKY_COOSYSTEM(coo), value, maxch)
+	default:
+	    call error (0, "SKY_GETSTR: Unknown coordinate system parameter")
+	}
+end
diff --git a/pkg/xtools/skywcs/sktransform.x b/pkg/xtools/skywcs/sktransform.x
new file mode 100644
index 00000000..a8cf87c3
--- /dev/null
+++ b/pkg/xtools/skywcs/sktransform.x
@@ -0,0 +1,577 @@
+include <math.h>
+include "skywcsdef.h"
+include "skywcs.h"
+
+# SK_ULTRAN -- Transform the sky coordinates from the input coordinate
+# system to the output coordinate system using the units conversions as
+# appropriate.
+
+procedure sk_ultran (cooin, cooout, ilng, ilat, olng, olat, npts) 
+
+pointer	cooin		#I pointer to the input coordinate system structure
+pointer	cooout		#I pointer to the output coordinate system structure
+double	ilng[ARB]	#I the input ra/longitude in radians
+double	ilat[ARB]	#I the input dec/latitude in radians
+double	olng[ARB]	#O the output ra/longitude in radians
+double	olat[ARB]	#O the output dec/latitude in radians
+int	npts		#I the number of points to be converted
+
+double	tilng, tilat, tolng, tolat
+int	i
+
+begin
+	do i = 1, npts {
+
+	    switch (SKY_NLNGUNITS(cooin)) {
+	    case SKY_HOURS:
+		tilng = DEGTORAD(15.0d0 * ilng[i])
+	    case SKY_DEGREES:
+		tilng = DEGTORAD(ilng[i])
+	    case SKY_RADIANS:
+		tilng = ilng[i]
+	    default:
+		tilng = ilng[i]
+	    }
+	    switch (SKY_NLATUNITS(cooin)) {
+	    case SKY_HOURS:
+		tilat = DEGTORAD(15.0d0 * ilat[i])
+	    case SKY_DEGREES:
+		tilat = DEGTORAD(ilat[i])
+	    case SKY_RADIANS:
+		tilat = ilat[i]
+	    default:
+		tilat = ilat[i]
+	    }
+
+	    call sk_lltran (cooin, cooout, tilng, tilat, INDEFD, INDEFD,
+		0.0d0, 0.0d0, tolng, tolat)
+
+	    switch (SKY_NLNGUNITS(cooout)) {
+	    case SKY_HOURS:
+		olng[i] = RADTODEG(tolng) / 15.0d0
+	    case SKY_DEGREES:
+		olng[i] = RADTODEG(tolng)
+	    case SKY_RADIANS:
+		olng[i] = tolng
+	    default:
+		olng[i] = tolng
+	    }
+	    switch (SKY_NLATUNITS(cooout)) {
+	    case SKY_HOURS:
+		olat[i] = RADTODEG(tolat) / 15.0d0
+	    case SKY_DEGREES:
+		olat[i] = RADTODEG(tolat)
+	    case SKY_RADIANS:
+		olat[i] = tolat
+	    default:
+		olat[i] = tolat
+	    }
+	}
+end
+
+
+# SK_LLTRAN -- Transform the sky coordinate from the input coordinate
+# system to the output coordinate system assuming that all the coordinate
+# are in radians.
+
+procedure sk_lltran (cooin, cooout, ilng, ilat, ipmlng, ipmlat, px, rv,
+	olng, olat)
+
+pointer	cooin		#I pointer to the input coordinate system structure
+pointer	cooout		#I pointer to the output coordinate system structure
+double	ilng		#I the input ra/longitude in radians
+double	ilat		#I the input dec/latitude in radians
+double	ipmlng		#I the input proper motion in ra in radians
+double	ipmlat		#I the input proper motion in dec in radians
+double	px		#I the input parallax in arcseconds
+double	rv		#I the input radial velocity in km / second
+double	olng		#O the output ra/longitude in radians
+double	olat		#O the output dec/latitude in radians
+
+int	pmflag
+double	pmr, pmd
+double	sl_epj(), sl_epb()
+
+begin
+	# Test for the case where the input coordinate system is the
+	# same as the output coordinate system.
+	if (SKY_CTYPE(cooin) == SKY_CTYPE(cooout)) {
+
+	    switch (SKY_CTYPE(cooin)) {
+
+	    case CTYPE_EQUATORIAL:
+		call sk_equatorial (cooin, cooout, ilng, ilat, ipmlng,
+		    ipmlat, px, rv, olng, olat)
+
+	    case CTYPE_ECLIPTIC:
+		if (SKY_EPOCH(cooin) == SKY_EPOCH(cooout)) {
+		    olng = ilng
+		    olat = ilat
+		} else {
+		    call sl_eceq (ilng, ilat, SKY_EPOCH(cooin), olng, olat)
+		    call sl_eqec (olng, olat, SKY_EPOCH(cooout), olng, olat)
+		}
+
+	    default:
+		olng = ilng
+		olat = ilat
+	    }
+
+	    return
+	}
+
+	# Compute proper motions ?
+	if (! IS_INDEFD(ipmlng) && ! IS_INDEFD(ipmlat))
+	    pmflag = YES
+	else
+	    pmflag = NO
+
+	# Cover the remaining cases.
+	switch (SKY_CTYPE(cooin)) {
+
+	# The input system is equatorial.
+	case CTYPE_EQUATORIAL:
+
+	    switch (SKY_RADECSYS(cooin)) {
+
+	    case EQTYPE_FK4, EQTYPE_FK4NOE:
+	        if (pmflag == YES) {
+		    call sl_pm (ilng, ilat, ipmlng, ipmlat, px, rv,
+		        sl_epb (SKY_EPOCH(cooin)), sl_epb (SKY_EPOCH(cooout)),
+			olng, olat)
+	        } else {
+		    olng = ilng
+		    olat = ilat
+	        }
+		if (SKY_RADECSYS(cooin) == EQTYPE_FK4)
+		    call sl_suet (olng, olat, SKY_EQUINOX(cooin), olng, olat)
+		if (SKY_EQUINOX(cooin) != 1950.0d0)
+		    call sl_prcs (1, SKY_EQUINOX(cooin), 1950.0d0, olng, olat) 
+		call sl_adet (olng, olat, 1950.0d0, olng, olat)
+		if (pmflag == YES)
+		    call sl_f45z (olng, olat, sl_epb(SKY_EPOCH(cooout)),
+		        olng, olat)
+		else
+		    call sl_f45z (olng, olat, sl_epb (SKY_EPOCH(cooin)),
+		        olng, olat)
+
+	    case EQTYPE_FK5:
+	        if (pmflag == YES) {
+		    call sl_pm (ilng, ilat, ipmlng, ipmlat, px, rv,
+		        sl_epj (SKY_EPOCH(cooin)), sl_epj(SKY_EPOCH(cooout)),
+			olng, olat)
+	        } else {
+	            olng = ilng
+	            olat = ilat
+		}
+		if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+
+	    case EQTYPE_ICRS:
+	        if (pmflag == YES) {
+		    call sl_pm (ilng, ilat, ipmlng, ipmlat, px, rv,
+		        sl_epj (SKY_EPOCH(cooin)), sl_epj(SKY_EPOCH(cooout)),
+			olng, olat)
+	        } else {
+	            olng = ilng
+	            olat = ilat
+		}
+		if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_hf5z (olng, olat, 2000.0d0, olng, olat, pmr, pmd)
+
+	    case EQTYPE_GAPPT:
+		call sl_amp (ilng, ilat, SKY_EPOCH(cooin), 2000.0d0, olng, olat)
+
+	    }
+
+	    switch (SKY_CTYPE(cooout)) {
+
+	    # The output coordinate system is ecliptic.
+	    case CTYPE_ECLIPTIC:
+		call sl_eqec (olng, olat, SKY_EPOCH(cooout), olng, olat)
+
+	    # The output coordinate system is galactic.
+	    case CTYPE_GALACTIC:
+		call sl_eqga (olng, olat, olng, olat)
+
+	    # The output coordinate system is supergalactic.
+	    case CTYPE_SUPERGALACTIC:
+		call sl_eqga (olng, olat, olng, olat)
+		call sl_gasu (olng, olat, olng, olat)
+
+	    default:
+	        olng = ilng
+	        olat = ilat
+	    }
+
+	# The input coordinate system is ecliptic.
+	case CTYPE_ECLIPTIC:
+
+	    call sl_eceq (ilng, ilat, SKY_EPOCH(cooin), olng, olat)
+	    switch (SKY_CTYPE(cooout)) {
+
+	    # The output coordinate system is equatorial.
+	    case CTYPE_EQUATORIAL:
+
+		switch (SKY_RADECSYS(cooout)) {
+		case EQTYPE_FK4, EQTYPE_FK4NOE:
+		    call sl_f54z (olng, olat, sl_epb(SKY_EPOCH(cooout)),
+		        olng, olat, pmr, pmd)
+		    call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 1950.0d0)
+			call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+		    if (SKY_RADECSYS(cooout) == EQTYPE_FK4)
+		        call sl_adet (olng, olat, SKY_EQUINOX(cooout),
+			    olng, olat)
+
+		case EQTYPE_FK5:
+		    if (SKY_EQUINOX(cooout) != 2000.0d0)
+			call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_ICRS:
+		    #call sl_f5hz (olng, olat, sl_epj(SKY_EPOCH(cooin)),
+		        #olng, olat)
+		    call sl_f5hz (olng, olat, 2000.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 2000.0d0)
+			call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_GAPPT:
+		    call sl_map (olng, olat, 0.0d0, 0.0d0, px, 0.0d0,
+			2000.0d0, SKY_EPOCH(cooout), olng, olat)
+		}
+
+	    # The output coordinate system is galactic.
+	    case CTYPE_GALACTIC:
+		call sl_eqga (olng, olat, olng, olat)
+
+	    # The output system is supergalactic.
+	    case CTYPE_SUPERGALACTIC:
+		call sl_eqga (olng, olat, olng, olat)
+		call sl_gasu (olng, olat, olng, olat)
+
+	    default:
+	        olng = ilng
+	        olat = ilat
+	    }
+
+	# The input coordinate system is galactic.
+	case CTYPE_GALACTIC:
+
+	    switch (SKY_CTYPE(cooout)) {
+
+	    # The output coordinate system is equatorial.
+	    case CTYPE_EQUATORIAL:
+	        call sl_gaeq (ilng, ilat, olng, olat)
+
+		switch (SKY_RADECSYS(cooout)) {
+		case EQTYPE_FK4, EQTYPE_FK4NOE:
+		    call sl_f54z (olng, olat, sl_epb(SKY_EPOCH(cooout)),
+		        olng, olat, pmr, pmd)
+		    call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 1950.0d0)
+			call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+		    if (SKY_RADECSYS(cooout) == EQTYPE_FK4)
+		        call sl_adet (olng, olat, SKY_EQUINOX(cooout),
+			    olng, olat)
+
+		case EQTYPE_FK5:
+		    if (SKY_EQUINOX(cooout) != 2000.0d0)
+			call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_ICRS:
+		    call sl_f5hz (olng, olat, 2000.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 2000.0d0)
+			call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_GAPPT:
+		    call sl_map (olng, olat, 0.0d0, 0.0d0, px, 0.0d0,
+			2000.0d0, SKY_EPOCH(cooout), olng, olat)
+		}
+
+	    # The output coordinate system is ecliptic.
+	    case CTYPE_ECLIPTIC:
+		call sl_gaeq (ilng, ilat, olng, olat)
+		call sl_eqec (olng, olat, SKY_EPOCH(cooout), olng, olat)
+
+	    # The output coordinate system is supergalactic.
+	    case CTYPE_SUPERGALACTIC:
+		call sl_gasu (ilng, ilat, olng, olat)
+
+	    default:
+	        olng = ilng
+	        olat = ilat
+	    }
+
+	# The input coordinates are supergalactic.
+	case CTYPE_SUPERGALACTIC:
+
+	    switch (SKY_CTYPE(cooout)) {
+
+	    case CTYPE_EQUATORIAL:
+		call sl_suga (ilng, ilat, olng, olat)
+
+		switch (SKY_RADECSYS(cooout)) {
+
+		case EQTYPE_FK4:
+		    call sl_gaeq (olng, olat, olng, olat)
+		    call sl_f54z (olng, olat, sl_epb (SKY_EPOCH(cooout)),
+		        olng, olat, pmr, pmd)
+		    call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 1950.0d0)
+			call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+		    call sl_adet (olng, olat, SKY_EQUINOX(cooout), olng, olat)
+
+		case EQTYPE_FK4NOE:
+		    call sl_gaeq (olng, olat, olng, olat)
+		    call sl_f54z (olng, olat, sl_epb (SKY_EPOCH(cooout)),
+		        olng, olat, pmr, pmd)
+		    call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 1950.0d0)
+			call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_FK5:
+		    call sl_gaeq (olng, olat, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 2000.0d0)
+			call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_ICRS:
+		    call sl_gaeq (olng, olat, olng, olat)
+		    call sl_f5hz (olng, olat, 2000.0d0, olng, olat)
+		    if (SKY_EQUINOX(cooout) != 2000.0d0)
+			call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout),
+			    olng, olat) 
+
+		case EQTYPE_GAPPT:
+		    call sl_gaeq (olng, olat, olng, olat)
+		    call sl_map (olng, olat, 0.0d0, 0.0d0, px, 0.0d0,
+			2000.0d0, SKY_EPOCH(cooout), olng, olat)
+		}
+
+	    case CTYPE_ECLIPTIC:
+		call sl_suga (ilng, ilat, olng, olat)
+		call sl_gaeq (olng, olat, olng, olat)
+		call sl_eqec (olng, olat, SKY_EPOCH(cooout), olng, olat)
+
+	    case CTYPE_GALACTIC:
+		call sl_suga (ilng, ilat, olng, olat)
+
+	    default:
+	        olng = ilng
+	        olat = ilat
+	    }
+
+	default:
+	    olng = ilng
+	    olat = ilat
+	}
+end
+
+
+# SK_EQUATORIAL -- Convert / precess equatorial coordinates.
+
+procedure sk_equatorial (cooin, cooout, ilng, ilat, ipmlng, ipmlat,
+	px, rv, olng, olat)
+
+pointer	cooin		#I the input coordinate system structure
+pointer	cooout		#I the output coordinate system structure
+double	ilng		#I the input ra in radians
+double	ilat		#I the input dec in radians
+double	ipmlng		#I the input proper motion in ra in radians
+double	ipmlat		#I the input proper motion in dec in radians
+double	px		#I the input parallax in arcseconds
+double	rv		#I the input radial valocity in km / second
+double	olng		#O the output ra in radians
+double	olat		#O the output dec in radians
+
+int	pmflag
+double	pmr, pmd
+double	sl_epb(), sl_epj()
+
+begin
+	# Check to see whether or not conversion / precession is necessary.
+	if ((SKY_RADECSYS(cooin) == SKY_RADECSYS(cooout)) &&
+	    (SKY_EQUINOX(cooin) == SKY_EQUINOX(cooout)) &&
+	    (SKY_EPOCH(cooin) == SKY_EPOCH(cooout))) {
+	    olng = ilng
+	    olat = ilat
+	    return
+	}
+
+	# Compute proper motions ?
+	if (! IS_INDEFD(ipmlng) && ! IS_INDEFD(ipmlat))
+	    pmflag = YES
+	else
+	    pmflag = NO
+
+	switch (SKY_RADECSYS(cooin)) {
+
+	# The input coordinate system is FK4 with or without the E terms.
+	case EQTYPE_FK4, EQTYPE_FK4NOE:
+
+	    if (pmflag == YES) {
+		call sl_pm (ilng, ilat, ipmlng, ipmlat, px, rv,
+		    sl_epb (SKY_EPOCH(cooin)), sl_epb (SKY_EPOCH(cooout)),
+		    olng, olat)
+	    } else {
+		olng = ilng
+		olat = ilat
+	    }
+	    if (SKY_RADECSYS(cooin) == EQTYPE_FK4)
+		call sl_suet (olng, olat, SKY_EQUINOX(cooin), olng, olat)
+	    if (SKY_EQUINOX(cooin) != 1950.0d0)
+	        call sl_prcs (1, SKY_EQUINOX(cooin), 1950.0d0, olng, olat) 
+	    call sl_adet (olng, olat, 1950.0d0, olng, olat)
+	    if (pmflag == YES)
+	        call sl_f45z (olng, olat, sl_epb (SKY_EPOCH(cooout)),
+		    olng, olat)
+	    else
+	        call sl_f45z (olng, olat, sl_epb (SKY_EPOCH(cooin)),
+		    olng, olat)
+
+	    switch (SKY_RADECSYS(cooout)) {
+
+	    # The output coordinate system is FK4 with and without the E terms.
+	    case EQTYPE_FK4, EQTYPE_FK4NOE:
+		call sl_f54z (olng, olat, sl_epb (SKY_EPOCH(cooout)),
+		    olng, olat, pmr, pmd)
+	        call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		if (SKY_EQUINOX(cooout) != 1950.0d0)
+	            call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout),
+		        olng, olat) 
+		if (SKY_RADECSYS(cooout) == EQTYPE_FK4)
+	            call sl_adet (olng, olat, SKY_EQUINOX(cooout), olng, olat)
+
+	    # The output coordinate system is FK5.
+	    case EQTYPE_FK5:
+		if (SKY_EQUINOX(cooout) != 2000.0d0)
+		    call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout), olng, olat) 
+
+	    # The output coordinate system is ICRS (Hipparcos).
+	    case EQTYPE_ICRS:
+		call sl_f5hz (olng, olat, 2000.0d0, olng, olat)
+		if (SKY_EQUINOX(cooout) != 2000.0d0)
+		    call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout), olng, olat) 
+
+	    # The output coordinate system is geocentric apparent.
+	    case EQTYPE_GAPPT:
+		call sl_map (olng, olat, 0.0d0, 0.0d0, px, 0.0d0, 2000.0d0,
+		    SKY_EPOCH(cooout), olng, olat)
+	    }
+
+	# The input coordinate system is FK5 or geocentric apparent.
+	case EQTYPE_FK5, EQTYPE_GAPPT:
+
+	    if (SKY_RADECSYS(cooin) == EQTYPE_FK5) {
+	        if (pmflag == YES) {
+		    call sl_pm (ilng, ilat, ipmlng, ipmlat, px, rv,
+		        sl_epj (SKY_EPOCH(cooin)), sl_epj (SKY_EPOCH(cooout)),
+			    olng, olat)
+	        } else {
+	            olng = ilng
+	            olat = ilat
+		}
+	    } else
+	        call sl_amp (ilng, ilat, SKY_EPOCH(cooin), 2000.0d0, olng, olat)
+
+	    switch (SKY_RADECSYS(cooout)) {
+
+	    # The output coordinate system is FK4 with or without the E terms.
+	    case EQTYPE_FK4, EQTYPE_FK4NOE:
+	        if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_f54z (olng, olat, sl_epb(SKY_EPOCH(cooout)),
+		    olng, olat, pmr, pmd)
+		call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		if (SKY_EQUINOX(cooout) != 1950.0d0)
+		    call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout), olng, olat) 
+		if (SKY_RADECSYS(cooout) == EQTYPE_FK4)
+	            call sl_adet (olng, olat, SKY_EQUINOX(cooout), olng, olat)
+
+	    # The output coordinate system is FK5.
+	    case EQTYPE_FK5:
+		if (SKY_EQUINOX(cooin) != SKY_EQUINOX(cooout))
+		    call sl_prcs (2, SKY_EQUINOX(cooin), SKY_EQUINOX(cooout),
+		        olng, olat) 
+
+	    # The output coordinate system is ICRS.
+	    case EQTYPE_ICRS:
+	        if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_f5hz (olng, olat, sl_epj(SKY_EPOCH(cooin)), olng, olat)
+	        if (SKY_EQUINOX(cooout) != 2000.0d0)
+		    call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout), olng, olat) 
+
+	    # The output coordinate system is geocentric apparent.
+	    case EQTYPE_GAPPT:
+	        if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_map (olng, olat, 0.0d0, 0.0d0, px, 0.0d0, 2000.0d0,
+		    SKY_EPOCH(cooout), olng, olat)
+	    }
+
+	# The input coordinate system is ICRS.
+	case EQTYPE_ICRS:
+
+	    if (pmflag == YES) {
+		call sl_pm (ilng, ilat, ipmlng, ipmlat, px, rv,
+		    sl_epj (SKY_EPOCH(cooin)), sl_epj (SKY_EPOCH(cooout)),
+		    olng, olat)
+	    } else {
+	        olng = ilng
+	        olat = ilat
+	    }
+
+	    switch (SKY_RADECSYS(cooout)) {
+
+	    # The output coordinate system is FK4 with or without the E terms.
+	    case EQTYPE_FK4, EQTYPE_FK4NOE:
+	        if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_hf5z (olng, olat, 2000.0d0, olng, olat,
+		    pmr, pmd)
+		call sl_f54z (olng, olat, sl_epb(SKY_EPOCH(cooout)), olng, olat,
+		    pmr, pmd)
+		call sl_suet (olng, olat, 1950.0d0, olng, olat)
+		if (SKY_EQUINOX(cooout) != 1950.0d0)
+		    call sl_prcs (1, 1950.0d0, SKY_EQUINOX(cooout), olng, olat) 
+		if (SKY_RADECSYS(cooout) == EQTYPE_FK4)
+	            call sl_adet (olng, olat, SKY_EQUINOX(cooout), olng, olat)
+
+	    # The output coordinate system is FK5.
+	    case EQTYPE_FK5:
+	        if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_hf5z (olng, olat, sl_epj(SKY_EPOCH(cooout)),
+		    olng, olat, pmr, pmd)
+	        if (SKY_EQUINOX(cooout) != 2000.0d0)
+		    call sl_prcs (2, 2000.0d0, SKY_EQUINOX(cooout), olng, olat) 
+
+	    # The output coordinate system is ICRS.
+	    case EQTYPE_ICRS:
+		if (SKY_EQUINOX(cooin) != SKY_EQUINOX(cooout))
+		    call sl_prcs (2, SKY_EQUINOX(cooin), SKY_EQUINOX(cooout),
+		        olng, olat) 
+
+	    # The output coordinate system is geocentric apparent.
+	    case EQTYPE_GAPPT:
+	        if (SKY_EQUINOX(cooin) != 2000.0d0)
+		    call sl_prcs (2, SKY_EQUINOX(cooin), 2000.0d0, olng, olat) 
+		call sl_hf5z (olng, olat, sl_epj(SKY_EPOCH(cooout)),
+		    olng, olat, pmr, pmd)
+		call sl_map (olng, olat, 0.0d0, 0.0d0, px, 0.0d0, 2000.0d0,
+		    SKY_EPOCH(cooout), olng, olat)
+
+	    }
+
+	}
+end
diff --git a/pkg/xtools/skywcs/skwrdstr.x b/pkg/xtools/skywcs/skwrdstr.x
new file mode 100644
index 00000000..a7c6b359
--- /dev/null
+++ b/pkg/xtools/skywcs/skwrdstr.x
@@ -0,0 +1,53 @@
+
+# SK_WRDSTR -- Search a dictionary string for a given string index number.
+# This is the opposite function of strdic(), that returns the index for
+# given string.  The entries in the dictionary string are separated by
+# a delimiter character which is the first character of the dictionary
+# string.  The index of the string found is returned as the function value.
+# Otherwise, if there is no string for that index, a zero is returned.
+
+int procedure sk_wrdstr (index, outstr, maxch, dict)
+
+int     index                   #I the string index
+char    outstr[ARB]             #O the output string as found in dictionary
+int     maxch                   #I the maximum length of output string
+char    dict[ARB]               #I the dictionary string
+
+int     i, len, start, count
+
+int     strlen()
+
+begin
+        # Clear output string
+        outstr[1] = EOS
+
+        # Return if the dictionary is not long enough
+        if (dict[1] == EOS)
+            return (0)
+
+        # Initialize counters
+        count = 1
+        len   = strlen (dict)
+
+        # Search the dictionary string. This loop only terminates
+        # successfully if the index is found. Otherwise the procedure
+        # returns with and error condition.
+        for (start = 2; count < index; start = start + 1) {
+            if (dict[start] == dict[1])
+                count = count + 1
+            if (start == len)
+                return (0)
+        }
+
+        # Extract the output string from the dictionary
+        for (i = start; dict[i] != EOS && dict[i] != dict[1]; i = i + 1) {
+            if (i - start + 1 > maxch)
+                break
+            outstr[i - start + 1] = dict[i]
+        }
+
+        outstr[i - start + 1] = EOS
+
+        # Return index for output string
+        return (count)
+end
diff --git a/pkg/xtools/skywcs/skwrite.x b/pkg/xtools/skywcs/skwrite.x
new file mode 100644
index 00000000..2e779b09
--- /dev/null
+++ b/pkg/xtools/skywcs/skwrite.x
@@ -0,0 +1,510 @@
+include "skywcsdef.h"
+include "skywcs.h"
+
+
+# SK_IIPRINT -- Print a summary of the input image or list coordinate system.
+
+procedure sk_iiprint (label, imagesys, mw, coo)
+
+char	label[ARB]		#I the input label
+char	imagesys[ARB]		#I the input image name and wcs
+pointer	mw			#I pointer to the image wcs
+pointer	coo			#I pointer to the coordinate system structure
+
+begin
+	if (mw == NULL)
+	    call sk_inprint (label, imagesys, SKY_CTYPE(coo),
+	        SKY_RADECSYS(coo), SKY_EQUINOX(coo), SKY_EPOCH(coo))
+	else
+	    call sk_imprint (label, imagesys, SKY_CTYPE(coo), SKY_PLNGAX(coo),
+	        SKY_PLATAX(coo), SKY_WTYPE(coo), SKY_PIXTYPE(coo),
+		SKY_RADECSYS(coo), SKY_EQUINOX(coo), SKY_EPOCH(coo))
+end
+
+
+# SK_IIWRITE -- Write a summary of the input image or list coordinate system
+# to the output file
+
+procedure sk_iiwrite (fd, label, imagesys, mw, coo)
+
+int	fd			#I the output file descriptor
+char	label[ARB]		#I the input label
+char	imagesys[ARB]		#I the input image name and wcs
+pointer	mw			#I pointer to the image wcs
+pointer	coo			#I pointer to the coordinate system structure
+
+begin
+	if (mw == NULL)
+	    call sk_inwrite (fd, label, imagesys, SKY_CTYPE(coo),
+	        SKY_RADECSYS(coo), SKY_EQUINOX(coo), SKY_EPOCH(coo))
+	else
+	    call sk_imwrite (fd, label, imagesys, SKY_CTYPE(coo),
+	        SKY_PLNGAX(coo), SKY_PLATAX(coo), SKY_WTYPE(coo),
+		SKY_PIXTYPE(coo), SKY_RADECSYS(coo), SKY_EQUINOX(coo),
+		SKY_EPOCH(coo))
+end
+
+
+# SK_INPRINT -- Print a summary of the input list coordinate system.
+# This should probably be a call to sk_inwrite with the file descriptor
+# set to STDOUT to avoid duplication of code. There was a reason for
+# having two routines at one point but I can't remember what it was ...
+
+procedure sk_inprint (label, system, ctype, radecsys, equinox, epoch)
+
+char	label[ARB]		#I the input label
+char	system[ARB]		#I the input system
+int	ctype			#I the input coordinate type
+int	radecsys		#I the input equatorial reference system
+double	equinox			#I the input equinox
+double	epoch			#I the input epoch of the observation
+
+pointer	sp, radecstr
+double	sl_epj(), sl_epb()
+int	sk_wrdstr()
+
+begin
+	call smark (sp)
+	call salloc (radecstr, SZ_FNAME, TY_CHAR)
+
+	switch (ctype) {
+
+	case CTYPE_EQUATORIAL:
+	    if (sk_wrdstr (radecsys, Memc[radecstr], SZ_FNAME,
+	        EQTYPE_LIST) <= 0)
+	        call strcpy ("FK5", Memc[radecstr], SZ_FNAME)
+	    call strupr (Memc[radecstr])
+	    call printf ("%s: %s  Coordinates: equatorial %s\n")
+	        call pargstr (label)
+		call pargstr (system)
+		call pargstr (Memc[radecstr])
+	    switch (radecsys) {
+	    case EQTYPE_GAPPT:
+		call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		    call pargd (epoch)
+		    if (IS_INDEFD(epoch)) {
+			call pargd (INDEFD)
+			call pargd (INDEFD)
+		    } else {
+		        call pargd (sl_epj (epoch))
+		        call pargd (sl_epb (epoch))
+		    }
+	    case EQTYPE_FK5, EQTYPE_ICRS:
+		call printf ("    Equinox: J%0.3f Epoch: J%0.8f MJD: %0.5f\n")
+		    call pargd (equinox)
+		    call pargd (sl_epj(epoch))
+		    call pargd (epoch)
+	    default:
+		call printf ("    Equinox: B%0.3f Epoch: B%0.8f MJD: %0.5f\n")
+		    call pargd (equinox)
+		    call pargd (sl_epb(epoch))
+		    call pargd (epoch)
+	    }
+
+	case CTYPE_ECLIPTIC:
+	    call printf ("%s: %s  Coordinates: ecliptic\n")
+		call pargstr (label)
+		call pargstr (system)
+	    call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		if (IS_INDEFD(epoch)) {
+		    call pargd (INDEFD)
+		    call pargd (INDEFD)
+		} else {
+		    call pargd (sl_epj(epoch))
+		    call pargd (sl_epb(epoch))
+		}
+
+	case CTYPE_GALACTIC:
+	    call printf ("%s: %s  Coordinates: galactic\n")
+		call pargstr (label)
+		call pargstr (system)
+	    call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj (epoch))
+		call pargd (sl_epb (epoch))
+
+	case CTYPE_SUPERGALACTIC:
+	    call printf ("%s: %s  Coordinates: supergalactic\n")
+		call pargstr (label)
+		call pargstr (system)
+	    call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj (epoch))
+		call pargd (sl_epb (epoch))
+
+	}
+
+	call sfree (sp)
+end
+
+
+# SK_INWRITE -- Write a summary of the input coordinate system.
+
+procedure sk_inwrite (fd, label, system, ctype, radecsys, equinox, epoch)
+
+int	fd			#I the output file descriptor
+char	label[ARB]		#I the input label
+char	system[ARB]		#I the input system
+int	ctype			#I the input coordinate type
+int	radecsys		#I the input equatorial reference system
+double	equinox			#I the input equinox
+double	epoch			#I the input epoch of the observation
+
+pointer	sp, radecstr
+double	sl_epj(), sl_epb()
+int	sk_wrdstr()
+
+begin
+	call smark (sp)
+	call salloc (radecstr, SZ_FNAME, TY_CHAR)
+
+	switch (ctype) {
+
+	case CTYPE_EQUATORIAL:
+	    if (sk_wrdstr (radecsys, Memc[radecstr], SZ_FNAME,
+	        EQTYPE_LIST) <= 0)
+	        call strcpy ("FK5", Memc[radecstr], SZ_FNAME)
+	    call strupr (Memc[radecstr])
+	    call fprintf (fd, "# %s: %s  Coordinates: equatorial %s\n")
+	        call pargstr (label)
+	        call pargstr (system)
+		call pargstr (Memc[radecstr])
+	    switch (radecsys) {
+	    case EQTYPE_GAPPT:
+		call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		    call pargd (epoch)
+		    if (IS_INDEFD(epoch)) {
+		        call pargd (INDEFD)
+		        call pargd (INDEFD)
+		    } else {
+		        call pargd (sl_epj(epoch))
+		        call pargd (sl_epb(epoch))
+		    }
+	    case EQTYPE_FK5, EQTYPE_ICRS:
+		call fprintf (fd,
+		    "#     Equinox: J%0.3f Epoch: J%0.8f MJD: %0.5f\n")
+		    call pargd (equinox)
+		    call pargd (sl_epj(epoch))
+		    call pargd (epoch)
+	    default:
+		call fprintf (fd,
+		    "#     Equinox: B%0.3f Epoch: B%0.8f MJD: %0.5f\n")
+		    call pargd (equinox)
+		    call pargd (sl_epb(epoch))
+		    call pargd (epoch)
+	    }
+
+	case CTYPE_ECLIPTIC:
+	    call fprintf (fd, "# %s: %s  Coordinates: ecliptic\n")
+		call pargstr (label)
+		call pargstr (system)
+	    call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		if (IS_INDEFD(epoch)) {
+		    call pargd (INDEFD)
+		    call pargd (INDEFD)
+		} else {
+		    call pargd (sl_epj(epoch))
+		    call pargd (sl_epb(epoch))
+		}
+
+	case CTYPE_GALACTIC:
+	    call fprintf (fd, "# %s: %s  Coordinates: galactic\n")
+		call pargstr (label)
+		call pargstr (system)
+	    call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj(epoch))
+		call pargd (sl_epb(epoch))
+
+	case CTYPE_SUPERGALACTIC:
+	    call fprintf (fd, "# %s: %s  Coordinates: supergalactic\n")
+		call pargstr (label)
+		call pargstr (system)
+	    call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj(epoch))
+		call pargd (sl_epb(epoch))
+
+	}
+
+	call sfree (sp)
+end
+
+
+# SK_IMPRINT -- Print a summary of the input image coordinate system.
+# This should probably be a call to sk_imwrite with the file descriptor
+# set to STDOUT to avoid duplication of code. There was a reason for
+# having two routines at one point but I can't remember what it was ...
+
+procedure sk_imprint (label, imagesys, ctype, lngax, latax, wtype, ptype,
+	radecsys, equinox, epoch)
+
+char	label[ARB]		#I input label
+char	imagesys[ARB]		#I the input image name and system
+int	ctype			#I the image coordinate type
+int	lngax			#I the image ra/glon/elon axis
+int	latax			#I the image dec/glat/elat axis
+int	wtype			#I the image projection type
+int	ptype			#I the image image wcs type
+int	radecsys		#I the image equatorial reference system
+double	equinox			#I the image equinox
+double	epoch			#I the image epoch of the observation
+
+pointer	sp, imname, projstr, wcsstr, radecstr
+double	sl_epj(), sl_epb()
+int	sk_wrdstr()
+
+begin
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+	call salloc (projstr, SZ_FNAME, TY_CHAR)
+	call salloc (wcsstr, SZ_FNAME, TY_CHAR)
+	call salloc (radecstr, SZ_FNAME, TY_CHAR)
+
+	call sscan (imagesys)
+	    call gargwrd (Memc[imname], SZ_FNAME)
+	if (sk_wrdstr (wtype, Memc[projstr], SZ_FNAME, WTYPE_LIST) <= 0)
+	    call strcpy ("linear", Memc[projstr], SZ_FNAME)
+	call strupr (Memc[projstr])
+	if (sk_wrdstr (ptype, Memc[wcsstr], SZ_FNAME, PIXTYPE_LIST) <= 0)
+	    call strcpy ("world", Memc[wcsstr], SZ_FNAME)
+	call strlwr (Memc[wcsstr])
+
+	switch (ctype) {
+
+	case CTYPE_EQUATORIAL:
+	    if (sk_wrdstr (radecsys, Memc[radecstr], SZ_FNAME,
+	        EQTYPE_LIST) <= 0)
+	        call strcpy ("FK5", Memc[radecstr], SZ_FNAME)
+	    call strupr (Memc[radecstr])
+	    call printf (
+	    "%s: %s %s  Projection: %s  Ra/Dec axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    switch (radecsys) {
+	    case EQTYPE_GAPPT:
+	        call printf ("    Coordinates: equatorial %s\n")
+		    call pargstr (Memc[radecstr])
+		call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		    call pargd (epoch)
+		    if (IS_INDEFD(epoch)) {
+			call pargd (INDEFD)
+			call pargd (INDEFD)
+		    } else {
+		        call pargd (sl_epj(epoch))
+		        call pargd (sl_epb(epoch))
+		    }
+	    case EQTYPE_FK5, EQTYPE_ICRS:
+	        call printf ("    Coordinates: equatorial %s Equinox: J%0.3f\n")
+		    call pargstr (Memc[radecstr])
+		    call pargd (equinox)
+	        call printf ("    Epoch: J%0.8f MJD: %0.5f\n")
+		    call pargd (sl_epj (epoch))
+		    call pargd (epoch)
+	    default:
+	        call printf ("    Coordinates: equatorial %s Equinox: B%0.3f\n")
+		    call pargstr (Memc[radecstr])
+		    call pargd (equinox)
+	        call printf ("    Epoch: B%0.8f MJD: %0.5f\n")
+		    call pargd (sl_epb (epoch))
+		    call pargd (epoch)
+	    }
+
+	case CTYPE_ECLIPTIC:
+	    call printf (
+	    "%s: %s %s  Projection: %s  Elong/Elat axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    call printf ("    Coordinates: ecliptic\n")
+	    call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		if (IS_INDEFD(epoch)) {
+		    call pargd (INDEFD)
+		    call pargd (INDEFD)
+		} else {
+		    call pargd (sl_epj(epoch))
+		    call pargd (sl_epb(epoch))
+		}
+
+	case CTYPE_GALACTIC:
+	    call printf (
+	    "%s: %s %s  Projection: %s  Glong/Glat axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    call printf ("    Coordinates: galactic\n")
+	    call printf ("    MJD: %0.5f  Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj (epoch))
+		call pargd (sl_epb (epoch))
+
+	case CTYPE_SUPERGALACTIC:
+	    call printf (
+	    "%s: %s %s  Projection: %s  Slong/Slat axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    call printf ("    Coordinates: supergalactic\n")
+	    call printf ("    MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj (epoch))
+		call pargd (sl_epb (epoch))
+	}
+
+	call sfree (sp)
+end
+
+
+# SK_IMWRITE -- Write a summary of the image coordinate system to the
+# output file.
+
+procedure sk_imwrite (fd, label, imagesys, ctype, lngax, latax, wtype, ptype,
+	radecsys, equinox, epoch)
+
+int	fd			#I the output file descriptor
+char	label[ARB]		#I input label
+char	imagesys[ARB]		#I the input image name and wcs
+int	ctype			#I the image coordinate type
+int	lngax			#I the image ra/glon/elon axis
+int	latax			#I the image dec/glat/elat axis
+int	wtype			#I the image projection type
+int	ptype			#I the image image wcs type
+int	radecsys		#I the image equatorial reference system
+double	equinox			#I the image equinox
+double	epoch			#I the image epoch of the observation
+
+pointer	sp, imname, projstr, wcsstr, radecstr
+double	sl_epj(), sl_epb()
+int	sk_wrdstr()
+
+begin
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+	call salloc (projstr, SZ_FNAME, TY_CHAR)
+	call salloc (wcsstr, SZ_FNAME, TY_CHAR)
+	call salloc (radecstr, SZ_FNAME, TY_CHAR)
+
+	call sscan (imagesys)
+	    call gargwrd (Memc[imname], SZ_FNAME)
+	if (sk_wrdstr (wtype, Memc[projstr], SZ_FNAME, WTYPE_LIST) <= 0)
+	    call strcpy ("linear", Memc[projstr], SZ_FNAME)
+	call strupr (Memc[projstr])
+	if (sk_wrdstr (ptype, Memc[wcsstr], SZ_FNAME, PIXTYPE_LIST) <= 0)
+	    call strcpy ("world", Memc[wcsstr], SZ_FNAME)
+	call strlwr (Memc[wcsstr])
+
+	switch (ctype) {
+
+	case CTYPE_EQUATORIAL:
+	    if (sk_wrdstr (radecsys, Memc[radecstr], SZ_FNAME,
+	        EQTYPE_LIST) <= 0)
+	        call strcpy ("FK5", Memc[radecstr], SZ_FNAME)
+	    call strupr (Memc[radecstr])
+	    call fprintf (fd,
+	    "# %s: %s %s  Projection: %s  Ra/Dec axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    switch (radecsys) {
+	    case EQTYPE_GAPPT:
+	        call fprintf (fd, "#     Coordinates: equatorial %s\n")
+		    call pargstr (Memc[radecstr])
+		call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		    call pargd (epoch)
+		    if (IS_INDEFD(epoch)) {
+		        call pargd (INDEFD)
+		        call pargd (INDEFD)
+		    } else {
+		        call pargd (sl_epj(epoch))
+		        call pargd (sl_epb(epoch))
+		    }
+	    case EQTYPE_FK5, EQTYPE_ICRS:
+	        call fprintf (fd,
+		    "#     Coordinates: equatorial %s Equinox: J%0.3f\n")
+		    call pargstr (Memc[radecstr])
+		    call pargd (equinox)
+	        call fprintf (fd, "#     Epoch: J%0.8f MJD: %0.5f\n")
+		    call pargd (sl_epj(epoch))
+		    call pargd (epoch)
+	    default:
+	        call fprintf (fd,
+		    "#     Coordinates: equatorial %s Equinox: B%0.3f\n")
+		    call pargstr (Memc[radecstr])
+		    call pargd (equinox)
+	        call fprintf (fd, "#     Epoch: B%0.8f MJD: %0.5f\n")
+		    call pargd (sl_epb (epoch))
+		    call pargd (epoch)
+	    }
+
+	case CTYPE_ECLIPTIC:
+	    call fprintf (fd,
+	    "# %s: %s %s  Projection: %s  Elong/Elat axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    call fprintf (fd, "#     Coordinates: ecliptic\n")
+	    call fprintf (fd, "#     MJD: %0.5f  Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		if (IS_INDEFD(epoch)) {
+		    call pargd (INDEFD)
+		    call pargd (INDEFD)
+		} else {
+		    call pargd (sl_epj(epoch))
+		    call pargd (sl_epb(epoch))
+		}
+
+	case CTYPE_GALACTIC:
+	    call fprintf (fd,
+	    "# %s: %s %s  Projection: %s  Glong/Glat axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    call fprintf (fd, "#     Coordinates: galactic\n")
+	    call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj(epoch))
+		call pargd (sl_epb(epoch))
+
+	case CTYPE_SUPERGALACTIC:
+	    call fprintf (fd,
+	    "# %s: %s %s  Projection: %s  Slong/Slat axes: %d/%d\n")
+		call pargstr (label)
+		call pargstr (Memc[imname])
+		call pargstr (Memc[wcsstr])
+	        call pargstr (Memc[projstr])
+	        call pargi (lngax)
+	        call pargi (latax)
+	    call fprintf (fd, "#     Coordinates: supergalactic\n")
+	    call fprintf (fd, "#     MJD: %0.5f Epoch: J%0.8f B%0.8f\n")
+		call pargd (epoch)
+		call pargd (sl_epj(epoch))
+		call pargd (sl_epb(epoch))
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/xtools/skywcs/skywcs.h b/pkg/xtools/skywcs/skywcs.h
new file mode 100644
index 00000000..85b664c0
--- /dev/null
+++ b/pkg/xtools/skywcs/skywcs.h
@@ -0,0 +1,133 @@
+# Public definitions file for the SKYWCS library.
+
+# Define the SKYWCS library parameters.
+
+define	S_VXOFF		1
+define	S_VYOFF		2
+define	S_VXSTEP	3
+define	S_VYSTEP	4
+define	S_EQUINOX	5
+define	S_EPOCH		6
+define	S_CTYPE		7
+define	S_RADECSYS	8
+define	S_WTYPE		9
+define	S_PLNGAX	10
+define	S_PLATAX	11
+define	S_XLAX		12
+define	S_YLAX		13
+define	S_PIXTYPE	14
+define	S_NLNGAX	15
+define	S_NLATAX	16
+define	S_NLNGUNITS	17
+define	S_NLATUNITS	18
+define	S_COOSYSTEM	19
+define	S_STATUS	20
+
+# Define the list of supported fundamental coordinate systems.
+
+define	FTYPE_LIST	"|fk4|noefk4|fk5|icrs|apparent|ecliptic|galactic|\
+supergalactic|"
+
+define	FTYPE_FK4		1
+define	FTYPE_FK4NOE		2
+define	FTYPE_FK5		3
+define	FTYPE_ICRS		4
+define	FTYPE_GAPPT		5
+define	FTYPE_ECLIPTIC		6
+define	FTYPE_GALACTIC		7
+define	FTYPE_SUPERGALACTIC	8
+
+# Define the list of supported coordinate systems.
+
+define	CTYPE_LIST	"|equatorial|ecliptic|galactic|supergalactic|"
+
+define	CTYPE_EQUATORIAL	1
+define	CTYPE_ECLIPTIC		2
+define	CTYPE_GALACTIC		3
+define	CTYPE_SUPERGALACTIC	4
+
+# Define the supported equatoral reference systems.
+
+define  EQTYPE_LIST     "|fk4|fk4-no-e|fk5|icrs|gappt|"
+
+define	EQTYPE_FK4	1
+define	EQTYPE_FK4NOE	2
+define	EQTYPE_FK5	3
+define	EQTYPE_ICRS	4
+define	EQTYPE_GAPPT	5
+
+# Define the input coordinate file longitude latitude units.
+
+define	SKY_LNG_UNITLIST	"|degrees|radians|hours|"
+define	SKY_LAT_UNITLIST	"|degrees|radians|"
+
+define	SKY_DEGREES	1
+define	SKY_RADIANS	2
+define	SKY_HOURS	3
+
+# Define the list of supported image sky projection types.
+
+define  WTYPE_LIST      "|lin|azp|tan|sin|stg|arc|zpn|zea|air|cyp|car|\
+mer|cea|cop|cod|coe|coo|bon|pco|gls|par|ait|mol|csc|qsc|tsc|tnx|zpx|tpv|"
+
+define  PTYPE_LIST      "|z|z|z|z|z|z|z|z|z|c|c|c|c|n|n|n|n|c|c|c|c|c|c|c|c|c|\
+x|x|z|"
+
+define  WTYPE_LIN       1
+define  WTYPE_AZP       2
+define  WTYPE_TAN       3
+define  WTYPE_SIN       4
+define  WTYPE_STG       5
+define  WTYPE_ARC       6
+define  WTYPE_ZPN       7
+define  WTYPE_ZEA       8
+define  WTYPE_AIR       9
+define  WTYPE_CYP       10
+define  WTYPE_CAR       11
+define  WTYPE_MER       12
+define  WTYPE_CEA       13
+define  WTYPE_COP       14
+define  WTYPE_COD       15
+define  WTYPE_COE       16
+define  WTYPE_COO       17
+define  WTYPE_BON       18
+define  WTYPE_PCO       19
+define  WTYPE_GLS       20
+define  WTYPE_PAR       21
+define  WTYPE_AIT       22
+define  WTYPE_MOL       23
+define  WTYPE_CSC       24
+define  WTYPE_QSC       25
+define  WTYPE_TSC       26
+define  WTYPE_TNX       27
+define  WTYPE_ZPX       28
+define  WTYPE_TPV       29
+
+define  PTYPE_NAMES     "|z|c|n|x|"
+
+define  PTYPE_ZEN       1
+define  PTYPE_CYL       2
+define  PTYPE_CON       3
+define  PTYPE_EXP       4
+
+# Define the supported image axis types.
+
+define  AXTYPE_LIST     "|ra|dec|glon|glat|elon|elat|slon|slat|"
+
+define  AXTYPE_RA       1
+define  AXTYPE_DEC      2
+define  AXTYPE_GLON     3
+define  AXTYPE_GLAT     4
+define  AXTYPE_ELON     5
+define  AXTYPE_ELAT     6
+define  AXTYPE_SLON     7
+define  AXTYPE_SLAT     8
+
+# Define the supported image pixel coordinate systems.
+
+define	PIXTYPE_LIST	"|logical|tv|physical|world|"
+
+define  PIXTYPE_LOGICAL		1
+define  PIXTYPE_TV		2
+define  PIXTYPE_PHYSICAL	3
+define  PIXTYPE_WORLD		4
diff --git a/pkg/xtools/skywcs/skywcsdef.h b/pkg/xtools/skywcs/skywcsdef.h
new file mode 100644
index 00000000..433247bd
--- /dev/null
+++ b/pkg/xtools/skywcs/skywcsdef.h
@@ -0,0 +1,24 @@
+# The SKYWCS library structure.
+
+define	LEN_SKYCOOSTRUCT	(30 + SZ_FNAME + 1)
+
+define	SKY_VXOFF	Memd[P2D($1)]	    # logical ra/longitude offset
+define	SKY_VYOFF	Memd[P2D($1+2)]	    # logical dec/tatitude offset
+define	SKY_VXSTEP	Memd[P2D($1+4)]	    # logical ra/longitude stepsize
+define	SKY_VYSTEP	Memd[P2D($1+6)]	    # logical dec/latitude stepsize
+define	SKY_EQUINOX	Memd[P2D($1+8)]	    # equinox of ra/dec system (B or J)
+define	SKY_EPOCH	Memd[P2D($1+10)]    # epoch of observation (MJD)
+define	SKY_CTYPE	Memi[$1+12]	    # celestial coordinate system code
+define	SKY_RADECSYS	Memi[$1+13]	    # ra/dec system code
+define	SKY_WTYPE	Memi[$1+14]	    # sky projection function code
+define	SKY_PLNGAX	Memi[$1+15]	    # physical ra/longitude axis
+define	SKY_PLATAX	Memi[$1+16]	    # physical dec/latitude axis
+define	SKY_XLAX	Memi[$1+17]	    # logical ra/longitude axis
+define	SKY_YLAX	Memi[$1+18]	    # logical dec/latitude axis
+define	SKY_PIXTYPE	Memi[$1+19]	    # iraf wcs system code
+define	SKY_NLNGAX	Memi[$1+20]	    # length of ra/longitude axis
+define	SKY_NLATAX	Memi[$1+21]	    # length of dec/latitude axis
+define	SKY_NLNGUNITS	Memi[$1+22]	    # the native ra/longitude units
+define	SKY_NLATUNITS	Memi[$1+23]	    # the native dec/latitude units
+define	SKY_STATUS	Memi[$1+24]	    # the status (OK or ERR)
+define	SKY_COOSYSTEM	Memc[P2C($1+25)]    # the coordinate system name