# include "ppSim.h"

psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) {

    pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename);
    PS_ASSERT_PTR_NON_NULL (readout, NULL);

    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    return sources;
}

bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe) {

    bool status;

    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    if (!status) return NULL;

    for (int i = 0; i < sources->n; i++) {
	pmSource *source = sources->data[i];

        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    }
    return true;
}

bool ppSimPhotomReadout(pmConfig *config, const pmFPAview *view) {

    psTimerStart ("psphotReadout");

    // select the current recipe
    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    if (!recipe) {
        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
        return false;
    }

# if 0    
    // set the photcode for this image
    if (!psphotAddPhotcode (recipe, config, view, "PPSIM.CHIP")) {
        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
        return false;
    }
# endif

    // find the currently selected readout. 
    // we always perform photometry on the mosaiced chip
    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.CHIP");
    PS_ASSERT_PTR_NON_NULL (readout, false);

    // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile
    psArray *realSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES");
    PS_ASSERT_PTR_NON_NULL (realSources, false);

    // the newly injected fake sources are saved on the PPSIM.OUTPUT pmFPAfile
    psArray *injectedSources = ppSimSelectSources (config, view, "PPSIM.OUTPUT");
    PS_ASSERT_PTR_NON_NULL (injectedSources, false);

    // XXX need to define the source pixels
    ppSimDefinePixels (realSources, readout, recipe);
    ppSimDefinePixels (injectedSources, readout, recipe);

    // make a copy of the fake-source parameters so these can be saved independently from the
    // measured fake-source parameters
    psArray *fakeSources = psArrayAlloc (injectedSources->n);
    for (int i = 0; i < injectedSources->n; i++) {
      fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]);
    }

    // load the forced source lists
    psArray *forceSources = ppSimLoadForceSources (config, view);
    psAssert (forceSources, "failed to load force photometry sources");

    // Generate the mask and weight images, including the user-defined analysis region of interest
    psphotSetMaskAndWeight (config, readout, recipe);

    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    // this function uses PSPHOT.PSF.LOAD as the pmFPAfile 
    pmPSF *psf = psphotLoadPSF (config, view, recipe);
    assert (psf);

    // remove all sources 
    psphotRemoveAllSources (realSources, recipe);

    // generate a background model (median, smoothed image)
    if (!psphotModelBackground (config, view, "PPSIM.CHIP")) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    }
    if (!psphotSubtractBackground (config, view, "PPSIM.CHIP")) {
        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    }

    // add noise for real (subtracted) objects
    psphotAddNoise (readout, realSources, recipe);

    // XXX fake sources should measure peak->x,y, force sources should not
    psMaskType maskVal = 0xff; 
    psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal);
    ppSimDetections (significance, recipe, fakeSources);
    ppSimDetections (significance, recipe, forceSources);
    psFree (significance);

    // remove noise for subtracted objects (ie, return to normal noise level)
    psphotSubNoise (readout, realSources, recipe);

    // replace all sources
    psphotReplaceAllSources (realSources, recipe);

    // construct an initial model for each object
    psphotGuessModels (readout, realSources, recipe, psf);
    psphotGuessModels (readout, fakeSources, recipe, psf);
    psphotGuessModels (readout, forceSources, recipe, psf);
    
    psArray *sources = NULL;

    // linear fit to real + fake sources
    sources = ppSimMergeSources (realSources, fakeSources);
    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract
    psphotReplaceAllSources (sources, recipe);
    psFree (sources); // only frees the merged references

    // linear fit to real + forced sources (if any were selected)
    if (forceSources->n) {
      sources = ppSimMergeSources (realSources, forceSources);
      psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract
      psphotReplaceAllSources (sources, recipe);
      psFree (sources); // only frees the merged references
    }

    // XXX do we need to measure aperture photometry corrections?
    // XXX do we store the pre and post correction magnitudes?
    // XXX for the fake sources, these must be identically zero
    // XXX for the force sources, need the apresid to put the mags on the correct system
    if (0) {
	if (!psphotApResid (readout, forceSources, recipe, psf)) {
	    psLogMsg ("psphot", 3, "failed on psphotApResid");
	    return psphotReadoutCleanup (config, readout, recipe, NULL, psf, NULL);
	}
    }

    // calculate source magnitudes (for which set??)
    pmReadout *background = psphotSelectBackground (config, view, false);
    psphotMagnitudes(fakeSources, recipe, psf, background);
    psphotMagnitudes(forceSources, recipe, psf, background);

    // drop the references to the image pixels held by each source
    psphotSourceFreePixels (realSources);
    psphotSourceFreePixels (fakeSources);
    psphotSourceFreePixels (forceSources);

    // create the exported-metadata and free local data
    // XXX this places the sources on readout->analysis as PSPHOT.SOURCES.  modify?
    // (or don't supply the sources, and do this with a different function)
    psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL);

    // add forceSources to the readout (real and fake already there)
    // these are outputs: we need to generate the fpa structure
    
    pmCell    *forceCell    = pmFPAfileThisCell (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceCell, "no cell?");
    pmChip    *forceChip    = pmFPAfileThisChip (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceChip, "no chip?");
    pmReadout *forceReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.SOURCES");
    if (!forceReadout) {
	forceReadout = pmReadoutAlloc (forceCell);
	psFree (forceReadout); // there is a copy on 'cell' as well
    }
    psAssert (forceReadout, "no forceReadout?");
    pmChipSetDataStatus (forceChip, true);
    psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "forced photometry ", forceSources);
    psFree (forceSources);

    pmCell    *fakeCell    = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?");
    pmChip    *fakeChip    = pmFPAfileThisChip (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeChip, "no chip?");
    pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES");
    if (!fakeReadout) {
	fakeReadout = pmReadoutAlloc (fakeCell);
	psFree (fakeReadout); // there is a copy on 'cell' as well
    }
    psAssert (fakeReadout, "no fakeReadout?");
    pmChipSetDataStatus (fakeChip, true);
    psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources);

    return true;
}
