Index: trunk/ppSim/configure.ac =================================================================== --- trunk/ppSim/configure.ac (revision 17915) +++ trunk/ppSim/configure.ac (revision 18011) @@ -17,7 +17,8 @@ AC_SYS_LARGEFILE -PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0]) -PKG_CHECK_MODULES([PSASTRO], [psastro >= 0.9.0]) +PKG_CHECK_MODULES([PSPHOT], [psphot >= 0.8.0]) +PKG_CHECK_MODULES([PSASTRO], [psastro >= 0.9.0]) IPP_STDOPTS Index: trunk/ppSim/notes.txt =================================================================== --- trunk/ppSim/notes.txt (revision 17915) +++ trunk/ppSim/notes.txt (revision 18011) @@ -1,2 +1,83 @@ + + +2008.05.17 + + pmFPAfiles used in ppSim: + + PPSIM.INPUT : an input image loaded by ppSim; the simulated features + may be added to this image. In fact, the features are + generated in the output image, and the two images are + added (as images) together. + + PPSIM.OUTPUT : the output image. this file is used for the main + ppSimLoop function + + PPSIM.CHIP : chip-mosaiced version of the output image. only + generated if photometry is requested. Is this + actually used as an output image? + + PPSIM.REAL.SOURCES : real sources loaded by PPSIM. these sources + are subtracted from the image before testing + for the detectability of the fake sources and + before measuring the flux of fake or forced + soures. They should be the sources measured in + a previous psphot run, though any input list + could be used. Results are very ill-defined if + the sources do not correspond to actual image + sources! + + PSPHOT.PSF.LOAD : input PSF model used to measure the input sources + + PPSIM.SOURCES : output fake sources (injected fake sources with real values) + PPSIM.FAKE.SOURCES : output fake photometry (measured photometry for all fake sources) + PPSIM.FORCE.SOURCES : output force photometry + + ** used within psphot functions: + PSPHOT.BACKMDL : model background generated by psphotModelBackground + PSPHOT.BACKMDL.STDEV : model background error generated by psphotModelBackground + PSPHOT.BACKGND : full-scale model background generated by psphotSubtractBackground + PSPHOT.BACKSUB : background-subtracted image + +2008.05.15 + + For fake and force photometry, we need to load the known existing + sources, then perform a complete linear fit solution to the complete + set of both real (known) sources plus the fake and/or forced + photometry positions (as PSFs). The fake and forced photometry + positions need to be measured independently. This provides the real + photometry for these sources ** in the presence of the other real + sources **. In order to test the detectability of the fake sources + (or the forced sources, for that matter), we need to generate the + smoothed detection image, and determine the peak value at the + positions of the sources. + + ppSimPhotomReadout outline: + + * we have three source lists (loaded before the function is called): + * realSources : these are loaded from a CMF-style file (or equiv) + * fakeSources : these are defined internally, and need a link from + the input fake source to the measured fake source + * forceSources : these are loaded from a DVO database via getstar + + * we require a supplied PSF + + * need to subtract the background (before or after subtracting the + sources?) + + * subtract the real sources + * model the background + * build the detection image + * measure the detectability of the fake and force sources + * replace the real sources + + * merge real + fake sources + * linear fit to merged source list + * replace the sources (make 'no-subtract' option?) + + * merge real + forced sources + * linear fit to merged source list + * replace the sources (make 'no-subtract' option?) + + ppSim development work still needed: Index: trunk/ppSim/src/Makefile.am =================================================================== --- trunk/ppSim/src/Makefile.am (revision 17915) +++ trunk/ppSim/src/Makefile.am (revision 18011) @@ -1,27 +1,38 @@ bin_PROGRAMS = ppSim ppSimSequence -ppSim_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) $(ppSim_CFLAGS) -ppSim_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS) +ppSim_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PSASTRO_CFLAGS) $(ppSim_CFLAGS) +ppSim_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_LIBS) $(PSASTRO_LIBS) ppSim_SOURCES = \ - ppSim.c \ - ppSimArguments.c \ - ppSimCreate.c \ - ppSimLoadStars.c \ - ppSimMakeStars.c \ - ppSimMakeGalaxies.c \ - ppSimMakeBiassec.c \ - ppSimMakeBias.c \ - ppSimMakeDark.c \ - ppSimMakeSky.c \ - ppSimInsertStars.c \ - ppSimInsertGalaxies.c \ - ppSimAddOverscan.c \ - ppSimAddNoise.c \ - ppSimSaturate.c \ - ppSimBounds.c \ - ppSimStars.c \ - ppSimSetPSF.c \ - ppSimUtils.c \ - ppSimLoop.c \ + ppSim.c \ + ppSimArguments.c \ + ppSimCreate.c \ + ppSimLoadStars.c \ + ppSimMakeStars.c \ + ppSimMakeGalaxies.c \ + ppSimMakeBiassec.c \ + ppSimMakeBias.c \ + ppSimMakeDark.c \ + ppSimMakeSky.c \ + ppSimInsertStars.c \ + ppSimInsertGalaxies.c \ + ppSimAddOverscan.c \ + ppSimAddNoise.c \ + ppSimSaturate.c \ + ppSimBounds.c \ + ppSimStars.c \ + ppSimSetPSF.c \ + ppSimUtils.c \ + ppSimLoop.c \ + ppSimLoadSpots.c \ + ppSimPhotom.c \ + ppSimPhotomReadoutFake.c \ + ppSimPhotomReadoutForce.c \ + ppSimPhotomFiles.c \ + ppSimLoadForceSources.c \ + ppSimMergeReadouts.c \ + ppSimDetections.c \ + ppSimMergeSources.c \ + ppSimMosaicChip.c \ + ppSimRandomGaussian.c \ ppSimBadPixels.c Index: trunk/ppSim/src/ppSim.c =================================================================== --- trunk/ppSim/src/ppSim.c (revision 17915) +++ trunk/ppSim/src/ppSim.c (revision 18011) @@ -1,3 +1,12 @@ # include "ppSim.h" + +int failure (pmConfig *config, psExit status, char *message) { + ppSimRandomGaussianFree (); + psErrorStackPrint(stderr, message); + psFree(config); + pmModelClassCleanup(); + psLibFinalize(); + exit (status); +} int main(int argc, char *argv[]) @@ -8,22 +17,20 @@ pmConfig *config = pmConfigRead(&argc, argv, PPSIM_RECIPE); // Configuration if (!config) { - psErrorStackPrint(stderr, "Unable to read configurations."); - exit(PS_EXIT_CONFIG_ERROR); + failure (config, PS_EXIT_CONFIG_ERROR, "Unable to read configurations."); } - ppSimArguments(argc, argv, config); + if (!ppSimArguments(argc, argv, config)) { + failure (config, PS_EXIT_CONFIG_ERROR, "Error parsing command-line arguments"); + } if (!ppSimCreate(config)) { - psErrorStackPrint(stderr, "Unable to create output file."); - psFree(config); - exit(PS_EXIT_CONFIG_ERROR); + failure (config, PS_EXIT_CONFIG_ERROR, "Unable to create output file."); } if (!ppSimLoop(config)) { - psErrorStackPrint(stderr, "Unable to generate data."); - psFree(config); - exit(PS_EXIT_SYS_ERROR); + failure (config, PS_EXIT_SYS_ERROR, "Unable to generate data."); } + ppSimRandomGaussianFree (); psFree(config); pmModelClassCleanup(); @@ -32,4 +39,4 @@ psLibFinalize(); - return PS_EXIT_SUCCESS; + exit (PS_EXIT_SUCCESS); } Index: trunk/ppSim/src/ppSim.h =================================================================== --- trunk/ppSim/src/ppSim.h (revision 17915) +++ trunk/ppSim/src/ppSim.h (revision 18011) @@ -13,12 +13,13 @@ #include #include +#include #define PPSIM_RECIPE "PPSIM" -#define OUTPUT_FILE "PPSIM.OUTPUT" +// #define OUTPUT_FILE "PPSIM.OUTPUT" // 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; } \ +#define COMPARE(VALUE,MIN,MAX) { \ + if (VALUE < MIN) { MIN = VALUE; } \ + if (VALUE > MAX) { MAX = VALUE; } \ } @@ -70,7 +71,7 @@ /// Parse command-line arguments -void ppSimArguments(int argc, char *argv[], ///< Command-line arguments +bool ppSimArguments(int argc, char **argv, ///< Command-line arguments pmConfig *config ///< Configuration - ); + ); /// Create output file @@ -78,23 +79,25 @@ /// Returns a borrowed pointer to the FPA file. pmFPAfile *ppSimCreate(pmConfig *config ///< Configuration - ); + ); // Return bounds of a chip, based on the concepts psRegion *ppSimChipBounds(const pmChip *chip, // Chip for which to determine size pmFPAview *view // View for chip - ); + ); // Return bounds of an FPA, based on the concepts psRegion *ppSimFPABounds(const pmFPA *fpa // FPA for which to determine size - ); + ); /// Loop over the output file, generating simulated data bool ppSimLoop(pmConfig *config ///< Configuration - ); + ); psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config); -psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) ; +psVector *ppSimMakeBias (bool *status, pmReadout *readout, pmConfig *config, const psRandom *rng) ; bool ppSimMakeDark (pmReadout *readout, pmConfig *config); bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config); + +bool ppSimLoadSpots (pmFPA *fpa, pmConfig *config); bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config); @@ -127,4 +130,10 @@ bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config); +bool ppSimMosaicChip(pmConfig *config, const psMaskType blankMask, const pmFPAview *view, + const char *outFile, const char *inFile); + +bool ppSimPhotom (pmConfig *config, pmFPAview *view); + + /// Add bad pixels to an image bool ppSimBadPixels(pmReadout *readout, ///< Readout for which to generate bad pixels @@ -139,3 +148,54 @@ float ppSimMagToFlux (float mag, float zp); +float ppSimArgToRecipeF32(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ); + +int ppSimArgToRecipeS32(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ); + +char *ppSimArgToRecipeStr(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ); + +bool ppSimArgToRecipeBool(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ); + +ppSimType ppSimTypeFromString (char *typeStr); +char *ppSimTypeToString (ppSimType type); + +float ppSimGetZeroPoint (psMetadata *recipe, char *filter); + +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view); + +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma); +double ppSimRandomGaussianNorm (const psRandom *rnd); +void ppSimRandomGaussianFree(); + +bool ppSimPhotomFiles (pmConfig *config, pmFPAfile *fakeFile, pmFPAfile *forceFile); + +bool ppSimPhotomReadoutFake(pmConfig *config, const pmFPAview *view); +bool ppSimPhotomReadoutForce(pmConfig *config, const pmFPAview *view); + +psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view); +bool ppSimDetections (psImage *significance, psMetadata *recipe, psArray *sources); +psArray *ppSimMergeSources (psArray *in1, psArray *in2); + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename); +bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe); + #endif Index: trunk/ppSim/src/ppSimAddNoise.c =================================================================== --- trunk/ppSim/src/ppSimAddNoise.c (revision 17915) +++ trunk/ppSim/src/ppSimAddNoise.c (revision 18011) @@ -32,5 +32,5 @@ // (generate a static array with the cumulative distribution, use the // random number to select a bin from the histogram) - signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * psRandomGaussian(rng); + signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * ppSimRandomGaussianNorm(rng); signal->data.F32[y][x] /= gain; // Converting to ADU } Index: trunk/ppSim/src/ppSimAddOverscan.c =================================================================== --- trunk/ppSim/src/ppSimAddOverscan.c (revision 17915) +++ trunk/ppSim/src/ppSimAddOverscan.c (revision 18011) @@ -4,4 +4,8 @@ bool mdok; + + // skip this step if we did not generate a bias level + if (biasCols == NULL) return true; + if (biasRows == NULL) return true; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe Index: trunk/ppSim/src/ppSimArguments.c =================================================================== --- trunk/ppSim/src/ppSimArguments.c (revision 17915) +++ trunk/ppSim/src/ppSimArguments.c (revision 18011) @@ -1,3 +1,3 @@ -#include "ppSim.h" +# include "ppSim.h" // Print usage information and die @@ -18,51 +18,40 @@ } -// Get a value from the command-line arguments and add it to recipe options -bool valueArgRecipe(psMetadata *options, // Target to which to add value - const char *recipeName, // Name for value in the recipe - psMetadata *arguments, // Command-line arguments - const char *argName // Argument name in the command-line arguments - ) +// this function supplements the RECIPE:OPTIONS folder with command-line options +bool ppSimArguments(int argc, char **argv, pmConfig *config) { - bool status; // Status of MD lookup - float value = psMetadataLookupF32(&status, arguments, argName); // Value of interest - if (isnan(value)) return true; - status = psMetadataAddF32(options, PS_LIST_TAIL, recipeName, 0, NULL, value); - return status; -} - -// this function supplements the RECIPE:OPTIONS folder with command-line options -void ppSimArguments(int argc, char *argv[], pmConfig *config) -{ - bool mdok; // Status of MD lookup + bool status; assert(config); - // save the following command-line options in the arguments structure + // save the following command-line options in the arguments structure. these will later be + // parsed and moved to the config->recipes:PPSIM_RECIPE folder + psMetadata *arguments = psMetadataAlloc(); // Command-line arguments - psMetadataAddStr(arguments, PS_LIST_TAIL, "-format", 0, "Camera format name", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Exposure type (BIAS|DARK|FLAT|OBJECT)", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-filter", 0, "Filter name", NULL); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-exptime", 0, "Exposure time (s)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-biaslevel", 0, "Bias level (e)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-biasrange", 0, "Bias range (e)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-darkrate", 0, "Dark rate (e/s)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-flatsigma", 0, "Flat sigma", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-flatrate", 0, "Flat rate (e/s)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-shuttertime", 0, "Shutter time (s)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-skyrate", 0, "Sky rate (e/s)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-skymags", 0, "Sky brightness in mags / square arcsec", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-ra", 0, "RA (degrees)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-dec", 0, "Dec (degrees)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-pa", 0, "Position angle (degrees)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-scale", 0, "Plate scale (arcsec/pixel)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Photometric zero point", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-seeing", 0, "Seeing FWHM (arcsec)", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-starslum", 0, "Fake star luminosity function slope", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsmag", 0, "Brightest magnitude for fake stars", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsdensity", 0, "Density of fake stars at magnitude", NAN); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-psfclass", 0, "Type of PSF model", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-galmodel", 0, "Type of Galaxy model", NULL); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-format", 0, "Camera format name", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Exposure type (BIAS|DARK|FLAT|OBJECT)", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-filter", 0, "Filter name", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-exptime", 0, "Exposure time (s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-biaslevel", 0, "Bias level (e)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-biasrange", 0, "Bias range (e)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-darkrate", 0, "Dark rate (e/s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-flatsigma", 0, "Flat sigma", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-flatrate", 0, "Flat rate (e/s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-shuttertime", 0, "Shutter time (s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-skyrate", 0, "Sky rate (e/s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-skymags", 0, "Sky brightness in mags / square arcsec", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-ra", 0, "RA (degrees)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-dec", 0, "Dec (degrees)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-pa", 0, "Position angle (degrees)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-scale", 0, "Plate scale (arcsec/pixel)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Photometric zero point", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-seeing", 0, "Seeing FWHM (arcsec)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-starslum", 0, "Fake star luminosity function slope", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsmag", 0, "Brightest magnitude for fake stars", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsdensity", 0, "Density of fake stars at magnitude", NAN); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-psfclass", 0, "Type of PSF model", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-galmodel", 0, "Type of Galaxy model", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1); + psMetadataAddBool(arguments, PS_LIST_TAIL, "+photom", 0, "Perform photometry on fake sources", false); if (psArgumentGet (argc, argv, "-h")) { usage(argv[0], arguments, config); } @@ -71,19 +60,77 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); - - // only one of -camera and -file is needed - bool status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); - if (!config->camera && !status) { - psErrorStackPrint(stderr, "A camera name (-camera NAME) or an image (-file NAME) must be specified"); - exit(PS_EXIT_CONFIG_ERROR); + if (psErrorCodeLast () != PS_ERR_NONE) { psAbort ("problem with -psf or -psflist option"); } + + pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT.SOURCES", "-cmf", "-cmflist"); + if (psErrorCodeLast () != PS_ERR_NONE) { psAbort ("problem with -cmf or -cmflist option"); } + + // Only one of -camera and -file is needed or allowed. The -camera option would have been + // already parsed by pmConfigRead in ppSim.c and resulted in a value for config->camera + bool loadImage = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-input", "-inputlist"); + if (!config->camera && !loadImage) { + psError(PS_ERR_IO, true, "A camera name (-camera NAME) or an input image (-input NAME) must be specified"); + psFree(arguments); + return false; + } + if (config->camera && loadImage) { + psError(PS_ERR_IO, true, "Only one of (-camera NAME) and (-file NAME) may be specified"); + psFree(arguments); + return false; } if (!psArgumentParse(arguments, &argc, argv)) { - exit(PS_EXIT_CONFIG_ERROR); + psError(PS_ERR_IO, false, "error in command-line arguments"); + psFree(arguments); + return false; } if (argc != 2) { - psErrorStackPrint(stderr, "Missing output filename"); - exit(PS_EXIT_CONFIG_ERROR); + psError(PS_ERR_IO, true, "Missing output filename"); + psFree(arguments); + return false; + } + + // output filename + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + // save the additional recipe values based on command-line options. These options override + // the PPSIM recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PPSIM_RECIPE); + if (!options) { + psError(PS_ERR_IO, false, "Unable to find recipe options for %s", PPSIM_RECIPE); + psFree(arguments); + return false; + } + + // these arguments may be used whether the input image is created or loaded + ppSimArgToRecipeF32(&status, options, "STARS.LUM", arguments, "-starslum"); + ppSimArgToRecipeF32(&status, options, "STARS.MAG", arguments, "-starsmag"); + ppSimArgToRecipeF32(&status, options, "STARS.DENSITY", arguments, "-starsdensity"); + ppSimArgToRecipeBool(&status, options, "PHOTOM", arguments, "+photom"); + + // if we are loading the input image (not creating it), then we can skip the remaining arguments + if (loadImage) { + // if we are supplying an input image, it is so we may supply fake stars; force it to + // be treated as an OBJECT image. + psMetadataAddStr(options, PS_LIST_TAIL, "IMAGE.TYPE", 0, "Exposure type", "OBJECT"); + + // check for these options as well + ppSimArgToRecipeStr(&status, options, "PSF.MODEL", arguments, "-psfclass"); // PSF model class + ppSimArgToRecipeStr(&status, options, "GALAXY.MODEL", arguments, "-galmodel"); // Galaxy model name + + // 'seeing' is not required: we can load a psf-model instead; but if not, it is allowed + ppSimArgToRecipeF32(&status, options, "SEEING", arguments, "-seeing"); // seeing (FWHM in arcsec) + + // 'scale' is not required: we can use the WCS instead; but if not, it is allowed + ppSimArgToRecipeF32(&status, options, "PIXEL.SCALE", arguments, "-scale"); // Plate scale + + // XXX we need to be more flexible: get this from the input image header or WCS + float ra0 = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight + float dec0 = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight + psMetadataAddF32(options, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0); + psMetadataAddF32(options, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0); + + psFree (arguments); + return true; } @@ -96,95 +143,48 @@ psMetadata *formats = psMetadataLookupMetadata(NULL, config->camera, "FORMATS"); // The camera formats if (!formats) { - psErrorStackPrint(stderr, "Unable to find FORMATS in camera configuration."); + psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration."); psFree(arguments); - psFree(config); - exit(PS_EXIT_CONFIG_ERROR); + return false; } psMetadata *format = psMetadataLookupMetadata(NULL, formats, formatName); // Format of interest if (!format) { - psErrorStackPrint(stderr, "Unable to find format %s in camera FORMATS.", formatName); + psError(PS_ERR_IO, false, "Unable to find format %s in camera FORMATS.", formatName); psFree(arguments); - psFree(config); - exit(PS_EXIT_CONFIG_ERROR); + return false; } config->format = psMemIncrRefCounter(format); } - // specify the type of simulated image to produce - // XXX this should not be required if we supplied INPUT - const char *typeStr = psMetadataLookupStr(NULL, arguments, "-type"); // Exposure type - if (!typeStr) { - psErrorStackPrint(stderr, "An exposure type must be specified using -type"); - exit(PS_EXIT_CONFIG_ERROR); - } - ppSimType type = PPSIM_TYPE_NONE; // Type to simulate - if (strcasecmp(typeStr, "BIAS") == 0) { - type = PPSIM_TYPE_BIAS; - } else if (strcasecmp(typeStr, "DARK") == 0) { - type = PPSIM_TYPE_DARK; - } else if (strcasecmp(typeStr, "FLAT") == 0) { - type = PPSIM_TYPE_FLAT; - } else if (strcasecmp(typeStr, "OBJECT") == 0) { - type = PPSIM_TYPE_OBJECT; + char *typeStr = ppSimArgToRecipeStr (&status, options, "IMAGE.TYPE", arguments, "-type"); // Requested exposure type + if (typeStr == NULL) { + psError(PS_ERR_IO, false, "An exposure type must be specified using -type"); + psFree(arguments); + return false; + } + ppSimType type = ppSimTypeFromString (typeStr); + + ppSimArgToRecipeStr(&status, options, "FILTER", arguments, "-filter"); // Filter name + + // set the exposure time + if (type == PPSIM_TYPE_BIAS) { + psMetadataAddF32(options, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", 0.0); } else { - psErrorStackPrint(stderr, "Unrecognised exposure type: %s", typeStr); - exit(PS_EXIT_CONFIG_ERROR); - } - assert(type != PPSIM_TYPE_NONE); - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "TYPE", 0, "Exposure type", type); - - const char *filter = psMetadataLookupStr(NULL, arguments, "-filter"); // Filter name - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter); - - // save the following additional recipe values based on command-line options - // these options override the PPSIM recipe values loaded from recipe files - psMetadata *options = pmConfigRecipeOptions (config, PPSIM_RECIPE); - if (!options) { - psErrorStackPrint(stderr, "Unable to find recipe options for %s", PPSIM_RECIPE); - psFree(arguments); - psFree(config); - exit(PS_EXIT_CONFIG_ERROR); - } - - float expTime; - if (type == PPSIM_TYPE_BIAS) { - expTime = 0.0; - } else { - expTime = psMetadataLookupF32(NULL, arguments, "-exptime"); // Exposure time - if (isnan(expTime)) { - psErrorStackPrint(stderr, "The exposure time must be specified using -exptime"); - exit(PS_EXIT_CONFIG_ERROR); - } - } - psMetadataAddF32(options, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", expTime); + ppSimArgToRecipeF32(&status, options, "EXPTIME", arguments, "-exptime"); + if (!status) { + psError(PS_ERR_IO, false, "The exposure time must be specified using -exptime"); + psFree(arguments); + return false; + } + } // these values all get moved from arguments to RECIPE:OPTIONS - valueArgRecipe(options, "BIAS.LEVEL", arguments, "-biaslevel"); - valueArgRecipe(options, "BIAS.RANGE", arguments, "-biasrange"); - valueArgRecipe(options, "DARK.RATE", arguments, "-darkrate"); - valueArgRecipe(options, "FLAT.SIGMA", arguments, "-flatsigma"); - valueArgRecipe(options, "FLAT.RATE", arguments, "-flatrate"); - valueArgRecipe(options, "SHUTTER.TIME", arguments, "-shuttertime"); - valueArgRecipe(options, "SKY.RATE", arguments, "-skyrate"); - valueArgRecipe(options, "STARS.LUM", arguments, "-starslum"); - valueArgRecipe(options, "STARS.MAG", arguments, "-starsmag"); - valueArgRecipe(options, "STARS.DENSITY", arguments, "-starsdensity"); - - psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe - - int biasOrder = psMetadataLookupS32(&mdok, recipe, "BIAS.ORDER"); // Overscan polynomial order - if (!mdok) { - psWarning("BIAS.ORDER(S32) is not set in the recipe %s --- assuming %d", PPSIM_RECIPE, biasOrder); - } - psMetadataAddS32(options, PS_LIST_TAIL, "BIAS.ORDER", 0, "Overscan polynomial order", biasOrder); - - int binning = psMetadataLookupS32(NULL, arguments, "-bin"); // Binning in x and y - if (binning <= 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Binning (%d) is non-positive.", binning); - exit(PS_EXIT_CONFIG_ERROR); - } - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BINNING", 0, "Binning in x and y", binning); - - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + ppSimArgToRecipeF32(&status, options, "BIAS.LEVEL", arguments, "-biaslevel"); + ppSimArgToRecipeF32(&status, options, "BIAS.RANGE", arguments, "-biasrange"); + ppSimArgToRecipeF32(&status, options, "DARK.RATE", arguments, "-darkrate"); + ppSimArgToRecipeF32(&status, options, "FLAT.SIGMA", arguments, "-flatsigma"); + ppSimArgToRecipeF32(&status, options, "FLAT.RATE", arguments, "-flatrate"); + ppSimArgToRecipeF32(&status, options, "SHUTTER.TIME", arguments, "-shuttertime"); + ppSimArgToRecipeF32(&status, options, "SKY.RATE", arguments, "-skyrate"); + ppSimArgToRecipeS32(&status, options, "BINNING", arguments, "-bin"); if (type == PPSIM_TYPE_OBJECT) { @@ -197,63 +197,25 @@ // XXX scale and zp should be supplied by the config file (allow override, but this is camera-dependent) if (isnan(ra0) || isnan(dec0) || isnan(pa) || isnan(seeing)) { - psErrorStackPrint(stderr, "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type"); - exit(PS_EXIT_CONFIG_ERROR); - } - - float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point - if (isnan(zp)) { - // use the filter to get the zeropoint from the recipe - psMetadataItem *zpItem = psMetadataLookup (recipe, "ZEROPTS"); - // check that item is multi... - - psArray *entries = psListToArray (zpItem->data.list); - - // search for matching filter - for (int i = 0; i < entries->n; i++) { - psMetadataItem *item = entries->data[i]; - psMetadata *entry = item->data.V; - - char *filterName = psMetadataLookupStr (&status, entry, "FILTER"); - assert (filterName); - - if (strcmp(filterName, filter)) continue; - - zp = psMetadataLookupF32 (&status, entry, "ZERO_PT"); - assert (status); - break; - } - psFree (entries); - } - - float scale = psMetadataLookupF32(NULL, arguments, "-scale"); // Plate scale - if (isnan(scale)) { - scale = psMetadataLookupF32 (&status, recipe, "PIXEL.SCALE"); - } - - float skymags = psMetadataLookupF32(NULL, arguments, "-skymags"); // Position angle - if (isnan(skymags)) { - skymags = psMetadataLookupF32 (&status, recipe, "SKY.MAGS"); - } - + psError(PS_ERR_IO, false, "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type"); + psFree(arguments); + return false; + } + + ppSimArgToRecipeF32(&status, options, "PIXEL.SCALE", arguments, "-scale"); // Plate scale (arcsec / pixel) + ppSimArgToRecipeF32(&status, options, "SKY.MAGS", arguments, "-skymags"); // Plate scale + ppSimArgToRecipeStr(&status, options, "PSF.MODEL", arguments, "-psfclass"); // PSF model class + ppSimArgToRecipeStr(&status, options, "GALAXY.MODEL", arguments, "-galmodel"); // Galaxy model name + + // the user supplies FWHM in arcsec; here we convert to Sigma in pixels + psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing FWHM (arcsec)", seeing); psMetadataAddF32(options, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0); psMetadataAddF32(options, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0); psMetadataAddF32(options, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",pa * M_PI / 180.0); - // the user supplies FWHM in arcsec; here we convert to Sigma in pixels - psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); - - psMetadataAddF32(options, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale); - psMetadataAddF32(options, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp); - psMetadataAddF32(options, PS_LIST_TAIL, "SKY.MAGS", 0, "sky surface brightness", skymags); - - const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // PSF model class - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "PSF.MODEL", 0, "PSF model class", psfClass); - - const char *galModel = psMetadataLookupStr(NULL, arguments, "-galmodel"); // Galaxy model name - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "GALAXY.MODEL", 0, "Galaxy model", galModel); + ppSimArgToRecipeF32(&status, options, "ZEROPOINT", arguments, "-zp"); // Zero point } psFree(arguments); - return; + return true; } @@ -262,12 +224,6 @@ PSPHOT.PSF INPUT - TYPE - FILTER - BIAS.ORDER - BINNING OUTPUT - PSF.MODEL - GALAXY.MODEL - - all othr values should come from the recipe + + all other values should come from the recipe */ Index: trunk/ppSim/src/ppSimCreate.c =================================================================== --- trunk/ppSim/src/ppSimCreate.c (revision 17915) +++ trunk/ppSim/src/ppSimCreate.c (revision 18011) @@ -1,127 +1,226 @@ # include "ppSim.h" - -// XXX this function forces us to define the camera (on the command line) and format (via the -// PPSIM.OUTPUT entry). In this case, we need to set config->format,formatName based on these -// values. This will be a problem when we want to load an input image (in order to add fake -// stars). We will need to add some logic in ppSimArguments to distinguish the cases of 1) -// input image and 2) specified camera pmFPAfile *ppSimCreate(pmConfig *config) { + PS_ASSERT_PTR_NON_NULL(config, NULL); + bool status; - bool simImage = false; - PS_ASSERT_PTR_NON_NULL(config, NULL); pmFPA *fpa = NULL; // the input image defines the camera. if it is not supplied, the user must have // supplied a camera and other metadata on the command line - pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT"); + pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPSIM.INPUT", "INPUT"); if (!input) { - // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes - if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { - psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName); - return NULL; - } - - simImage = true; - fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the observation - if (!fpa) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); - return NULL; - } - + // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes + if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { + psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName); + return NULL; + } } else { - simImage = false; - if (input->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_IO, true, "PPIMAGE.INPUT is not of type IMAGE"); - return NULL; - } - fpa = input->fpa; - } - - // define the output image file - pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE); - if (!file) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s. " - "Did you forget to specify the format?", OUTPUT_FILE); - return NULL; - } - 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; - } - file->save = true; - - config->format = psMemIncrRefCounter (file->format); - config->formatName = psStringCopy (file->formatName); - - // have we supplied a psf model? - if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { - bool status = false; - - // tie the psf file to the chipMosaic - pmFPAfileBindFromArgs(&status, file, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); - psFree(fpa); - psFree(file); - return NULL; - } - } - - // XXX only invoke this code for OBJECT types of images + // If an image is supplied, we still generate a fake image and merge them together downstream + // (otherwise, we get the variance wrong). + if (input->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.INPUT is not of type IMAGE"); + return NULL; + } + } + + // generate the fpa structure used by the output camera (determined from INPUT or specified) + assert (config->camera); + fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the observation + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + return NULL; + } + + // define the output image file -- this is the basis for the ppSimLoop + pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSIM.OUTPUT"); + if (!output) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from PPSIM.OUTPUT. Did you forget to specify the format?"); + return NULL; + } + if (output->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "PPSIM.OUTPUT type is not IMAGE"); + psFree(fpa); + return NULL; + } + // XXX we should not require the output image to be written + output->save = true; + + config->format = psMemIncrRefCounter (output->format); + config->formatName = psStringCopy (output->formatName); + + // the recipe is now fully realized for the desired camera + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + + if (type == PPSIM_TYPE_OBJECT) { + // adjust the seeing by the scale + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); + if (isnan(seeing)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "seeing is not defined"); + psFree(fpa); + return NULL; + } + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); + if (isnan(scale)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "pixel scale is not defined"); + psFree(fpa); + return NULL; + } + psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + + // if we have been supplied an input image, but no ra & dec, use the input image values + if (input) { + float ra = psMetadataLookupF32(&status, recipe, "RA"); + if (isnan(ra)) { + ra = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.RA"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "RA", PS_META_REPLACE, "ra (radians)", ra); + } + float dec = psMetadataLookupF32(&status, recipe, "DEC"); + if (isnan(dec)) { + dec = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.DEC"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "DEC", PS_META_REPLACE, "dec (radians)", dec); + } + } + } + + if ((type == PPSIM_TYPE_OBJECT) || (type == PPSIM_TYPE_FLAT)) { + // determine the zeropoint from the filter + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); + if (isnan(zp)) { + char *filter = psMetadataLookupStr(&status, recipe, "FILTER"); + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + } + + // For photometry, we operate on the chip-mosaicked image. we create a copy of the mosaicked + // image for psphot so we can write out a clean image + bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes + if (doPhotom) { + + // XXX at the moment, we can perform photometry on the fake sources and the forced + // photometry positions, but not one or the other. Also, we only support photometry in + // the context of an image previously analysed by psphot. Add support for: + // * photometry of a pure fake image (requires peak detection and psf creation) + // * photometry of forced source positions without fake image (this might work, it just + // requires not generating any fake features). + + if (!input) { + psError(PS_ERR_UNKNOWN, false, "input image not found; currently required for photometry"); + return NULL; + } + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *fakeImage = pmFPAfileDefineChipMosaic(config, output->fpa, "PPSIM.FAKE.CHIP"); + if (!fakeImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FAKE.CHIP")); + psFree(fpa); + return NULL; + } + if (fakeImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.FAKE.CHIP is not of type IMAGE"); + psFree(fpa); + return NULL; + } + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *forceImage = NULL; + if (input) { + forceImage = pmFPAfileDefineChipMosaic(config, input->fpa, "PPSIM.FORCE.CHIP"); + if (!forceImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FORCE.CHIP")); + psFree(fpa); + return NULL; + } + if (forceImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.FORCE.CHIP is not of type IMAGE"); + psFree(fpa); + return NULL; + } + } + + // define associated psphot input/output files + if (!ppSimPhotomFiles (config, fakeImage, forceImage)) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files for psphot"); + psFree(fpa); + return NULL; + } + } else { + // have we supplied a psf model? this happens in ppSimPhotomFiles if we request a photometry + // analysis. however, even if we do not, a psf model may be used to generate the fake + // sources. + if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { + // tie the psf file to the chipMosaic + pmFPAfileBindFromArgs(&status, output, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); + psFree(fpa); + return NULL; + } + } + } + // PPSIM.SOURCES carries the constructed, fake sources with their true parameters - pmFPAfile *simSources = pmFPAfileDefineOutput (config, file->fpa, "PPSIM.SOURCES"); + // XXX only invoke this code for OBJECT types of images? + pmFPAfile *simSources = pmFPAfileDefineOutput (config, output->fpa, "PPSIM.SOURCES"); if (!simSources) { - psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES"); - return false; + psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES"); + return false; } simSources->save = true; - // if we have loaded an input image, we do not need to populate the fpa - if (!simImage) { - return file; - } - - pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes + // if we have loaded an input image, we derive certain values from the image, if possible + if (input) { + // we need to extract certain metadata from the image and populate the recipe. + // or else we need to set the fpa concepts based on the recipe options... + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + ppSimArgToRecipeF32(&status, recipe, "EXPTIME", input->fpa->concepts, "FPA.EXPOSURE"); + char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", input->fpa->concepts, "FPA.FILTER"); + + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + + pmFPALevel phuLevel = pmFPAPHULevel(output->format); // Level at which PHU goes pmFPAview *view = pmFPAviewAlloc(0);// View for current level if (phuLevel == PM_FPA_LEVEL_FPA) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - psFree(file); - return NULL; - } + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } } pmChip *chip; // Chip from FPA while ((chip = pmFPAviewNextChip(view, fpa, 1))) { - if (phuLevel == PM_FPA_LEVEL_CHIP) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - psFree(file); - return NULL; - } - } - - pmCell *cell; // Cell from chip - while ((cell = pmFPAviewNextCell(view, fpa, 1))) { - if (phuLevel == PM_FPA_LEVEL_CELL) { - if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); - psFree(fpa); - psFree(view); - psFree(file); - return NULL; - } - } - } + if (phuLevel == PM_FPA_LEVEL_CHIP) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + + pmCell *cell; // Cell from chip + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CELL) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + } } @@ -129,4 +228,4 @@ psFree(view); - return file; + return output; } Index: trunk/ppSim/src/ppSimDetectionLimits.c =================================================================== --- trunk/ppSim/src/ppSimDetectionLimits.c (revision 18011) +++ trunk/ppSim/src/ppSimDetectionLimits.c (revision 18011) @@ -0,0 +1,47 @@ +# include "ppSim.h" + +// compare injected sources (PPSIM.SOURCES) and measured fake sources (PPSIM.FAKE.SOURCES) +// to determine detection limits and recovered magnitude errors +bool ppSimDetectionLimits (pmConfig *config, pmFPAview *view) { + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // XXX do we need to ask if the analysis was performed to 1st or 2nd detection limit? + float NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + psArray *injectedSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.SOURCES"); + psArray *measuredSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.FAKE.SOURCES"); + psAssert (injectedSources->n == measuredSources->n, "mis-match between injected and measured sources"); + + psVector *mag = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *dmag = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *SN = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *detect = psVectorAlloc (injectedSources->n, PS_TYPE_BOOL); + + // first measure, for each source, if it was detected (SN > limit) + // and the magnitude offset (dM = M_inject - M_measure) + for (int i = 0; i < injectedSources->n; i++) { + pmSource *injectSource = injectedSources->data[i]; + pmSource *measureSource = measuredSources->data[i]; + + SN->data.F32[i] = measureSource->peak->SN; + dmag->data.F32[i] = injectSource->psfMag - measureSource->psfMag; + mag->data.F32[i] = injectSource->psfMag; + detect->data.Bool[i] = (measureSource->peak->SN >= NSIGMA_PEAK); + } + + // generate a histogram consisting of the instrumental magnitude bin, the + // XXX what is resolution? (user parameter?) + + return true; +} Index: trunk/ppSim/src/ppSimDetections.c =================================================================== --- trunk/ppSim/src/ppSimDetections.c (revision 18011) +++ trunk/ppSim/src/ppSimDetections.c (revision 18011) @@ -0,0 +1,30 @@ +# include "ppSim.h" + +bool ppSimDetections (psImage *significance, psMetadata *recipe, psArray *sources) { + psAssert (sources, "programming error: ppSimDetections passed NULL sources"); + + bool status; + + // for each source, measure the significance of the peak at the given coordinate + // where does this get stored? + + // XXX need to get the effective Area from the PSF sigma + float SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "SIGMA_SMOOTH"); + psAssert (status, "SIGMA_SMOOTH missing: call psphotSignificanceImage first"); + float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH); + + int row0 = significance->row0; + int col0 = significance->col0; + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + psAssert (peak, "peak is not defined for the source"); + + peak->value = significance->data.F32[peak->y-row0][peak->x-col0]; + peak->SN = sqrt(peak->value*effArea); + + } + return true; +} Index: trunk/ppSim/src/ppSimInsertGalaxies.c =================================================================== --- trunk/ppSim/src/ppSimInsertGalaxies.c (revision 17915) +++ trunk/ppSim/src/ppSimInsertGalaxies.c (revision 18011) @@ -10,7 +10,4 @@ if (!galaxies->n) { return true; } - - // XXX is this needed? - // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources pmCell *cell = readout->parent; @@ -38,5 +35,5 @@ if (isnan(skyRate)) { float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok); - float scale = psMetadataLookupF32(&mdok, recipe, "SCALE"); assert (mdok); + float scale = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE"); assert (mdok); float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); @@ -56,8 +53,8 @@ int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y // determine the galaxy model - char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // galaxy model name + char *modelName = psMetadataLookupStr(&mdok, recipe, "GALAXY.MODEL"); // galaxy model name if (modelName == NULL) { modelName = defaultModel; Index: trunk/ppSim/src/ppSimInsertStars.c =================================================================== --- trunk/ppSim/src/ppSimInsertStars.c (revision 17915) +++ trunk/ppSim/src/ppSimInsertStars.c (revision 18011) @@ -15,7 +15,4 @@ if (!stars->n) { return true; } - - // XXX is this needed? - // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources pmCell *cell = readout->parent; @@ -43,5 +40,5 @@ if (isnan(skyRate)) { float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok); - float scale = psMetadataLookupF32(&mdok, recipe, "SCALE"); assert (mdok); + float scale = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE"); assert (mdok); float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); @@ -61,5 +58,5 @@ int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); @@ -117,10 +114,7 @@ source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / 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); - psFree(source); // Drop reference + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY); // Blow away the image parts of the source, which makes the memory explode @@ -132,8 +126,12 @@ RESET(source->psfFlux); RESET(source->blends); + + // add the sources to the source array + psArrayAdd (sources, 100,source); + psFree(source); // Drop reference } - // 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); + // NOTE: the pmFPAfile "PPSIM.OUTPUT" points at these sources + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "fake sources", sources); psFree(sources); @@ -141,3 +139,2 @@ return true; } - Index: trunk/ppSim/src/ppSimLoadForceSources.c =================================================================== --- trunk/ppSim/src/ppSimLoadForceSources.c (revision 18011) +++ trunk/ppSim/src/ppSimLoadForceSources.c (revision 18011) @@ -0,0 +1,113 @@ +# include "ppSim.h" + +// Reset a pointer: free and set to NULL +#define RESET(PTR) \ + psFree(PTR); \ + PTR = NULL; + +psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view) { + + bool status; + + pmFPAfile *outFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); + psAssert(outFile, "missing PPSIM.OUTPUT"); + psAssert(outFile->fpa, "missing fpr for PPSIM.OUTPUT"); + + pmChip *outChip = pmFPAviewThisChip (view, outFile->fpa); + pmPSF *psf = psMetadataLookupPtr (&status, outChip->analysis, "PSPHOT.PSF"); + assert (psf); + + psArray *spots = psMetadataLookupPtr(&status, outFile->fpa->analysis, "FORCED.SPOTS"); + PS_ASSERT_PTR_NON_NULL (spots, NULL); + + // in this program, we are looping over the output image, rather than the input as in ppImage + pmFPAfile *srcFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.REAL.SOURCES"); + psAssert(srcFile, "missing PPSIM.REAL.SOURCES"); + + // select the fpa and chip from PPSIM.REAL.SOURCES, which carries the astrometry in the headers + pmFPA *fpa = srcFile->fpa; + pmChip *chip = pmFPAviewThisChip (view, fpa); + + // we use the readout for the PPSIM.FORCE.CHIP, which contains the chip-mosaicked image(s) + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.CHIP"); + + // XXX we should have only one cell and readout per chip, right? + + // read WCS from existing output from psphot & psastro + // XXX can we be more flexible? ie, use the header of PPSIM.FORCE.CHIP? + bool bilevelAstrometry = false; + status = psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); + psAssert (status, "failed to set astrometry"); + + // the pixel scale is only used to set the focal-plane coordinate system + // setting this value to 1 makes this system have the units of pixels + float pixelScale = 1.0; + status = psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry); + psAssert (status, "failed to set astrometry"); + + float minX = 0.0; + float maxX = readout->image->numCols; + float minY = 0.0; + float maxY = readout->image->numRows; + + // the sources is a subset within range of this chip + psArray *sources = psArrayAllocEmpty (100); + + // Select the spots within range of this readout. Project the spots to this chip + for (int i = 0; i < spots->n; i++) { + pmAstromObj *obj = spots->data[i]; + + psProject (obj->TP, obj->sky, fpa->toSky); + psPlaneTransformApply (obj->FP, fpa->fromTPA, obj->TP); + psPlaneTransformApply (obj->chip, chip->fromFPA, obj->FP); + + // limit the X,Y range of the objs to the selected chip + if (obj->chip->x < minX) continue; + if (obj->chip->x > maxX) continue; + if (obj->chip->y < minY) continue; + if (obj->chip->y > maxY) continue; + + // convert the pmAstromObj to pmSource + + // instantiate a model for the PSF at this location, Io = 1.0 + pmModel *model = pmModelFromPSFforXY (psf, obj->chip->x, obj->chip->y, 1.0); + + // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) + // XXX use a fixed radius?? + // float radius = model->modelRadius (model->params, roughNoise); + // radius = PS_MAX (radius, 1.0); + float radius = 5.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; + + // these are not really needed since we will be fitting for them + source->psfMag = NAN; + source->errMag = NAN; + + // insert the source flux in the image + // XXX not sure the offset is really 0,0 + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, 0.0, 0.0); + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, 0.0, 0.0); + + // Blow away the image parts of the source, which makes the memory explode + RESET(source->pixels); + RESET(source->weight); + RESET(source->maskObj); + RESET(source->maskView); + RESET(source->modelFlux); + RESET(source->psfFlux); + RESET(source->blends); + + psArrayAdd (sources, 100, source); + psFree(source); // Drop local reference + } + psTrace ("psastro", 4, "Added %ld sources\n", sources->n); + + return sources; +} Index: trunk/ppSim/src/ppSimLoadSpots.c =================================================================== --- trunk/ppSim/src/ppSimLoadSpots.c (revision 17915) +++ trunk/ppSim/src/ppSimLoadSpots.c (revision 18011) @@ -1,21 +1,40 @@ -#include "ppSim.h" +# include "ppSim.h" +// Load the relevant forced-photometry positions for this field. This function does not determine +// the pixel coordinates, merely the celestial coordinates for sources in the general vicinity. The +// sources are saved on the fpa->analysis metadata as FORCED.SPOTS. We always create an entry; it +// will be empty if no sources were selected or forced photometry is not desired. bool ppSimLoadSpots (pmFPA *fpa, pmConfig *config) { - bool mdok; - assert (stars); + bool status; - psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe - bool forcedPhot = psMetadataLookupBool(&mdok, recipe, "FORCED.PHOT"); // Density of fakes - if (!forcedPhot) return true; + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - // Read catalogue stars using psastro - // XXX probably need to modify this... + bool forcedPhot = psMetadataLookupBool(&status, recipe, "FORCED.PHOT"); // Density of fakes + if (!forcedPhot) { + psArray *spots = psArrayAllocEmpty (1); + psMetadataAddArray (fpa->analysis, PS_LIST_TAIL, "FORCED.SPOTS", PS_META_REPLACE, "forced photometry positions", spots); + psFree (spots); // free the extra (local) reference; a copy remains on fpa->analysis + return true; + } + + // We read the catalogue stars using psastroLoadRefstars psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE); if (!astroRecipe) { psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE); - return NULL; + return false; } + + float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) + float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + + if (isnan(ra0) || isnan(dec0)) { + psError(PS_ERR_UNKNOWN, false, "image boresite coords not defined."); + return false; + } + + char *catdir = psMetadataLookupStr(NULL, recipe, "FORCED.CATDIR"); // Size of FPA @@ -24,45 +43,29 @@ psFree(bounds); - float x0fpa = 0.5*(bounds->x0 + bounds->x1); - float y0fpa = 0.5*(bounds->y0 + bounds->y1); + // modify the PSASTRO recipe to use the values desired for FORCED.PHOT. we can use a view on + // the PSASTRO recipe because we are not calling psastro elsewhere in this program. XXX need to + // use the WCS to define the overlap region + psMetadataAddStr(astroRecipe, PS_LIST_TAIL, "DVO.CATDIR", PS_META_REPLACE, "", catdir); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", ra0 - radius); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "RA_MAX", PS_META_REPLACE, "", ra0 + radius); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", dec0 - radius); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", dec0 + radius); + + // tell psastroLoadRefstars to ignore photcode and maglim + psMetadataAddStr(astroRecipe, PS_LIST_TAIL, "DVO.GETSTAR.PHOTCODE", PS_META_REPLACE, "", "NONE"); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DVO.GETSTAR.MAG.MAX", PS_META_REPLACE, "", NAN); - // XXX need to use the WCS to define the overlap region - 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) { + // load refstars from the catalog + psArray *spots = psastroLoadRefstars(config); + if (!spots || spots->n == 0) { psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars."); - psFree(refStars); - return NULL; + psFree(spots); + return false; } - psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n); + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", spots->n); - // convert the pmAstromObj sources to the storage format used by the CMF files + psMetadataAddArray (fpa->analysis, PS_LIST_TAIL, "FORCED.SPOTS", PS_META_REPLACE, "forced photometry positions", spots); + psFree (spots); // free the extra (local) reference; a copy remains on fpa->analysis - long oldSize = stars->n; - stars = psArrayRealloc (stars, refStars->n); - - // Conversion loop - for (long i = 0; i < refStars->n; i++) { - ppSimStar *star = ppSimStarAlloc (); - - pmAstromObj *ref = refStars->data[i]; // Reference star - star->ra = ref->sky->r; // RA of star - star->dec = ref->sky->d; // Dec of star - star->mag = ref->Mag; // Magnitude of star - - // Apply rotation, make FPA center of boresite - // convert the ra,dec to x,y using examples in psastro. - // star->x = cos(pa) * xi - sin(pa) * eta + x0fpa; - // star->y = sin(pa) * xi + cos(pa) * eta + y0fpa; - - stars->data[oldSize + i] = star; - - psTrace("ppSim", 10, "Adding catalogue star: %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux); - } - - // XXX these need to be saved on PSPHOT.INPUT.CMF return true; } Index: trunk/ppSim/src/ppSimLoadStars.c =================================================================== --- trunk/ppSim/src/ppSimLoadStars.c (revision 17915) +++ trunk/ppSim/src/ppSimLoadStars.c (revision 18011) @@ -1,3 +1,3 @@ -#include "ppSim.h" +# include "ppSim.h" bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config) { @@ -18,5 +18,5 @@ } - // XXX push these into the recipe in ppSimArguments() + // relevant metadata float zp = psMetadataLookupF32(NULL, recipe, "ZEROPOINT"); // Photometric zero point float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) @@ -24,5 +24,5 @@ float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) float seeing = psMetadataLookupF32(NULL, recipe, "SEEING"); // Seeing SIGMA (pixels) - float scale = psMetadataLookupF32(NULL, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time (sec) @@ -40,5 +40,5 @@ psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius); psArray *refStars = psastroLoadRefstars(config); - if (!refStars || refStars->n == 0) { + if (!refStars) { psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars."); psFree(refStars); Index: trunk/ppSim/src/ppSimLoop.c =================================================================== --- trunk/ppSim/src/ppSimLoop.c (revision 17915) +++ trunk/ppSim/src/ppSimLoop.c (revision 18011) @@ -1,34 +1,48 @@ #include "ppSim.h" + +# define ESCAPE(CODE,MSG) { \ + psError(CODE, false, MSG); \ + psFree (rng); \ + return false; } bool ppSimLoop(pmConfig *config) { + bool status; + + // XXX if we are supplying a PSF, then we should use that to specify the seeing. + // we will need to force the psf to be loaded here (deactivate everyone, activate psf, load + // it, then calculate seeing as appropriate). + PS_ASSERT_PTR_NON_NULL(config, PS_EXIT_PROG_ERROR); - pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, OUTPUT_FILE); // Output file + // in this program, we are looping over the output image, rather than the input as in ppImage + pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); // Output file assert(file); + pmFPA *fpa = file->fpa; // FPA for file assert(fpa); - 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 - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y - - // Load catalogue stars + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y + psArray *stars = psArrayAllocEmpty (1); - if (type == PPSIM_TYPE_OBJECT) { - ppSimLoadStars (stars, fpa, config); - } - - // Add random stars - if (type == PPSIM_TYPE_OBJECT) { - ppSimMakeStars (stars, fpa, config, rng); - } - - // Add random galaxies psArray *galaxies = psArrayAllocEmpty (1); if (type == PPSIM_TYPE_OBJECT) { - ppSimMakeGalaxies (galaxies, fpa, config, rng); + // Load forced-photometry positions (these are placed on fpa->analysis for use in ppSimPhotomReadout) + if (!ppSimLoadSpots (fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load forced-photometry spots"); + + // Load catalogue stars + if (!ppSimLoadStars (stars, fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load catalog stars"); + + // Add random stars + if (!ppSimMakeStars (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random stars"); + + // Add random galaxies + if (!ppSimMakeGalaxies (galaxies, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random galaxies"); } @@ -46,4 +60,7 @@ ppSimUpdateConceptsFPA (fpa, config); + if (fpa->hdu) { // XXX only do this if there is no INPUT image + if (!ppSimInitHeader(config, fpa, NULL, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); + } pmChip *chip; // Chip from FPA @@ -98,8 +115,6 @@ psVector *biasCols = ppSimMakeBiassec (cell, config); - for (int i = 0; i < cell->readouts->n; i++) { - - pmReadout *readout = cell->readouts->data[i]; - assert (readout); + pmReadout *readout; + while ((readout = pmFPAviewNextReadout (view, fpa, 1))) { // if we have not read in a weight or generated a fake image above, we need to @@ -117,19 +132,20 @@ } - psVector *biasRows = ppSimMakeBias (readout, config, rng); + psVector *biasRows = ppSimMakeBias (&status, readout, config, rng); + if (!status) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure"); if (type == PPSIM_TYPE_BIAS) goto done; - ppSimMakeDark (readout, config); + if (!ppSimMakeDark (readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure"); if (type == PPSIM_TYPE_DARK) goto done; - ppSimMakeSky (readout, expCorr, type, config); + if (!ppSimMakeSky (readout, expCorr, type, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating sky background"); if (type == PPSIM_TYPE_FLAT) goto done; if (type == PPSIM_TYPE_OBJECT) { - ppSimInsertStars (readout, expCorr, stars, config); + if (!ppSimInsertStars (readout, expCorr, stars, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inserting stars"); } if (type == PPSIM_TYPE_OBJECT) { - ppSimInsertGalaxies (readout, expCorr, galaxies, config); + if (!ppSimInsertGalaxies (readout, expCorr, galaxies, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inserting galaxies"); } @@ -137,10 +153,10 @@ done: - ppSimAddNoise(readout->image, readout->weight, cell, config, rng); - ppSimSaturate(readout, config); - - ppSimBadPixels(readout, config, rng); - - ppSimAddOverscan (readout, config, biasCols, biasRows, rng); + if (!ppSimAddNoise(readout->image, readout->weight, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise"); + if (!ppSimSaturate(readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem setting saturation levels"); + + if (!ppSimBadPixels(readout, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding bad pixels"); + + if (!ppSimAddOverscan (readout, config, biasCols, biasRows, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding overscan region"); psFree(biasRows); @@ -148,11 +164,15 @@ readout->parent->data_exists = true; readout->parent->parent->data_exists = true; + + // if there is an input image, merge it with the simulated image + if (!ppSimMergeReadouts (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); } psFree(biasCols); - ppSimUpdateConceptsCell (cell, config); + if (!ppSimUpdateConceptsCell (cell, config)) ESCAPE (PS_ERR_UNKNOWN, "problem updating cell concepts"); if (cell->hdu) { - ppSimInitHeader(config, NULL, NULL, cell); + // XXX only do this if there is no INPUT image? + if (!ppSimInitHeader(config, NULL, NULL, cell)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); } @@ -166,7 +186,13 @@ } + // XXX why no UpdateConceptsChip?? + if (chip->hdu) { - ppSimInitHeader(config, NULL, chip, NULL); - } + // XXX only do this if there is no INPUT image + if (!ppSimInitHeader(config, NULL, chip, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); + } + + // we perform photometry on the readouts of this chip in the output + if (!ppSimPhotom (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem performing photometry"); if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { @@ -174,16 +200,10 @@ psFree(rng); psFree(view); - // return PS_EXIT_SYS_ERROR; return false; } - } psFree(stars); psFree(galaxies); - - if (fpa->hdu) { - ppSimInitHeader(config, fpa, NULL, NULL); - } if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { @@ -191,5 +211,4 @@ psFree(rng); psFree(view); - // return PS_EXIT_SYS_ERROR; return false; } Index: trunk/ppSim/src/ppSimMakeBias.c =================================================================== --- trunk/ppSim/src/ppSimMakeBias.c (revision 17915) +++ trunk/ppSim/src/ppSimMakeBias.c (revision 18011) @@ -1,11 +1,16 @@ # include "ppSim.h" -psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) { +psVector *ppSimMakeBias (bool *status, pmReadout *readout, pmConfig *config, const psRandom *rng) { bool mdok; + + if (status) *status = true; pmCell *cell = readout->parent; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + bool bias = psMetadataLookupBool(&mdok, recipe, "BIAS"); // Generate a Bias? + if (!bias) return NULL; float biasLevel = psMetadataLookupF32(NULL, recipe, "BIAS.LEVEL"); // Bias level @@ -19,4 +24,5 @@ if (!mdok) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + *status = false; return NULL; } Index: trunk/ppSim/src/ppSimMakeBiassec.c =================================================================== --- trunk/ppSim/src/ppSimMakeBiassec.c (revision 17915) +++ trunk/ppSim/src/ppSimMakeBiassec.c (revision 18011) @@ -6,4 +6,7 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + bool bias = psMetadataLookupBool(&mdok, recipe, "BIAS"); // Generate a Bias? + if (!bias) return NULL; psList *biassec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // Bias regions Index: trunk/ppSim/src/ppSimMakeDark.c =================================================================== --- trunk/ppSim/src/ppSimMakeDark.c (revision 17915) +++ trunk/ppSim/src/ppSimMakeDark.c (revision 18011) @@ -11,6 +11,20 @@ psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe - float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate - float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + bool dark = psMetadataLookupBool(&mdok, recipe, "DARK"); // Generate a DARK? + if (!dark) return true; + + float darkRate = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate + if (isnan(darkRate)) { + psError(PS_ERR_UNKNOWN, false, "missing DARK.RATE\n"); + return false; + } + + float expTime = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time + if (isnan(expTime)) { + psError(PS_ERR_UNKNOWN, false, "missing EXPTIME\n"); + return false; + } + + psTrace("ppSim", 6, "darkRate: %f, expTime: %f\n", darkRate, expTime); for (int y = 0; y < signal->numRows; y++) { Index: trunk/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- trunk/ppSim/src/ppSimMakeGalaxies.c (revision 17915) +++ trunk/ppSim/src/ppSimMakeGalaxies.c (revision 18011) @@ -9,5 +9,5 @@ bool galaxyFake = psMetadataLookupBool(&mdok, recipe, "GALAXY.FAKE"); // Density of fakes - if (!galaxyFake) return NULL; + if (!galaxyFake) return true; float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope @@ -35,5 +35,5 @@ float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point float seeing = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(&mdok, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) float skyRate = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { @@ -42,5 +42,5 @@ } - if (galaxyDensity <= 0) return NULL; + if (galaxyDensity <= 0) return true; // Size of FPA @@ -66,5 +66,5 @@ psLogMsg("ppSim", PS_LOG_INFO, "Image noise is above brightest random galaxy --- no random galaxy added."); - return NULL; + return true; } @@ -138,4 +138,4 @@ } } - return galaxies; + return true; } Index: trunk/ppSim/src/ppSimMakeSky.c =================================================================== --- trunk/ppSim/src/ppSimMakeSky.c (revision 17915) +++ trunk/ppSim/src/ppSimMakeSky.c (revision 18011) @@ -1,5 +1,7 @@ # include "ppSim.h" -// XXX add bounds to the inputs? +// this function sets the skyRate to a value for the night sky (SKY.RATE or SKY.MAGS) or for a +// flat-field image (FLAT.RATE). Include a shutter correction and a scattered light source + bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config) { @@ -15,14 +17,22 @@ psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time - float flatSigma = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient - float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate - float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time - float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate - if (isnan(skyRate)) { - float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); - float scale = psMetadataLookupF32(&status, recipe, "SCALE"); assert (status); - float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + bool sky = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux? + bool flat = psMetadataLookupBool(&status, recipe, "FLAT"); // Apply flat-field term? + + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + + float flatSigma = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient + float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate + float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time + float scatterFrac = psMetadataLookupF32(&status, recipe, "SCATTER.FRAC"); // scattered light fraction (max) + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + if (!isnan(skyMags)) { + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + if (type == PPSIM_TYPE_FLAT) { + skyRate = flatRate; } @@ -37,12 +47,14 @@ int yParityCell = psMetadataLookupS32(&status, cell->concepts, "CELL.YPARITY"); - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y // Size of FPA psRegion *bounds = ppSimFPABounds (fpa); + int dXfpa = bounds->x1 - bounds->x0; + int dYfpa = bounds->y1 - bounds->y0; // 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); + x0Chip -= 0.5 * dXfpa; + y0Chip -= 0.5 * dYfpa; for (int y = 0; y < signal->numRows; y++) { @@ -58,21 +70,30 @@ float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0; - // Gaussian flat-field over the FPA - float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) / - flatSigma / sqrtf(2.0 * M_PI); + // Gaussian flat-field over the FPA with flatValue = 1.0 at the field center + float flatValue = 1.0; + if (flat) { + // we make the flat-field have a response of 1.0 at the field center (like a vignetting) + flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))); + } - 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; - } + float scatterRate = 0.0; - expCorr->data.F32[y][x] = realExpTime / expTime; + if (sky) { + // add a scattered light term to the flat-field images (no + if (type == PPSIM_TYPE_FLAT) { + float xF = 2.0*(xFPA / dXfpa) - 1.0; + scatterRate = scatterFrac * PS_SQR(xF); + } - // Sky background - float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky - signal->data.F32[y][x] += skyFlux; - variance->data.F32[y][x] += skyFlux; + // Sky background + float skyFlux = (skyRate * (flatValue + scatterRate)) * realExpTime; // Flux from sky + signal->data.F32[y][x] += skyFlux; + variance->data.F32[y][x] += skyFlux; + } + + // used later to modify the star and galaxy photometry + if (expCorr) { + expCorr->data.F32[y][x] = realExpTime / expTime; + } // TO DO: Add fringes Index: trunk/ppSim/src/ppSimMakeStars.c =================================================================== --- trunk/ppSim/src/ppSimMakeStars.c (revision 17915) +++ trunk/ppSim/src/ppSimMakeStars.c (revision 18011) @@ -1,3 +1,7 @@ # include "ppSim.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, MSG); \ + return false; } bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { @@ -23,5 +27,11 @@ float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) - float scale = psMetadataLookupF32(&status, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + + if (isnan(darkRate)) darkRate = 0.0; + if (isnan(expTime)) ESCAPE("EXPTIME is not defined"); + if (isnan(zp)) ESCAPE("ZEROPOINT is not defined"); + if (isnan(seeing)) ESCAPE("SEEING is not defined"); + if (isnan(scale)) ESCAPE("PIXEL.SCALE is not defined"); bool flatLum = psMetadataLookupBool(&status,recipe, "STARS.FLAT.LUM"); // were real stars generated? @@ -47,5 +57,10 @@ refMag = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.MAG.PEAK"); // Star brightest magnitude refSum = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.SUM.PEAK"); // Star brightest magnitude - assert (status); + + // if we tried and failed to load reference stars, set more artificial limits + if (!status) { + refMag = brightMag; + refSum = 1; + } } else { refMag = brightMag; Index: trunk/ppSim/src/ppSimMergeReadouts.c =================================================================== --- trunk/ppSim/src/ppSimMergeReadouts.c (revision 18011) +++ trunk/ppSim/src/ppSimMergeReadouts.c (revision 18011) @@ -0,0 +1,48 @@ +# include "ppSim.h" + +// XXX add bounds to the inputs? +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view) { + + // bool mdok; + + // if we have an input image, we need to add the synthetic data on top of it below + // XXX we may potentially use this input file to skip missing elements + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSIM.INPUT"); + if (!input) return true; + + pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa); + if (!inReadout) return true; + + if (!inReadout->weight) { + if (!pmReadoutGenerateWeight(inReadout, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating weight"); + return false; + } + } + + // output must exist or we made a programming error + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); // Output file + assert(output); + + // XXX require outReadout? + pmReadout *outReadout = pmFPAviewThisReadout (view, output->fpa); + if (!outReadout) return true; + + psImage *inSignal = inReadout->image; + psImage *inVariance = inReadout->weight; + + psImage *outSignal = outReadout->image; + psImage *outVariance = outReadout->weight; + + assert (inSignal->numRows == outSignal->numRows); + assert (inSignal->numCols == outSignal->numCols); + + for (int y = 0; y < inSignal->numRows; y++) { + for (int x = 0; x < inSignal->numCols; x++) { + outSignal->data.F32[y][x] += inSignal->data.F32[y][x]; + outVariance->data.F32[y][x] += inVariance->data.F32[y][x]; + } + } + return true; +} + Index: trunk/ppSim/src/ppSimMergeSources.c =================================================================== --- trunk/ppSim/src/ppSimMergeSources.c (revision 18011) +++ trunk/ppSim/src/ppSimMergeSources.c (revision 18011) @@ -0,0 +1,17 @@ +# include "ppSim.h" + +psArray *ppSimMergeSources (psArray *in1, psArray *in2) { + + psArray *out = psArrayAlloc (in1->n + in2->n); + + for (int i = 0; i < in1->n; i++) { + out->data[i] = psMemIncrRefCounter (in1->data[i]); + } + + int nOff = in1->n; + for (int i = nOff; i < out->n; i++) { + out->data[i] = psMemIncrRefCounter (in2->data[i-nOff]); + } + + return out; +} Index: trunk/ppSim/src/ppSimPhotom.c =================================================================== --- trunk/ppSim/src/ppSimPhotom.c (revision 17915) +++ trunk/ppSim/src/ppSimPhotom.c (revision 18011) @@ -8,18 +8,32 @@ pmReadout *readout; + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes + if (!doPhotom) return true; + psphotModelClassInit (); - // find or define a pmFPAfile PSPHOT.INPUT - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + int blankMask = 0; // XXX not sure what this should be set to... + ppSimMosaicChip(config, blankMask, view, "PPSIM.FORCE.CHIP", "PPSIM.INPUT"); + ppSimMosaicChip(config, blankMask, view, "PPSIM.FAKE.CHIP", "PPSIM.OUTPUT"); + + // use PPSIM.FAKE.CHIP as the base since it is guaranteed to exist (derived from PPSIM.OUTPUT) + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PPSIM.FAKE.CHIP"); if (!status) { - psError(PSPHOT_ERR_CONFIG, false, "PSPHOT.INPUT I/O file is not defined"); + psError(PSPHOT_ERR_CONFIG, false, "PPSIM.FAKE.CHIP I/O file is not defined"); return false; } - // we make a new copy of the output chip to keep psphot from modifying the output image - // XXX is this needed for ppSim? (yes, unless we do not do photometry and make a simulated image) - pmChip *oldChip = pmFPAviewThisChip (view, input->src); - pmChip *newChip = pmFPAviewThisChip (view, input->fpa); - pmChipCopy (newChip, oldChip); + // XXX If we want to be able to write out both the positive and residual images, we will + // need to define an additional pmFPAfile container, and then make a copy as is done for + // PPIMAGE.CHIP -> PSPHOT.INPUT in ppImage. At the moment, PPSIM.CHIP will write out the + // mosaicked image with the sources subtracted. + + // XXX if input -> PSPHOT.INPUT, which is tied to PPSIM.CHIP with pmFPAfileDefineFromFile, + // then the following code creates the chip copy: + // pmChip *oldChip = pmFPAviewThisChip (view, input->src); + // pmChip *newChip = pmFPAviewThisChip (view, input->fpa); + // pmChipCopy (newChip, oldChip); // iterate over the cells and readout for this chip @@ -33,8 +47,13 @@ // run the actual photometry analysis - if (!psphotReadout (config, view)) { + if (!ppSimPhotomReadoutFake (config, view)) { psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); return false; } + if (!ppSimPhotomReadoutForce (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + return false; + } + // ppSimDetectionLimits (config, view); } } Index: trunk/ppSim/src/ppSimPhotomFiles.c =================================================================== --- trunk/ppSim/src/ppSimPhotomFiles.c (revision 18011) +++ trunk/ppSim/src/ppSimPhotomFiles.c (revision 18011) @@ -0,0 +1,83 @@ +# include "ppSim.h" + +// define the needed / desired I/O files +bool ppSimPhotomFiles (pmConfig *config, pmFPAfile *fakeFile, pmFPAfile *forceFile) { + + bool status = false; + + // select recipe options supplied on command line + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + + // optionally save the background model (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, DX, DY, "PSPHOT.BACKMDL"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKMDL"); + return false; + } + output->save = true; + } + // optionally save the background model's standard deviation (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL.STDEV")) { + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, DX, DY, "PSPHOT.BACKMDL.STDEV"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKMDL.STDEV"); + return false; + } + output->save = true; + } + // optionally save the full background image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKGND")) { + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, 1, 1, "PSPHOT.BACKGND"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKGND"); + return false; + } + output->save = true; + } + // optionally save the background-subtracted image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKSUB")) { + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, 1, 1, "PSPHOT.BACKSUB"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKSUB"); + return false; + } + output->save = true; + } + + // the input sources are required for sensible measurements of the fake or force sources + pmFPAfile *cmffile = pmFPAfileDefineFromArgs (&status, config, "PPSIM.REAL.SOURCES", "INPUT.SOURCES"); + if (!cmffile) { + psError(PS_ERR_UNKNOWN, false, "input sources not found; required for photometry (or error with PPSIM.REAL.SOURCES filerule)"); + return false; + } + + // the fake sources are carried on the fakeFile->fpa structures (PPSIM.FAKE.CHIP) + pmFPAfile *outsources = pmFPAfileDefineOutput (config, fakeFile->fpa, "PPSIM.FAKE.SOURCES"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PPSIM.FAKE.SOURCES"); + return false; + } + outsources->save = true; + + // the forced sources are carried on the forceFile->fpa structures (PPSIM.FORCE.CHIP) + pmFPAfile *outforced = pmFPAfileDefineOutput (config, forceFile->fpa, "PPSIM.FORCE.SOURCES"); + if (!outforced) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PPSIM.FORCE.SOURCES"); + return false; + } + outforced->save = true; + + // we bind the psffile to the chip-mosaicked image so we have a reference to the image + pmFPAfile *psffile = pmFPAfileBindFromArgs(&status, fakeFile, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!psffile) { + psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.PSF.LOAD"); + return false; + } + + return true; +} Index: trunk/ppSim/src/ppSimPhotomReadout.c =================================================================== --- trunk/ppSim/src/ppSimPhotomReadout.c (revision 18011) +++ trunk/ppSim/src/ppSimPhotomReadout.c (revision 18011) @@ -0,0 +1,184 @@ +# include "ppSim.h" + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) { + + pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename); + PS_ASSERT_PTR_NON_NULL (readout, NULL); + + psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + return sources; +} + +bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe) { + + bool status; + + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + if (!status) return NULL; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + } + return true; +} + +bool ppSimPhotomReadout(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + +# if 0 + // set the photcode for this image + if (!psphotAddPhotcode (recipe, config, view, "PPSIM.CHIP")) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } +# endif + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile + psArray *realSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES"); + PS_ASSERT_PTR_NON_NULL (realSources, false); + + // the newly injected fake sources are saved on the PPSIM.OUTPUT pmFPAfile + psArray *injectedSources = ppSimSelectSources (config, view, "PPSIM.OUTPUT"); + PS_ASSERT_PTR_NON_NULL (injectedSources, false); + + // XXX need to define the source pixels + ppSimDefinePixels (realSources, readout, recipe); + ppSimDefinePixels (injectedSources, readout, recipe); + + // make a copy of the fake-source parameters so these can be saved independently from the + // measured fake-source parameters + psArray *fakeSources = psArrayAlloc (injectedSources->n); + for (int i = 0; i < injectedSources->n; i++) { + fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]); + } + + // load the forced source lists + psArray *forceSources = ppSimLoadForceSources (config, view); + psAssert (forceSources, "failed to load force photometry sources"); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // remove all sources + psphotRemoveAllSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PPSIM.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PPSIM.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, fakeSources); + ppSimDetections (significance, recipe, forceSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAllSources (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, fakeSources, recipe, psf); + psphotGuessModels (readout, forceSources, recipe, psf); + + psArray *sources = NULL; + + // linear fit to real + fake sources + sources = ppSimMergeSources (realSources, fakeSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + + // linear fit to real + forced sources (if any were selected) + if (forceSources->n) { + sources = ppSimMergeSources (realSources, forceSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + } + + // XXX do we need to measure aperture photometry corrections? + // XXX do we store the pre and post correction magnitudes? + // XXX for the fake sources, these must be identically zero + // XXX for the force sources, need the apresid to put the mags on the correct system + if (0) { + if (!psphotApResid (readout, forceSources, recipe, psf)) { + psLogMsg ("psphot", 3, "failed on psphotApResid"); + return psphotReadoutCleanup (config, readout, recipe, NULL, psf, NULL); + } + } + + // calculate source magnitudes (for which set??) + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(fakeSources, recipe, psf, background); + psphotMagnitudes(forceSources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (fakeSources); + psphotSourceFreePixels (forceSources); + + // create the exported-metadata and free local data + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + // add forceSources to the readout (real and fake already there) + // these are outputs: we need to generate the fpa structure + + pmCell *forceCell = pmFPAfileThisCell (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceCell, "no cell?"); + pmChip *forceChip = pmFPAfileThisChip (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceChip, "no chip?"); + pmReadout *forceReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.SOURCES"); + if (!forceReadout) { + forceReadout = pmReadoutAlloc (forceCell); + psFree (forceReadout); // there is a copy on 'cell' as well + } + psAssert (forceReadout, "no forceReadout?"); + pmChipSetDataStatus (forceChip, true); + psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "forced photometry ", forceSources); + psFree (forceSources); + + pmCell *fakeCell = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?"); + pmChip *fakeChip = pmFPAfileThisChip (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeChip, "no chip?"); + pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES"); + if (!fakeReadout) { + fakeReadout = pmReadoutAlloc (fakeCell); + psFree (fakeReadout); // there is a copy on 'cell' as well + } + psAssert (fakeReadout, "no fakeReadout?"); + pmChipSetDataStatus (fakeChip, true); + psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + + return true; +} Index: trunk/ppSim/src/ppSimPhotomReadoutFake.c =================================================================== --- trunk/ppSim/src/ppSimPhotomReadoutFake.c (revision 18011) +++ trunk/ppSim/src/ppSimPhotomReadoutFake.c (revision 18011) @@ -0,0 +1,124 @@ +# include "ppSim.h" + +bool ppSimPhotomReadoutFake(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + +# if 0 + // set the photcode for this image + if (!psphotAddPhotcode (recipe, config, view, "PPSIM.FAKE.CHIP")) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } +# endif + + // *** in this section, perform the photometry for real + fake sources on PPSIM.FAKE.CHIP *** + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile + psArray *realMeasuredSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES"); + PS_ASSERT_PTR_NON_NULL (realMeasuredSources, false); + + // the newly injected fake sources are saved on the PPSIM.OUTPUT pmFPAfile + psArray *injectedSources = ppSimSelectSources (config, view, "PPSIM.OUTPUT"); + PS_ASSERT_PTR_NON_NULL (injectedSources, false); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // XXX need to define the source pixels + ppSimDefinePixels (realMeasuredSources, readout, recipe); + ppSimDefinePixels (injectedSources, readout, recipe); + + // make a copy of the injected fake-source parameters so these can be saved independently + // from the measured fake-source parameters + psArray *fakeSources = psArrayAlloc (injectedSources->n); + for (int i = 0; i < injectedSources->n; i++) { + fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]); + } + + // make a copy of the measured real-source parameters so these can be fitted here and in + // ppSimPhotomReadoutForce, if desired. + psArray *realSources = psArrayAlloc (realMeasuredSources->n); + for (int i = 0; i < realMeasuredSources->n; i++) { + realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]); + } + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // remove all sources + psphotRemoveAllSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PPSIM.FAKE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PPSIM.FAKE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, fakeSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAllSources (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, fakeSources, recipe, psf); + + // linear fit to real + fake sources + psArray *sources = ppSimMergeSources (realSources, fakeSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + + // calculate source magnitudes (for which set??) + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(fakeSources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (fakeSources); + + // create the exported-metadata and free local data + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + pmCell *fakeCell = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?"); + pmChip *fakeChip = pmFPAfileThisChip (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeChip, "no chip?"); + pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES"); + if (!fakeReadout) { + fakeReadout = pmReadoutAlloc (fakeCell); + psFree (fakeReadout); // there is a copy on 'cell' as well + } + psAssert (fakeReadout, "no fakeReadout?"); + pmChipSetDataStatus (fakeChip, true); + psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + + return true; +} Index: trunk/ppSim/src/ppSimPhotomReadoutForce.c =================================================================== --- trunk/ppSim/src/ppSimPhotomReadoutForce.c (revision 18011) +++ trunk/ppSim/src/ppSimPhotomReadoutForce.c (revision 18011) @@ -0,0 +1,117 @@ +# include "ppSim.h" + +bool ppSimPhotomReadoutForce(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + +# if 0 + // set the photcode for this image + if (!psphotAddPhotcode (recipe, config, view, "PPSIM.FAKE.CHIP")) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } +# endif + + // *** in this section, perform the photometry for real + force sources on PPSIM.FORCE.CHIP *** + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile + psArray *realMeasuredSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES"); + PS_ASSERT_PTR_NON_NULL (realMeasuredSources, false); + + // load the forced source lists + psArray *forceSources = ppSimLoadForceSources (config, view); + psAssert (forceSources, "failed to load force photometry sources"); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // XXX need to define the source pixels + ppSimDefinePixels (realMeasuredSources, readout, recipe); + ppSimDefinePixels (forceSources, readout, recipe); + + // make a copy of the measured real-source parameters so these can be fitted here and in + // ppSimPhotomReadoutForce, if desired. + psArray *realSources = psArrayAlloc (realMeasuredSources->n); + for (int i = 0; i < realMeasuredSources->n; i++) { + realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]); + } + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // remove all sources + psphotRemoveAllSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PPSIM.FORCE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PPSIM.FORCE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, forceSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAllSources (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, forceSources, recipe, psf); + + // linear fit to real + force sources + psArray *sources = ppSimMergeSources (realSources, forceSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + + // calculate source magnitudes (for which set??) + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(forceSources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (forceSources); + + // create the exported-metadata and free local data + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + pmCell *forceCell = pmFPAfileThisCell (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceCell, "no cell?"); + pmChip *forceChip = pmFPAfileThisChip (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceChip, "no chip?"); + pmReadout *forceReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.SOURCES"); + if (!forceReadout) { + forceReadout = pmReadoutAlloc (forceCell); + psFree (forceReadout); // there is a copy on 'cell' as well + } + psAssert (forceReadout, "no forceReadout?"); + pmChipSetDataStatus (forceChip, true); + psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "force photometry ", forceSources); + + return true; +} Index: trunk/ppSim/src/ppSimRandomGaussian.c =================================================================== --- trunk/ppSim/src/ppSimRandomGaussian.c (revision 18011) +++ trunk/ppSim/src/ppSimRandomGaussian.c (revision 18011) @@ -0,0 +1,102 @@ +# include "ppSim.h" + +static int Ngaussint = 0; +static double *gaussint = NULL; + +extern double drand48(); + +double p_ppSimGaussian (double x, double mean, double sigma) { + + double f; + + f = exp (-0.5 * PS_SQR(x - mean) / PS_SQR(sigma)) / sqrt(2 * M_PI * PS_SQR(sigma)); + + return (f); + +} + +void ppSimRandomGaussianFree() +{ + psFree (gaussint); + return; +} + +void ppSimRandomGaussianAlloc (int Nbin) { + + gaussint = (double *) psAlloc(Nbin*sizeof(double)); + return; +} + +/* integrate a gaussian from -5 sigma to +5 sigma */ +void p_ppSimRandomGaussianInit () { + + int i; + long A, B; + double val, x, dx, dx1, dx2, dx3, df; + double mean, sigma; + + /* no need to generate this if it already exists */ + if (gaussint) return; + + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + + Ngaussint = 0x1000; + ppSimRandomGaussianAlloc (Ngaussint + 1); + + val = 0; + dx = 1.0 / Ngaussint; + dx1 = dx / 3.0; + dx2 = 2.0*dx/3.0; + dx3 = dx; + mean = 0.0; + sigma = 1.0; + + for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx) { + df = (3.0*p_ppSimGaussian(x , mean, sigma) + + 9.0*p_ppSimGaussian(x+dx1, mean, sigma) + + 9.0*p_ppSimGaussian(x+dx2, mean, sigma) + + 3.0*p_ppSimGaussian(x+dx3, mean, sigma)) * (dx1/8.0); + val += df; + if (val > (i + 0.5) / (double) Ngaussint) { + gaussint[i] = x + dx / 2.0; + i++; + } + } +} + +// XXX we are using drand48() rather than the random var supplied by rnd +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma) { + + int i; + double y; + + if (gaussint == NULL) { + p_ppSimRandomGaussianInit (); + } + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]*sigma + mean; + + return (y); + +} + +// XXX we are using drand48() rather than the random var supplied by rnd +double ppSimRandomGaussianNorm (const psRandom *rnd) { + + int i; + double y; + + if (gaussint == NULL) { + p_ppSimRandomGaussianInit (); + } + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]; + + return (y); +} Index: trunk/ppSim/src/ppSimSetPSF.c =================================================================== --- trunk/ppSim/src/ppSimSetPSF.c (revision 17915) +++ trunk/ppSim/src/ppSimSetPSF.c (revision 18011) @@ -21,5 +21,5 @@ float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) - char *psfModelName = psMetadataLookupStr(&status, config->arguments, "PSF.MODEL"); // Name of PSF model + char *psfModelName = psMetadataLookupStr(&status, recipe, "PSF.MODEL"); // Name of PSF model if (psfModelName == NULL) { psfModelName = defaultModel; Index: trunk/ppSim/src/ppSimUtils.c =================================================================== --- trunk/ppSim/src/ppSimUtils.c (revision 17915) +++ trunk/ppSim/src/ppSimUtils.c (revision 18011) @@ -15,8 +15,8 @@ float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) - float scale = psMetadataLookupF32(NULL, recipe, "SCALE"); // plate scale in arcsec / pixel + float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE"); // plate scale in arcsec / pixel scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y float x0 = 0.0, y0 = 0.0; @@ -95,4 +95,28 @@ } +char *ppSimTypeToString (ppSimType type) { + + char *typeStr; + + switch (type) { + case PPSIM_TYPE_BIAS: typeStr = psStringCopy ("BIAS"); break; + case PPSIM_TYPE_DARK: typeStr = psStringCopy ("DARK"); break; + case PPSIM_TYPE_FLAT: typeStr = psStringCopy ("FLAT"); break; + case PPSIM_TYPE_OBJECT: typeStr = psStringCopy ("OBJECT"); break; + default: + psAbort("Should never get here."); + } + return (typeStr); +} + +ppSimType ppSimTypeFromString (char *typeStr) { + + if (!strcasecmp (typeStr, "BIAS")) return PPSIM_TYPE_BIAS; + if (!strcasecmp (typeStr, "DARK")) return PPSIM_TYPE_DARK; + if (!strcasecmp (typeStr, "FLAT")) return PPSIM_TYPE_FLAT; + if (!strcasecmp (typeStr, "OBJECT")) return PPSIM_TYPE_OBJECT; + psAbort("Should never get here."); +} + bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config) { @@ -103,28 +127,14 @@ float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time - const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name + const char *filter = psMetadataLookupStr(NULL, recipe, "FILTER"); // Filter name if (!filter) { filter = "NONE"; } - ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate - - // 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); + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + + 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.EXPOSURE", 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); @@ -139,5 +149,5 @@ psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe - int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time @@ -153,2 +163,139 @@ return true; } + +bool ppSimRecipeValidation (pmConfig *config) { + + bool status; + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + int binning = psMetadataLookupS32(&status, recipe, "BINNING"); // Binning in x and y + if (binning <= 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Binning (%d) is non-positive.", binning); + exit(PS_EXIT_CONFIG_ERROR); + } + return true; +} + +// Get a value from the command-line arguments and add it to recipe options +float ppSimArgToRecipeF32(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ) +{ + bool myStatus; + float value = psMetadataLookupF32(&myStatus, arguments, argName); // Value of interest + if (status) { *status = myStatus; } + if (isnan(value)) return value; + + psMetadataAddF32(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +// Get a value from the command-line arguments and add it to recipe options +int ppSimArgToRecipeS32(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ) +{ + bool myStatus; + int value = psMetadataLookupS32(&myStatus, arguments, argName); // Value of interest + if (status) { *status = myStatus; } + + psMetadataAddS32(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +// Get a value from the command-line arguments and add it to recipe options +bool ppSimArgToRecipeBool(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ) +{ + bool myStatus; + bool value = psMetadataLookupS32(&myStatus, arguments, argName); // Value of interest + if (status) { *status = myStatus; } + + psMetadataAddBool(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +// Get a value from the command-line arguments and add it to recipe options +char *ppSimArgToRecipeStr(bool *status, + psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ) +{ + bool myStatus; + + char *value = psMetadataLookupStr(&myStatus, arguments, argName); // Value of interest + if (status) { + *status = myStatus; + } + psMetadataAddStr(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +float ppSimGetZeroPoint (psMetadata *recipe, char *filter) { + + bool mdok; + float zp; + + // use the filter to get the zeropoint from the recipe + psMetadataItem *zpItem = psMetadataLookup (recipe, "ZEROPTS"); + // check that item is multi... + + psArray *entries = psListToArray (zpItem->data.list); + + // search for matching filter + for (int i = 0; i < entries->n; i++) { + psMetadataItem *item = entries->data[i]; + psMetadata *entry = item->data.V; + + char *filterName = psMetadataLookupStr (&mdok, entry, "FILTER"); + assert (filterName); + + if (strcmp(filterName, filter)) continue; + + zp = psMetadataLookupF32 (&mdok, entry, "ZERO_PT"); + assert (mdok); + psFree (entries); + return zp; + } + psFree (entries); + return NAN; +} + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) { + + pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename); + PS_ASSERT_PTR_NON_NULL (readout, NULL); + + psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + return sources; +} + +bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe) { + + bool status; + + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + if (!status) return NULL; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + } + return true; +} +