Changeset 18011

Timestamp:

Jun 8, 2008, 4:03:31 PM (18 years ago)

Author:

eugene

Message:

update from changes on eam_branch_20080511 : adding photometry of fake sources, force photometry; major cleanups

Location:

trunk/ppSim

Files:

• configure.ac (modified) (1 diff)
• notes.txt (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/ppSim.c (modified) (3 diffs)
• src/ppSim.h (modified) (5 diffs)
• src/ppSimAddNoise.c (modified) (1 diff)
• src/ppSimAddOverscan.c (modified) (1 diff)
• src/ppSimArguments.c (modified) (6 diffs)
• src/ppSimCreate.c (modified) (2 diffs)
• src/ppSimDetectionLimits.c (added)
• src/ppSimDetections.c (added)
• src/ppSimInsertGalaxies.c (modified) (3 diffs)
• src/ppSimInsertStars.c (modified) (6 diffs)
• src/ppSimLoadForceSources.c (added)
• src/ppSimLoadSpots.c (modified) (2 diffs)
• src/ppSimLoadStars.c (modified) (4 diffs)
• src/ppSimLoop.c (modified) (9 diffs)
• src/ppSimMakeBias.c (modified) (2 diffs)
• src/ppSimMakeBiassec.c (modified) (1 diff)
• src/ppSimMakeDark.c (modified) (1 diff)
• src/ppSimMakeGalaxies.c (modified) (5 diffs)
• src/ppSimMakeSky.c (modified) (4 diffs)
• src/ppSimMakeStars.c (modified) (3 diffs)
• src/ppSimMergeReadouts.c (added)
• src/ppSimMergeSources.c (added)
• src/ppSimPhotom.c (modified) (2 diffs)
• src/ppSimPhotomFiles.c (added)
• src/ppSimPhotomReadout.c (added)
• src/ppSimPhotomReadoutFake.c (added)
• src/ppSimPhotomReadoutForce.c (added)
• src/ppSimRandomGaussian.c (added)
• src/ppSimSetPSF.c (modified) (1 diff)
• src/ppSimUtils.c (modified) (5 diffs)

trunk/ppSim/configure.ac

@@ -17 +17 @@
 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

trunk/ppSim/notes.txt

@@ +1 @@
+
+
+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:

trunk/ppSim/src/Makefile.am

@@ -1 +1 @@
 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
 

trunk/ppSim/src/ppSim.c

@@ -1 +1 @@
 # 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 +17 @@
     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 +39 @@
     psLibFinalize();
 
-    return PS_EXIT_SUCCESS;
+    exit (PS_EXIT_SUCCESS);
 }

trunk/ppSim/src/ppSim.h

@@ -13 +13 @@
 #include <psmodules.h>
 #include <psastro.h>
+#include <psphot.h>
 
 #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 +71 @@
 
 /// 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 +79 @@
 /// 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 +130 @@
 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 +148 @@
 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

trunk/ppSim/src/ppSimAddNoise.c

@@ -32 +32 @@
             // (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
         }

trunk/ppSim/src/ppSimAddOverscan.c

@@ -4 +4 @@
 
     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

trunk/ppSim/src/ppSimArguments.c

@@ -1 @@
-#include "ppSim.h"
+# include "ppSim.h"
 
 // Print usage information and die
@@ -18 +18 @@
 }
 
-// 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 +60 @@
 
     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 +143 @@
         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 +197 @@
         // 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 +224 @@
    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
 */

trunk/ppSim/src/ppSimCreate.c

@@ -1 +1 @@
 # 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 +228 @@
     psFree(view);
 
-    return file;
+    return output;
 }

trunk/ppSim/src/ppSimInsertGalaxies.c

@@ -10 +10 @@
     
     if (!galaxies->n) { return true; }
-
-    // XXX is this needed?
-    // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
 
     pmCell *cell = readout->parent;
@@ -38 +35 @@
     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 +53 @@
     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;

trunk/ppSim/src/ppSimInsertStars.c

@@ -15 +15 @@
 
     if (!stars->n) { return true; }
-
-    // XXX is this needed?
-    // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
 
     pmCell *cell = readout->parent;
@@ -43 +40 @@
     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 +58 @@
     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 +114 @@
         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 +126 @@
         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 +139 @@
     return true;
 }
-

trunk/ppSim/src/ppSimLoadSpots.c

@@ -1 @@
-#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 +43 @@
     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;
 }

trunk/ppSim/src/ppSimLoadStars.c

@@ -1 @@
-#include "ppSim.h"
+# include "ppSim.h"
 
 bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config) {
@@ -18 +18 @@
     }
 
-    // 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 +24 @@
     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 +40 @@
     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);

trunk/ppSim/src/ppSimLoop.c

@@ -1 +1 @@
 #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 +60 @@
 
     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 +115 @@
             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 +132 @@
                 }
 
-                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 +153 @@
 
             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 +164 @@
                 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 +186 @@
         }
 
+        // 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 +200 @@
             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 +211 @@
         psFree(rng);
         psFree(view);
-        // return PS_EXIT_SYS_ERROR;
         return false;
     }

trunk/ppSim/src/ppSimMakeBias.c

@@ -1 +1 @@
 # 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 +24 @@
         if (!mdok) {
             psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
+            *status = false;
             return NULL;
         }

trunk/ppSim/src/ppSimMakeBiassec.c

@@ -6 +6 @@
 
     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

trunk/ppSim/src/ppSimMakeDark.c

@@ -11 +11 @@
     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++) {

trunk/ppSim/src/ppSimMakeGalaxies.c

@@ -9 +9 @@
 
     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 +35 @@
     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 +42 @@
     }
 
-    if (galaxyDensity <= 0) return NULL;
+    if (galaxyDensity <= 0) return true;
 
     // Size of FPA
@@ -66 +66 @@
         psLogMsg("ppSim", PS_LOG_INFO,
                  "Image noise is above brightest random galaxy --- no random galaxy added.");
-        return NULL;
+        return true;
     }
 
@@ -138 +138 @@
         }
     }
-    return galaxies;
+    return true;
 }

trunk/ppSim/src/ppSimMakeSky.c

@@ -1 +1 @@
 # 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 +17 @@
     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 +47 @@
     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 +70 @@
             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

trunk/ppSim/src/ppSimMakeStars.c

@@ -1 +1 @@
 # 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 +27 @@
     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 +57 @@
         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;

trunk/ppSim/src/ppSimPhotom.c

@@ -8 +8 @@
     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 +47 @@
 
             // 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);
         }
     }

trunk/ppSim/src/ppSimSetPSF.c

@@ -21 +21 @@
     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;

trunk/ppSim/src/ppSimUtils.c

@@ -15 +15 @@
     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 +95 @@
 }
 
+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 +127 @@
     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 +149 @@
     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 +163 @@
     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;
+}
+