# include "ppSim.h"

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


bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {

    bool mdok;

    assert (readout);
    assert (stars);

    if (!stars->n) { return true; }

    pmCell *cell = readout->parent;
    pmChip *chip = cell->parent;

    // XXX this is an estimate of the sky noise based on the inputs to the image simulation.
    // XXX update this to allow the estimate based on the measured sky background
    // XXX this is missing the gain.
    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe

    float expTime   = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
    float darkRate  = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate

    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
    if (isnan(readnoise)) {
        psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
        readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
        if (!mdok) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
            return false;
        }
    }

    float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate
    if (isnan(skyRate)) {
	float zp      = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok);
	float scale   = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE");     assert (mdok);
	float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    }

    // Rough noise estimate, appropriate for entire cell (use for source radius?)
    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);

    int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0");
    int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0");
    int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY");
    int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY");

    int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");

    int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y

    pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF");
    assert (psf);

    int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning);
    int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning);

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

    // add sources to the readout image & weight
    psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n);
    for (long i = 0; i < stars->n; i++) {
        ppSimStar *star = stars->data[i];
        psTrace("ppSim", 10, "Inserting star at %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux);

        // star->x,y are in fpa coordinates

        // Position on the cell and peak flux
        float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
        float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);

        // Position on the cell and peak flux
        float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
        float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);

        // XXX Note, the below does not put the edges of stars on the readout if they fall slightly off
        // This will be visible as cut-off stars at amplifier boundaries (e.g., Megacam)
        if (xCell < 0) continue;
        if (yCell < 0) continue;
        if (xCell > readout->image->numCols) continue;
        if (yCell > readout->image->numRows) continue;
        // XXX need to apply col0, row0 if readout is a subarray

        // Apply the expCorr to the star->flux before setting the model flux
        float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell];

        // instantiate a model for the PSF at this location, set desired flux
        pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
        pmModelSetFlux (model, flux);

        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
        float radius = model->modelRadius (model->params, roughNoise);
        radius = PS_MAX (radius, 1.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;

        source->psfMag = -2.5*log10(flux);
        source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;

        // insert the source flux in the image
        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY);

        // 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);

        // add the sources to the source array
        psArrayAdd (sources, 100,source);
        psFree(source);                 // Drop reference
    }

    // NOTE: the pmFPAfile "PPSIM.OUTPUT" points at these sources
    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "fake sources", sources);
    psFree(sources);

    // XXX many leaks in here, i think
    return true;
}
