# include "psastroInternal.h"

bool psastroChooseRefstars (pmConfig *config, psArray *refs) {

    bool status;
    pmChip *chip = NULL;
    pmCell *cell = NULL;
    pmReadout *readout = NULL;

    // select the current recipe
    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    if (!recipe) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n");
        return false;
    }

    // select the input data sources
    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    if (!input) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n");
        return false;
    }

    float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA");
    if (!status) fieldExtra = 0.0;

    bool matchLumFunc = psMetadataLookupBool (&status, recipe, "PSASTRO.MATCH.LUMFUNC");

    pmFPAview *view = pmFPAviewAlloc (0);
    pmFPA *fpa = input->fpa;

    // this loop selects the matched stars for all chips
    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
            // XXX there can only be one readout per chip in astrometry, right?
            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
                if (! readout->data_exists) { continue; }

                // read WCS data from the corresponding header
                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);

                int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1");
                int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2");

                float minX = -fieldExtra*Nx;
                float maxX = (1+fieldExtra)*Nx;
                float minY = -fieldExtra*Ny;
                float maxY = (1+fieldExtra)*Ny;

                // the refstars is a subset within range of this chip
                psArray *refstars = psArrayAllocEmpty (100);

                // select the reference objects within range of this readout
                // project the reference objects to this chip
                for (int i = 0; i < refs->n; i++) {
                    pmAstromObj *ref = pmAstromObjCopy(refs->data[i]);

                    psProject (ref->TP, ref->sky, fpa->toSky);
                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);

                    // limit the X,Y range of the refs to the selected chip
                    if (ref->chip->x < minX) goto skip;
                    if (ref->chip->x > maxX) goto skip;
                    if (ref->chip->y < minY) goto skip;
                    if (ref->chip->y > maxY) goto skip;

                    psArrayAdd (refstars, 100, ref);
                skip:
                    psFree (ref);
                }
                psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);

		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
		psFree (refstars);

		if (matchLumFunc) {
		    // in this case, no PSASTRO.REFSTARS is added to readout->analysis
		    if (!psastroRefstarSubset (readout)) {
			psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n");
			return false;
		    }
		}
            }
        }
    }

    psFree (view);
    return true;
}
