/** @file psastroLoadGhosts.c
 *
 *  @brief calculate ghost FPA and Chip positions for the stars loaded on the FPA based on the model
 *
 *  @ingroup psastro
 *
 *  @author IfA
 *  @version $Revision: 1.7 $
 *  @date $Date: 2009-02-07 02:03:34 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

# include "psastroInternal.h"

# define GET_2D_POLY(NAME,OUT) \
    md = psMetadataLookupMetadata (&status, ghostModel, NAME); \
    if (!md) { \
	psError(PSASTRO_ERR_CONFIG, true, "Missing %s in model file %s", NAME, ghostFile); \
	goto escape; \
    } \
    OUT = psPolynomial2DfromMetadata(md); \
    if (!OUT) { \
	psError(PSASTRO_ERR_CONFIG, true, "Trouble interpretting %s in model file %s", NAME, ghostFile); \
	goto escape; \
    }

# define GET_1D_POLY(NAME,OUT) \
    md = psMetadataLookupMetadata (&status, ghostModel, NAME); \
    if (!md) { \
	psError(PSASTRO_ERR_CONFIG, true, "Missing %s in model file %s", NAME, ghostFile); \
	goto escape; \
    } \
    OUT = psPolynomial1DfromMetadata(md);	\
    if (!OUT) { \
	psError(PSASTRO_ERR_CONFIG, true, "Trouble interpretting %s in model file %s", NAME, ghostFile); \
	goto escape; \
    }

/**
 * calculate ghost FPA and Chip positions for the stars loaded on the FPA
 */
bool psastroLoadGhosts (pmConfig *config) {

    bool status;
    pmChip *chip = NULL;
    pmCell *cell = NULL;
    pmReadout *readout = NULL;
    float zeropt, exptime;
    psMetadata *md = NULL;
    psPolynomial2D *centerX = NULL;
    psPolynomial2D *centerY = NULL;
    psPolynomial1D *outerMajor = NULL;
    psPolynomial1D *outerMinor = NULL;
    psPolynomial1D *innerMajor = NULL;
    psPolynomial1D *innerMinor = NULL;
    psMetadata *ghostModel = NULL;

    psLogMsg ("psastro", PS_LOG_INFO, "determine ghost positions");

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

    bool REFSTAR_MASK_GHOST = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GHOST");
    if (!REFSTAR_MASK_GHOST) return true;

    char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL");
    if (!strcasecmp(ghostFile, "NONE")) return true;

    double MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG");

    if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) {
	psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model");
        return false;
    }

    pmFPAview *view = pmFPAviewAlloc (0);

    GET_2D_POLY ("GHOST.CENTER.X", centerX);
    GET_2D_POLY ("GHOST.CENTER.Y", centerY);

    GET_1D_POLY ("GHOST.OUTER.MAJOR", outerMajor);
    GET_1D_POLY ("GHOST.OUTER.MINOR", outerMinor);
    GET_1D_POLY ("GHOST.INNER.MAJOR", innerMajor);
    GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor);

    // select the input astrometry data (also carries the refstars)
    pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    if (!astrom) {
        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
	goto escape;
    }
    pmFPA *fpa = astrom->fpa;

    // really error-out here?  or just skip?
    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
        psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
	goto escape;
    }

    // recipe values are given in instrumental magnitudes
    // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime)
    float MagOffset = zeropt + 2.5*log10(exptime);
    MAX_MAG += MagOffset;

    // 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; }
        if (!chip->fromFPA) { 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; }

                // select the raw objects for this readout (loaded in psastroExtract.c)
                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
                if (refstars == NULL) { continue; }

                // identify the bright stars of interest
                for (int i = 0; i < refstars->n; i++) {
                    pmAstromObj *ref = refstars->data[i];
                    if (ref->Mag > MAX_MAG) continue;

		    // Ghost model:
		    // X(ghost) = -X + dX(X,Y) --> GHOST.CENTER.X
		    // Y(ghost) = -Y + dY(X,Y) --> GHOST.CENTER.Y
		    // R1inner  = R1inner_0 + r dR1inner
		    // R2inner  = R2inner_0 + r dR1inner
		    // R1outer  = R1outer_0 + r dR1inner
		    // R2outer  = R2outer_0 + r dR1inner

		    psastroGhost *ghost = psastroGhostAlloc ();
		    ghost->srcFP->x = ref->FP->x; 
		    ghost->srcFP->y = ref->FP->y;

		    // XXX it is stupid that this takes -X_fpa,-Y_fpa --> the analysis script is reporting the guess FPA position, and the fit is using that...
		    ghost->FP->x = -ref->FP->x + psPolynomial2DEval(centerX, -ref->FP->x, -ref->FP->y);
		    ghost->FP->y = -ref->FP->y + psPolynomial2DEval(centerY, -ref->FP->x, -ref->FP->y);

		    float rSrc = hypot (ref->FP->x, ref->FP->y);

		    ghost->inner.major = psPolynomial1DEval (innerMajor, rSrc);
		    ghost->inner.minor = psPolynomial1DEval (innerMinor, rSrc);
		    ghost->inner.theta = atan2(ref->FP->y, ref->FP->x);

		    ghost->outer.major = psPolynomial1DEval (outerMajor, rSrc);
		    ghost->outer.minor = psPolynomial1DEval (outerMinor, rSrc);
		    ghost->outer.theta = atan2(ref->FP->y, ref->FP->x);

		    // report instrumental ghost star mags
		    ghost->Mag = ref->Mag - MagOffset;
		    
		    // XXX this code yields a single chip: we need to provide results for any chips
		    // which encompass the full size of the ghost
		    pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, -ghost->srcFP->x, -ghost->srcFP->y);
		    fprintf (stderr, "raw chip position: %f, %f ", ghost->chip->x, ghost->chip->y);

		    ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y);
		    fprintf (stderr, "-> model chip position: %f, %f\n", ghost->chip->x, ghost->chip->y);

		    if (!ghostChip) goto skip;
		    if (!ghostChip->cells) goto skip;
		    if (!ghostChip->cells->n) goto skip;
		    pmCell *ghostCell = ghostChip->cells->data[0];
		    if (!ghostCell) goto skip;
		    if (!ghostCell->readouts) goto skip;
		    if (!ghostCell->readouts->n) goto skip;
		    pmReadout *ghostReadout = ghostCell->readouts->data[0];
		    if (!ghostReadout) goto skip;

		    psArray *ghosts = psMetadataLookupPtr (&status, ghostReadout->analysis, "PSASTRO.GHOSTS");
		    if (ghosts == NULL) { 
			ghosts = psArrayAllocEmpty (100);
			if (!psMetadataAdd (ghostReadout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts)) {
			  psError(PSASTRO_ERR_CONFIG, false, "failure to add ghosts to readout");
			  goto escape;
			}
			psFree (ghosts);
		    }

		    psArrayAdd (ghosts, 100, ghost);

		skip:
		    psFree (ghost);
                }
            }
        }
    }

    psastroExtractFreeChipBounds();

    psFree (centerX);
    psFree (centerY);
    psFree (innerMajor);
    psFree (innerMinor);
    psFree (outerMajor);
    psFree (outerMinor);
    psFree (ghostModel);
    psFree (view);
    return true;

escape:
    psFree (centerX);
    psFree (centerY);
    psFree (innerMajor);
    psFree (innerMinor);
    psFree (outerMajor);
    psFree (outerMinor);
    psFree (ghostModel);
    psFree (view);
    return false;
}
