Index: /tags/ipp-20170121/ippconfig/gpc1/pswarp.config
===================================================================
--- /tags/ipp-20170121/ippconfig/gpc1/pswarp.config	(revision 39981)
+++ /tags/ipp-20170121/ippconfig/gpc1/pswarp.config	(revision 39982)
@@ -3,2 +3,3 @@
 BKG.XGRID		S32	50
 BKG.YGRID		S32	50
+ADDITIONAL_WCS_ORDERS   S32     4               # By default, use a bilevel astrometric model with 4 additional orders.
Index: /tags/ipp-20170121/ippconfig/recipes/pswarp.config
===================================================================
--- /tags/ipp-20170121/ippconfig/recipes/pswarp.config	(revision 39981)
+++ /tags/ipp-20170121/ippconfig/recipes/pswarp.config	(revision 39982)
@@ -25,4 +25,6 @@
 BKG.YGRID              S32    4
 
+ADDITIONAL_WCS_ORDERS S32     4              # By default, use a bilevel astrometric model with 4 additional orders.
+
 # Default recipe for warping
 WARP	METADATA
Index: /tags/ipp-20170121/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- /tags/ipp-20170121/psModules/src/astrom/pmAstrometryWCS.c	(revision 39981)
+++ /tags/ipp-20170121/psModules/src/astrom/pmAstrometryWCS.c	(revision 39982)
@@ -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;
+        }
     }
 
@@ -522,78 +522,78 @@
 
     if (fpa->toSky == NULL) {
-	psFree(fpa->toTPA);
-	psFree(fpa->fromTPA);
-	fpa->toTPA = psPlaneTransformIdentity (1);
-	fpa->fromTPA = psPlaneTransformIdentity (1);
-	fpa->toSky = toSky;
+        psFree(fpa->toTPA);
+        psFree(fpa->fromTPA);
+        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;
@@ -649,13 +649,17 @@
 
     psFree (chip->toFPA);
-    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;
+    if ((fabs(wcs->crpix1) > 0.01) || (fabs(wcs->crpix2) > 0.01)) {
+      chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
+    } else {
+      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
@@ -682,8 +686,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;
+        }
     }
 
@@ -691,5 +695,13 @@
     // the region defines the FPA pixels covered by the tranformation
     psFree (fpa->fromTPA);
-    psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4);
+
+    int additional_orders = 4;  // This is the number of orders that should be added.
+    bool status = false;
+    int config_additional_orders = psMetadataLookupS32(&status,fpa->analysis, "ADDITIONAL_WCS_ORDERS");
+    if (status) { // If this keyword is found, then use that number of orders instead
+      additional_orders = config_additional_orders; 
+    }
+    
+    psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+additional_orders, fpa->toTPA->x->nY+additional_orders);
     fpa->fromTPA = psPlaneTransformInvert(myPT, fpa->toTPA, region, 50);
     psFree (myPT);
@@ -729,13 +741,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];
@@ -761,9 +773,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;
     }
 
@@ -798,8 +810,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;
+        }
     }
 
@@ -885,7 +897,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;
     }
 
@@ -905,8 +917,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;
+        }
     }
 
@@ -931,16 +943,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;
+        }
     }
 
@@ -948,6 +960,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;
     }
 
@@ -958,14 +970,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);
@@ -985,11 +997,11 @@
 
     if (outFPA == NULL) {
-	outFPA = inFPA;
+        outFPA = inFPA;
     }
     if (outChip == NULL) {
-	outChip = inChip;
+        outChip = inChip;
     }
     if (outputBounds == NULL) {
-	outputBounds = pmChipPixels(outChip);
+        outputBounds = pmChipPixels(outChip);
     }
 
@@ -997,6 +1009,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;
     }
 
@@ -1005,6 +1017,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;
     }
 
@@ -1012,15 +1024,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;
     }
 
@@ -1074,30 +1086,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);
+        }
     }
 
@@ -1107,8 +1119,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;
     }
 
@@ -1116,15 +1128,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);
 
     }
@@ -1137,7 +1149,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;
     }
 
@@ -1162,5 +1174,5 @@
 
     if (!wcs)
-	return;
+        return;
     psFree (wcs->trans);
     psFree (wcs->toSky);
Index: /tags/ipp-20170121/pswarp/src/pswarpLoadAstrometry.c
===================================================================
--- /tags/ipp-20170121/pswarp/src/pswarpLoadAstrometry.c	(revision 39981)
+++ /tags/ipp-20170121/pswarp/src/pswarpLoadAstrometry.c	(revision 39982)
@@ -96,4 +96,14 @@
     // apply the bilevel astrometry elements to the target
     if (bilevelAstrometry) {
+      int additional_orders = 4;
+      bool status = false;
+      int config_additional_orders = psMetadataLookupS32(&status,config->arguments,"ADDITIONAL_WCS_ORDERS");
+      if (status) {
+	additional_orders = config_additional_orders;
+      }
+      psMetadataAddS32(astrom->fpa->analysis,PS_LIST_TAIL, "ADDITIONAL_WCS_ORDERS", PS_META_REPLACE,
+		       "Additional orders for bilevel fit.", additional_orders);
+      
+	
 	if (!pmAstromReadBilevelMosaic(astrom->fpa, phu->header)) {
 	    psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell.");
Index: /tags/ipp-20170121/pswarp/src/pswarpOptions.c
===================================================================
--- /tags/ipp-20170121/pswarp/src/pswarpOptions.c	(revision 39981)
+++ /tags/ipp-20170121/pswarp/src/pswarpOptions.c	(revision 39982)
@@ -92,4 +92,12 @@
     }
 
+    // See if we should use a more accurate WCS model.
+    int config_additional_orders = psMetadataLookupS32(&status,recipe, "ADDITIONAL_WCS_ORDERS");
+    if (!status) {
+      // We did not find this recipe option.  We are likely updating an old config.
+      config_additional_orders = 0;
+    }
+    
+
     
     // Set recipe values in the recipe (since we've possibly altered some)
@@ -114,4 +122,5 @@
     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid);
     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid);
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "ADDITIONAL_WCS_ORDERS", PS_META_REPLACE, "Additional orders for bilevel fit.", config_additional_orders);
 
     psTrace("pswarp", 1, "Done with pswarpOptions...\n");