Index: trunk/psastro/src/psastro.c
===================================================================
--- trunk/psastro/src/psastro.c	(revision 41284)
+++ trunk/psastro/src/psastro.c	(revision 41285)
@@ -51,4 +51,6 @@
     }
 
+    psLogMsg("psastro", 3, "TIMEMARK: psastroDataLoad: %f sec\n", psTimerMark ("complete"));
+
     psMetadata *stats = psMetadataAlloc(); // Statistics, for output
     psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
@@ -61,4 +63,6 @@
         exit(PS_EXIT_SYS_ERROR);
     }
+
+    psLogMsg("psastro", 3, "TIMEMARK: psastroAnalysis: %f sec\n", psTimerMark ("complete"));
 
     // write out the results
Index: trunk/psastro/src/psastroAnalysis.c
===================================================================
--- trunk/psastro/src/psastroAnalysis.c	(revision 41284)
+++ trunk/psastro/src/psastroAnalysis.c	(revision 41285)
@@ -135,4 +135,7 @@
         }
     }
+
+    psLogMsg("psastro", 3, "TIMEMARK: psastroChipAstrom: %f sec\n", psTimerMark ("complete"));
+
     if (mosastro) {
         if (!psastroMosaicAstrom (config)) {
@@ -150,4 +153,6 @@
         }
     }
+
+    psLogMsg("psastro", 3, "TIMEMARK: psastroMosaicAstrom: %f sec\n", psTimerMark ("complete"));
 
     if (!skipastro) {
Index: trunk/psastro/src/psastroChipAstrom.c
===================================================================
--- trunk/psastro/src/psastroChipAstrom.c	(revision 41284)
+++ trunk/psastro/src/psastroChipAstrom.c	(revision 41285)
@@ -28,4 +28,13 @@
     }
 
+    // In this function, we are fitting the chip->FPA model, keeping the toTPA static.  
+    // If the toTPA model has order > 1, then the WCS terms need to use the BiLevel
+    // model.  It is possible for this conversion to fail.  However, if we are going to try
+    // again with the mosaic model (fitting both), then we need to react to failure differently
+    bool mosastro  = psMetadataLookupBool (&status, config->arguments, "PSASTRO.MOSAIC.MODE");
+    if (!status) {
+        mosastro  = psMetadataLookupBool (&status, recipe, "PSASTRO.MOSAIC.MODE");
+    }
+
     // select the input data sources
     pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT");
@@ -39,4 +48,6 @@
 
     int numGoodChips = 0;               // Number of chips for which astrometry succeeds
+
+    bool useBilevelWCS = (fpa->toTPA->x->nX > 1) || (fpa->toTPA->x->nY > 1);
 
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
@@ -145,12 +156,23 @@
                 psU64 astrom_chip_val = 1;
                 astrom_chip_val <<= view->chip;
-                psMetadataAddU64 (updates, PS_LIST_TAIL, "ASTROM_CHIPS", PS_META_REPLACE, 
-                    "chips that passed astrometry", astrom_chip_val);
+                psMetadataAddU64 (updates, PS_LIST_TAIL, "ASTROM_CHIPS", PS_META_REPLACE, "chips that passed astrometry", astrom_chip_val);
 
                 // write the elapsed time here; this will be updated in psastroMosaicAstrometry, if called
                 psMetadataAddF32 (updates, PS_LIST_TAIL, "DT_ASTR", PS_META_REPLACE, "elapsed psastro time", psTimerMark ("psastroAnalysis"));
 		
-		fpa->wcsCDkeys = psMetadataLookupBool(&status, recipe , "PSASTRO.WCS.USECDKEYS");
-                pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL);
+		// if the toTPA distortion model uses higher-order terms, then we need to use BiLevel astrometry 
+		if (useBilevelWCS) {
+		    // create the header keywords to descripe the results
+		    if (!pmAstromWriteBilevelChip (updates, chip, NONLIN_TOL)) {
+			if (!mosastro) { readout->data_exists = false; } // unless we are going to try again, give up on this chip
+			psError(PS_ERR_UNKNOWN, false, "invalid solution for %d,%d,%d\n", view->chip, view->cell, view->readout);
+			psErrorStackPrint(stderr, "failure to generate WCS keywords for one chip\n");
+			psErrorClear();
+			continue;
+		    }
+		} else {
+		    fpa->wcsCDkeys = psMetadataLookupBool(&status, recipe , "PSASTRO.WCS.USECDKEYS");
+		    pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL);
+		}
 
                 if (psTraceGetLevel("psastro.dump") > 0) {
@@ -179,4 +201,19 @@
     }
 
+    if (useBilevelWCS) {
+	// save WCS and analysis metadata in update header.
+	// (pull or create local view to entry on readout->analysis)
+	psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
+	if (!updates) {
+	    updates = psMetadataAlloc ();
+	    psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
+	    psFree (updates);
+	}
+	if (!pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL)) {
+	    psError(psErrorCodeLast(), false, "Failed to save header terms");
+	    return false;
+	}
+    }
+
     if (fpa->chips->n == 1 && numGoodChips == 0) {
         psError(PSASTRO_ERR_UNKNOWN, false, "Failed to fit single chip.");
Index: trunk/psastro/src/psastroChooseRefstars.c
===================================================================
--- trunk/psastro/src/psastroChooseRefstars.c	(revision 41284)
+++ trunk/psastro/src/psastroChooseRefstars.c	(revision 41285)
@@ -46,4 +46,10 @@
 
     bool matchLumFunc = psMetadataLookupBool (&status, recipe, "PSASTRO.MATCH.LUMFUNC");
+
+    int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.NITER");
+    if (!status) nIter = 3;
+
+    psF32 clumpScale = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.SCALE");
+    if (!status) clumpScale = 150;
 
     pmFPAview *view = pmFPAviewAlloc (0);
@@ -178,5 +184,5 @@
 		// generate a reduced subset excluding the clumps
 		// XXX do we need both REFSTARS and SUBSET? 
-		psArray *subset = psastroRemoveClumpsIterate(refstars, 150, 3);
+		psArray *subset = psastroRemoveClumpsIterate(refstars, clumpScale, nIter);
 		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset);
 
Index: trunk/psastro/src/psastroOneChipFit.c
===================================================================
--- trunk/psastro/src/psastroOneChipFit.c	(revision 41284)
+++ trunk/psastro/src/psastroOneChipFit.c	(revision 41285)
@@ -36,5 +36,5 @@
 
     // select the desired chip order
-    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32);
+    REQUIRED_RECIPE_VALUE (int defaultOrder, "PSASTRO.CHIP.ORDER", S32);
 
     // allowed limits for valid solutions
@@ -57,4 +57,9 @@
         }
 
+        sprintf (name, "PSASTRO.ONE.CHIP.ORDER.N%d", iter);
+        int order = psMetadataLookupS32 (&status, recipe, name);
+        if (!status) {
+            order = defaultOrder;
+        }
 
         // use small radius to match stars
@@ -76,4 +81,15 @@
 	    match = unique;
 	}
+
+	// XXX check if we correctly applied the new transformation:
+	if (psTraceGetLevel("psastro.dump") > 0) {
+	  char *filename = NULL;
+	  char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
+	  psStringAppend (&filename, "match.pref.%s.%d.dat", chipname, iter);
+	  psastroDumpMatchedStars (filename, rawstars, refstars, match);
+	  psFree (filename);
+	  filename = NULL;
+	}
+
 
         // modify the order to correspond to the actual number of matched stars:
@@ -137,4 +153,14 @@
         psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS
 
+	// XXX check if we correctly applied the new transformation:
+	if (psTraceGetLevel("psastro.dump") > 0) {
+	  char *filename = NULL;
+	  char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
+	  psStringAppend (&filename, "match.post.%s.%d.dat", chipname, iter);
+	  psastroDumpMatchedStars (filename, rawstars, refstars, match);
+	  psFree (filename);
+	  filename = NULL;
+	}
+
         // toSky converts from FPA & TPA units (microns) to sky units (radians)
         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
Index: trunk/psastro/src/psastroRemoveClumps.c
===================================================================
--- trunk/psastro/src/psastroRemoveClumps.c	(revision 41284)
+++ trunk/psastro/src/psastroRemoveClumps.c	(revision 41285)
@@ -35,4 +35,10 @@
     }
 
+    int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.NITER");
+    if (!status) nIter = 3;
+
+    psF32 clumpScale = psMetadataLookupS32 (&status, recipe, "PSASTRO.REFSTAR.CLUMP.SCALE");
+    if (!status) clumpScale = 150;
+
     pmFPAview *view = pmFPAviewAlloc (0);
     pmFPA *fpa = input->fpa;
@@ -58,5 +64,5 @@
 		// XXX do we need both RAWSTARS and SUBSET? 
 		// XXX put these parameters in the recipe, please
-		psArray *subset = psastroRemoveClumpsIterate(rawstars, 150, 3);
+		psArray *subset = psastroRemoveClumpsIterate(rawstars, clumpScale, nIter);
 		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset);
 		psFree (subset);
@@ -65,5 +71,5 @@
                 if ((gridstars == rawstars) || (gridstars == NULL)) { continue; }
 
-		psArray *gridstars_subset = psastroRemoveClumpsIterate(gridstars, 150, 3);
+		psArray *gridstars_subset = psastroRemoveClumpsIterate(gridstars, clumpScale, nIter);
 		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.GRID.RAWSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", gridstars_subset);
 		psFree (gridstars_subset);
@@ -79,5 +85,5 @@
     psArray *newset = psMemIncrRefCounter (input);
     for (int i = 0; i < nIter; i++) {
-	psArray *subset = psastroRemoveClumps (newset, 150);
+	psArray *subset = psastroRemoveClumps (newset, scale);
 	psFree (newset);
 	newset = subset;