# 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 psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {

    bool status;
    pmAstromStats *stats = NULL;

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

    // do we need to get a rough initial match?
    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 *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);
    }

    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->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 (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;

    // 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 (stats);
    psFree (results);
    psFree (fitStats);
    return validSolution;
}

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

