Index: trunk/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 10873)
+++ trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 10921)
@@ -7,6 +7,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-01-01 21:05:46 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-01-05 20:20:06 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -236,4 +236,5 @@
     wcs->crpix2 = psMetadataLookupF64 (&status, header, "CRPIX2");
     wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type);
+    // XXX I think this is wrong for linear proj
 
     // test the CDELTi varient
@@ -253,5 +254,6 @@
         }
 
-        // test the PC00i00j varient:
+        // FITS WCS PCi,j has units of unity
+        // wcs->trans has units of degrees/pixel
         wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC001001"); // == PC1_1
         wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC001002"); // == PC1_2
@@ -283,10 +285,10 @@
     // test the CDi_j varient
     if (cdKeys) {
-        wcs->cdelt1 = 1.0;
-        wcs->cdelt2 = 1.0;
         wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1
         wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2
         wcs->trans->y->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD2_1"); // == PC2_1
         wcs->trans->y->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD2_2"); // == PC2_2
+        wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
+        wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
         return wcs;
     }
@@ -321,7 +323,7 @@
 
     // XXX make it optional to write out CDi_j terms, or other versions
-    // solve for CDELT1,2 (degrees / pixel)
-    double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->y->coeff[1][0]);
-    double cdelt2 = hypot (wcs->trans->x->coeff[0][1], wcs->trans->y->coeff[0][1]);
+    // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
+    double cdelt1 = wcs->cdelt1;
+    double cdelt2 = wcs->cdelt2;
     psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
     psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
@@ -336,4 +338,5 @@
     // XXX currently, Elixir/DVO cannot accept mixed orders
     // XXX need to respect the masks
+    // XXX is wcs->cdelt1,2 always consistent?
     int fitOrder = wcs->trans->x->nX;
     if (fitOrder > 1) {
@@ -345,7 +348,7 @@
                     continue;
                 sprintf (name, "PCA1X%1dY%1d", i, j);
-                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
+                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
                 sprintf (name, "PCA2X%1dY%1d", i, j);
-                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
+                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
             }
         }
@@ -357,40 +360,24 @@
 
 // interpret header WCS (only handles traditional WCS for the moment)
-// pixelScale is the pixel size in microns per pixel
+// pixelScale is the pixel size in microns/pixel
 bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double pixelScale)
 {
     psPlaneTransform *toFPA;
 
-    /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
-     * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following 
-     * elements defined:
-     * type (projection type)
-     * crval1,2 (in RA,DEC degrees)
-     * crpix1,2 
-     * cdelt1,2 (in degrees / pixel)
-     * pixelScale (microns / pixel)
-     * 
-     * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
-     * 1) the output is in microns (not degrees): divide by cdelt1,2
-     * 2) X,Y are applied directly, without an applied Xo,Yo offset
-     * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
-     */
-
-    // create transformation with 0,0 reference pixel
+    // create transformation with 0,0 reference pixel and units of degrees/pixel
     toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
 
-    // modify scale of toFPA to yield L,M in microns
-    double cdelt1 = hypot (toFPA->x->coeff[1][0], toFPA->x->coeff[0][1]);
-    double cdelt2 = hypot (toFPA->y->coeff[1][0], toFPA->y->coeff[0][1]);
+    // modify scale of toFPA to have units of microns/pixel
+    // cdelt1,2 has units of degree/pixel
     for (int i = 0; i <= toFPA->x->nX; i++) {
         for (int j = 0; j <= toFPA->x->nX; j++) {
-            toFPA->x->coeff[i][j] *= pixelScale/cdelt1;
-            toFPA->y->coeff[i][j] *= pixelScale/cdelt2;
-        }
-    }
-
-    // scale from TPA (microns) to degrees
-    double pdelt1 = cdelt1 / pixelScale;
-    double pdelt2 = cdelt2 / pixelScale;
+            toFPA->x->coeff[i][j] *= pixelScale/wcs->cdelt1;
+            toFPA->y->coeff[i][j] *= pixelScale/wcs->cdelt2;
+        }
+    }
+
+    // pdelt1,2 has units of degree/micron
+    double pdelt1 = wcs->cdelt1 / pixelScale;
+    double pdelt2 = wcs->cdelt2 / pixelScale;
 
     // projection from TPA (linear microns) to SKY (radians)
@@ -485,5 +472,5 @@
      */
 
-    // create transformation with 0,0 reference pixel
+    // create transformation with 0,0 reference pixel and units of microns/pixel
     chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
 
@@ -502,9 +489,18 @@
 bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs)
 {
-    // projection from TPA to SKY
-    fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->toSky->Xs, wcs->toSky->Ys, wcs->toSky->type);
+    // projection from TPA (microns) to SKY (radians)
+    // cdelt1,2 has units of degrees/micron
+    fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->cdelt1*PM_RAD_DEG, wcs->cdelt2*PM_RAD_DEG, wcs->toSky->type);
 
     // create transformation with 0,0 reference pixel
     fpa->toTPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
+
+    // convert fpa->toTPA to units of unity (microns/micron)
+    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
+        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
+            fpa->toTPA->x->coeff[i][j] /= wcs->cdelt1;
+            fpa->toTPA->y->coeff[i][j] /= wcs->cdelt2;
+        }
+    }
 
     // XXX this needs to perform the full (non-linear) inversion
@@ -525,4 +521,5 @@
 
     // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
+    // XXX not really: toTPA needs to have unity scale for distortion fitting function
     if (!psPlaneTransformIsDiagonal (fpa->toTPA))
         psAbort ("psastro", "invalid TPA transformation");
@@ -535,5 +532,5 @@
     pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
 
-    // convert projection from TPA to SKY into wcs projection (degrees to radians)
+    // convert projection from FPA to SKY into wcs projection (degrees to radians)
     wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type);
     wcs->crval1 = fpa->toSky->R*PM_DEG_RAD;
@@ -543,5 +540,5 @@
     psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
 
-    // create wcs transform from toFPA, result converts pixels to microns
+    // create wcs transform from toFPA, resulting transformation has units of microns/pixel
     // adjust wcs transform to use center as reference coordinate
     psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
@@ -552,19 +549,20 @@
     psFree (center);
 
-    // pdelt1,2 converts from microns->degrees
+    // pdelt1,2 has units of degrees/micron
     double pdelt1 = fpa->toSky->Xs * PM_DEG_RAD;
     double pdelt2 = fpa->toSky->Ys * PM_DEG_RAD;
 
-    // cdelt1,2 converts from pixels->degrees
-    wcs->cdelt1 = pdelt1 * hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
-    wcs->cdelt2 = pdelt2 * hypot (chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
-
-    // convert wcs->trans to a matrix which yields L,M in degrees
+    // convert wcs->trans to a matrix with units of degrees/pixel
     for (int i = 0; i <= wcs->trans->x->nX; i++) {
-        for (int j = 0; j <= wcs->trans->x->nX; j++) {
+        for (int j = 0; j <= wcs->trans->x->nY; j++) {
             wcs->trans->x->coeff[i][j] *= pdelt1;
             wcs->trans->y->coeff[i][j] *= pdelt2;
         }
     }
+
+    // cdelt1,2 has units of degrees/pixel
+    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
+    wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
+
     return wcs;
 }
@@ -588,7 +586,6 @@
 
     // Chip to FPA transformation is a Cartesian 'projection'
+    // reference pixel for FPA is 0.0, 0.0
     wcs->toSky = psProjectionAlloc (0.0, 0.0, 1.0, 1.0, PS_PROJ_WRP);
-
-    // reference pixel for FPA is 0.0, 0.0
     wcs->crval1 = 0.0;
     wcs->crval2 = 0.0;
@@ -598,4 +595,5 @@
 
     // adjust wcs transform to use center as reference coordinate
+    // resulting transformation has units of microns/pixel
     psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
 
@@ -605,7 +603,7 @@
     psFree (center);
 
-    // output coordinates are in FPA pixels
-    wcs->cdelt1 = 1.0;
-    wcs->cdelt2 = 1.0;
+    // output coordinates are in microns : CDELT1,2 has units of microns/pixel
+    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
+    wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
 
     return wcs;
@@ -633,4 +631,5 @@
 
     // adjust wcs transform to use center as reference coordinate
+    // resulting transformation has units of unity (microns/micron)
     psPlaneTransformSetCenter (wcs->trans, fpa->toTPA, center->x, center->y);
 
@@ -640,7 +639,19 @@
     psFree (center);
 
-    // output coordinates are in FPA pixels
-    wcs->cdelt1 = 1.0;
-    wcs->cdelt2 = 1.0;
+    // pdelt1,2 has units of degrees/micron
+    double pdelt1 = fpa->toSky->Xs * PM_DEG_RAD;
+    double pdelt2 = fpa->toSky->Ys * PM_DEG_RAD;
+
+    // convert wcs->trans to units of degree/micron
+    for (int i = 0; i <= wcs->trans->x->nX; i++) {
+        for (int j = 0; j <= wcs->trans->x->nY; j++) {
+            wcs->trans->x->coeff[i][j] *= pdelt1;
+            wcs->trans->y->coeff[i][j] *= pdelt2;
+        }
+    }
+
+    // cdelt1,2 has units of degrees/micron
+    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
+    wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
 
     return wcs;
@@ -699,2 +710,17 @@
 */
 
+/* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
+ * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following 
+ * elements defined:
+ * type (projection type)
+ * crval1,2 (in RA,DEC degrees)
+ * crpix1,2 
+ * cdelt1,2 (in degrees / pixel)
+ * pixelScale (microns / pixel)
+ * 
+ * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
+ * 1) the output is in microns (not degrees): divide by cdelt1,2
+ * 2) X,Y are applied directly, without an applied Xo,Yo offset
+ * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
+ */
+