Changeset 39697 for branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c

Timestamp:

Sep 15, 2016, 3:10:46 PM (10 years ago)

Author:

eugene

Message:

rather than attempt to force the WRP chip-level CRVAL1,2 entries to be (0,0), we re-work the keywords to use CRPIX1,2 = (0,0)

File:

• branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c (modified) (40 diffs)

branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c

@@ -43 +43 @@
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     if (!wcs) {
-        return false;
+        return false;
     }
 
@@ -70 +70 @@
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     if (!wcs) {
-        return false;
+        return false;
     }
 
@@ -86 +86 @@
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     if (!wcs) {
-        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from header");
-        return false;
+        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from header");
+        return false;
     }
 
@@ -96 +96 @@
 
     if (!status1 || !status2) {
-        Nx = psMetadataLookupS32 (&status1, header, "IMNAXIS1");
-        Ny = psMetadataLookupS32 (&status2, header, "IMNAXIS2");
+        Nx = psMetadataLookupS32 (&status1, header, "IMNAXIS1");
+        Ny = psMetadataLookupS32 (&status2, header, "IMNAXIS2");
     }
 
     if (!status1 || !status2) {
-        Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
-        Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
+        Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
+        Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
     }
 
     if (!status1 || !status2) {
-        psFree (wcs);
-        psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header");
-        return false;
+        psFree (wcs);
+        psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header");
+        return false;
     }
 
@@ -123 +123 @@
     pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol);
     if (!wcs) {
-        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
-        return false;
+        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
+        return false;
     }
 
@@ -139 +139 @@
     pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol);
     if (!wcs) {
-        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
-        return false;
+        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
+        return false;
     }
 
@@ -163 +163 @@
 
     if (chip == NULL)
-        return false;
+        return false;
     if (sky == NULL)
-        return false;
+        return false;
     if (wcs == NULL)
-        return false;
+        return false;
 
     psPlane *Chip = psPlaneAlloc();
@@ -188 +188 @@
 
     if (chip == NULL)
-        return false;
+        return false;
     if (sky == NULL)
-        return false;
+        return false;
     if (wcs == NULL)
-        return false;
+        return false;
 
     psError(PS_ERR_UNKNOWN, true, "not yet implemented: needs to invert the transformation");
@@ -223 +223 @@
     char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
     if (!status) {
-        psLogMsg ("psastro", 5, "warning: no WCS metadata in header\n");
-        return NULL;
+        psLogMsg ("psastro", 5, "warning: no WCS metadata in header\n");
+        return NULL;
     }
 
@@ -232 +232 @@
     type = psProjectTypeFromString (ctype);
     if (type == PS_PROJ_NTYPE) {
-        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
-        return NULL;
+        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
+        return NULL;
     }
 
@@ -243 +243 @@
 
     if (cdKeys && pcKeys) {
-        // XXX make this an option
-        psLogMsg ("psastro", 5, "warning: both CDi_j and PC00i00j defined in headers, using PC00i00j terms\n");
+        // XXX make this an option
+        psLogMsg ("psastro", 5, "warning: both CDi_j and PC00i00j defined in headers, using PC00i00j terms\n");
     }
     if (!cdKeys && !pcKeys) {
-        psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms");
-        // XXX we could default here to RA, DEC, ROTANGLE
-        return NULL;
+        psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms");
+        // XXX we could default here to RA, DEC, ROTANGLE
+        return NULL;
     }
     if (isPoly) {
-        if (!pcKeys) {
-            psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
-            return NULL;
-        }
-        if (fitOrder == 0)
-            fitOrder = 1;
-        if ((fitOrder > 3) || (fitOrder < 1)) {
-            psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
-            return NULL;
-        }
+        if (!pcKeys) {
+            psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
+            return NULL;
+        }
+        if (fitOrder == 0)
+            fitOrder = 1;
+        if ((fitOrder > 3) || (fitOrder < 1)) {
+            psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
+            return NULL;
+        }
     } else {
-        fitOrder = 1;
+        fitOrder = 1;
     }
 
@@ -280 +280 @@
     // XXX if type == ZPN, look for PV2_%d elements:
     if (type == PS_PROJ_ZPN) {
-      psVector *maxRadial = psVectorAlloc (21, PS_TYPE_F64);
-      for (int i = 0; i <= 20; i++) {
-        char name[64];
-        snprintf (name, 64, "PV2_%d", i);
-
-        maxRadial->data.F64[i] = 0.0;
-        double value = psMetadataLookupF64 (&status, header, name);
-
-        if (status) {
-          maxRadial->data.F64[i] = value;
-          maxRadial->n = i;
-        }
-
-        // PV2_1 is implicit if not present
-        if ((i == 1) && !status) {
-          maxRadial->data.F64[i] = 1.0;
-          continue;
-        }
-      }
-      maxRadial->n ++;
-      wcs->toSky->radial = maxRadial;
+        psVector *maxRadial = psVectorAlloc (21, PS_TYPE_F64);
+        for (int i = 0; i <= 20; i++) {
+            char name[64];
+            snprintf (name, 64, "PV2_%d", i);
+
+            maxRadial->data.F64[i] = 0.0;
+            double value = psMetadataLookupF64 (&status, header, name);
+
+            if (status) {
+                maxRadial->data.F64[i] = value;
+                maxRadial->n = i;
+            }
+
+            // PV2_1 is implicit if not present
+            if ((i == 1) && !status) {
+                maxRadial->data.F64[i] = 1.0;
+                continue;
+            }
+        }
+        maxRadial->n ++;
+        wcs->toSky->radial = maxRadial;
     }
 
@@ -314 +314 @@
     // test the CDELTi varient
     if (pcKeys) {
-        wcs->wcsCDkeys = 0;
-        wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
-        wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
-
-        // test the CROTAi varient:
-        // XXX double check lambda..
-        double rotate = psMetadataLookupF64 (&status, header, "CROTA2");
-        if (status) {
-            wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1
-            wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
-            wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
-            wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
-            return wcs;
-        }
-
-        // 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
-        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) {
-            // Elixir-style polynomial terms
-            // XXX currently, Elixir/DVO cannot accept mixed orders
-            for (int i = 0; i <= fitOrder; i++) {
-                for (int j = 0; j <= fitOrder; j++) {
-                    if (i + j < 2)
-                        continue;
-                    if (i + j > fitOrder) {
-                        wcs->trans->x->coeffMask[i][j] = PS_POLY_MASK_SET;
-                        wcs->trans->y->coeffMask[i][j] = PS_POLY_MASK_SET;
-                        continue;
-                    }
-                    sprintf (name, "PCA1X%1dY%1d", i, j);
-                    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) * psMetadataLookupF64 (&status, header, name);
-                }
-            }
-        }
-        return wcs;
+        wcs->wcsCDkeys = 0;
+        wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
+        wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
+
+        // test the CROTAi varient:
+        // XXX double check lambda..
+        double rotate = psMetadataLookupF64 (&status, header, "CROTA2");
+        if (status) {
+            wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1
+            wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
+            wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
+            wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
+            return wcs;
+        }
+
+        // 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
+        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) {
+            // Elixir-style polynomial terms
+            // XXX currently, Elixir/DVO cannot accept mixed orders
+            for (int i = 0; i <= fitOrder; i++) {
+                for (int j = 0; j <= fitOrder; j++) {
+                    if (i + j < 2)
+                        continue;
+                    if (i + j > fitOrder) {
+                        wcs->trans->x->coeffMask[i][j] = PS_POLY_MASK_SET;
+                        wcs->trans->y->coeffMask[i][j] = PS_POLY_MASK_SET;
+                        continue;
+                    }
+                    sprintf (name, "PCA1X%1dY%1d", i, j);
+                    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) * psMetadataLookupF64 (&status, header, name);
+                }
+            }
+        }
+        return wcs;
     }
 
     // test the CDi_j varient
     if (cdKeys) {
-        wcs->wcsCDkeys = 1;
-
-        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;
+        wcs->wcsCDkeys = 1;
+
+        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;
     }
     psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");
@@ -402 +402 @@
 
     if (wcs->toSky->type == PS_PROJ_ZPN) {
-      psAssert (wcs->toSky->radial, "missing radial vector");
-      for (int i = 0; i < wcs->toSky->radial->n; i++) {
-        if (wcs->toSky->radial->data.F64[i] == 0.0) continue;
-        snprintf (name, 16, "PV2_%d", i);
-        psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->toSky->radial->data.F64[i]);
-      }
+        psAssert (wcs->toSky->radial, "missing radial vector");
+        for (int i = 0; i < wcs->toSky->radial->n; i++) {
+            if (wcs->toSky->radial->data.F64[i] == 0.0) continue;
+            snprintf (name, 16, "PV2_%d", i);
+            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->toSky->radial->data.F64[i]);
+        }
     }
 
@@ -414 +414 @@
     if (!wcs->wcsCDkeys) {
 
-      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);
-
-      // test the PC00i00j varient:
-      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
-      // 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) {
-        for (int i = 0; i <= fitOrder; i++) {
-          for (int j = 0; j <= fitOrder; j++) {
-            if (i + j < 2)
-              continue;
-            if (i + j > fitOrder)
-              continue;
-            sprintf (name, "PCA1X%1dY%1d", i, 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(cdelt1, i) / pow(cdelt2, j));
-          }
-        }
-        psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
-      }
-
-      // remove any existing 'CDi_j style' wcs keywords
-      if (psMetadataLookup(header, "CD1_1")) {
-        psMetadataRemoveKey(header, "CD1_1");
-        psMetadataRemoveKey(header, "CD1_2");
-        psMetadataRemoveKey(header, "CD2_1");
-        psMetadataRemoveKey(header, "CD2_2");
-      }
-
-      // Remove 'CDi_jX' WCS keywords
-      psString cd11 = psStringCopy("CD1_1 ");
-      psString cd12 = psStringCopy("CD1_2 ");
-      psString cd21 = psStringCopy("CD2_1 ");
-      psString cd22 = psStringCopy("CD2_2 ");
-      for (char extra = 'A'; extra <= 'Z'; extra++) {
-          cd11[strlen(cd11)-1] = extra;
-          if (psMetadataLookup(header, cd11)) {
-              cd12[strlen(cd12)-1] = extra;
-              cd21[strlen(cd21)-1] = extra;
-              cd22[strlen(cd22)-1] = extra;
-              psMetadataRemoveKey(header, cd11);
-              psMetadataRemoveKey(header, cd12);
-              psMetadataRemoveKey(header, cd21);
-              psMetadataRemoveKey(header, cd22);
-          }
-      }
-      psFree(cd11);
-      psFree(cd12);
-      psFree(cd21);
-      psFree(cd22);
+        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);
+
+        // test the PC00i00j varient:
+        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
+        // 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) {
+            for (int i = 0; i <= fitOrder; i++) {
+                for (int j = 0; j <= fitOrder; j++) {
+                    if (i + j < 2)
+                        continue;
+                    if (i + j > fitOrder)
+                        continue;
+                    sprintf (name, "PCA1X%1dY%1d", i, 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(cdelt1, i) / pow(cdelt2, j));
+                }
+            }
+            psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
+        }
+
+        // remove any existing 'CDi_j style' wcs keywords
+        if (psMetadataLookup(header, "CD1_1")) {
+            psMetadataRemoveKey(header, "CD1_1");
+            psMetadataRemoveKey(header, "CD1_2");
+            psMetadataRemoveKey(header, "CD2_1");
+            psMetadataRemoveKey(header, "CD2_2");
+        }
+
+        // Remove 'CDi_jX' WCS keywords
+        psString cd11 = psStringCopy("CD1_1 ");
+        psString cd12 = psStringCopy("CD1_2 ");
+        psString cd21 = psStringCopy("CD2_1 ");
+        psString cd22 = psStringCopy("CD2_2 ");
+        for (char extra = 'A'; extra <= 'Z'; extra++) {
+            cd11[strlen(cd11)-1] = extra;
+            if (psMetadataLookup(header, cd11)) {
+                cd12[strlen(cd12)-1] = extra;
+                cd21[strlen(cd21)-1] = extra;
+                cd22[strlen(cd22)-1] = extra;
+                psMetadataRemoveKey(header, cd11);
+                psMetadataRemoveKey(header, cd12);
+                psMetadataRemoveKey(header, cd21);
+                psMetadataRemoveKey(header, cd22);
+            }
+        }
+        psFree(cd11);
+        psFree(cd12);
+        psFree(cd21);
+        psFree(cd22);
 
 
     } else {
 
-      psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
-      psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]);
-      psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]);
-      psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]);
-
-      if (psMetadataLookup(header, "PC001001")) {
-        psMetadataRemoveKey(header, "PC001001");
-        psMetadataRemoveKey(header, "PC001002");
-        psMetadataRemoveKey(header, "PC002001");
-        psMetadataRemoveKey(header, "PC002002");
-      }
+        psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
+        psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]);
+        psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]);
+        psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]);
+
+        if (psMetadataLookup(header, "PC001001")) {
+            psMetadataRemoveKey(header, "PC001001");
+            psMetadataRemoveKey(header, "PC001002");
+            psMetadataRemoveKey(header, "PC002001");
+            psMetadataRemoveKey(header, "PC002002");
+        }
     }
 
@@ -507 +507 @@
     // 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/wcs->cdelt1;
-            toFPA->y->coeff[i][j] *= pixelScale/wcs->cdelt2;
-        }
+        for (int j = 0; j <= toFPA->x->nX; j++) {
+            toFPA->x->coeff[i][j] *= pixelScale/wcs->cdelt1;
+            toFPA->y->coeff[i][j] *= pixelScale/wcs->cdelt2;
+        }
     }
 
@@ -524 +524 @@
         psFree(fpa->toTPA);
         psFree(fpa->fromTPA);
-        fpa->toTPA = psPlaneTransformIdentity (1);
-        fpa->fromTPA = psPlaneTransformIdentity (1);
-        fpa->toSky = toSky;
+        fpa->toTPA = psPlaneTransformIdentity (1);
+        fpa->fromTPA = psPlaneTransformIdentity (1);
+        fpa->toSky = toSky;
     } else {
 
-        // this section allows the loaded chip to be included in an fpa structure in which
-        // other chips have already been loaded (ie, the fpa->toTPA, fpa->toSky components have
-        // already been defined).  we have to adjust to match the existing transformation.
-
-        if (fpa->toTPA == NULL)
-            psAbort("projection defined, tangent-plane not defined");
-        if (fpa->fromTPA == NULL)
-            psAbort("projection defined, tangent-plane not defined");
-
-        // convert from pixels on this chip to pixels on reference chip
-        // rX has units of refpixels / pixel
-        double rX = toSky->Xs / fpa->toSky->Xs;
-        double rY = toSky->Ys / fpa->toSky->Ys;
-
-        for (int i = 0; i <= toFPA->x->nX; i++) {
-            for (int j = 0; j <= toFPA->x->nY; j++) {
-                toFPA->x->coeff[i][j] *= rX;
-                toFPA->y->coeff[i][j] *= rY;
-            }
-        }
-
-        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
-        // XXX this only works if toTPA is at most a linear transformation
-        psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY);
-        for (int i = 0; i <= toFPA->x->nX; i++) {
-          for (int j = 0; j <= toFPA->x->nY; j++) {
-            double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j];
-            double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j];
-            toFPAnew->x->coeff[i][j] = f1 + f2;
-
-            double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j];
-            double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j];
-            toFPAnew->y->coeff[i][j] = g1 + g2;
-          }
-        }
-        toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0];
-        toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0];
-
-        psFree (toFPA);
-        toFPA = toFPAnew;
-
-        // adjust reference pixel for new toSky reference coordinate
-        // find the FPA coordinate of 0,0 for this chip.
-        psPlane *fpOld = psPlaneAlloc();
-        psPlane *fpNew = psPlaneAlloc();
-        psPlane *tp = psPlaneAlloc();
-        psSphere *sky = psSphereAlloc();
-
-        sky->r = toSky->R;
-        sky->d = toSky->D;
-        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
-        psPlaneTransformApply (fpOld, fpa->fromTPA, tp);
-
-        sky->r = fpa->toSky->R;
-        sky->d = fpa->toSky->D;
-        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
-        psPlaneTransformApply (fpNew, fpa->fromTPA, tp);
-
-        toFPA->x->coeff[0][0] -= fpNew->x - fpOld->x;
-        toFPA->y->coeff[0][0] -= fpNew->y - fpOld->y;
-
-        psFree (sky);
-        psFree (tp);
-        psFree (fpOld);
-        psFree (fpNew);
-
-        psFree (toSky);
+        // this section allows the loaded chip to be included in an fpa structure in which
+        // other chips have already been loaded (ie, the fpa->toTPA, fpa->toSky components have
+        // already been defined).  we have to adjust to match the existing transformation.
+
+        if (fpa->toTPA == NULL)
+            psAbort("projection defined, tangent-plane not defined");
+        if (fpa->fromTPA == NULL)
+            psAbort("projection defined, tangent-plane not defined");
+
+        // convert from pixels on this chip to pixels on reference chip
+        // rX has units of refpixels / pixel
+        double rX = toSky->Xs / fpa->toSky->Xs;
+        double rY = toSky->Ys / fpa->toSky->Ys;
+
+        for (int i = 0; i <= toFPA->x->nX; i++) {
+            for (int j = 0; j <= toFPA->x->nY; j++) {
+                toFPA->x->coeff[i][j] *= rX;
+                toFPA->y->coeff[i][j] *= rY;
+            }
+        }
+
+        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
+        // XXX this only works if toTPA is at most a linear transformation
+        psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY);
+        for (int i = 0; i <= toFPA->x->nX; i++) {
+            for (int j = 0; j <= toFPA->x->nY; j++) {
+                double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j];
+                double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j];
+                toFPAnew->x->coeff[i][j] = f1 + f2;
+
+                double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j];
+                double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j];
+                toFPAnew->y->coeff[i][j] = g1 + g2;
+            }
+        }
+        toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0];
+        toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0];
+
+        psFree (toFPA);
+        toFPA = toFPAnew;
+
+        // adjust reference pixel for new toSky reference coordinate
+        // find the FPA coordinate of 0,0 for this chip.
+        psPlane *fpOld = psPlaneAlloc();
+        psPlane *fpNew = psPlaneAlloc();
+        psPlane *tp = psPlaneAlloc();
+        psSphere *sky = psSphereAlloc();
+
+        sky->r = toSky->R;
+        sky->d = toSky->D;
+        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+        psPlaneTransformApply (fpOld, fpa->fromTPA, tp);
+
+        sky->r = fpa->toSky->R;
+        sky->d = fpa->toSky->D;
+        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+        psPlaneTransformApply (fpNew, fpa->fromTPA, tp);
+
+        toFPA->x->coeff[0][0] -= fpNew->x - fpOld->x;
+        toFPA->y->coeff[0][0] -= fpNew->y - fpOld->y;
+
+        psFree (sky);
+        psFree (tp);
+        psFree (fpOld);
+        psFree (fpNew);
+
+        psFree (toSky);
     }
 
@@ -613 +613 @@
     // XXX if the inversion fails, we probably do not have a valid transform anyway
     if (!chip->fromFPA) {
-      psWarning ("failed to find a valid transformation");
-      psFree (chip->toFPA);
-      return false;
+        psWarning ("failed to find a valid transformation");
+        psFree (chip->toFPA);
+        return false;
     }
 
     // this can take a very long time...
     while (fpa->toSky->R < 0)
-        fpa->toSky->R += 2.0*M_PI;
+        fpa->toSky->R += 2.0*M_PI;
     while (fpa->toSky->R > 2.0*M_PI)
-        fpa->toSky->R -= 2.0*M_PI;
+        fpa->toSky->R -= 2.0*M_PI;
 
     fpa->wcsCDkeys = wcs->wcsCDkeys;
 
     psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
-             chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
-             chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
-             chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
+             chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
+             chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
+             chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
 
     psTrace ("psastro", 5, "frFPA: %f %f  (%f,%f),(%f,%f)\n",
-             chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
-             chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
-             chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
+             chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
+             chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
+             chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
 
     return true;
@@ -648 +648 @@
      */
 
-    // create transformation with 0,0 reference pixel and units of microns/pixel
     psFree (chip->toFPA);
-    chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
+    chip->toFPA = psPlaneTransformAlloc(wcs->trans->x->nX, wcs->trans->x->nY);
+
+    // copy the toFPA x,y, transformations to the wcs version
+    chip->toFPA->x = psPolynomial2DCopy (chip->toFPA->x, wcs->trans->x);
+    chip->toFPA->y = psPolynomial2DCopy (chip->toFPA->y, wcs->trans->y);
+
+    // these need to be set based on crval1,2
+    chip->toFPA->x->coeff[0][0] = wcs->crval1;
+    chip->toFPA->y->coeff[0][0] = wcs->crval2;
 
     // determine the inverse transformation: we need the chip pixels covered by this transform
@@ -675 +682 @@
     // 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;
-        }
+        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;
+        }
     }
 
@@ -720 +727 @@
 
     for (int i = 0; i <= toTPA->x->nX; i++) {
-        for (int j = 0; j <= toTPA->x->nY; j++) {
-            double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
-            double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
-            toTPA->x->coeff[i][j] = f1 + f2;
-
-            double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
-            double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
-            toTPA->y->coeff[i][j] = g1 + g2;
-        }
+        for (int j = 0; j <= toTPA->x->nY; j++) {
+            double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
+            double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
+            toTPA->x->coeff[i][j] = f1 + f2;
+
+            double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
+            double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
+            toTPA->y->coeff[i][j] = g1 + g2;
+        }
     }
     toTPA->x->coeff[0][0] += fpa->toTPA->x->coeff[0][0];
@@ -752 +759 @@
     psPlane *center = psPlaneTransformGetCenter (tpa1, tol);
     if (!center) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
-        psFree (toTPA);
-        psFree (tpa1);
-        psFree (wcs);
-        return NULL;
+        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
+        psFree (toTPA);
+        psFree (tpa1);
+        psFree (wcs);
+        return NULL;
     }
 
@@ -789 +796 @@
     // 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->nY; j++) {
-            wcs->trans->x->coeff[i][j] *= pdelt1;
-            wcs->trans->y->coeff[i][j] *= pdelt2;
-        }
+        for (int j = 0; j <= wcs->trans->x->nY; j++) {
+            wcs->trans->x->coeff[i][j] *= pdelt1;
+            wcs->trans->y->coeff[i][j] *= pdelt2;
+        }
     }
 
@@ -811 +818 @@
 */
 
-// convert the chip-level toFPA to a wcs polynomial transformation
+// convert the chip-level toFPA to a wcs polynomial transformation.  the pmAstromWCS
+// structure represents a single layer transformation (e.g., RA-TAN, RA-WRP).  Here we are
+// converting the chip-level to a WRP projection in the structure.  Later, this will be
+// converted to the WCS keywords
+
 pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, double tol)
 {
@@ -822 +833 @@
     pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
 
+    // copy the toFPA x,y, transformations to the wcs version
+    wcs->trans->x = psPolynomial2DCopy (wcs->trans->x, chip->toFPA->x);
+    wcs->trans->y = psPolynomial2DCopy (wcs->trans->y, chip->toFPA->y);
+
     // 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);
-    wcs->crval1 = 0.0;
-    wcs->crval2 = 0.0;
-
-    // given transformation, solve for coordinates which yields output coordinates of 0,0
-    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
-    if (!center) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
-        psFree (wcs);
-        return NULL;
-    }
-
-    // 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);
-
-    // calculated center is crpix1,2
-    wcs->crpix1 = center->x;
-    wcs->crpix2 = center->y;
-    psFree (center);
-
+
+    // reference pixel (CRPIX1,2) is (0.0, 0.0):
+    wcs->crpix1 = 0.0;
+    wcs->crpix2 = 0.0;
+
+    // we need to set CRVAL1,2 for the 0,0 pixel:
+    wcs->crval1 = psPolynomial2DEval (chip->toFPA->x, 0.0, 0.0);
+    wcs->crval2 = psPolynomial2DEval (chip->toFPA->y, 0.0, 0.0);
+
+    wcs->toSky->R = wcs->crval1*PM_RAD_DEG;
+    wcs->toSky->D = wcs->crval2*PM_RAD_DEG;
+
+    // these need to be set to 0.0 since they have been moved to crpix1,crpix2
+    wcs->trans->x->coeff[0][0] = 0.0;
+    wcs->trans->y->coeff[0][0] = 0.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]);
@@ -872 +883 @@
     psPlane *center = psPlaneTransformGetCenter (fpa->toTPA, tol);
     if (!center) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
-        psFree (wcs);
-        return NULL;
+        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
+        psFree (wcs);
+        return NULL;
     }
 
@@ -892 +903 @@
     // 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;
-        }
+        for (int j = 0; j <= wcs->trans->x->nY; j++) {
+            wcs->trans->x->coeff[i][j] *= pdelt1;
+            wcs->trans->y->coeff[i][j] *= pdelt2;
+        }
     }
 
@@ -918 +929 @@
     int k=0;
     for (int j=0; j<nSamples; j++) {
-        double y = bounds->y0 + (j * deltaY / nSamples);
-        for (int i=0; i<nSamples; i++) {
-            psPlane *s = psPlaneAlloc();
-            s->x = bounds->x0 + (i * deltaX / nSamples);
-            s->y = y;
-            psArraySet(src, k, s);
-            psPlane *d = psPlaneTransformApply(NULL, trans, s);
-            psArraySet(dst, k, d);
-            psFree(s);  // drop our refs to s and d
-            psFree(d);
-            ++k;
-        }
+        double y = bounds->y0 + (j * deltaY / nSamples);
+        for (int i=0; i<nSamples; i++) {
+            psPlane *s = psPlaneAlloc();
+            s->x = bounds->x0 + (i * deltaX / nSamples);
+            s->y = y;
+            psArraySet(src, k, s);
+            psPlane *d = psPlaneTransformApply(NULL, trans, s);
+            psArraySet(dst, k, d);
+            psFree(s);  // drop our refs to s and d
+            psFree(d);
+            ++k;
+        }
     }
 
@@ -935 +946 @@
 
     if (!psPlaneTransformFit(newTrans, src, dst, 0, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
-        return NULL;
+        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
+        return NULL;
     }
 
@@ -945 +956 @@
     printf("   i     chip_x  tpa_x     tpa_x_fit     dx         chip_y    tpa_y     tpa_y_fit     dy     dx > 0.5 || dy > 0.5\n");
     for (int i=0; i<psArrayLength(dst); i++) {
-        psPlane *d = (psPlane *) psArrayGet(dst, i);
-        psPlane *s = (psPlane *) psArrayGet(src, i);
-
-        new = psPlaneTransformApply(new, newTrans, s);
-
-        double xerr = new->x - d->x;
-        double yerr = new->y - d->y;
-        bool bigerr = (fabs(xerr) > .5) || (fabs(yerr) > .5);
-        printf("%4d %9.2f %9.2f %9.2f %9.4f     %9.2f %9.2f %9.2f %9.4f   %s\n"
-        , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : "");
+        psPlane *d = (psPlane *) psArrayGet(dst, i);
+        psPlane *s = (psPlane *) psArrayGet(src, i);
+
+        new = psPlaneTransformApply(new, newTrans, s);
+
+        double xerr = new->x - d->x;
+        double yerr = new->y - d->y;
+        bool bigerr = (fabs(xerr) > .5) || (fabs(yerr) > .5);
+        printf("%4d %9.2f %9.2f %9.2f %9.4f     %9.2f %9.2f %9.2f %9.4f   %s\n"
+               , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : "");
     }
     psFree(new);
@@ -972 +983 @@
 
     if (outFPA == NULL) {
-        outFPA = inFPA;
+        outFPA = inFPA;
     }
     if (outChip == NULL) {
-        outChip = inChip;
+        outChip = inChip;
     }
     if (outputBounds == NULL) {
-        outputBounds = pmChipPixels(outChip);
+        outputBounds = pmChipPixels(outChip);
     }
 
@@ -984 +995 @@
     psPlaneTransform *chipToTPA = psPlaneTransformCombine(NULL, inChip->toFPA, inFPA->toTPA, *outputBounds, 50);
     if (!chipToTPA) {
-        psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA");
-        return false;
+        psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA");
+        return false;
     }
 
@@ -992 +1003 @@
     psFree(chipToTPA);
     if (!chipToFPA) {
-        psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed");
-        return false;
+        psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed");
+        return false;
     }
 
@@ -999 +1010 @@
     psPlaneTransform *outToFPA;
     if (offset_x != 0. && offset_y != 0.) {
-        outToFPA = psPlaneTransformSetCenter(NULL, chipToFPA, offset_x, offset_y);
-        psFree(chipToFPA);
+        outToFPA = psPlaneTransformSetCenter(NULL, chipToFPA, offset_x, offset_y);
+        psFree(chipToFPA);
     } else {
-        outToFPA = chipToFPA;
+        outToFPA = chipToFPA;
     }
 
     psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 50);
     if (!outFromFPA) {
-        psFree(outToFPA);
-        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
-        return false;
+        psFree(outToFPA);
+        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
+        return false;
     }
 
@@ -1061 +1072 @@
 
     for (int j = 0; j < nSamples; j++) {
-        double y = bounds->y0 + (j * deltaY / nSamples);
-        for (int i =  0; i < nSamples; i++) {
-
-            psSphere srcSky;
-            psPlane *srcChip = psPlaneAlloc();
-            psPlane *dstTP = psPlaneAlloc();
-
-            srcChip->x = bounds->x0 + (i * deltaX / nSamples);
-            srcChip->y = y;
-
-            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
-            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
-            psDeproject (&srcSky, &srcTP, inFPA->toSky);
-
-            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
-
-            psProject (dstTP, &srcSky, outFPA->toSky);
-
-            srcChip->x -= bounds->x0;
-            srcChip->y -= bounds->y0;
-            psArrayAdd (src, 100, srcChip);
-            psArrayAdd (dst, 100, dstTP);
-
-            psFree(srcChip);  // drop our refs to s and d
-            psFree(dstTP);
-        }
+        double y = bounds->y0 + (j * deltaY / nSamples);
+        for (int i =  0; i < nSamples; i++) {
+
+            psSphere srcSky;
+            psPlane *srcChip = psPlaneAlloc();
+            psPlane *dstTP = psPlaneAlloc();
+
+            srcChip->x = bounds->x0 + (i * deltaX / nSamples);
+            srcChip->y = y;
+
+            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
+            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
+            psDeproject (&srcSky, &srcTP, inFPA->toSky);
+
+            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
+
+            psProject (dstTP, &srcSky, outFPA->toSky);
+
+            srcChip->x -= bounds->x0;
+            srcChip->y -= bounds->y0;
+            psArrayAdd (src, 100, srcChip);
+            psArrayAdd (dst, 100, dstTP);
+
+            psFree(srcChip);  // drop our refs to s and d
+            psFree(dstTP);
+        }
     }
 
@@ -1094 +1105 @@
 
     if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
-        psFree(src);
-        psFree(dst);
-        return NULL;
+        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
+        psFree(src);
+        psFree(dst);
+        return NULL;
     }
 
@@ -1103 +1114 @@
     for (int i = 0; i < src->n; i++) {
 
-        psSphere srcSky, dstSky;
-        psPlane *srcChip = src->data[i];
-        psPlane *dstTP   = dst->data[i];
-
-        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
-        psDeproject (&srcSky, &srcFP, outFPA->toSky);
-        psDeproject (&dstSky, dstTP, outFPA->toSky);
-
-        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
-        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
-        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
+        psSphere srcSky, dstSky;
+        psPlane *srcChip = src->data[i];
+        psPlane *dstTP   = dst->data[i];
+
+        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
+        psDeproject (&srcSky, &srcFP, outFPA->toSky);
+        psDeproject (&dstSky, dstTP, outFPA->toSky);
+
+        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
+        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
+        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
 
     }
@@ -1124 +1135 @@
     psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *bounds, 1);
     if (!newFromFPA) {
-        psFree(newToFPA);
-        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
-        return false;
+        psFree(newToFPA);
+        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
+        return false;
     }
 
@@ -1149 +1160 @@
 
     if (!wcs)
-        return;
+        return;
     psFree (wcs->trans);
     psFree (wcs->toSky);
@@ -1171 +1182 @@
 /*****
 
-For mosaic astrometry, we need to have a starting set of projection terms in which the
-chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
-microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
-result in substantial errors between the observed and predicted star positions on the focal
-plane: this is the measurement of the optical distortion in the camera.  At the same time,
-we need to carry around the transformations which allow us to make an accurate calculation
-of the position of the stars based on the input (per-chip) astrometry.  These
-transformations will allow us to match the raw and ref stars robustly.  To convert the
-per-chip astrometry (which may have been calculated with a different plate scale for each
-chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
-the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
-each chip (eg, cdelt1).  Thus, we need to carry around both the
+      For mosaic astrometry, we need to have a starting set of projection terms in which the
+      chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
+      microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
+      result in substantial errors between the observed and predicted star positions on the focal
+      plane: this is the measurement of the optical distortion in the camera.  At the same time,
+      we need to carry around the transformations which allow us to make an accurate calculation
+      of the position of the stars based on the input (per-chip) astrometry.  These
+      transformations will allow us to match the raw and ref stars robustly.  To convert the
+      per-chip astrometry (which may have been calculated with a different plate scale for each
+      chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
+      the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
+      each chip (eg, cdelt1).  Thus, we need to carry around both the
 
 *****/