Index: /branches/eam_branch_20070817/ppSim/src/Makefile.am =================================================================== --- /branches/eam_branch_20070817/ppSim/src/Makefile.am (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/Makefile.am (revision 14547) @@ -13,5 +13,5 @@ ppSimMakeDark.c \ ppSimMakeSky.c \ - ppSimInsertStars.c \ + ppSimInsertSources.c \ ppSimAddOverscan.c \ ppSimAddNoise.c \ Index: /branches/eam_branch_20070817/ppSim/src/ppSim.h =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSim.h (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSim.h (revision 14547) @@ -89,8 +89,9 @@ psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config); -psVector *ppSimMakeBias (psImage *signal, psImage *variance, pmCell *cell, pmConfig *config, const psRandom *rng) ; -bool ppSimMakeDark (psImage *signal, psImage *variance, pmConfig *config); -bool ppSimMakeSky (psImage *signal, psImage *variance, psImage *expCorr, ppSimType type, pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell); +psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) ; +bool ppSimMakeDark (pmReadout *readout, pmConfig *config); +bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config); bool ppSimInsertStars (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell); +bool ppSimInsertSources (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config); bool ppSimInitHeader(pmConfig *config, @@ -99,7 +100,6 @@ pmCell *cell); -bool ppSimSaturate(psImage *image, // Image to apply saturation - const pmConfig *config, - const pmCell *cell); // Saturation level +bool ppSimSaturate(pmReadout *readout, // Image to apply saturation + const pmConfig *config); // Saturation level bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config); @@ -108,10 +108,10 @@ bool ppSimAddOverscan (pmReadout *readout, pmConfig *config, psVector *biasCols, psVector *biasRows, psRandom *rng); -psImage *ppSimAddNoise(psImage *signal, // Signal image, modified and returned - const psImage *variance, // Variance image - const pmConfig *config, - const pmCell *cell, - const psRandom *rng // Random number generator - ); +bool ppSimAddNoise(psImage *signal, + psImage *variance, + const pmCell *cell, + const pmConfig *config, + const psRandom *rng // Random number generator + ); #endif Index: /branches/eam_branch_20070817/ppSim/src/ppSimAddNoise.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimAddNoise.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimAddNoise.c (revision 14547) @@ -2,9 +2,9 @@ // Add noise to an image -psImage *ppSimAddNoise(psImage *signal, // Signal image, modified and returned - const psImage *variance, // Variance image - const pmConfig *config, - const pmCell *cell, - const psRandom *rng // Random number generator +bool ppSimAddNoise(psImage *signal, // Signal image, modified and returned + psImage *variance, + const pmCell *cell, + const pmConfig *config, + const psRandom *rng // Random number generator ) { @@ -22,5 +22,5 @@ if (!mdok) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find GAIN in recipe."); - return NULL; + return false; } } @@ -38,5 +38,5 @@ // XXX why return this?? - return signal; + return true; } Index: /branches/eam_branch_20070817/ppSim/src/ppSimAddOverscan.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimAddOverscan.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimAddOverscan.c (revision 14547) @@ -29,5 +29,5 @@ psImageInit(variance, PS_SQR(readnoise)); - ppSimAddNoise(signal, variance, config, readout->parent, rng); + ppSimAddNoise(signal, variance, readout->parent, config, rng); psListAdd(readout->bias, PS_LIST_TAIL, signal); Index: /branches/eam_branch_20070817/ppSim/src/ppSimCreate.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimCreate.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimCreate.c (revision 14547) @@ -12,5 +12,4 @@ pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE); - psFree(fpa); if (!file) { psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s. " @@ -20,4 +19,5 @@ if (file->type != PM_FPA_FILE_IMAGE) { psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE); + psFree(fpa); psFree(file); return NULL; @@ -31,5 +31,7 @@ if (!status) { psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); - return status; + psFree(fpa); + psFree(file); + return NULL; } } @@ -44,4 +46,5 @@ psFree(fpa); psFree(view); + psFree(file); return NULL; } @@ -55,4 +58,5 @@ psFree(fpa); psFree(view); + psFree(file); return NULL; } @@ -66,4 +70,5 @@ psFree(fpa); psFree(view); + psFree(file); return NULL; } @@ -71,4 +76,6 @@ } } + + psFree(fpa); psFree(view); Index: /branches/eam_branch_20070817/ppSim/src/ppSimInsertSource.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimInsertSource.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimInsertSource.c (revision 14547) @@ -1,5 +1,5 @@ # include "ppSim.h" -bool ppSimInsertSources (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) { +bool ppSimInsertSource (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) { bool mdok; Index: /branches/eam_branch_20070817/ppSim/src/ppSimInsertSources.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimInsertSources.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimInsertSources.c (revision 14547) @@ -1,13 +1,17 @@ # include "ppSim.h" -bool ppSimInsertSources (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) { +bool ppSimInsertSources (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) { - bool mdok; assert (stars); + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + + // XXX we probably need this to get the size right +# if (0) + bool mdok; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate @@ -23,4 +27,8 @@ } + // Rough noise estimate, appropriate for entire cell (use for source radius?) + float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); +# endif + int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); @@ -35,25 +43,45 @@ int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y - // Rough noise estimate, appropriate for entire cell - float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); - pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSF"); + int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning); + int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); + + // add sources to the readout image & weight for (long i = 0; i < stars->n; i++) { ppSimStar *star = stars->data[i]; + // star->x,y are in fpa coordinates + // Position on the cell and peak flux - float x = PPSIM_FPA_TO_CELL(star->x, x0Cell, xParityCell, binning, x0Chip, xParityChip); - float y = PPSIM_FPA_TO_CELL(star->y, y0Cell, yParityCell, binning, y0Chip, yParityChip); + float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip); + float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip); - // instantiate a model for the PSF at this location (normalized or not?) - // XXX - pmModel *model = pmModelFromPSFforXY (psf, star->x, star->y, star->flux); + // Position on the cell and peak flux + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning); + + if (xCell < 0) continue; + if (yCell < 0) continue; + if (xCell > readout->image->numCols) continue; + if (yCell > readout->image->numRows) continue; + // XXX need to apply col0, row0 if readout is a subarray + + // XXX apply the expCorr to the star->flux before setting the model flux + + // instantiate a model for the PSF at this location, set desired flux + pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0); + pmModelSetFlux (model, star->flux); // construct a source, with model flux pixels set, based on the model - pmSource *source = pmSourceFromModel (model); + pmSource *source = pmSourceFromModel (model, readout, PM_SOURCE_TYPE_STAR); - ppSimInsertSource(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr); + // XXX add the sources to a source array + + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); } + + // XXX add the source array to the readout analysis + return true; } Index: /branches/eam_branch_20070817/ppSim/src/ppSimLoop.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimLoop.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimLoop.c (revision 14547) @@ -31,5 +31,5 @@ // load any needed files (eg, input image, PSF) if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(PSPHOT_ERR_DATA, false, "failed IO for fpa in ppSim\n"); + psError(PS_ERR_UNKNOWN, false, "failed IO for fpa in ppSim\n"); psFree(view); return false; @@ -43,12 +43,12 @@ // load any needed files (eg, input image, PSF) if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(PSPHOT_ERR_DATA, false, "failed IO for chip %d in ppSim\n", view->chip); + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in ppSim\n", view->chip); psFree (view); return false; } - pmPSF *psf = NULL; + // XXX check for a supplied PSF, otherwise generate a basic one if (type == PPSIM_TYPE_OBJECT) { - psf = ppSimSetPSF (config); + // XXX TBD : ppSimSetPSF (chip, config); } @@ -74,6 +74,6 @@ // 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 + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels + readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects @@ -82,24 +82,24 @@ } - psVector *biasRows = ppSimMakeBias (signal, variance, cell, config, rng); + psVector *biasRows = ppSimMakeBias (readout, config, rng); if (type == PPSIM_TYPE_BIAS) goto done; - ppSimMakeDark (signal, variance, config); + ppSimMakeDark (readout, config); if (type == PPSIM_TYPE_DARK) goto done; - ppSimMakeSky (signal, variance, expCorr, type, config, fpa, chip, cell); + ppSimMakeSky (readout, expCorr, type, config); if (type == PPSIM_TYPE_FLAT) goto done; + // if (type == PPSIM_TYPE_OBJECT) { + // ppSimInsertStars (readout, expCorr, stars, config, chip, cell); + // } + if (type == PPSIM_TYPE_OBJECT) { - ppSimInsertStars (signal, variance, expCorr, stars, config, chip, cell); + ppSimInsertSources (readout, expCorr, stars, config); } done: - ppSimAddNoise(signal, variance, config, cell, rng); - ppSimSaturate(signal, config, cell); - - readout->image = signal; - psFree(variance); - psFree(expCorr); + ppSimAddNoise(readout->image, readout->weight, cell, config, rng); + ppSimSaturate(readout, config); ppSimAddOverscan (readout, config, biasCols, biasRows, rng); Index: /branches/eam_branch_20070817/ppSim/src/ppSimMakeBias.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimMakeBias.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimMakeBias.c (revision 14547) @@ -1,7 +1,9 @@ # include "ppSim.h" -psVector *ppSimMakeBias (psImage *signal, psImage *variance, pmCell *cell, pmConfig *config, const psRandom *rng) { +psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) { bool mdok; + + pmCell *cell = readout->parent; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe @@ -20,4 +22,7 @@ } } + + psImage *signal = readout->image; + psImage *variance = readout->weight; int numRows = signal->numRows; Index: /branches/eam_branch_20070817/ppSim/src/ppSimMakeDark.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimMakeDark.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimMakeDark.c (revision 14547) @@ -2,5 +2,8 @@ // XXX add bounds to the inputs? -bool ppSimMakeDark (psImage *signal, psImage *variance, pmConfig *config) { +bool ppSimMakeDark (pmReadout *readout, pmConfig *config) { + + psImage *signal = readout->image; + psImage *variance = readout->weight; float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate Index: /branches/eam_branch_20070817/ppSim/src/ppSimMakeSky.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimMakeSky.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimMakeSky.c (revision 14547) @@ -2,5 +2,12 @@ // XXX add bounds to the inputs? -bool ppSimMakeSky (psImage *signal, psImage *variance, psImage *expCorr, ppSimType type, pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell) { +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 Index: /branches/eam_branch_20070817/ppSim/src/ppSimMakeSources.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimMakeSources.c (revision 14547) +++ /branches/eam_branch_20070817/ppSim/src/ppSimMakeSources.c (revision 14547) @@ -0,0 +1,75 @@ +# include "ppSim.h" + +bool ppSimMakeSources(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { + + assert (stars); + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + 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 skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate + float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + + float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point + 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) + + if (starsDensity <= 0) return true; + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + // Size of focal plane (can't I get this from the config?) + int xSize = bounds->x1 - bounds->x0; + int ySize = bounds->y1 - bounds->y0; + psFree(bounds); + + // 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); + return false; + } + + // 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 { + + // 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); + + long oldSize = stars->n; + psArrayRealloc (stars, stars->n + num); + + for (long i = 0; i < num; i++) { + ppSimStar *star = ppSimStarAlloc (); + + // make fpa center of distribution + star->x = psRandomUniform(rng) * xSize; // x position + star->y = psRandomUniform(rng) * ySize; // y position + star->peak = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux + star->flux = star->peak * sqrt(2.0*M_PI)*seeing; + stars->data[oldSize + i] = star; + } + stars->n = stars->nalloc; + } + return true; +} Index: /branches/eam_branch_20070817/ppSim/src/ppSimSaturate.c =================================================================== --- /branches/eam_branch_20070817/ppSim/src/ppSimSaturate.c (revision 14546) +++ /branches/eam_branch_20070817/ppSim/src/ppSimSaturate.c (revision 14547) @@ -2,10 +2,12 @@ // Apply saturation limit to image -bool ppSimSaturate(psImage *image, // Image to apply saturation - const pmConfig *config, - const pmCell *cell // Saturation level +bool ppSimSaturate(pmReadout *readout, // Image to apply saturation + const pmConfig *config // configuration data ) { bool mdok; + + psImage *image = readout->image; + pmCell *cell = readout->parent; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe