# include "ppSim.h"

// Reset a pointer: free and set to NULL
#define RESET(PTR) \
    psFree(PTR); \
    PTR = NULL;

psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view) {

    bool status;

    pmFPAfile *outFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT");
    psAssert(outFile, "missing PPSIM.OUTPUT");
    psAssert(outFile->fpa, "missing fpr for PPSIM.OUTPUT");

    pmChip *outChip = pmFPAviewThisChip (view, outFile->fpa);
    pmPSF *psf = psMetadataLookupPtr (&status, outChip->analysis, "PSPHOT.PSF");
    assert (psf);

    psArray *spots = psMetadataLookupPtr(&status, outFile->fpa->analysis, "FORCED.SPOTS");
    PS_ASSERT_PTR_NON_NULL (spots, NULL);

    // in this program, we are looping over the output image, rather than the input as in ppImage
    pmFPAfile *srcFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.REAL.SOURCES");
    psAssert(srcFile, "missing PPSIM.REAL.SOURCES");

    // select the fpa and chip from PPSIM.REAL.SOURCES, which carries the astrometry in the headers
    pmFPA     *fpa     = srcFile->fpa;
    pmChip    *chip    = pmFPAviewThisChip (view, fpa);

    // we use the readout for the PPSIM.FORCE.CHIP, which contains the chip-mosaicked image(s)
    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.CHIP");

    // XXX we should have only one cell and readout per chip, right?

    // read WCS from existing output from psphot & psastro
    // XXX can we be more flexible?  ie, use the header of PPSIM.FORCE.CHIP?
    bool bilevelAstrometry = false;
    status = psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
    psAssert (status, "failed to set astrometry");

    // the pixel scale is only used to set the focal-plane coordinate system
    // setting this value to 1 makes this system have the units of pixels
    float pixelScale = 1.0;
    status = psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry);
    psAssert (status, "failed to set astrometry");

    float minX = 0.0;
    float maxX = readout->image->numCols;
    float minY = 0.0;
    float maxY = readout->image->numRows;

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

    // Select the spots within range of this readout.  Project the spots to this chip
    for (int i = 0; i < spots->n; i++) {
	pmAstromObj *obj = spots->data[i];

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

	// limit the X,Y range of the objs to the selected chip
	if (obj->chip->x < minX) continue;
	if (obj->chip->x > maxX) continue;
	if (obj->chip->y < minY) continue;
	if (obj->chip->y > maxY) continue;

	// convert the pmAstromObj to pmSource

        // instantiate a model for the PSF at this location, Io = 1.0
        pmModel *model = pmModelFromPSFforXY (psf, obj->chip->x, obj->chip->y, 1.0);

        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
	// XXX use a fixed radius??
        // float radius = model->modelRadius (model->params, roughNoise);
        // radius = PS_MAX (radius, 1.0);
	float radius = 5.0;

        // construct a source, with model flux pixels set, based on the model
        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);

        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
        // XXX i should be applying the gain and the correct effective area
        // psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
        // psF64 Area = 2.0 * M_PI * axes.major * axes.minor;

	// these are not really needed since we will be fitting for them
        source->psfMag = NAN;
        source->errMag = NAN;

        // insert the source flux in the image
	// XXX not sure the offset is really 0,0
        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, 0.0, 0.0);
        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, 0.0, 0.0);

        // Blow away the image parts of the source, which makes the memory explode
        RESET(source->pixels);
        RESET(source->weight);
        RESET(source->maskObj);
        RESET(source->maskView);
        RESET(source->modelFlux);
        RESET(source->psfFlux);
        RESET(source->blends);

	psArrayAdd (sources, 100, source);
        psFree(source);                 // Drop local reference
    }
    psTrace ("psastro", 4, "Added %ld sources\n", sources->n);

    return sources;
}
