#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>
#include <psastro.h>

#include "ppSim.h"

#define STAR_RANGE_MIN 4.5              // Minimum pixel range for adding stars, sigma
#define STAR_BG_FRAC 0.1                // Fraction of background error to go out to when adding stars

#define LOG10(X) (log(X)/log(10))       // Base 10 logarithm

// Return FPA position, given a cell position; calculations are all done in pixel units
static inline float cell2fpa(float pos, // Cell position
                             int cell0, // Cell offset
                             int cellParity, // Cell parity
                             int binning, // Binning factor
                             int chip0, // Chip offset
                             int chipParity // Chip parity
    )
{
    return (chip0 + binning * chipParity * (cell0 + cellParity * pos));
}

// Return cell position, given an FPA position; calculations are all done in pixel units
static inline float fpa2cell(float pos, // FPA position
                             int cell0, // Cell offset
                             int cellParity, // Cell parity
                             int binning, // Binning factor
                             int chip0, // Chip offset
                             int chipParity // Chip parity
    )
{
    return ((pos - chip0) * chipParity - cell0) * cellParity / binning;
}

// Compare a value with minimum and maximum values, replacing where required.
#define COMPARE(VALUE,MIN,MAX) { \
        if (VALUE < MIN) { \
            MIN = VALUE; \
        } \
        if (VALUE > MAX) { \
            MAX = VALUE; \
        } \
    }

// Return bounds of a chip, based on the concepts
static psRegion *chipBounds(const pmChip *chip, // Chip for which to determine size
                            pmFPAview *view // View for chip
    )
{
    assert(chip);
    assert(view);

    // Bounds of chip
    int xMin = INT_MAX;
    int xMax = - INT_MAX;
    int yMin = INT_MAX;
    int yMax = - INT_MAX;

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

    if ((xParityChip != 1 && xParityChip != -1) || (yParityChip != 1 && yParityChip != -1)) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Chip parities are not set.");
        psFree(view);
        return NULL;
    }

    pmCell *cell;                   // Cell from chip
    while ((cell = pmFPAviewNextCell(view, chip->parent, 1))) {
        int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");
        int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");

        if (xSize == 0 || ySize == 0) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell sizes are not set.");
            psFree(view);
            return NULL;
        }

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

        if ((xParityCell != 1 && xParityCell != -1) || (yParityCell != 1 && yParityCell != -1)) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell parities are not set.");
            psFree(view);
            return NULL;
        }

        // "Left", "Right", "Bottom" and "Top" don't take into account the parity
        int cellLeft = cell2fpa(0, x0Cell, xParityCell, 1, x0Chip, xParityChip);
        int cellRight = cell2fpa(xSize, x0Cell, xParityCell, 1, x0Chip, xParityChip);
        int cellBottom = cell2fpa(0, y0Cell, yParityCell, 1, y0Chip, yParityChip);
        int cellTop = cell2fpa(ySize, y0Cell, yParityCell, 1, y0Chip, yParityChip);

        COMPARE(cellLeft, xMin, xMax);
        COMPARE(cellRight, xMin, xMax);
        COMPARE(cellBottom, yMin, yMax);
        COMPARE(cellTop, yMin, yMax);
    }

    return psRegionAlloc(xMin, xMax, yMin, yMax);
}

// Return bounds of an FPA, based on the concepts
static psRegion *fpaBounds(const pmFPA *fpa       // FPA for which to determine size
    )
{
    assert(fpa);

    pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA

    // Bounds of focal plane
    int xMin = INT_MAX;
    int xMax = - INT_MAX;
    int yMin = INT_MAX;
    int yMax = - INT_MAX;

    pmChip *chip;                       // Chip from FPA
    while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
        psRegion *bounds = chipBounds(chip, view); // Bounds for chip
        COMPARE(bounds->x0, xMin, xMax);
        COMPARE(bounds->x1, xMin, xMax);
        COMPARE(bounds->y0, yMin, yMax);
        COMPARE(bounds->y1, yMin, yMax);
        psFree(bounds);
    }

    psFree(view);
    return psRegionAlloc(xMin, xMax, yMin, yMax);
}


// Add noise to an image
static psImage *addNoise(psImage *signal, // Signal image, modified and returned
                         const psImage *variance, // Variance image
                         psRandom *rng, // Random number generator
                         float gain     // CCD gain (e/ADU)
    )
{
    assert(signal->type.type == PS_TYPE_F32);
    assert(signal->numCols == variance->numCols && signal->numRows == variance->numRows);

    // Add the noise into the image
    for (int y = 0; y < signal->numRows; y++) {
        for (int x = 0; x < signal->numCols; x++) {
            signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * psRandomGaussian(rng);
            signal->data.F32[y][x] /= gain; // Converting to ADU
        }
    }

    return signal;
}

// Apply saturation limit to image
static void saturate(psImage *image, // Image to apply saturation
                     float saturation // Saturation level
    )
{
    for (int y = 0; y < image->numRows; y++) {
        for (int x = 0; x < image->numCols; x++) {
            if (image->data.F32[y][x] > saturation) {
                image->data.F32[y][x] = saturation;
            }
        }
    }
    return;
}

// Return star flux at a position
static inline float starFlux(float x, float y, // Position of interest
                             float x0, float y0, // Centre of star
                             float seeing // Seeing FWHM (pixels)
    )
{
#ifdef STAR_GAUSSIAN
    // Gaussian star
    return exp(-0.5 * (PS_SQR(x - x0) + PS_SQR(y - y0)) / PS_SQR(seeing));
#else
    // Waussian star
    float z = (PS_SQR(x - x0) + PS_SQR(y - y0)) / (2.0 * PS_SQR(seeing));
    return 1.0 / (1.0 + z + PS_SQR(z));
#endif
}

// Return size of star in pixels
static inline int starSize(float noise, float peak, float seeing)
{
#ifdef STAR_GAUSSIAN
    // Gaussian star (solving Gaussian)
    float target = STAR_BG_FRAC * seeing * sqrt(M_2_PI) * noise / peak;
    return target < 1.0 ? 2.0 * seeing * sqrtf(-logf(target)) + 0.5 : seeing * STAR_RANGE_MIN;
#else
    // Waussian star (solving Waussian where z >> 1 and peak/sqrt(bg) >> 1)
    float target = STAR_BG_FRAC * noise / peak;
    return PS_MAX(sqrtf(1.0 / target) + 0.5, seeing * STAR_RANGE_MIN);
#endif
}

// Add a star into the signal and variance images
static void star(psImage *signal,       // Signal image, to which to add star
                 psImage *variance,     // Variance image, to which to add star
                 float x0, float y0,    // Position of star
                 float peak,            // Peak flux of star
                 float noise,           // Rough noise estimate
                 float seeing,          // Seeing for star
                 const psImage *correction // Exposure correction as a function of position
    )
{
    int size = starSize(noise, peak, seeing); // Size of star

    // Range in image pixels on which to add star
    int xMin = PS_MAX(0, x0 - size);
    int xMax = PS_MIN(signal->numCols - 1, x0 + size);
    int yMin = PS_MAX(0, y0 - size);
    int yMax = PS_MIN(signal->numRows - 1, y0 + size);

    for (int y = yMin; y <= yMax; y++) {
        for (int x = xMin; x <= xMax; x++) {
            float star = peak * correction->data.F32[y][x] * starFlux(x, y, x0, y0, seeing); // Star flux
            signal->data.F32[y][x] += star;
            variance->data.F32[y][x] += star;
        }
    }
    return;
}

// Generate a header containing WCS keywords
static psMetadata *wcsHeader(float ra0, float dec0, // Boresight (radians)
                             float pa,// Position angle (radians)
                             float scale, // Pixel scale (radians/pix)
                             float x0, float y0, // Position of 0,0 on the FPA
                             int xParity, int yParity // Parity in x and y
                             )
{
    psMetadata *header = psMetadataAlloc(); // Header, to return
    pmAstromWCS *wcs = pmAstromWCSAlloc(1, 1); // WCS structure
    wcs->toSky = psProjectionAlloc(ra0, dec0, scale * xParity, scale * yParity, PS_PROJ_TAN);
    wcs->cdelt1 = scale * PM_DEG_RAD * xParity;
    wcs->cdelt2 = scale * PM_DEG_RAD * yParity;
    wcs->crpix1 = x0;
    wcs->crpix2 = y0;
    wcs->trans->x->coeff[1][0] = cos(pa) * wcs->cdelt1;
    wcs->trans->x->coeff[0][1] = -sin(pa) * wcs->cdelt1;
    wcs->trans->y->coeff[1][0] = sin(pa) * wcs->cdelt2;
    wcs->trans->y->coeff[0][1] = cos(pa) * wcs->cdelt2;
    // These aren't used by pmAstromWCStoHeader, but set them anyway
    wcs->trans->x->coeff[0][0] = ra0;
    wcs->trans->y->coeff[0][0] = dec0;
    wcs->trans->x->coeff[1][1] = 0.0;
    wcs->trans->y->coeff[1][1] = 0.0;

    pmAstromWCStoHeader(header, wcs);
    psFree(wcs);

    return header;
}


psExit ppSimLoop(pmConfig *config)
{
    PS_ASSERT_PTR_NON_NULL(config, PS_EXIT_PROG_ERROR);

    pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, OUTPUT_FILE); // Output file
    assert(file);
    pmFPA *fpa = file->fpa;             // FPA for file
    assert(fpa);

    bool mdok;                          // Status of MD lookups

    psRegion *bounds = fpaBounds(fpa);  // Bounds of FPA
    if (!bounds || (bounds->x0 == 0.0 && bounds->x1 == 0.0) || (bounds->y0 == 0.0 && bounds->y1 == 0.0)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to determine bounds of FPA");
        return PS_EXIT_CONFIG_ERROR;
    }

    float biasLevel = psMetadataLookupF32(NULL, config->arguments, "BIAS.LEVEL"); // Bias level
    float biasRange = psMetadataLookupF32(NULL, config->arguments, "BIAS.RANGE"); // Bias range
    int biasOrder = psMetadataLookupS32(NULL, config->arguments, "BIAS.ORDER"); // Bias order
    float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
    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
    float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
    float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude
    float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes
    int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y

    float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
    const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name
    if (!filter) {
        filter = "NONE";
    }
    ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate

    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator

    psVector *ra = NULL;                // Star RAs
    psVector *dec = NULL;               // Star Declinations
    psVector *mag = NULL;               // Star magnitudes
    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") *
        M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
    float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
    float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
    float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
    float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)

    // Add catalogue stars
    if (type == PPSIM_TYPE_OBJECT) {
        // Read catalogue stars using psastro
        psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE);
        if (!astroRecipe) {
            psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE);
            psFree(rng);
            psFree(bounds);
            return PS_EXIT_CONFIG_ERROR;
        }

        // Size of FPA
        float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale;

        psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius);
        psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius);
        psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius);
        psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius);
        psArray *refStars = psastroLoadRefstars(config);
        if (!refStars || refStars->n == 0) {
            psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars.");
            psFree(refStars);
            psFree(rng);
            psFree(bounds);
            return PS_EXIT_CONFIG_ERROR;
        }
        psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);

        ra = psVectorAlloc(refStars->n, PS_TYPE_F32);
        dec = psVectorAlloc(refStars->n, PS_TYPE_F32);
        mag = psVectorAlloc(refStars->n, PS_TYPE_F32);

        // Conversion loop
        for (long i = 0; i < refStars->n; i++) {
            pmAstromObj *ref = refStars->data[i]; // Reference star
            float raStar = ref->sky->r; // RA of star
            float decStar = ref->sky->d; // Dec of star
            float magStar = ref->Mag;       // Magnitude of star


            // Convert to x,y position on tangent plane, in pixels
            float div = sin(decStar) * sin(dec0) + cos(decStar) * cos(dec0) * cos(raStar - ra0); // Divisor
            float xi = cos(decStar) * sin(raStar - ra0) / div / scale;
            float eta = (sin(decStar) * cos(dec0) - cos(decStar) * sin(dec0) * cos(raStar - ra0)) /
                div / scale;

            // Apply rotation
            ra->data.F32[i] = cos(pa) * xi - sin(pa) * eta;
            dec->data.F32[i] = sin(pa) * xi + cos(pa) * eta;

            // Convert magnitude to peak flux
            mag->data.F32[i] = powf(10.0, -0.4 * (magStar - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
        }
        psFree(refStars);

        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0);
        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0);
    }

    // Add random stars
    if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) {

        // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it
        // shouldn't be too bad.
        float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
        if (!mdok) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
            psFree(bounds);
            psFree(rng);
            return PS_EXIT_CONFIG_ERROR;
        }

        // Peak fluxes: faintest and brightest levels for random stars
        float faint = 0.1 * 10.0 * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
        float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
        if (bright < faint) {
            psLogMsg("ppSim", PS_LOG_INFO,
                     "Image noise is above brightest random star --- no random stars added.");
        } else {
            // Size of focal plane
            int xSize = bounds->x1 - bounds->x0;
            int ySize = bounds->y1 - bounds->y0;

            // Normalisation, set by the specified stellar density at the specified bright magnitude
            float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
                powf(bright, starsLum);

            // Total number of stars down to the faint flux end
            long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;

            psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
                     -2.5 * LOG10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);

            // Extend the vectors
            long oldSize;               // Old size of ra, dec, mag vectors
            if (ra && dec && mag) {
                oldSize = ra->n;
                ra = psVectorRealloc(ra, oldSize + num);
                dec = psVectorRealloc(dec, oldSize + num);
                mag = psVectorRealloc(mag, oldSize + num);
                ra->n = dec->n = mag->n = oldSize + num;
            } else {
                oldSize = 0;
                ra = psVectorAlloc(num, PS_TYPE_F32);
                dec = psVectorAlloc(num, PS_TYPE_F32);
                mag = psVectorAlloc(num, PS_TYPE_F32);
            }

            for (long i = 0; i < num; i++) {
                ra->data.F32[oldSize + i] = (psRandomUniform(rng) - 0.5) * xSize ; // x position
                dec->data.F32[oldSize + i] = (psRandomUniform(rng) - 0.5) * ySize; // y position
                mag->data.F32[oldSize + i] = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
            }
        }
    }

    pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA

    // Update FPA concepts
    const char *typeStr;                // Exposure type String
    switch (type) {
      case PPSIM_TYPE_BIAS:   typeStr = "BIAS";   break;
      case PPSIM_TYPE_DARK:   typeStr = "DARK";   break;
      case PPSIM_TYPE_FLAT:   typeStr = "FLAT";   break;
      case PPSIM_TYPE_OBJECT: typeStr = "OBJECT"; break;
      default:
        psAbort("Should never get here.");
    }
    psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBSTYPE", PS_META_REPLACE,
                     "Observation type", typeStr);
    psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBJECT", PS_META_REPLACE,
                     "Observation name", typeStr);
    psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "FPA.EXPTIME", PS_META_REPLACE,
                     "Exposure time (sec)", expTime);
    psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.FILTERID", PS_META_REPLACE, "Filter name", filter);
    psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.FILTER", PS_META_REPLACE, "Filter name", filter);


    pmChip *chip;                       // Chip from FPA
    while ((chip = pmFPAviewNextChip(view, fpa, 1))) {

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

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

        pmCell *cell;                   // Cell from chip
        while ((cell = pmFPAviewNextCell(view, fpa, 1))) {

            // Size, position and orientation of cell
            int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning;
            int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning;
            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");

            float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // CCD gain, e/ADU
            if (isnan(gain)) {
                psWarning("CELL.GAIN is not set; reverting to recipe value GAIN.");
                gain = psMetadataLookupF32(&mdok, recipe, "GAIN");
                if (!mdok) {
                    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find GAIN in recipe.");
                    psFree(bounds);
                    psFree(rng);
                    return PS_EXIT_CONFIG_ERROR;
                }
            }
            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.");
                    psFree(bounds);
                    psFree(rng);
                    return PS_EXIT_CONFIG_ERROR;
                }
            }

            float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
            if (isnan(saturation)) {
                psWarning("CELL.SATURATION is not set; reverting to recipe value SATURATION.");
                readnoise = psMetadataLookupF32(&mdok, recipe, "SATURATION");
                if (!mdok) {
                    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find SATURATION in recipe.");
                    psFree(bounds);
                    psFree(rng);
                    return PS_EXIT_CONFIG_ERROR;
                }
            }

            int readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); // Read direction
            if (readdir != 1) {
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "CELL.READDIR = 1 is the only supported mode.");
                psFree(bounds);
                psFree(rng);
                return PS_EXIT_CONFIG_ERROR;
            }

            psList *biassec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // Bias regions
            int numBias = psListLength(biassec); // Number of bias sections
            psVector *biasCols;
            if (numBias > 0) {
                biasCols = psVectorAlloc(numBias, PS_TYPE_S32);
                psListIterator *iter = psListIteratorAlloc(biassec, PS_LIST_HEAD, false); // Iterator
                psRegion *bias;         // Bias region, from iteration
                int i = 0;              // Counter
                while ((bias = psListGetAndIncrement(iter))) {
                    biasCols->data.S32[i++] = bias->x1 = bias->x0;
                }
            } else {
                biasCols = psVectorAlloc(1, PS_TYPE_S32);
                biasCols->data.S32[0] = psMetadataLookupS32(&mdok, recipe, "OVERSCAN.SIZE");
                if (!mdok) {
                    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find OVERSCAN.SIZE in recipe.");
                    psFree(biasCols);
                    psFree(bounds);
                    psFree(rng);
                    return PS_EXIT_CONFIG_ERROR;
                }
                psRegion *newBias = psRegionAlloc(NAN, NAN, NAN, NAN); // New bias region, to be set later
                biassec = psListAlloc(newBias);
                psFree(newBias);
                psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE,
                              "Bias sections", biassec);
                psFree(biassec);
                numBias = 1;
            }

            // TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately

            // TO DO: Decide if cell is to be video readout
            int numReadouts = 1;        // Number of readouts in cell

            for (int i = 0; i < numReadouts; i++) {
                pmReadout *readout = pmReadoutAlloc(cell); // Readout within cell

                psImage *signal = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels
                psImage *variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels

                // Polynomial for bias
                psPolynomial1D *biasPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, biasOrder);
                for (int j = 0; j < biasOrder + 1; j++) {
                    biasPoly->coeff[j] = biasRange * psRandomGaussian(rng);
                }
                psVector *rowBias = psVectorAlloc(numRows, PS_TYPE_F32); // Bias value, per row
                int biasOffset = 0.5 * numRows * psRandomUniform(rng); // Offset to prevent common pattern

                psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects
                if (type == PPSIM_TYPE_OBJECT) {
                    expCorr = psImageAlloc(numCols, numRows, PS_TYPE_F32);
                }

                for (int y = 0; y < numRows; y++) {
                    // Adjust bias level for this row
                    rowBias->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) /
                                                              (float)numRows - 0.5) + biasLevel;
                    float yFPA = cell2fpa(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
                        (bounds->y1 - bounds->y0); // Relative y position in FPA

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

                        // Bias level
                        signal->data.F32[y][x] = rowBias->data.F32[y];
                        variance->data.F32[y][x] = PS_SQR(readnoise);

                        if (type == PPSIM_TYPE_BIAS) {
                            continue;
                        }

                        // Dark current
                        float darkCurrent = darkRate * expTime; // Dark current accumulated
                        signal->data.F32[y][x] += darkCurrent;
                        variance->data.F32[y][x] += darkCurrent;

                        if (type == PPSIM_TYPE_DARK) {
                            continue;
                        }

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

                        // Gaussian flat-field over the FPA
                        float flatValue = exp(-0.5 / PS_SQR(flatSigma) *
                                              (PS_SQR(yFPA) + PS_SQR(xFPA))) / flatSigma / sqrt(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] = flatValue * realExpTime / expTime;
                        // Sky background
                        float skyFlux = skyRate * realExpTime * flatValue; // Flux from sky
                        signal->data.F32[y][x] += skyFlux;
                        variance->data.F32[y][x] += skyFlux;

                        // TO DO: Add fringes

                    }
                }
                psFree(biasPoly);

                // Add stars
                if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) {
                    // Rough noise estimate, appropriate for entire cell
                    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);

                    for (long j = 0; j < ra->n; j++) {
                        // Position on the cell and peak flux
                        float x = fpa2cell(ra->data.F32[j], x0Cell, xParityCell, binning,
                                           x0Chip, xParityChip);
                        float y = fpa2cell(dec->data.F32[j], y0Cell, yParityCell, binning,
                                           y0Chip, yParityChip);
                        float flux = mag->data.F32[j];
                        star(signal, variance, x, y, flux, roughNoise, seeing, expCorr);
                    }
                }
                readout->image = addNoise(signal, variance, rng, gain);
                saturate(readout->image, saturation);
                psFree(variance);

                psFree(expCorr);

                // Add overscan
                for (int j = 0; j < numBias; j++) {
                    psImage *signal = psImageAlloc(biasCols->data.S32[j], numRows, PS_TYPE_F32); // Overscan
                    for (int y = 0; y < signal->numRows; y++) {
                        for (int x = 0; x < signal->numCols; x++) {
                            signal->data.F32[y][x] = rowBias->data.F32[y];
                        }
                    }
                    psImage *variance = psImageAlloc(biasCols->data.S32[j], numRows, PS_TYPE_F32); // Variance
                    psImageInit(variance, PS_SQR(readnoise));
                    addNoise(signal, variance, rng, gain);
                    psFree(variance);

                    psListAdd(readout->bias, PS_LIST_TAIL, signal);
                    psFree(signal);     // Drop reference
                }

                psFree(rowBias);

                readout->data_exists = true;
                psFree(readout);        // Drop reference
            }
            psFree(biasCols);


            psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,
                             "Exposure time (sec)", expTime);
            psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.DARKTIME", PS_META_REPLACE,
                             "Dark time (sec)", expTime);
            psMetadataAddS32(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE,
                             "Binning in x", binning);
            psMetadataAddS32(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE,
                             "Binning in y", binning);

            if (cell->hdu) {
                cell->hdu->header = wcsHeader(ra0, dec0, pa, scale * binning,
                                              fpa2cell(0.0, x0Cell, xParityCell, binning,
                                                       x0Chip, xParityChip),
                                              fpa2cell(0.0, y0Cell, yParityCell, binning,
                                                       y0Chip, yParityChip),
                                              xParityCell * xParityChip, yParityCell * yParityChip);
            }
            cell->data_exists = true;

            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                psError(PS_ERR_IO, false, "Unable to write file.");
                psFree(rng);
                psFree(view);
                psFree(bounds);
                return PS_EXIT_SYS_ERROR;
            }
        }

        if (chip->hdu) {
            chip->hdu->header = wcsHeader(ra0, dec0, pa, scale * binning,
                                          fpa2cell(0.0, 0, 1, binning, x0Chip, xParityChip),
                                          fpa2cell(0.0, 0, 1, binning, y0Chip, yParityChip),
                                          xParityChip, yParityChip);
        }
        chip->data_exists = true;

        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
            psError(PS_ERR_IO, false, "Unable to write file.");
            psFree(rng);
            psFree(view);
            psFree(bounds);
            return PS_EXIT_SYS_ERROR;
        }

    }

    if (fpa->hdu) {
        fpa->hdu->header = wcsHeader(ra0, dec0, pa, scale * binning, 0.5 * (bounds->x1 - bounds->x0),
                                     0.5 * (bounds->y1 - bounds->y0), 1, 1);
    }

    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
        psError(PS_ERR_IO, false, "Unable to write file.");
        psFree(rng);
        psFree(view);
        psFree(bounds);
        return PS_EXIT_SYS_ERROR;
    }

    psFree(ra);
    psFree(dec);
    psFree(mag);
    psFree(rng);
    psFree(view);
    psFree(bounds);

    return 0;
}
