Index: trunk/ppSim/src/ppSimLoop.c =================================================================== --- trunk/ppSim/src/ppSimLoop.c (revision 12834) +++ trunk/ppSim/src/ppSimLoop.c (revision 12873) @@ -9,11 +9,102 @@ #include "ppSim.h" -#define COMPARE(VALUE,MIN,MAX) \ - if (VALUE < MIN) { \ - MIN = VALUE; \ - } \ - if (VALUE > MAX) { \ - MAX = VALUE; \ - } +#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 chip0, // Chip offset + int chipParity // Chip parity + ) +{ + return (chip0 + 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 chip0, // Chip offset + int chipParity // Chip parity + ) +{ + return ((pos - chip0) * chipParity - cell0) * cellParity; +} + +// 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, x0Chip, xParityChip); + int cellRight = cell2fpa(xSize, x0Cell, xParityCell, x0Chip, xParityChip); + int cellBottom = cell2fpa(0, y0Cell, yParityCell, y0Chip, yParityChip); + int cellTop = cell2fpa(ySize, y0Cell, yParityCell, 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 @@ -33,52 +124,110 @@ 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"); - - 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, fpa, 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 = x0Chip + xParityChip * x0Cell; - int cellRight = x0Chip + xParityChip * (x0Cell + xParityCell * xSize); - int cellBottom = y0Chip + yParityChip * y0Cell; - int cellTop = y0Chip + yParityChip * (y0Cell + yParityCell * ySize); - - COMPARE(cellLeft, xMin, xMax); - COMPARE(cellRight, xMin, xMax); - COMPARE(cellBottom, yMin, yMax); - COMPARE(cellTop, yMin, yMax); - } + 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); } 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; +} + @@ -92,4 +241,6 @@ 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)) { @@ -98,6 +249,42 @@ } + psVector *ra = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.RA"); // Star RAs + psVector *dec = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.DEC"); // Star Declinations + psVector *mag = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.MAG"); // Star magnitudes + float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) + float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // Plate scale (arcsec/pixel) + scale *= M_PI / 3600.0 / 180.0; // Convert to radians/pixel + float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point + + if (ra && dec && mag) { + // Convert star ra,dec to pixels from FPA centre + 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) + + // Conversion loop + for (long i = 0; i < ra->n; i++) { + float div = sin(dec->data.F32[i]) * sin(dec0) + + cos(dec->data.F32[i]) * cos(dec0) * cos(ra->data.F32[i] - ra0); // Divisor + + // Convert to x,y position on tangent plane, in pixels + float xi = cos(dec->data.F32[i]) * sin(ra->data.F32[i] - ra0) / div / scale; + float eta = (sin(dec->data.F32[i]) * cos(dec0) - + cos(dec->data.F32[i]) * sin(dec0) * cos(ra->data.F32[i] - 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; + ra->data.F32[i] = xi; + dec->data.F32[i] = eta; + + // Convert magnitude to peak flux + mag->data.F32[i] = powf(10.0, -0.4 * (mag->data.F32[i] - zp)) / sqrt(2.0*M_PI) / seeing; + } + } + 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 @@ -105,8 +292,13 @@ 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 float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name 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 @@ -139,4 +331,43 @@ 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); + +#if 0 + // Add random stars + if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) { + // Peak fluxes: faintest and brightest levels for random stars + float faint = roughNoise; + 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 { + // Normalisation, set by the specified stellar density at the specified bright + // magnitude + float norm = starsDensity * numCols * numRows * 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; + + printf("Adding %ld stars to %f\n", num, faint); + + for (long j = 0; j < num; j++) { + float flux = expf((logf(j + 1) - logf(norm)) / starsLum); // Flux of star + + // TO DO: put stars on chips, rather than cells, so they can overlap cell + // boundaries + + // Position on cell + float x = psRandomUniform(rng) * numCols; + float y = psRandomUniform(rng) * numRows; + star(signal, variance, x, y, flux, roughNoise, seeing, expCorr); + } + } + } +#endif + pmCell *cell; // Cell from chip while ((cell = pmFPAviewNextCell(view, fpa, 1))) { @@ -152,11 +383,73 @@ 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 = psMetadataLookupF32(NULL, config->arguments, "GAIN"); - } - float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.RDNOISE"); // CCD read noise, e + 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.RDNOISE is not set; reverting to recipe value."); - readnoise = psMetadataLookupF32(NULL, config->arguments, "RDNOISE"); + 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; } @@ -172,16 +465,30 @@ 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 - float rowBias = biasLevel += psRandomGaussian(rng) * biasRange; - float yFPA = (y0Chip + yParityChip * (y0Cell + yParityCell * y) - 2 * bounds->y0) / + rowBias->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) / + (float)numRows - 0.5) + biasLevel; + float yFPA = cell2fpa(y, y0Cell, yParityCell, y0Chip, yParityChip) * 2.0 / (bounds->y1 - bounds->y0); // Relative y position in FPA for (int x = 0; x < numCols; x++) { - float xFPA = (x0Chip + xParityChip * (x0Cell + xParityCell * x) - 2 * bounds->x0) / + float xFPA = cell2fpa(x, x0Cell, xParityCell, x0Chip, xParityChip) * 2.0 / (bounds->x1 - bounds->x0); // Relative x position in FPA // Bias level - signal->data.F32[y][x] = rowBias; + signal->data.F32[y][x] = rowBias->data.F32[y]; variance->data.F32[y][x] = PS_SQR(readnoise); @@ -211,5 +518,5 @@ continue; } - + expCorr->data.F32[y][x] = flatValue * realExpTime / expTime; // Sky background float skyFlux = skyRate * realExpTime * flatValue; // Flux from sky @@ -219,25 +526,85 @@ // TO DO: Add fringes - // TO DO: Add stars } } - // TO DO: Add overscan - -#if 0 - // Add the noise into the image - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - signal->data.F32[y][x] /= gain; // Converting to ADU - signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x] / gain) * - psRandomGaussian(rng); + // Rough noise estimate, appropriate for entire cell + float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + + + // XXX Not sure I've got the scaling right on the peak flux: it's multiplied by the expTime + // here, and the exposure up above. + + // Add specified (catalogue) stars + if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) { + for (long j = 0; j < ra->n; j++) { + // Position on the cell + float x = fpa2cell(ra->data.F32[j], x0Cell, xParityCell, x0Chip, xParityChip); + float y = fpa2cell(dec->data.F32[j], y0Cell, yParityCell, y0Chip, yParityChip); + star(signal, variance, x, y, mag->data.F32[j] * expTime, + roughNoise, seeing, expCorr); } } -#endif - readout->image = signal; + + // Add random stars + if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) { + // Peak fluxes: faintest and brightest levels for random stars + float faint = roughNoise; + 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 { + // Normalisation, set by the specified stellar density at the specified bright + // magnitude + float norm = starsDensity * numCols * numRows * 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; + + printf("Adding %ld stars to %f\n", num, faint); + + for (long j = 0; j < num; j++) { + float flux = expf((logf(j + 1) - logf(norm)) / starsLum); // Flux of star + + // TO DO: put stars on chips, rather than cells, so they can overlap cell + // boundaries + + // Position on cell + float x = psRandomUniform(rng) * numCols; + float y = psRandomUniform(rng) * numRows; + 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); psFree(readout); // Drop reference } + psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,