Index: /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c	(revision 39696)
+++ /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c	(revision 39697)
@@ -43,5 +43,5 @@
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     if (!wcs) {
-        return false;
+	return false;
     }
 
@@ -70,5 +70,5 @@
     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     if (!wcs) {
-        return false;
+	return false;
     }
 
@@ -86,6 +86,6 @@
     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,17 +96,17 @@
 
     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,6 +123,6 @@
     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,6 +139,6 @@
     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,9 +163,9 @@
 
     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,9 +188,9 @@
 
     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,6 +223,6 @@
     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,6 +232,6 @@
     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,25 +243,25 @@
 
     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,25 +280,25 @@
     // 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,59 +314,59 @@
     // 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,10 +402,10 @@
 
     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,80 +414,80 @@
     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,8 +507,8 @@
     // 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,76 +524,76 @@
 	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,26 +613,26 @@
     // 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,7 +648,14 @@
      */
 
-    // 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,8 +682,8 @@
     // 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,13 +727,13 @@
 
     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,9 +759,9 @@
     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,8 +796,8 @@
     // 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,5 +818,9 @@
 */
 
-// 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,27 +833,27 @@
     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,7 +883,7 @@
     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,8 +903,8 @@
     // 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,16 +929,16 @@
     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,6 +946,6 @@
 
     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,14 +956,14 @@
     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,11 +983,11 @@
 
     if (outFPA == NULL) {
-        outFPA = inFPA;
+	outFPA = inFPA;
     }
     if (outChip == NULL) {
-        outChip = inChip;
+	outChip = inChip;
     }
     if (outputBounds == NULL) {
-        outputBounds = pmChipPixels(outChip);
+	outputBounds = pmChipPixels(outChip);
     }
 
@@ -984,6 +995,6 @@
     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,6 +1003,6 @@
     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,15 +1010,15 @@
     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,30 +1072,30 @@
 
     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,8 +1105,8 @@
 
     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,15 +1114,15 @@
     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,7 +1135,7 @@
     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,5 +1160,5 @@
 
     if (!wcs)
-        return;
+	return;
     psFree (wcs->trans);
     psFree (wcs->toSky);
@@ -1171,16 +1182,16 @@
 /*****
 
-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
 
 *****/