Index: trunk/psastro/src/psastroMosaicAstrom.c
===================================================================
--- trunk/psastro/src/psastroMosaicAstrom.c	(revision 10815)
+++ trunk/psastro/src/psastroMosaicAstrom.c	(revision 10830)
@@ -1,8 +1,10 @@
 # include "psastro.h"
+# define NONLIN_TOL 0.001 /* tolerance in pixels */
 
+// XXX require this fpa to have multiple chip extensions and a PHU?
 bool psastroMosaicAstrom (pmConfig *config, psArray *refs) {
 
     bool status;
-    psArray *grads = NULL;
+    psArray *gradients = NULL;
 
     // select the current recipe
@@ -20,44 +22,71 @@
     }
 
-    int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.DISTORTION.ORDER");
-    if (!status) order = 3.0;
-
     pmFPA *fpa = input->fpa;
 
-    psastroMosaicSetMatch (fpa, recipe);
+    // XXX before we do object matches, we need to fix failed chips
+    // compare chips with supplied mosaic model
+    // adjust significant outliers to match model
 
-    // XXX no input distortion model yet; existing fpa distortion is identity: replace it
-    // XXX make this a test if these are NULL; the load step can NULL them if nothing is loaded
-    // XXX if we have an input distortion model, the gradient is measuring the error rel to that model
-    // XXX how does this couple to the individual chip fits?
+    // given the existing per-chip astrometry, determine matches between raw and ref stars
+    psastroMosaicSetMatch (fpa, recipe, 0);
+
+    // fitted chips will follow the local plate-scale, hiding the distortion
+    // modify the chip->toFPA scaling to match knowledge about pixel scale
+    psastroMosaicCommonScale (fpa, recipe);
+
+    gradients = psastroMosaicGradients (fpa, recipe);
+
+    // allocate mosaic-level polynomial transformation and set masks needed by DVO
+    int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
+    if (!status) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
+        return false;
+    }
     psFree (fpa->toTPA);
-    psFree (fpa->fromTPA);
-    fpa->toTPA   = psPlaneDistortIdentity (order);
-    fpa->fromTPA = psPlaneDistortIdentity (order);
+    fpa->toTPA = psPlaneTransformAlloc (order, order);
+    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
+        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
+            if (i + j > order) {
+		fpa->toTPA->x->mask[i][j] = 1;
+		fpa->toTPA->y->mask[i][j] = 1;
+            }
+        }
+    }
 
-    grads = psastroMosaicGetGrads (fpa, recipe);
-
-    // fit the measured gradients with the telescope distortion model (3rd order polynomial)
-    pmAstromFitDistortion (fpa, grads, recipe);
-    psastroMosaicRescaleChips (fpa);
+    // fit the measured gradients with the telescope distortion model (polynomial order based on toTPA)
+    pmAstromFitDistortion (fpa, gradients, recipe);
+    psFree (gradients);
 
     // apply the new distortion terms up and down
     // re-perform the match with a slightly tighter circle
+    // XXX modify match radius
+    // XXX set chip.order to 1
     psastroMosaicSetAstrom (fpa);
-    psastroMosaicSetMatch (fpa, recipe);
+    psastroMosaicSetMatch (fpa, recipe, 1);
+    psastroMosaicChipAstrom (fpa, recipe, 1);
 
-    // fit the chips with linear fits
-    // re-match
-    psastroMosaicChipAstrom (fpa, recipe);
-    psastroMosaicSetMatch (fpa, recipe);
+    // do a second pass on the distortion with improved chip positions and rotations
+    // XXX do we need to iterate with this step until the distortions don't change? 
+    psastroMosaicCommonScale (fpa, recipe);
+    gradients = psastroMosaicGradients (fpa, recipe);
+    pmAstromFitDistortion (fpa, gradients, recipe);
+    psFree (gradients);
+    
+    // now fit the chips under the common distortion with higher-order terms
+    // XXX modify match radius
+    // XXX set chip.order to NORDER
+    psastroMosaicSetAstrom (fpa);
+    psastroMosaicSetMatch (fpa, recipe, 2);
+    psastroMosaicChipAstrom (fpa, recipe, 2);
 
-    // re-fit the chips with higher-order fits
-    // re-match
-    psastroMosaicChipAstrom (fpa, recipe);
-    psastroMosaicSetMatch (fpa, recipe);
+    // save WCS and analysis metadata in update header
+    // XXX need to add global summary statistics
+    psMetadata *updates = psMetadataAlloc();
+    pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL);
+    psMetadataAdd (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
+    psFree (updates);
 
-    // re-fit the chips with higher-order fits
-    psastroMosaicChipAstrom (fpa, recipe);
-    psastroMosaicHeaders (config);
+    // update the headers based on the results
+    // psastroMosaicHeaders (config);
 
     return true;
@@ -65,9 +94,9 @@
 
 /* coordinate frame hierachy
-   pixels (on a given readout)
-   cell
-   chip
-   FP (focal plane)
-   TP (tangent plane)
-   sky (ra, dec)
-*/
+ * pixels (on a given readout)
+ * cell
+ * chip
+ * FP (focal plane)
+ * TP (tangent plane)
+ * sky (ra, dec)
+ */