Index: trunk/ppSim/src/ppSimLoop.c =================================================================== --- trunk/ppSim/src/ppSimLoop.c (revision 14367) +++ trunk/ppSim/src/ppSimLoop.c (revision 14463) @@ -1,268 +1,3 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - #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) @@ -275,195 +10,28 @@ assert(fpa); - bool mdok; // Status of MD lookups + ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate - 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; - } + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator - 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 + // Load catalogue stars + psArray *stars = NULL; 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); + stars = ppSimLoadStars (fpa, config); } // 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 - } - } + // XXX put this in a wrapper (add to array of stars) + if (type == PPSIM_TYPE_OBJECT) { + ppSimMakeStars (stars, fpa, config, rng); } 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); - + ppSimUpdateConceptsFPA (fpa, config); 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 @@ -473,100 +41,21 @@ 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; - } + psVector *biasCols = ppSimMakeBiassec (cell, config); - // 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 - + int numReadouts = 1; for (int i = 0; i < numReadouts; i++) { pmReadout *readout = pmReadoutAlloc(cell); // Readout within cell + // TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately 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 @@ -575,95 +64,27 @@ } - 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 + psVector *biasRows = ppSimMakeBias (signal, variance, cell, config, rng); + if (type == PPSIM_TYPE_BIAS) goto done; + + ppSimMakeDark (signal, variance, config); + if (type == PPSIM_TYPE_DARK) goto done; - 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 + ppSimMakeSky (signal, variance, expCorr, type, config, fpa, chip, cell); + if (type == PPSIM_TYPE_FLAT) goto done; - // Bias level - signal->data.F32[y][x] = rowBias->data.F32[y]; - variance->data.F32[y][x] = PS_SQR(readnoise); + if (type == PPSIM_TYPE_OBJECT) { + ppSimInsertStars (signal, variance, expCorr, stars, config, chip, cell); + } - if (type == PPSIM_TYPE_BIAS) { - continue; - } + done: + ppSimAddNoise(signal, variance, config, cell, rng); + ppSimSaturate(signal, config, cell); - // 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); + readout->image = signal; 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); + ppSimAddOverscan (readout, config, biasCols, biasRows, rng); + psFree(biasRows); readout->data_exists = true; @@ -674,23 +95,9 @@ 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); + ppSimUpdateConceptsCell (cell, config); 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); + ppSimInitHeader(config, NULL, NULL, cell); } - cell->data_exists = true; if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { @@ -698,5 +105,4 @@ psFree(rng); psFree(view); - psFree(bounds); return PS_EXIT_SYS_ERROR; } @@ -704,10 +110,6 @@ 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); + ppSimInitHeader(config, NULL, chip, NULL); } - chip->data_exists = true; if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { @@ -715,5 +117,4 @@ psFree(rng); psFree(view); - psFree(bounds); return PS_EXIT_SYS_ERROR; } @@ -722,6 +123,5 @@ 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); + ppSimInitHeader(config, fpa, NULL, NULL); } @@ -730,14 +130,9 @@ 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;