# include "psastroInternal.h"

# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\
  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
  if (!status) { \
   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
   return false; } 

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"); 

    // 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"); 
    RADIUS *= pixelScale;

    // select the desired chip order
    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32, "failed to find single-chip fit order\n");

    // 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");

    // 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");

    if (gridSearch) {
	// find initial offset / rotation
	gridStats = pmAstromGridMatch (rawstars, refstars, recipe);
	if (gridStats == NULL) {
	    psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
	    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 (rawstars, refstars, recipe, gridStats);
	if (stats == NULL) {
	    psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
	    psFree (gridStats);
	    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);
    }

    // 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;
    }
    
    // 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

    // 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) {
	psastroDumpStars (refstars, "refstars.cal.dat");
    }

    if (psTraceGetLevel("psastro.plot") > 0) {
	psastroPlotOneChipFit (rawstars, refstars, match, results);
    }

    psFree (match);
    psFree (stats);
    psFree (results);
    psFree (fitStats);
    psFree (gridStats);
    return validSolution;
}

// psastroWriteStars ("raw.1.dat", rawstars);
// psastroWriteStars ("ref.1.dat", refstars);
// psastroWriteTransform (chip->toFPA);

