# include "ppSim.h"

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

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

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

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

    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, 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;
}

