[fitswcs] I've got this FITS WCS matrix

[Doug Mink]

Malcolm Currie asked:
> I've been given some Gemini NIRI data.  The WCS headers are listed
> below.
> 
> CTYPE1  = 'RA---TAN'           / the coordinate type for the first axis
> CRPIX1  =     1281.62863446163 / x-coordinate of reference pixel
> CRVAL1  =    0.698405727933413 / first axis value at ref pixel
> CTYPE2  = 'DEC--TAN'           / the coordinate type for the second axis
> CRPIX2  =     433.151678256953 / y-coordinate of reference pixel
> CRVAL2  =     35.8161778878783 / second axis value at ref pixel
> CD1_1   = 6.02715759044019E-06 / partial of first axis coord w.r.t. x
> CD1_2   = 1.27956486107373E-08 / partial of first axis coord w.r.t. y
> CD2_1   = -1.6056794631393E-08 / partial of second axis coord w.r.t. x
> CD2_2   = -5.96793770760736E-06 / partial of second axis coord w.r.t. y
> 
> Now for the purpose of registration in a pipeline I need to know the
> rotation from this CD matrix to determine the approximate locations of
> the source after dithering.  Using Equation 191 from WCS Paper II, I
> find that rho_a is 179.84 degrees and rho_b is 0.12 degrees to two
> decimal places.  What should one do in such cases?

Bill Joye and I wrestled with some rotation problems for a month several years
ago when ds9 needed to get a reliable rotation angle from any CD matrix.  I wrote
an 80-line subroutine in WCSTools, wcsrot() in libwcs/wcs.c, which has worked
in all of the situations we have encountered so far.  I think that most of the
variable names are clear.  Code for the entire package is at

http://tdc-www.harvard.edu/software/wcstools/wcstools.tar.gz

-Doug Mink

/* Compute image rotation */

void
wcsrotset (wcs)

struct WorldCoor *wcs;  /* World coordinate system structure */
{
     int off;
     double cra, cdec, xc, xn, xe, yc, yn, ye;

     /* If image is one-dimensional, leave rotation angle alone */
     if (wcs->nxpix < 1.5 || wcs->nypix < 1.5) {
         wcs->imrot = wcs->rot;
         wcs->pa_north = wcs->rot + 90.0;
         wcs->pa_east = wcs->rot + 180.0;
         return;
         }

     /* Do not try anything if image is LINEAR (not Cartesian projection) */
     if (wcs->syswcs == WCS_LINEAR)
         return;

     wcs->xinc = fabs (wcs->xinc);
     wcs->yinc = fabs (wcs->yinc);

     /* Compute position angles of North and East in image */
     xc = wcs->xrefpix;
     yc = wcs->yrefpix;
     pix2wcs (wcs, xc, yc, &cra, &cdec);
     if (wcs->coorflip) {
         wcs2pix (wcs, cra+wcs->yinc, cdec, &xe, &ye, &off);
         wcs2pix (wcs, cra, cdec+wcs->xinc, &xn, &yn, &off);
         }
     else {
         wcs2pix (wcs, cra+wcs->xinc, cdec, &xe, &ye, &off);
         wcs2pix (wcs, cra, cdec+wcs->yinc, &xn, &yn, &off);
         }
     wcs->pa_north = raddeg (atan2 (yn-yc, xn-xc));
     if (wcs->pa_north < -90.0)
         wcs->pa_north = wcs->pa_north + 360.0;
     wcs->pa_east = raddeg (atan2 (ye-yc, xe-xc));
     if (wcs->pa_east < -90.0)
         wcs->pa_east = wcs->pa_east + 360.0;

     /* Compute image rotation angle from North */
     if (wcs->pa_north < -90.0)
         wcs->imrot = 270.0 + wcs->pa_north;
     else
         wcs->imrot = wcs->pa_north - 90.0;

     /* Compute CROTA */
     if (wcs->coorflip) {
         wcs->rot = wcs->imrot + 90.0;
         if (wcs->rot > 180)
             wcs->rot = wcs->rot - 360.0;
         }
     else
         wcs->rot = wcs->imrot;

     /* Set image mirror flag based on axis orientation */
     wcs->imflip = 0;
     if (wcs->pa_east - wcs->pa_north < -80.0 &&
         wcs->pa_east - wcs->pa_north > -100.0)
         wcs->imflip = 1;
     if (wcs->pa_east - wcs->pa_north < 280.0 &&
         wcs->pa_east - wcs->pa_north > 260.0)
         wcs->imflip = 1;
     if (wcs->pa_north - wcs->pa_east > 80.0 &&
         wcs->pa_north - wcs->pa_east < 100.0)
         wcs->imflip = 1;
     if (wcs->coorflip) {
         if (wcs->imflip)
             wcs->yinc = -wcs->yinc;
         }
     else {
         if (!wcs->imflip)
             wcs->xinc = -wcs->xinc;
         }

     return;
}