/** @file psastroMosaicOneChip.c
 *
 *  @brief 
 *
 *  @ingroup libpsastro
 *
 *  @author IfA
 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-09 21:25:34 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

# 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];
    char orderWord[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
    snprintf (errorWord, 64, "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)
    REQUIRED_RECIPE_VALUE (int defaultOrder, "PSASTRO.MOSAIC.CHIP.ORDER", S32, "failed to find mosaic chip-level fit default order\n");

    snprintf (orderWord, 64, "PSASTRO.MOSAIC.CHIP.ORDER.N%d", iteration);
    int order = psMetadataLookupS32 (&status, recipe, orderWord);
    if (!status || (order == -1)) {
        order = defaultOrder;
    }

    // modify the order to correspond to the actual number of matched stars:
    int Ndof_min = 3;
    int order_max = 0.5*(sqrt(4*match->n - 4*Ndof_min + 1) - 3);
    order = PS_MIN (order, order_max);

    // if ((match->n < 17) && (order >= 3)) order = 2;
    // if ((match->n < 13) && (order >= 2)) order = 1;
    // if ((match->n <  9) && (order >= 1)) order = 0;

    if (order < 0) {
        psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order);
        return false;
    }

    psLogMsg ("psastro", PS_LOG_DETAIL, "mosaic fit chip order %d", order);

    // create output toFPA; set masks appropriate to the Elixir DVO astrometry format if we are
    // fitting 0th order, use the current polynomial of whatever order, with higher order
    // coefficients frozen to the current values
    if (order == 0) {
        // set FIT mask for all higher order terms of the existing solution
        // any existing SET masks will be retained.
        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
                if (i + j > 0) {
                    chip->toFPA->x->coeffMask[i][j] |= PS_POLY_MASK_FIT;
                    chip->toFPA->y->coeffMask[i][j] |= PS_POLY_MASK_FIT;
                }
            }
        }
    } else {
        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
    // only clip if we are fitting the chip parameters.
    psStats *fitStats = NULL;
//    if (order == 0) {
//      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");
//    } else {
        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 ((maxError > 0) && (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);

    //plot results
    pmAstromVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);

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