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

# include "psphot.h"

// given a pmSource which has been fitted using modelPSF, evaluate the
// resulting fit: did the fit succeed? is this object PSF-like? is this object 
// extended?  return status is TRUE for a valid PSF, FALSE otherwise.

// identify objects consistent with PSF shape/magnitude distribution
// we expect dparams[4],dparams[5] to have a scatter of:
// sigma_x / (S/N) * sqrt(2)
// 1 / (params[4],params[5])*(S/N)

// sigma_x : 1 / SX
// dsx : 1 / (SX * SN)
// dsx_o = hypot (1/(SX*SN), MIN_DSX)

// any objects which is consistent with the PSF should have 
// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag)
// this includes a minimum buffer (DS) for the brighter objects

// saturated stars should fall outside (larger), but have peaks above SATURATION
// extended sources should be larger, cosmic rays smaller
// we also reject objects with S/N too low or ChiSquare to high

// floor for DS value 
// XXX EAM : add to configuration?
# define MIN_DS 0.01

static float SATURATION;
static float PSF_MIN_SN;
static float PSF_MAX_CHI;
static float PSF_SHAPE_NSIGMA;

bool psphotInitLimitsPSF (psMetadata *recipe) {

    bool status;

    // we may set this differently here from the value used to mark likely saturated stars
    SATURATION       = psMetadataLookupF32 (&status, recipe, "SATURATION");
    PSF_MIN_SN       = psMetadataLookupF32 (&status, recipe, "PSF_MIN_SN");
    PSF_MAX_CHI      = psMetadataLookupF32 (&status, recipe, "PSF_MAX_CHI");
    PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_SHAPE_NSIGMA");

    return true;
}

// examine the model->status, fit parameters, etc and decide if the model succeeded
// set the source->type and source->mode appropriately
bool psphotEvalPSF (pmSource *source, pmModel *model) { 

    int keep;
    float dSX, dSY, SX, SY, SN;
    float nSx, nSy, Chi;

    // do we actually have a valid PSF model?
    if (model == NULL) {
	source->mode &= ~PM_SOURCE_MODE_FITTED;
	return false;
    }

    // did the model fit fail for one or another reason?
    switch (model->status) {
      case PM_MODEL_SUCCESS:
	break;
      case PM_MODEL_UNTRIED:
	source->mode &= ~PM_SOURCE_MODE_FITTED; 
	return false;
      case PM_MODEL_BADARGS:
      case PM_MODEL_NONCONVERGE:
      case PM_MODEL_OFFIMAGE:
      default:
	source->mode |= PM_SOURCE_MODE_FAIL;
	return false;
    }

    // unless we prove otherwise, this object is a star.
    source->type = PM_SOURCE_TYPE_STAR;

    // the following source->mode information pertains to modelPSF:
    source->mode |= PM_SOURCE_MODE_PSFMODEL;

    // if the object has fitted peak above saturation, label as SATSTAR
    // this is a valid PSF object, but ignore the other quality tests
    // remember: fit does not use saturated pixels (masked)
    // XXX no extended object can saturate and stay extended...
    if (model->params->data.F32[1] >= SATURATION) {
	if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	    psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
	}
	source->mode |=  PM_SOURCE_MODE_SATSTAR;
	return true;
    } 

    // if the object has a fitted peak below 0, the fit did not converge cleanly
    if (model->params->data.F32[1] <= 0) {
	source->mode |= PM_SOURCE_MODE_FAIL;
	return false;
    } 

    // if the source was predicted to be a SATSTAR, but it fitted below saturation, 
    // make a note to the user
    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
	source->mode &= ~PM_SOURCE_MODE_SATSTAR;
    }

    SN  = model->params->data.F32[1]/model->dparams->data.F32[1];
    SX  = model->params->data.F32[4];
    SY  = model->params->data.F32[5];
    dSX = model->dparams->data.F32[4];
    dSY = model->dparams->data.F32[5];
    Chi = model->chisqNorm / model->nDOF;

    // swing of sigma_x,y in sigmas
    nSx = dSX / hypot (MIN_DS, 1 / (SX * SN));
    nSy = dSY / hypot (MIN_DS, 1 / (SY * SN));

    // assign PM_SOURCE_MODE_GOODSTAR to bright objects within PSF region of dparams[]
    keep = TRUE;
    keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA);
    keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA);
    keep &= (SN > PSF_MIN_SN);
    keep &= (Chi < PSF_MAX_CHI);
    if (keep) return true;

    // this source is not a star, warn if it was a PSFSTAR
    if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality   (%f, %f  :  %f %f %f %f)\n", 
		  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
    } else {
	psLogMsg ("psphot", 5, "fails PSF fit (%f, %f  :  %f %f %f %f)\n", 
		  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
    }

    // object appears to be small, suspected defect
    if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) {
	source->type = PM_SOURCE_TYPE_DEFECT;
	return false;
    }

    // object appears to be large, suspected extended source
    if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) {
	source->type = PM_SOURCE_TYPE_EXTENDED;
	return false;
    }

    // poor-quality fit; only keep if nothing else works...
    source->mode |= PM_SOURCE_MODE_POOR;
    return false;
}	

// examine the model->status, fit parameters, etc and decide if the model succeeded
// set the source->type and source->mode appropriately
bool psphotEvalDBL (pmSource *source, pmModel *model) { 

    // do we actually have a valid PSF model?
    if (model == NULL) {
	source->mode &= ~PM_SOURCE_MODE_FITTED;
	return false;
    }

    // did the model fit fail for one or another reason?
    switch (model->status) {
      case PM_MODEL_SUCCESS:
	break;
      case PM_MODEL_UNTRIED:
	source->mode &= ~PM_SOURCE_MODE_FITTED; 
	return false;
      case PM_MODEL_BADARGS:
      case PM_MODEL_NONCONVERGE:
      case PM_MODEL_OFFIMAGE:
      default:
	source->mode |= PM_SOURCE_MODE_FAIL;
	return false;
    }

    // unless we prove otherwise, this object is a star.
    source->type = PM_SOURCE_TYPE_STAR;

    // the following source->mode information pertains to modelPSF:
    source->mode |= PM_SOURCE_MODE_PSFMODEL;

    // if the object has fitted peak above saturation, label as SATSTAR
    // this is a valid PSF object, but ignore the other quality tests
    // remember: fit does not use saturated pixels (masked)
    // XXX no extended object can saturate and stay extended...
    if (model->params->data.F32[1] >= SATURATION) {
	if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
	    psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
	}
	source->mode |=  PM_SOURCE_MODE_SATSTAR;
	return true;
    } 

    // if the object has a fitted peak below 0, the fit did not converge cleanly
    if (model->params->data.F32[1] <= 0) {
	source->mode |= PM_SOURCE_MODE_FAIL;
	return false;
    } 

    // if the source was predicted to be a SATSTAR, but it fitted below saturation, 
    // make a note to the user
    if (source->mode & PM_SOURCE_MODE_SATSTAR) {
	psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
	source->mode &= ~PM_SOURCE_MODE_SATSTAR;
    }
    return true;
}	
