# 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 psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration) {

    bool status;
    char errorWord[64];

    PS_ASSERT_PTR_NON_NULL(chip,    false);
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(recipe,  false);
    PS_ASSERT_PTR_NON_NULL(updates, false);

    // select the raw objects for this readout
    psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    if (rawstars == NULL) return false;

    psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    if (refstars == NULL) return false;

    psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
    if (match == NULL) return false;

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

    // allowed limits for valid solutions
    sprintf (errorWord, "PSASTRO.MOSAIC.MAX.ERROR.N%d", iteration);
    REQUIRED_RECIPE_VALUE (float maxError, errorWord, F32, "failed to find single-chip max allowed error\n");
    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MOSAIC.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");

    // set the order of the per-chip fit (higher order only if iteration > 0)
    int order = 1;
    if (iteration > 0) {
	// select the desired chip order
	REQUIRED_RECIPE_VALUE (order, "PSASTRO.MOSAIC.CHIP.ORDER", S32, "failed to find mosaic chip-level fit order\n");

	// modify the order to correspond to the actual number of matched stars:
	if ((match->n < 15) && (order >= 3)) order = 2;
	if ((match->n < 11) && (order >= 2)) order = 1;
	if ((match->n <  8) && (order >= 1)) order = 0;
	if (order < 1) {
	    psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 
	    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);
    fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
    fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");

    // need to pass in an update header, sent in from above
    pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
    if (!results) {
	psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n");
	return false;
    }

    // toSky converts from FPA & TPA units (microns) to sky units (radians)
    pmFPA *fpa = chip->parent;
    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;

    // 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;
    int astNstar = results->yStats->clippedNvalues;

    // if we clip away too many stars, the order may be invalid
    if (order == 3) { minNstar = PS_MAX (15, minNstar); }
    if (order == 2) { minNstar = PS_MAX (11, minNstar); }
    if (order == 1) { minNstar = PS_MAX ( 8, minNstar); }

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

    psFree (fitStats);
    psFree (results);
    return validSolution;
}
