Index: trunk/ppSim/src/Makefile.am =================================================================== --- trunk/ppSim/src/Makefile.am (revision 14667) +++ trunk/ppSim/src/Makefile.am (revision 14668) @@ -14,5 +14,5 @@ ppSimMakeDark.c \ ppSimMakeSky.c \ - ppSimInsertSources.c \ + ppSimInsertStars.c \ ppSimInsertGalaxies.c \ ppSimAddOverscan.c \ Index: trunk/ppSim/src/ppSim.h =================================================================== --- trunk/ppSim/src/ppSim.h (revision 14667) +++ trunk/ppSim/src/ppSim.h (revision 14668) @@ -90,23 +90,12 @@ ); -// Add a star into the signal and variance images -bool ppSimInsertStar(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 - ); - -bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config); -bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); - psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config); 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 ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config); +bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); +bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config); bool ppSimInitHeader(pmConfig *config, Index: trunk/ppSim/src/ppSimInsertSource.c =================================================================== --- trunk/ppSim/src/ppSimInsertSource.c (revision 14667) +++ (revision ) @@ -1,55 +1,0 @@ -# include "ppSim.h" - -bool ppSimInsertSource (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) { - - bool mdok; - - assert (stars); - - 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 - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate - 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."); - return false; - } - } - - int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); - int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); - int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY"); - int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY"); - - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y - - // Rough noise estimate, appropriate for entire cell - float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); - - pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSF"); - - for (long i = 0; i < stars->n; i++) { - ppSimStar *star = stars->data[i]; - - // 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); - - pmModel *model = pmModelFromPSFforXY (psf, x, y, star->flux); - - ppSimInsertSource(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr); - } - return true; -} - Index: trunk/ppSim/src/ppSimInsertSources.c =================================================================== --- trunk/ppSim/src/ppSimInsertSources.c (revision 14667) +++ (revision ) @@ -1,113 +1,0 @@ -# include "ppSim.h" - -bool ppSimInsertSources (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) { - - bool mdok; - - assert (readout); - assert (stars); - - if (!stars->n) { return true; } - - // XXX is this needed? - // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources - - pmCell *cell = readout->parent; - pmChip *chip = cell->parent; - - // XXX this is an estimate of the sky noise based on the inputs to the image simulation. - // XXX update this to allow the estimate based on the measured sky background - // XXX this is missing the gain. - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time - float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate - 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."); - return false; - } - } - - // Rough noise estimate, appropriate for entire cell (use for source radius?) - float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); - - int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); - int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); - int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY"); - int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY"); - - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y - - pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); - assert (psf); - - int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning); - int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); - - psArray *sources = psArrayAllocEmpty (stars->n); - - // 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 xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip); - float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip); - - // 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); - - // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) - float radius = model->modelRadius (model->params, roughNoise); - radius = PS_MAX (radius, 1.0); - - // construct a source, with model flux pixels set, based on the model - pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR); - - // XXX set the mag & err values (should this be done in pmSourceFromModel?) - // XXX i should be applying the gain and the correct effective area - psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - source->psfMag = -2.5*log10(star->flux); - source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux; - - // XXX add the sources to a source array - - // insert the source flux in the image - pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); - psArrayAdd (sources, 100,source); - } - - // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT" - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); - - // XXX many leaks in here, i think - return true; -} - Index: trunk/ppSim/src/ppSimInsertStars.c =================================================================== --- trunk/ppSim/src/ppSimInsertStars.c (revision 14667) +++ trunk/ppSim/src/ppSimInsertStars.c (revision 14668) @@ -1,13 +1,24 @@ # include "ppSim.h" -bool ppSimInsertStars (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) { +bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) { bool mdok; + assert (readout); assert (stars); - + + if (!stars->n) { return true; } + + // XXX is this needed? + // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources + + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + + // XXX this is an estimate of the sky noise based on the inputs to the image simulation. + // XXX update this to allow the estimate based on the measured sky background + // XXX this is missing the gain. 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 +34,7 @@ } + // Rough noise estimate, appropriate for entire cell (use for source radius?) + float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); @@ -35,15 +49,65 @@ 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 (&mdok, chip->analysis, "PSPHOT.PSF"); + assert (psf); + int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning); + int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); + + psArray *sources = psArrayAllocEmpty (stars->n); + + // 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); - ppSimInsertStar(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr); + float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip); + float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip); + + // 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); + + // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) + float radius = model->modelRadius (model->params, roughNoise); + radius = PS_MAX (radius, 1.0); + + // construct a source, with model flux pixels set, based on the model + pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR); + + // XXX set the mag & err values (should this be done in pmSourceFromModel?) + // XXX i should be applying the gain and the correct effective area + psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + source->psfMag = -2.5*log10(star->flux); + source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux; + + // XXX add the sources to a source array + + // insert the source flux in the image + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + psArrayAdd (sources, 100,source); } + + // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT" + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); + + // XXX many leaks in here, i think return true; } + Index: trunk/ppSim/src/ppSimLoop.c =================================================================== --- trunk/ppSim/src/ppSimLoop.c (revision 14667) +++ trunk/ppSim/src/ppSimLoop.c (revision 14668) @@ -125,5 +125,5 @@ if (type == PPSIM_TYPE_OBJECT) { - ppSimInsertSources (readout, expCorr, stars, config); + ppSimInsertStars (readout, expCorr, stars, config); } Index: trunk/ppSim/src/ppSimMakeSources.c =================================================================== --- trunk/ppSim/src/ppSimMakeSources.c (revision 14667) +++ (revision ) @@ -1,75 +1,0 @@ -# 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; -}