diff options
Diffstat (limited to 'pkg/xtools/skywcs/doc/README')
| -rw-r--r-- | pkg/xtools/skywcs/doc/README | 301 |
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) |