# include "ppSim.h"

// XXX add bounds to the inputs?
bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config) {

    bool status;

    psImage *signal = readout->image;
    psImage *variance = readout->weight;

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

    float expTime     = psMetadataLookupF32(&status, config->arguments, "EXPTIME"); // Exposure time
    float flatSigma   = psMetadataLookupF32(&status, config->arguments, "FLAT.SIGMA"); // Flat-field coefficient
    float flatRate    = psMetadataLookupF32(&status, config->arguments, "FLAT.RATE"); // Flat-field rate
    float shutterTime = psMetadataLookupF32(&status, config->arguments, "SHUTTER.TIME"); // Shutter time

    float skyRate     = psMetadataLookupF32(&status, config->arguments, "SKY.RATE"); // Sky rate
    if (isnan(skyRate)) {
	float zp      = psMetadataLookupF32(&status, config->arguments, "ZEROPOINT"); assert (status);
	float scale   = psMetadataLookupF32(&status, config->arguments, "SCALE");     assert (status);
	float skyMags = psMetadataLookupF32(&status, config->arguments, "SKY.MAGS");  assert (status);
	skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
    }

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

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

    int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y

    // Size of FPA
    psRegion *bounds = ppSimFPABounds (fpa);

    // Correct chip offsets so that boresight is in the middle of the FPA
    x0Chip -= 0.5 * (bounds->x1 - bounds->x0);
    y0Chip -= 0.5 * (bounds->y1 - bounds->y0);

    for (int y = 0; y < signal->numRows; y++) {

        float yFPA = PPSIM_CELL_TO_FPA(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
            (bounds->y1 - bounds->y0); // Relative y position in FPA

        for (int x = 0; x < signal->numCols; x++) {
            float xFPA = PPSIM_CELL_TO_FPA(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 /
                (bounds->x1 - bounds->x0); // Relative x position in FPA

            // Shutter: adjust exposure time
            float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;

            // Gaussian flat-field over the FPA
            float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) /
                flatSigma / sqrtf(2.0 * M_PI);

            if (type == PPSIM_TYPE_FLAT) {
                float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field
                signal->data.F32[y][x] += flatFlux;
                variance->data.F32[y][x] += flatFlux;
                continue;
            }

            expCorr->data.F32[y][x] = realExpTime / expTime;

            // Sky background
            float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky
            signal->data.F32[y][x] += skyFlux;
            variance->data.F32[y][x] += skyFlux;

            // TO DO: Add fringes

        }
    }
    psFree(bounds);
    return true;
}

