Index: /trunk/psastro/src/psastroOneChipFit.c
===================================================================
--- /trunk/psastro/src/psastroOneChipFit.c	(revision 16070)
+++ /trunk/psastro/src/psastroOneChipFit.c	(revision 16070)
@@ -0,0 +1,166 @@
+# include "psastroInternal.h"
+
+# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE)\
+  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
+  if (!status) { \
+   psAbort ("Failed to find %s in recipe", NAME); }
+
+bool psastroOneChipFit (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {
+
+    bool status;
+
+    // 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); 
+    RADIUS *= pixelScale;
+
+    // select the desired chip order
+    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32);
+
+    // allowed limits for valid solutions
+    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32);
+    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32);
+
+    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");
+	    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);
+	    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->coeffMask[i][j] = PS_POLY_MASK_SET;
+		    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
+		}
+	    }
+	}
+
+	// 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 (fitStats);
+	    return false;
+	}
+    
+	// determine fromFPA transformation and apply new transformation to raw & ref stars
+	psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    
+	// 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);
+
+	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;
+
+    // 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;
+
+    // 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;
+    int astNstar = results->yStats->clippedNvalues;
+
+    bool validSolution = true;
+
+    // XXX should these result in errors or be handled another way?
+    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
+    if (astError > maxError) {
+        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
+	validSolution = false;
+    }
+    if (astNstar < minNstar) {
+        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
+	validSolution = false;
+    }
+
+    // DVO expects NASTRO = 0 if we fail to find a solution
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
+    if (validSolution) {
+	psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
+	psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
+    } else {
+	psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
+	psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
+    }
+    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
+
+    // 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);
+	psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
+	psastroDumpStars (refstars, "refstars.cal.dat");
+    }
+
+    if (psTraceGetLevel("psastro.plot") > 0) {
+	psastroPlotOneChipFit (rawstars, refstars, match, recipe);
+    }
+
+    psFree (match);
+    psFree (results);
+    psFree (fitStats);
+    return validSolution;
+}
+
+// psastroWriteStars ("raw.1.dat", rawstars);
+// psastroWriteStars ("ref.1.dat", refstars);
+// psastroWriteTransform (chip->toFPA);
+
Index: /trunk/psastro/src/psastroOneChipGrid.c
===================================================================
--- /trunk/psastro/src/psastroOneChipGrid.c	(revision 16070)
+++ /trunk/psastro/src/psastroOneChipGrid.c	(revision 16070)
@@ -0,0 +1,64 @@
+# include "psastroInternal.h"
+
+# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE)\
+  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
+  if (!status) { \
+   psAbort ("Failed to find %s in recipe", NAME); }
+
+bool psastroOneChipGrid (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {
+
+    bool status;
+    pmAstromStats *stats = NULL;
+
+    // do we need to get a rough initial match?
+    REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool);
+    if (!gridSearch) return true;
+
+    // do we need to get a rough initial match?
+    REQUIRED_RECIPE_VALUE (int maxNstar, "PSASTRO.GRID.NSTAR.MAX", S32);
+
+    // generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N)
+    psArray *subset = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
+    for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) {
+	subset->data[i] = psMemIncrRefCounter (rawstars->data[i]);
+    }
+
+    // XXX set clump scale from recipe
+    psArray *gridStars = psastroRemoveClumps (subset, 150);
+    psFree (subset);
+
+    psArray *refSubset = psastroRemoveClumps (refstars, 150);
+
+    psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refSubset->n);
+
+    // find initial offset / rotation / scale
+    pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refSubset, recipe);
+    if (gridStats == NULL) {
+	psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
+	psFree (gridStars);
+	psFree (refSubset);
+	return false;
+    }
+    psLogMsg ("psastro", 3, "basic grid search result - offset: %f,%f pixels, rotation: %f deg\n", gridStats->offset.x, gridStats->offset.y, DEG_RAD*gridStats->angle);
+
+    // tweak the position by finding peak of matches stars
+    stats = pmAstromGridTweak (gridStars, refSubset, recipe, gridStats);
+    if (stats == NULL) {
+	psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
+	psFree (gridStats);
+	psFree (gridStars);
+	psFree (refSubset);
+	return false;
+    }
+    psLogMsg ("psastro", 3, "tweak grid search result - offset: %f,%f pixels, rotation: %f deg\n", stats->offset.x, stats->offset.y, DEG_RAD*stats->angle);
+
+    // adjust the chip.toFPA terms only
+    pmAstromGridApply (chip->toFPA, stats);
+    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    psFree (gridStats);
+    psFree (gridStars);
+    psFree (refSubset);
+    psFree (stats);
+
+    return true;
+}