# 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)
    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    if (sources == NULL)
        return false;

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

    // sort in ascending magnitude order
    psArraySort (inStars, psastroSortByMag);

    // we are going to select the brighter Nmax subset for astrometry 
    int nMax = psMetadataLookupS32 (&status, recipe, "PSASTRO.MAX.NSTAR");
    if (!status || (nMax < 10)) nMax = 300; // 10 is really somewhat absurd as a lower limit

    // choose the first nMax sources
    psArray *rawStars = psArrayAlloc (PS_MIN (nMax, inStars->n));
    for (int i = 0; i < rawStars->n; i++) {
	rawStars->data[i] = psMemIncrRefCounter (inStars->data[i]);
    }

    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\n", sources->n, rawStars->n, inStars->n);

    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;

        pmAstromObj *obj = pmAstromObjAlloc ();

        // is the source magnitude calibrated in any sense?
        obj->pix->x = PAR[2];
        obj->pix->y = PAR[3];
        obj->Mag = source->psfMag;

        // 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)
{
    pmAstromObj *A = *(pmAstromObj **) a;
    pmAstromObj *B = *(pmAstromObj **) 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);
}

