Repatch (from linux) of OSX IRAF

1 files changed, 301 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)