Index: trunk/psastro/src/psastroOneChip.c
===================================================================
--- trunk/psastro/src/psastroOneChip.c	(revision 13165)
+++ trunk/psastro/src/psastroOneChip.c	(revision 15201)
@@ -1,38 +1,49 @@
 # include "psastroInternal.h"
 
-# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\
+# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE)\
   VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
   if (!status) { \
-   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
-   return false; } 
+   psAbort ("Failed to find %s in recipe", NAME); }
 
 bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {
 
     bool status;
-    pmAstromStats *gridStats = NULL;
     pmAstromStats *stats = NULL;
 
-    // supplied radius is in pixels
-    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32, "Failed to lookup matching radius"); 
+    // default value for match/fit : radius is in pixels
+    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
+
+    // run the match/fit sequence NITER times
+    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
 
     // correct radius to FP units (physical pixel scale in microns per pixel)
-    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); 
+    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
     RADIUS *= pixelScale;
 
     // select the desired chip order
-    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32, "failed to find single-chip fit order\n");
+    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32);
 
     // allowed limits for valid solutions
-    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32, "failed to find single-chip max allowed error\n");
-    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
+    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32);
+    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32);
 
     // do we need to get a rough initial match?
-    REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool, "failed to find chip grid-search option\n");
+    REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool);
+
+    // do we need to get a rough initial match?
+    REQUIRED_RECIPE_VALUE (int maxNstar, "PSASTRO.GRID.NSTAR.MAX", S32);
 
     if (gridSearch) {
+	// generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N)
+	psArray *gridStars = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
+	for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) {
+	    gridStars->data[i] = psMemIncrRefCounter (rawstars->data[i]);
+	}
+
 	// find initial offset / rotation
-	gridStats = pmAstromGridMatch (rawstars, refstars, recipe);
+	pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refstars, recipe);
 	if (gridStats == NULL) {
 	    psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
+	    psFree (gridStars);
 	    return false;
 	}
@@ -40,8 +51,9 @@
 
 	// tweak the position by finding peak of matches stars
-	stats = pmAstromGridTweak (rawstars, refstars, recipe, gridStats);
+	stats = pmAstromGridTweak (gridStars, refstars, recipe, gridStats);
 	if (stats == NULL) {
 	    psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
 	    psFree (gridStats);
+	    psFree (gridStars);
 	    return false;
 	}
@@ -51,66 +63,100 @@
 	pmAstromGridApply (chip->toFPA, stats);
 	psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
-    }
-
-    // use small radius to match stars
-    psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS);
-    if (match == NULL) {
-	psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
-	psFree (stats);
 	psFree (gridStats);
-        return false;
-    }
-
-    // modify the order to correspond to the actual number of matched stars:
-    if ((match->n < 11) && (order >= 3)) order = 2;
-    if ((match->n <  7) && (order >= 2)) order = 1;
-    if ((match->n <  4) && (order >= 1)) order = 0;
-    if (order < 1) {
-	psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 
-	psFree (match);
-	psFree (stats);
-	psFree (gridStats);
-	return false; 
-    } 
-
-    // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
-    psFree (chip->toFPA);
-    chip->toFPA = psPlaneTransformAlloc (order, order);
-    for (int i = 0; i <= chip->toFPA->x->nX; i++) {
-        for (int j = 0; j <= chip->toFPA->x->nY; j++) {
-            if (i + j > order) {
-		chip->toFPA->x->mask[i][j] = 1;
-		chip->toFPA->y->mask[i][j] = 1;
-            }
-        }
-    }
-
-    // XXX allow statitic to be set by the user
-    psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    // psStats *fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
-    fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
-
-    // improved fit for astrometric terms
-    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
-    if (!results) {
-	psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
-	psFree (match);
-	psFree (stats);
-	psFree (fitStats);
-	psFree (gridStats);
-        return false;
-    }
+	psFree (gridStars);
+    }
+
+    psArray *match = NULL;
+    psStats *fitStats = NULL;
+    pmAstromFitResults *results = NULL;
+
+    for (int iter = 0; iter < nIter; iter++) {
+	
+	char name[128];
+
+	sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
+	float radius = psMetadataLookupF32 (&status, recipe, name);
+	radius *= pixelScale;
+	if (!status || (radius == 0.0)) {
+	    radius = RADIUS;
+	}
+
+
+	// use small radius to match stars
+	match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
+	if (match == NULL) {
+	    psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
+	    psFree (stats);
+	    return false;
+	}
+
+	// modify the order to correspond to the actual number of matched stars:
+	if ((match->n < 11) && (order >= 3)) order = 2;
+	if ((match->n <  7) && (order >= 2)) order = 1;
+	if ((match->n <  4) && (order >= 1)) order = 0;
+	if (order < 1) {
+	    psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 
+	    psFree (match);
+	    psFree (stats);
+	    return false; 
+	} 
+
+	// create output toFPA; set masks appropriate to the Elixir DVO astrometry format
+	psFree (chip->toFPA);
+	chip->toFPA = psPlaneTransformAlloc (order, order);
+	for (int i = 0; i <= chip->toFPA->x->nX; i++) {
+	    for (int j = 0; j <= chip->toFPA->x->nY; j++) {
+		if (i + j > order) {
+		    chip->toFPA->x->mask[i][j] = 1;
+		    chip->toFPA->y->mask[i][j] = 1;
+		}
+	    }
+	}
+
+	// XXX allow statitic to be set by the user
+	// fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+	fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+	fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
+	fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
+
+	// improved fit for astrometric terms
+	results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
+	if (!results) {
+	    psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
+	    psFree (match);
+	    psFree (stats);
+	    psFree (fitStats);
+	    return false;
+	}
     
+	// determine fromFPA transformation and apply new transformation to raw & ref stars
+	psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    
+	if (iter < nIter - 1) {
+	    psFree (fitStats);
+	    psFree (results);
+	    psFree (match);
+	}
+
+	// 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;
+	// float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
+	float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
+	int astNstar = results->yStats->clippedNvalues;
+	psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
+    }
+
+
+
     // 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;
 
     // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
-    float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
-    // float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale;
+    // float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
+    float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale;
 
     // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
-    float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
-    // float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
+    // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
+    float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     int astNstar = results->yStats->clippedNvalues;
 
@@ -140,17 +186,16 @@
     psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
 
-    // determine fromFPA transformation and apply new transformation to raw & ref stars
-    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
+    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     
     // XXX check if we correctly applied the new transformation:
     if (psTraceGetLevel("psastro.dump") > 0) {
 	psastroDumpRawstars (rawstars, fpa, chip);
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) {
+	psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
 	psastroDumpStars (refstars, "refstars.cal.dat");
     }
 
     if (psTraceGetLevel("psastro.plot") > 0) {
-	psastroPlotOneChipFit (rawstars, refstars, match, results);
+	psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     }
 
@@ -159,5 +204,4 @@
     psFree (results);
     psFree (fitStats);
-    psFree (gridStats);
     return validSolution;
 }