Index: /trunk/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- /trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 10872)
+++ /trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 10873)
@@ -7,6 +7,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-30 03:27:19 $
+ *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-01-01 21:05:46 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -48,5 +48,5 @@
 // convert toFPA / toSky components to pmAstromWCS
 // tolerance is convergence for inversion of non-linear terms in pixels
-bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol)
+bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, double tol)
 {
     // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
@@ -95,5 +95,5 @@
 
 // convert chip->toFPA components to bilevel WCS
-bool pmAstromWriteBilevelChip (psMetadata *header, const pmChip *chip, float tol)
+bool pmAstromWriteBilevelChip (psMetadata *header, const pmChip *chip, double tol)
 {
     pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol);
@@ -107,5 +107,5 @@
 
 // convert fpa->toTPA, fpa->toSky components to bilevel WCS
-bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol)
+bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, double tol)
 {
     pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol);
@@ -198,6 +198,6 @@
     // XXX add check for CROTA2
     int fitOrder = psMetadataLookupS32 (&isPoly, header, "NPLYTERM");
-    psMetadataLookupF32 (&pcKeys, header, "PC001001");
-    psMetadataLookupF32 (&cdKeys, header, "CD1_1");
+    psMetadataLookupF64 (&pcKeys, header, "PC001001");
+    psMetadataLookupF64 (&cdKeys, header, "CD1_1");
 
     if (cdKeys && pcKeys) {
@@ -231,18 +231,18 @@
     // and then define a transformation from degrees to degrees
 
-    wcs->crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");
-    wcs->crval2 = psMetadataLookupF32 (&status, header, "CRVAL2");
-    wcs->crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1");
-    wcs->crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2");
+    wcs->crval1 = psMetadataLookupF64 (&status, header, "CRVAL1");
+    wcs->crval2 = psMetadataLookupF64 (&status, header, "CRVAL2");
+    wcs->crpix1 = psMetadataLookupF64 (&status, header, "CRPIX1");
+    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);
 
     // test the CDELTi varient
     if (pcKeys) {
-        wcs->cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1");
-        wcs->cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2");
+        wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
+        wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
 
         // test the CROTAi varient:
         // XXX double check lambda..
-        double rotate = psMetadataLookupF32 (&status, header, "CROTA2");
+        double rotate = psMetadataLookupF64 (&status, header, "CROTA2");
         if (status) {
             wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1
@@ -254,8 +254,8 @@
 
         // test the PC00i00j varient:
-        wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1
-        wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2
-        wcs->trans->y->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1
-        wcs->trans->y->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2
+        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
+        wcs->trans->y->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC002001"); // == PC2_1
+        wcs->trans->y->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC002002"); // == PC2_2
 
         if (isPoly) {
@@ -272,7 +272,7 @@
                     }
                     sprintf (name, "PCA1X%1dY%1d", i, j);
-                    wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name);
+                    wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF64 (&status, header, name);
                     sprintf (name, "PCA2X%1dY%1d", i, j);
-                    wcs->trans->y->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name);
+                    wcs->trans->y->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF64 (&status, header, name);
                 }
             }
@@ -285,8 +285,8 @@
         wcs->cdelt1 = 1.0;
         wcs->cdelt2 = 1.0;
-        wcs->trans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1
-        wcs->trans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2
-        wcs->trans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1
-        wcs->trans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2
+        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
         return wcs;
     }
@@ -328,8 +328,8 @@
 
     // test the PC00i00j varient:
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
+    psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
+    psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
+    psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
+    psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
 
     // Elixir-style polynomial terms
@@ -345,7 +345,7 @@
                     continue;
                 sprintf (name, "PCA1X%1dY%1d", i, j);
-                psMetadataAddF32 (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(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
                 sprintf (name, "PCA2X%1dY%1d", i, j);
-                psMetadataAddF32 (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(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
             }
         }
@@ -390,20 +390,14 @@
     }
 
-    // scale from FPA to TPA (degrees / micron)
+    // scale from TPA (microns) to degrees
     double pdelt1 = cdelt1 / pixelScale;
     double pdelt2 = cdelt2 / pixelScale;
-    float rX = 1.0;
-    float rY = 1.0;
-
-    // projection from TPA ("linear" degrees) to SKY (radians)
-    psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, PM_RAD_DEG, PM_RAD_DEG, wcs->toSky->type);
+
+    // projection from TPA (linear microns) to SKY (radians)
+    psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, PM_RAD_DEG*pdelt1, PM_RAD_DEG*pdelt2, wcs->toSky->type);
 
     if (fpa->toSky == NULL) {
         fpa->toTPA = psPlaneTransformIdentity (1);
         fpa->fromTPA = psPlaneTransformIdentity (1);
-        fpa->toTPA->x->coeff[1][0] = pdelt1;
-        fpa->toTPA->y->coeff[0][1] = pdelt2;
-        fpa->fromTPA->x->coeff[1][0] = 1.0 / pdelt1;
-        fpa->fromTPA->y->coeff[0][1] = 1.0 / pdelt2;
         fpa->toSky = toSky;
     } else {
@@ -420,6 +414,6 @@
         // 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];
-        rY = pdelt2 / fpa->toTPA->y->coeff[0][1];
+        double rX = toSky->Xs / fpa->toSky->Xs;
+        double rY = toSky->Ys / fpa->toSky->Ys;
         for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
             for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
@@ -527,5 +521,5 @@
 // convert toFPA / toSky components to pmAstromWCS
 // tolerance is allowed error in center solution in pixels
-pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, float tol)
+pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, double tol)
 {
 
@@ -549,4 +543,5 @@
     psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
 
+    // create wcs transform from toFPA, result converts pixels to microns
     // adjust wcs transform to use center as reference coordinate
     psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
@@ -557,15 +552,17 @@
     psFree (center);
 
-    // cdelt1,2 convert from pixels->degrees
-    double cdelt1 = fpa->toTPA->x->coeff[1][0];
-    double cdelt2 = fpa->toTPA->y->coeff[0][1];
-    wcs->cdelt1 = cdelt1;
-    wcs->cdelt2 = cdelt2;
-
-    // convert wcs->trans to a matrix which yields L,M in pixels
+    // pdelt1,2 converts from microns->degrees
+    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
     for (int i = 0; i <= wcs->trans->x->nX; i++) {
         for (int j = 0; j <= wcs->trans->x->nX; j++) {
-            wcs->trans->x->coeff[i][j] *= cdelt1;
-            wcs->trans->y->coeff[i][j] *= cdelt2;
+            wcs->trans->x->coeff[i][j] *= pdelt1;
+            wcs->trans->y->coeff[i][j] *= pdelt2;
         }
     }
@@ -580,5 +577,5 @@
 
 // convert the chip-level toFPA to a wcs polynomial transformation
-pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, float tol)
+pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, double tol)
 {
     // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
@@ -616,5 +613,5 @@
 
 // convert the fpa-level toTPA, toSky to a wcs polynomial transformation
-pmAstromWCS *pmAstromWCSBilevelMosaicFromFPA (const pmFPA *fpa, float tol)
+pmAstromWCS *pmAstromWCSBilevelMosaicFromFPA (const pmFPA *fpa, double tol)
 {
     // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY