Changeset 10603

Timestamp:

Dec 9, 2006, 4:06:47 PM (20 years ago)

Author:

magnier

Message:

moving WCS to modules, adding high-order terms

Location:

trunk/psModules

Files:

• src/astrom/pmAstrometryDistortion.c (modified) (5 diffs)
• src/astrom/pmAstrometryWCS.c (modified) (15 diffs)
• src/astrom/pmAstrometryWCS.h (modified) (2 diffs)
• src/camera/pmFPA.c (modified) (2 diffs)
• src/camera/pmFPA.h (modified) (2 diffs)
• src/psmodules.h (modified) (1 diff)
• test/astrom/Makefile.am (modified) (1 diff)

trunk/psModules/src/astrom/pmAstrometryDistortion.c

@@ -7 +7 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-10-24 22:55:04 $
+*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-12-10 02:06:47 $
 *
 *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -181 +181 @@
     // the order of the gradient fits need to be 1 less than the distortion term
     // determine the gradient order(s) from the fpa->toTP structure
-    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->x->nX-1, fpa->toTangentPlane->x->nY);
-    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->x->nX,   fpa->toTangentPlane->x->nY-1);
+    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTPA->x->nX-1, fpa->toTPA->x->nY);
+    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTPA->x->nX,   fpa->toTPA->x->nY-1);
 
     // set masks based on fpa->toTP
-    for (int i = 0; i < fpa->toTangentPlane->x->nX; i++) {
-        for (int j = 0; j < fpa->toTangentPlane->x->nY; j++) {
-            if (fpa->toTangentPlane->x->mask[i][j][0][0]) {
+    for (int i = 0; i < fpa->toTPA->x->nX; i++) {
+        for (int j = 0; j < fpa->toTPA->x->nY; j++) {
+            if (fpa->toTPA->x->mask[i][j][0][0]) {
                 localX->mask[i-1][j] = 1;
                 localY->mask[i][j-1] = 1;
@@ -203 +203 @@
 
     // update fpa->toTP distortion terms
-    fpa->toTangentPlane->x->coeff[0][0][0][0] = 0;
-    for (int i = 1; i < fpa->toTangentPlane->x->nX; i++) {
-        if (!fpa->toTangentPlane->x->mask[i][0][0][0]) {
-            fpa->toTangentPlane->x->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
-        }
-    }
-    for (int j = 1; j < fpa->toTangentPlane->x->nY; j++) {
-        if (!fpa->toTangentPlane->x->mask[0][j][0][0]) {
-            fpa->toTangentPlane->x->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
-        }
-    }
-    for (int i = 1; i < fpa->toTangentPlane->x->nX; i++) {
-        for (int j = 1; j < fpa->toTangentPlane->x->nY; j++) {
-            if (!fpa->toTangentPlane->x->mask[i][j][0][0]) {
-                fpa->toTangentPlane->x->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
+    fpa->toTPA->x->coeff[0][0][0][0] = 0;
+    for (int i = 1; i < fpa->toTPA->x->nX; i++) {
+        if (!fpa->toTPA->x->mask[i][0][0][0]) {
+            fpa->toTPA->x->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
+        }
+    }
+    for (int j = 1; j < fpa->toTPA->x->nY; j++) {
+        if (!fpa->toTPA->x->mask[0][j][0][0]) {
+            fpa->toTPA->x->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
+        }
+    }
+    for (int i = 1; i < fpa->toTPA->x->nX; i++) {
+        for (int j = 1; j < fpa->toTPA->x->nY; j++) {
+            if (!fpa->toTPA->x->mask[i][j][0][0]) {
+                fpa->toTPA->x->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
             }
         }
@@ -226 +226 @@
     // the order of the gradient fits need to be 1 less than the distortion term
     // determine the gradient order(s) from the fpa->toTP structure
-    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->y->nX-1, fpa->toTangentPlane->y->nY);
-    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->y->nX,   fpa->toTangentPlane->y->nY-1);
+    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTPA->y->nX-1, fpa->toTPA->y->nY);
+    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTPA->y->nX,   fpa->toTPA->y->nY-1);
 
     // set masks based on fpa->toTP
-    for (int i = 0; i < fpa->toTangentPlane->y->nX; i++) {
-        for (int j = 0; j < fpa->toTangentPlane->y->nY; j++) {
-            if (fpa->toTangentPlane->y->mask[i][j][0][0]) {
+    for (int i = 0; i < fpa->toTPA->y->nX; i++) {
+        for (int j = 0; j < fpa->toTPA->y->nY; j++) {
+            if (fpa->toTPA->y->mask[i][j][0][0]) {
                 localX->mask[i-1][j] = 1;
                 localY->mask[i][j-1] = 1;
@@ -259 +259 @@
 
     // update fpa->toTP distortion terms
-    fpa->toTangentPlane->y->coeff[0][0][0][0] = 0;
-    for (int i = 1; i < fpa->toTangentPlane->y->nX; i++) {
-        if (!fpa->toTangentPlane->y->mask[i][0][0][0]) {
-            fpa->toTangentPlane->y->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
-        }
-    }
-    for (int j = 1; j < fpa->toTangentPlane->y->nY; j++) {
-        if (!fpa->toTangentPlane->y->mask[0][j][0][0]) {
-            fpa->toTangentPlane->y->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
-        }
-    }
-    for (int i = 1; i < fpa->toTangentPlane->y->nX; i++) {
-        for (int j = 1; j < fpa->toTangentPlane->y->nY; j++) {
-            if (!fpa->toTangentPlane->y->mask[i][j][0][0]) {
-                fpa->toTangentPlane->y->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
+    fpa->toTPA->y->coeff[0][0][0][0] = 0;
+    for (int i = 1; i < fpa->toTPA->y->nX; i++) {
+        if (!fpa->toTPA->y->mask[i][0][0][0]) {
+            fpa->toTPA->y->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
+        }
+    }
+    for (int j = 1; j < fpa->toTPA->y->nY; j++) {
+        if (!fpa->toTPA->y->mask[0][j][0][0]) {
+            fpa->toTPA->y->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
+        }
+    }
+    for (int i = 1; i < fpa->toTPA->y->nX; i++) {
+        for (int j = 1; j < fpa->toTPA->y->nY; j++) {
+            if (!fpa->toTPA->y->mask[i][j][0][0]) {
+                fpa->toTPA->y->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
             }
         }

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -7 +7 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-12-09 21:06:44 $
+*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-12-10 02:06:47 $
 *
 *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -22 +22 @@
 
 // interpret header WCS (only handles traditional WCS for the moment)
-// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
+// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
 bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic)
 {
 
-    psProjectionType type;
+    # if (0)
+        psProjectionType type;
     bool status, pcKeys, cdKeys;
     float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
@@ -55 +56 @@
 
     // what type of WCS keywords are available?
+    // XXX add check for CROTA2
+    int fitOrder = psMetadataLookupS32 (&isPoly, header, "NPLYTERM");
     psMetadataLookupF32 (&pcKeys, header, "PC001001");
     psMetadataLookupF32 (&cdKeys, header, "CD1_1");
 
     if (cdKeys && pcKeys) {
-        psLogMsg ("psastro", 2, "warning: both CDi_j and PC00i00j defined in headers, using CDi_j terms\n");
+        // XXX make this an option
+        psLogMsg ("psastro", 2, "warning: both CDi_j and PC00i00j defined in headers, using PC00i00j terms\n");
     }
     if (!cdKeys && !pcKeys) {
@@ -66 +70 @@
         return false;
     }
+    if (isPoly && !pcKeys) {
+        psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
+        return false;
+        if (fitOrder = 0)
+            fitOrder = 1;
+        if ((fitOrder > 3) || (fitOrder < 1)) {
+            psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
+            return false;
+        }
+    } else {
+        fitOrder = 1;
+    }
+
+    // construct a transformation from X,Y in pixels to L,M in pixels
+    // NOTE that the WCS keywords convert X,Y to degrees first (using cdelt1,2)
+    // and then define a transformation from degrees to degrees
+    psPlaneTransform *wcsTrans = psPlaneTransformAlloc (fitOrder, fitOrder);
 
     crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");
@@ -78 +99 @@
 
         // test the CROTAi varient:
+        // XXX double check lambda..
         double rotate = psMetadataLookupF32 (&status, header, "CROTA2");
         if (status) {
             double Lambda = cdelt2 / cdelt1;
-            pc1_1 =  cos(rotate*RAD_DEG);
-            pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
-            pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
-            pc2_2 =  cos(rotate*RAD_DEG);
+            wcsTrans->x->coeff[1][0] = +cos(rotate*PS_RAD_DEG); // == PC1_1
+            wcsTrans->x->coeff[0][1] = -sin(rotate*PS_RAD_DEG) * Lambda; // == PC1_2
+            wcsTrans->y->coeff[1][0] = +sin(rotate*PS_RAD_DEG) / Lambda; // == PC2_1
+            wcsTrans->y->coeff[1][0] = +cos(rotate*PS_RAD_DEG); // == PC2_2
             goto got_matrix;
         }
 
         // test the PC00i00j varient:
-        pc1_1 = psMetadataLookupF32 (&status, header, "PC001001");
-        if (status) {
-            pc1_2 = psMetadataLookupF32 (&status, header, "PC001002");
-            pc2_1 = psMetadataLookupF32 (&status, header, "PC002001");
-            pc2_2 = psMetadataLookupF32 (&status, header, "PC002002");
-
-            // XXX EAM : add Elixir polynomial terms here eventually
+        wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1
+        wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2
+        wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1
+        wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2
+
+        if (isPoly) {
+            // Elixir-style polynomial terms
+            for (int i = 0; i < wcsTrans->x->nX; i++) {
+                for (int j = 0; j < wcsTrans->x->nX; j++) {
+                    if (i + j < 2)
+                        continue;
+                    if (i + j > fitOrder)
+                        continue;
+                    sprintf (name, "PCA1dX%1dY%1d", i, j);
+                    wcsTrans->x->coeff[i][j] = psMetadataLookupF32 (&status, header, name);
+                }
+            }
+            // XXX currently, Elixir/DVO cannot accept mixed orders
+            for (int i = 0; i < wcsTrans->y->nX; i++) {
+                for (int j = 0; j < wcsTrans->y->nX; j++) {
+                    if (i + j < 2)
+                        continue;
+                    if (i + j > fitOrder)
+                        continue;
+                    sprintf (name, "PCA2dX%1dY%1d", i, j);
+                    wcsTrans->y->coeff[i][j] = psMetadataLookupF32 (&status, header, name);
+                }
+            }
             goto got_matrix;
         }
@@ -104 +147 @@
     // test the CDi_j varient
     if (cdKeys) {
-        pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1");
-        pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2");
-        pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1");
-        pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2");
-
-        // renormalize to cdelt1, cdelt2, etc
-        double scale = hypot (pc1_1, pc1_2);
+        wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1
+        wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2
+        wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1
+        wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2
+
+        // normalize rotation matrix, generate cdelt1, cdelt2
+        double scale = hypot (wcsTrans->x->coeff[1][0], wcsTrans->y->coeff[0][1]);
         cdelt1 = cdelt2 = scale;
-        pc1_1 /= scale;
-        pc1_2 /= scale;
-        pc2_1 /= scale;
-        pc2_2 /= scale;
+        wcsTrans->x->coeff[1][0] /= scale;
+        wcsTrans->x->coeff[0][1] /= scale;
+        wcsTrans->y->coeff[1][0] /= scale;
+        wcsTrans->y->coeff[0][1] /= scale;
         goto got_matrix;
     }
@@ -140 +183 @@
     *****/
 
+    /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
+       wcsTrans, which will con
+       along with the addition
+       PC00i00j representation:
+     * type (CTYPE)
+     * crval1,2 (in RA,DEC degrees)
+     * crpix1,2 
+     * cdelt1,2 (in degrees / pixel)
+     * pci,j (normalized)
+     * 
+     * we also have plateScale (radians / physical TPA units)
+     * now we convert to pmFPA terms
+     */
+
+    /*** XXXX need to extend these formulae to higher-order terms ***/
+
     {
-        double rX = 1.0;
-        double rY = 1.0;
-
         // XXX free an existing toFPA?
         psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1);
 
-        // basic transformation from chip to FPA
+        //
+
+        // basic transformation from chip to FPA (FPA in pixels)
         toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
         toFPA->x->coeff[1][0] = pc1_1;
@@ -158 +216 @@
         toFPA->y->mask[1][1]  = 1;
 
+        // scale from FPA to TPA (microns / pixel)
+        double pdelt1 = cdelt1*PS_RAD_DEG / plateScale;
+        double pdelt2 = cdelt2*PS_RAD_DEG / plateScale;
+        float rX = 1.0;
+        float rY = 1.0;
+
         // projection from TPA to SKY
-        psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
+        psProjection *toSky = psProjectionAlloc (crval1*PS_RAD_DEG, crval2*PS_RAD_DEG, plateScale, plateScale, type);
 
         if (fpa->toSky == NULL)
         {
-            // XXX for now, use the identity for TPA <--> FPA
-            fpa->toTangentPlane = psPlaneDistortIdentity (1);
-            fpa->fromTangentPlane = psPlaneDistortIdentity (1);
+            fpa->toTPA = psPlaneDistortIdentity (1);
+            fpa->fromTPA = psPlaneDistortIdentity (1);
+            fpa->toTPA->x->coeff[1][0][0][0] = pdelt1;
+            fpa->toTPA->y->coeff[0][1][0][0] = pdelt2;
+            fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1;
+            fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2;
             fpa->toSky = toSky;
         } else
         {
-            if (fpa->toTangentPlane == NULL)
+            if (fpa->toTPA == NULL)
                 psAbort ("wcs", "projection defined, tangent-plane not defined");
-            if (fpa->fromTangentPlane == NULL)
+            if (fpa->fromTPA == NULL)
                 psAbort ("wcs", "projection defined, tangent-plane not defined");
 
-            // adjust for common toSky for mosaic:
-            // ignore the TPA since toTPA is identity
+            // convert from pixels on this chip to pixels on reference chip
+            // rX has units of refpixels / pixel
+            rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0];
+            rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0];
+
+            // correct to common plate scale (output is in refpixel units)
+            toFPA->x->coeff[0][0] *= rX;
+            toFPA->x->coeff[1][0] *= rX;
+            toFPA->x->coeff[0][1] *= rX;
+            toFPA->y->coeff[0][0] *= rY;
+            toFPA->y->coeff[1][0] *= rY;
+            toFPA->y->coeff[0][1] *= rY;
+
+            // adjust for common toSky, toTPA for mosaic:
             // find the FPA coordinate of 0,0 for this chip.
+            psPlane *chip = psPlaneAlloc();
             psPlane *fp = psPlaneAlloc();
-            psPlane *chip = psPlaneAlloc();
+            psPlane *tp = psPlaneAlloc();
             psSphere *sky = psSphereAlloc();
             chip->x = chip->y = 0;
 
             psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
-            p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate
-            p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+            psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0);
+            p_psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
+            p_psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+            psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
 
             toFPA->x->coeff[0][0] = fp->x;
             toFPA->y->coeff[0][0] = fp->y;
-
-            rX = toSky->Xs/fpa->toSky->Xs;
-            rY = toSky->Ys/fpa->toSky->Ys;
-
-            // correct to common plate scale
-            toFPA->x->coeff[1][0] *= rX;
-            toFPA->x->coeff[0][1] *= rX;
-            toFPA->y->coeff[1][0] *= rY;
-            toFPA->y->coeff[0][1] *= rY;
 
             psFree (fp);
@@ -207 +280 @@
         chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA);
 
-        while (fpa->toSky->R <        0)
+        // this can take a very long time...
+        while (fpa->toSky->R < 0)
             fpa->toSky->R += 2.0*M_PI;
         while (fpa->toSky->R > 2.0*M_PI)
@@ -213 +287 @@
 
         // remove the correction to the common plate scale
+        // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry
+        // and 2) we are going to measure the mosaic distortion in the next step.
+        // XXX perhaps make this its own function? (I'll need to store rX somewhere).
         if (isMosaic)
         {
+            chip->toFPA->x->coeff[0][0] /= rX;
             chip->toFPA->x->coeff[1][0] /= rX;
             chip->toFPA->x->coeff[0][1] /= rX;
+            chip->toFPA->y->coeff[0][0] /= rY;
             chip->toFPA->y->coeff[1][0] /= rY;
             chip->toFPA->y->coeff[0][1] /= rY;
@@ -232 +311 @@
 
         psLogMsg ("psastro", 3, "field center: %f,%f, plate scale: %f,%f (arcsec/pixel)\n",
-                  DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D,
-                  3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys);
-    }
+                  PS_DEG_RAD*fpa->toSky->R, PS_DEG_RAD*fpa->toSky->D,
+                  3600*PS_DEG_RAD*fpa->toSky->Xs, 3600*PS_DEG_RAD*fpa->toSky->Ys);
+    }
+    # endif
     return true;
 }
@@ -293 +373 @@
 
         // solve for CDELT1,2 (degrees / pixel)
-        cdelt1 = DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
-    cdelt2 = DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
+        cdelt1 = PS_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
+    cdelt2 = PS_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
 
     // L,M = toFPA(X,Y)
@@ -365 +445 @@
         coords[0].crval1 -= 360.0;
 
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*DEG_RAD);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", Xo);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", Yo);
@@ -397 +477 @@
     # else
 
-        psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*DEG_RAD);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*DEG_RAD);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale);
 
     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
@@ -486 +566 @@
 
     // XXX need to handle the plateScale??
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*DEG_RAD);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", 0.0);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", 0.0);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale);
 
     if (toTP->x->nX != toTP->x->nY)

trunk/psModules/src/astrom/pmAstrometryWCS.h

@@ -7 +7 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-12-09 21:06:44 $
+*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-12-10 02:06:47 $
 *
 *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -25 +25 @@
 bool psPlaneDistortIsIdentity (psPlaneDistort *distort);
 
+# define PS_DEG_RAD 57.295779513082322
+# define PS_RAD_DEG  0.017453292519943
+
 #endif // PM_ASTROMETRY_WCS_H

trunk/psModules/src/camera/pmFPA.c

@@ -109 +109 @@
     # if FPA_ASTROM
 
-    psFree(fpa->fromTangentPlane);
-    psFree(fpa->toTangentPlane);
-    psFree(fpa->projection);
+    psFree(fpa->fromTPA);
+    psFree(fpa->toTPA);
+    psFree(fpa->toSky);
     # endif
 }
@@ -326 +326 @@
     #if FPA_ASTROM
 
-    tmpFPA->fromTangentPlane = NULL;
-    tmpFPA->toTangentPlane = NULL;
-    tmpFPA->projection = NULL;
+    tmpFPA->fromTPA = NULL;
+    tmpFPA->toTPA = NULL;
+    tmpFPA->toSky = NULL;
     #endif
 

trunk/psModules/src/camera/pmFPA.h

@@ -9 +9 @@
 /// @author Eugene Magnier, IfA
 ///
-/// @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2006-10-21 03:30:05 $
+/// @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2006-12-10 02:06:47 $
 ///
 /// Copyright 2005-2006 Institute for Astronomy, University of Hawaii
@@ -34 +34 @@
     #if FPA_ASTROM
     // Astrometric transformations
-    psPlaneDistort *fromTangentPlane;   ///< Transformation from tangent plane to focal plane, or NULL
-    psPlaneDistort *toTangentPlane;     ///< Transformation from focal plane to tangent plane, or NULL
-    psProjection *projection;           ///< Projection from tangent plane to sky, or NULL
+    psPlaneDistort *fromTPA;  ///< Transformation from tangent plane to focal plane, or NULL
+    psPlaneDistort *toTPA;  ///< Transformation from focal plane to tangent plane, or NULL
+    psProjection *toSky;         ///< Projection from tangent plane to sky, or NULL
     #endif
     // Information

trunk/psModules/src/psmodules.h

@@ -52 +52 @@
 
 // the following headers are from psModule:astrom
+#include <pmAstrometryWCS.h>
 #include <pmAstrometryObjects.h>
 #include <pmAstrometryDistortion.h>

trunk/psModules/test/astrom/Makefile.am

@@ -1 @@
-AM_LDFLAGS = -L$(top_builddir)/src -lpsmodules $(PSMODULES_LIBS)
-AM_CFLAGS  = @AM_CFLAGS@ $(PSMODULES_CFLAGS) $(SRCINC)
-AM_CPPFLAGS = $(SRCINC) $(PSLIB_CFLAGS)
 
-check_PROGRAMS =
+
+AM_CPPFLAGS = \
+        $(SRCINC) \
+        -I$(top_srcdir)/test/tap/src \
+        -I$(top_srcdir)/test/pstap/src \
+        $(PSMODULES_CFLAGS)
+
+AM_LDFLAGS = \
+        $(top_builddir)/src/libpsmodules.la  \
+        $(top_builddir)/test/tap/src/libtap.la \
+        $(top_builddir)/test/pstap/src/libpstap.la \
+        $(PSMODULES_LIBS)
+
+check_PROGRAMS = \
+        tap_pmAstrometryWCS
 
 test: check