#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

# include "psphot.h"

// try PSF models and select best option
pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) {

    bool            status;
    char           *modelName;
    pmPSF          *psf = NULL;
    pmPSFtry       *try = NULL;
    psArray        *stars = NULL;

    psTimerStart ("psphot");

    // check if a PSF model is supplied by the user
    psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    if (psf != NULL) return psf;

    // examine PSF sources in S/N order (brightest first)
    sources = psArraySort (sources, psphotSortBySN);

    // array to store candidate PSF stars
    int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS");
    if (!status) NSTARS = PS_MIN (sources->n, 200);

    // use poissonian errors or local-sky errors
    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
    if (!status) POISSON_ERRORS = true;
    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);

    stars = psArrayAlloc (sources->n);

    // select the candidate PSF stars (pointers to original sources)
    for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
        pmSource *source = sources->data[i];
        if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source);
    }
    psLogMsg ("psphot.pspsf", 4, "selected candidate %ld PSF objects\n", stars->n);

    if (stars->n == 0) {
	psError(PSPHOT_ERR_PSF, true, "Failed to find any PSF candidates");
	return NULL;
    }

    // get the fixed PSF fit radius
    // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
    float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");

    // get the list pointers for the PSF_MODEL entries
    psList *list = NULL;
    psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
    if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");

    if (mdi->type == PS_DATA_STRING) {
        list = psListAlloc(NULL);
        psListAdd (list, PS_LIST_HEAD, mdi);
    } else {
        if (mdi->type != PS_DATA_METADATA_MULTI) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
        list = psMemIncrRefCounter(mdi->data.list);
    }

    // set up an array to store the results
    psArray *models = psArrayAlloc (list->n);
    models->n = list->n;

    // try each model option listed in config
    psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE);
    for (int i = 0; i < models->n; i++) {
        psMetadataItem *item = psListGetAndIncrement (iter);
        modelName = item->data.V;
        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS);
    }
    psFree (iter);
    psFree (list);
    psFree (stars);

    // select the best of the models
    // here we are using the clippedStdev on the metric as the indicator
    try = models->data[0];
    int bestN = 0;
    float bestM = try->psf->dApResid;
    for (int i = 1; i < models->n; i++) {
        try = models->data[i];
        float M = try->psf->dApResid;
        if (M < bestM) {
            bestM = M;
            bestN = i;
        }
    }

    // use the best model:
    try = models->data[bestN];

    // unset the PSFSTAR flag for stars not used for PSF model
    for (int i = 0; i < try->sources->n; i++) {
        pmSource *source = try->sources->data[i];
        if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) {
            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
        }
    }

    // save only the best model;
    psf = psMemIncrRefCounter(try->psf);
    psFree (models);

    psMetadataAdd (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "number of stars used to make PSF", psf->nPSFstars);

    modelName = pmModelGetType (psf->type);
    psLogMsg ("psphot.pspsf", 3, "select psf model: %f sec\n", psTimerMark ("psphot"));
    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);

    return (psf);
}

bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {

    psEllipseShape shape;
    psEllipseAxes axes;

    psImage *image = readout->image;

    pmModel *modelEXT = pmModelAlloc (psf->type);

    // make a model with unit central intensity at the image center
    modelEXT->params->data.F32[0] = 0;
    modelEXT->params->data.F32[1] = 1;
    modelEXT->params->data.F32[2] = 0.5*image->numCols;
    modelEXT->params->data.F32[3] = 0.5*image->numRows;

    // construct a PSF model at this coordinate
    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);

    // get the correct model-radius function
    pmModelRadius modelRadiusFunc = pmModelRadius_GetFunction (psf->type);

    // get the model full-width at half-max
    psF64 FWHM_X = 2*modelRadiusFunc (modelPSF->params, 0.5);

    shape.sx  = modelPSF->params->data.F32[4];
    shape.sy  = modelPSF->params->data.F32[5];
    shape.sxy = modelPSF->params->data.F32[6];
    axes = psEllipseShapeToAxes (shape);

    psF64 FWHM_Y = FWHM_X * (axes.minor / axes.major);

    psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
    psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
    psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE",  PS_DATA_F32 | PS_META_REPLACE, "PSF angle",           axes.theta);

    psFree (modelEXT);
    psFree (modelPSF);

    return true;
}

bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) {

    // without the PSF model, we can only rely on the source->moments
    // as a measure of the FWHM values

    int FWHM_N = 0;
    double FWHM_X = 0.0;
    double FWHM_Y = 0.0;
    double FWHM_T = 0.0;

    psEllipseMoments moments;
    psEllipseAxes axes;

    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];
        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;

        // moments->Sx,Sy are roughly sigma_x,y
        moments.x2 = PS_SQR(source->moments->Sx);
        moments.y2 = PS_SQR(source->moments->Sy);
        moments.xy = source->moments->Sxy;

        axes = psEllipseMomentsToAxes (moments);

        FWHM_X += axes.major * 2.35;
        FWHM_Y += axes.minor * 2.35;
        FWHM_T += axes.theta;
        FWHM_N ++;
    }

    FWHM_X /= FWHM_N;
    FWHM_Y /= FWHM_N;
    FWHM_T /= FWHM_N;

    psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
    psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
    psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE",  PS_DATA_F32 | PS_META_REPLACE, "PSF angle",           FWHM_T);

    return true;
}
