/** @file psastroConvert.c
 *
 *  @brief
 *
 *  @ingroup libpsastro
 *
 *  @author IfA
 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-07 02:03:34 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

# include "psastroInternal.h"
// leak free 2006.04.27

bool psastroConvertFPA (pmFPA *fpa, psMetadata *recipe) {

    pmChip *chip;
    pmCell *cell;
    pmReadout *readout;
    pmFPAview *view = pmFPAviewAlloc (0);

    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
        if (!chip->process || !chip->file_exists) { continue; }

        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) { continue; }

            // process each of the readouts
            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
                if (! readout->data_exists) { continue; }

                psastroConvertReadout (readout, recipe);
            }
        }
    }
    psFree (view);
    return true;
}

// pass/apply the WCS information?
bool psastroConvertReadout (pmReadout *readout, psMetadata *recipe) {

    bool status;

    // PSPHOT.SOURCES carries the pmSource objects (from psphot analysis or loaded externally)
    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    psAssert (detections, "missing detections?");

    psArray *sources = detections->allSources;
    psAssert (sources, "missing sources?");

    // convert the pmSource objects into pmAstromObj objects (drop !STAR and SATSTAR?)
    psArray *inStars = pmSourceToAstromObj (sources);

    // sort in ascending magnitude order
    // psArraySort (inStars, psastroSortByMag);
    // psVector *index = psArraySortIndex (sources, pmSourceSortBySN);
    psVector *index = psArraySortIndex (NULL, inStars, psastroSortByMag);

    // XXX need to exit gracefully is inStars->n is 0 (or 1?)

    // we are going to select the brighter Nmax subset for astrometry
    int nMax = psMetadataLookupS32 (&status, recipe, "PSASTRO.MAX.NRAW");
    if (!status || !nMax) nMax = inStars->n;

    // we are going to select the brighter Nmax subset for astrometry
    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.MAX.INST.MAG.RAW");
    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.MIN.INST.MAG.RAW");

    // we are going to select the brighter Nmax subset for astrometry
    pmSourceMode skip = PM_SOURCE_MODE_DEFAULT;
    char *ignoreList = psMetadataLookupStr (&status, recipe, "PSASTRO.IGNORE");
    if (ignoreList != NULL) {
      psArray *list = psStringSplitArray (ignoreList, ",", false);
      for (int i = 0; i < list->n; i++) {
        pmSourceMode mode = pmSourceModeFromString (list->data[i]);
        if (mode == PM_SOURCE_MODE_DEFAULT) {
          psWarning ("unknown source mode in PSASTRO.IGNORE, skipping");
          continue;
        }
        skip |= mode;
      }
      psFree (list);
    }

    // choose the first nMax sources
    int j = 0;
    psArray *rawStars = psArrayAlloc (PS_MIN (nMax, inStars->n));

    float mMin = +100.0;
    float mMax = -100.0;
    int nModeSkip = 0;
    int nFaintSkip = 0;
    int nBrightSkip = 0;
    int nInfSkip = 0;

    for (int i = 0; (i < inStars->n) && (j < rawStars->n); i++) {
        int n = index->data.S32[i];
        pmSource *source = sources->data[n];

        psTrace ("psastro", 6, "mag: %f +/- %f, mode: %x, skip: %x\n", source->psfMag, source->errMag, source->mode, skip);

        if (source->mode & skip) {
          nModeSkip ++;
          continue;
        }

        if ((iMagMax != 0.0) && (source->psfMag > iMagMax)) {
          nFaintSkip ++;
          continue;
        }
        if ((iMagMin != 0.0) && (source->psfMag < iMagMin)) {
          nBrightSkip ++;
          continue;
        }
        if (!isfinite(source->psfMag)) {
          nInfSkip ++;
          continue;
        }
        mMin = PS_MIN (mMin, source->psfMag);
        mMax = PS_MAX (mMax, source->psfMag);
        rawStars->data[j] = psMemIncrRefCounter (inStars->data[n]);
        j++;
    }
    rawStars->n = j;

    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawStars);

    psLogMsg ("psastro", 4, "loaded %ld sources, using %ld of %ld good sources (inst mag: %f to %f)\n", sources->n, rawStars->n, inStars->n, mMin, mMax);
    psLogMsg ("psastro", 4, "skip reasons: mode: %d, faint: %d, bright: %d, inf: %d\n", nModeSkip, nFaintSkip, nBrightSkip, nInfSkip);

    psFree (index);
    psFree (inStars);
    psFree (rawStars);

    return true;
}

// select a magnitude range?
psArray *pmSourceToAstromObj (psArray *sources) {

    psArray *objects = psArrayAllocEmpty (sources->n);

    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];

        // only accept the PSF sources?
        // XXX drop SATSTAR?
        // if (source->type != PM_SOURCE_TYPE_STAR) continue;

        pmModel *model = source->modelPSF;
        if (model == NULL) continue;

        psF32 *PAR = model->params->data.F32;
        psF32 *dPAR = model->dparams->data.F32;

        pmAstromObj *obj = pmAstromObjAlloc ();

        // is the source magnitude calibrated in any sense?
        obj->pix->x    = PAR[PM_PAR_XPOS];
        obj->pix->y    = PAR[PM_PAR_YPOS];
        obj->pix->xErr = dPAR[PM_PAR_XPOS];
        obj->pix->yErr = dPAR[PM_PAR_YPOS];
        obj->Mag = source->psfMag;
        obj->dMag = source->errMag;

        // XXX do we have the information giving the readout and cell offset?
        // for the moment, assume chip == cell == readout
        *obj->cell = *obj->pix;
        *obj->chip = *obj->cell;

        psArrayAdd (objects, 100, obj);
        psFree (obj);
    }
    return objects;
}

/**
 * sort by Mag (ascending)
 */
int psastroSortByMag (const void *a, const void *b)
{
    const pmAstromObj *A = a;
    const pmAstromObj *B = b;

    psF32 mA = (isfinite(A->Mag)) ? A->Mag : FLT_MAX;
    psF32 mB = (isfinite(B->Mag)) ? B->Mag : FLT_MAX;

    psF32 diff = mA - mB;
    if (diff > +FLT_EPSILON) return (+1);
    if (diff < -FLT_EPSILON) return (-1);
    return (0);
}

