Index: /trunk/ppSim/notes.txt
===================================================================
--- /trunk/ppSim/notes.txt	(revision 14657)
+++ /trunk/ppSim/notes.txt	(revision 14657)
@@ -0,0 +1,47 @@
+
+ppSim development work still needed:
+
+* save sources in output files
+* load input image
+
+* optionally apply detrend data (bias, dark, etc).  when adding data
+  to a loaded image, we only need to inject fake stars.
+
+* examine the range of star fluxes.  (too many are being generated).  
+* check the overlap between the real and random star fluxes
+
+* generate stars with / without known exposure correction (ie, stars
+  could be generated with known fluxes or know magnitudes)
+
+* generate galaxies (objects with random values for major, minor, theta)
+* generate convolved galaxies
+
+* should not have to supply ZP, Seeing, PA, Scale on command line for all ppSim runs
+* DVO.CATDIR path://SIMTEST is not being interpolated?
+* -D PSASTRO:DVO.CATDIR is not being interpretted?
+
+* fake stars do not have random magnitudes?
+* add generated stars to an output file (pmSource list?)
+* make star source an option (random / dvo database)
+* add the gain to the metadata
+
+ppSimLoop:
+
+  * load catalogs stars (ra, dec, mag, x, y)
+  * generate random stars (ra, dec, mag, x, y)
+  * loop over chip,cell,readout
+  ** generate output images (signal, variance)
+  ** add the bias
+  ** add the dark
+  ** add the sky signal 
+  ** add the stars (modified by the flat, shutter)
+  ** add the galaxies (modified by the flat, shutter)
+  ** add the poisson noise
+  ** add the overscan + readnoise
+  ** save the stars
+  ** save the galaxies
+
+  ** update the concept info
+  ** update the header wcs info
+  ** write out the image data
+
Index: /trunk/ppSim/psf.txt
===================================================================
--- /trunk/ppSim/psf.txt	(revision 14657)
+++ /trunk/ppSim/psf.txt	(revision 14657)
@@ -0,0 +1,59 @@
+
+Code to load a psf model and construct psf images.  This assumes there
+is an analysis loop which is performing pmFPAfileIOChecks at the
+appropriate levels. 
+
+** to supply a psf model from the command, define the command-line arguments (ppSimArguments.c)
+
+   pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist");
+
+** if a psf model is supplied, generate a pmFPAfile to carry it. bind
+   it to the pmFPAfile in your main loop (ppSimCreate.c):
+
+    // have we supplied a psf model?
+    if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) {
+	bool status = false;
+
+	// tie the psf file to the chipMosaic 
+	// file is the pmFPAfile used for the analysis loop
+        pmFPAfileBindFromArgs(&status, file, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF");
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD");
+	    // XXX free things to be freed...
+            return NULL;
+        }
+    }
+
+** generate source models from the psf and insert into a readout:
+
+   // load the psf model from the corresponding chip.  note: you may
+   // need to select the chip from the current readout, eg: (ppSimInsertSources.c)
+   // pmCell *cell = readout->parent;
+   // pmChip *chip = cell->parent;
+   pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF");
+
+
+   // instantiate a model for the PSF at this location, set desired flux
+
+   // note that pmModelFromPSFforXY takes the chip coordinate.  you many
+   // need to convert the x,y coordinate in the readout to/from the x,y
+   // coordinate. (ppSimInsertSources.c)
+   pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
+   pmModelSetFlux (model, flux);
+
+   // define the radius of valid pixels
+   float radius = model->modelRadius (model->params, noise);
+   radius = PS_MAX (radius, 1.0);
+
+   // construct a source, with model flux pixels set, based on the model
+   pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
+
+   // insert the source flux in the image.  dX,dY is the coordinate of
+   // the readout 0,0 pixel in the chip frame.
+   pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+
+** given a pmSource, generate the its cached model flux:
+
+   // this fills in the pixels in the image source->modelFlux
+   pmSourceCacheModel (source, maskVal);
+
Index: /trunk/ppSim/src/Makefile.am
===================================================================
--- /trunk/ppSim/src/Makefile.am	(revision 14656)
+++ /trunk/ppSim/src/Makefile.am	(revision 14657)
@@ -13,5 +13,5 @@
 	ppSimMakeDark.c         \
 	ppSimMakeSky.c          \
-	ppSimInsertStars.c      \
+	ppSimInsertSources.c    \
 	ppSimAddOverscan.c      \
 	ppSimAddNoise.c         \
@@ -19,4 +19,5 @@
 	ppSimBounds.c           \
 	ppSimStars.c            \
+	ppSimSetPSF.c           \
 	ppSimUtils.c            \
 	ppSimLoop.c
Index: /trunk/ppSim/src/ppSim.c
===================================================================
--- /trunk/ppSim/src/ppSim.c	(revision 14656)
+++ /trunk/ppSim/src/ppSim.c	(revision 14657)
@@ -4,4 +4,6 @@
 {
     psLibInit(NULL);
+    if (!pmModelClassInit ()) abort();
+
     pmConfig *config = pmConfigRead(&argc, argv, NULL); // Configuration
     if (!config) {
@@ -18,9 +20,8 @@
     }
 
-    psExit loopError;
-    if ((loopError = ppSimLoop(config))) {
+    if (!ppSimLoop(config)) {
         psErrorStackPrint(stderr, "Unable to generate data.");
         psFree(config);
-        exit(loopError);
+        exit(PS_EXIT_SYS_ERROR);
     }
 
Index: /trunk/ppSim/src/ppSim.h
===================================================================
--- /trunk/ppSim/src/ppSim.h	(revision 14656)
+++ /trunk/ppSim/src/ppSim.h	(revision 14657)
@@ -72,5 +72,5 @@
 
 /// Loop over the output file, generating simulated data
-psExit ppSimLoop(pmConfig *config ///< Configuration
+bool ppSimLoop(pmConfig *config ///< Configuration
 		 );
 
@@ -89,8 +89,9 @@
 
 psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config);
-psVector *ppSimMakeBias (psImage *signal, psImage *variance, pmCell *cell, pmConfig *config, const psRandom *rng) ;
-bool ppSimMakeDark (psImage *signal, psImage *variance, pmConfig *config);
-bool ppSimMakeSky (psImage *signal, psImage *variance, psImage *expCorr, ppSimType type, pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell);
+psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) ;
+bool ppSimMakeDark (pmReadout *readout, pmConfig *config);
+bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config);
 bool ppSimInsertStars (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell);
+bool ppSimInsertSources (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config);
 
 bool ppSimInitHeader(pmConfig *config,
@@ -99,7 +100,6 @@
 		     pmCell *cell);
 
-bool ppSimSaturate(psImage *image,	// Image to apply saturation
-		   const pmConfig *config,		       
-		   const pmCell *cell); // Saturation level
+bool ppSimSaturate(pmReadout *readout,	// Image to apply saturation
+		   const pmConfig *config); // Saturation level
 
 bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config);
@@ -108,10 +108,12 @@
 bool ppSimAddOverscan (pmReadout *readout, pmConfig *config, psVector *biasCols, psVector *biasRows, psRandom *rng);
 
-psImage *ppSimAddNoise(psImage *signal, // Signal image, modified and returned
-		       const psImage *variance, // Variance image
-		       const pmConfig *config,		       
-		       const pmCell *cell,
-		       const psRandom *rng // Random number generator
-		       );
+bool ppSimAddNoise(psImage *signal,
+		   psImage *variance,
+		   const pmCell *cell,		       
+		   const pmConfig *config,		       
+		   const psRandom *rng // Random number generator
+    );
+
+bool ppSimSetPSF (pmChip *chip, pmConfig *config);
 
 #endif
Index: /trunk/ppSim/src/ppSimAddNoise.c
===================================================================
--- /trunk/ppSim/src/ppSimAddNoise.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimAddNoise.c	(revision 14657)
@@ -2,9 +2,9 @@
 
 // Add noise to an image
-psImage *ppSimAddNoise(psImage *signal, // Signal image, modified and returned
-		       const psImage *variance, // Variance image
-		       const pmConfig *config,		       
-		       const pmCell *cell,
-		       const psRandom *rng // Random number generator
+bool ppSimAddNoise(psImage *signal, // Signal image, modified and returned
+		   psImage *variance,
+		   const pmCell *cell,
+		   const pmConfig *config,		       
+		   const psRandom *rng // Random number generator
     )
 {
@@ -22,5 +22,5 @@
 	if (!mdok) {
 	    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find GAIN in recipe.");
-	    return NULL;
+	    return false;
 	}
     }
@@ -38,5 +38,5 @@
 
     // XXX why return this??
-    return signal;
+    return true;
 }
 
Index: /trunk/ppSim/src/ppSimAddOverscan.c
===================================================================
--- /trunk/ppSim/src/ppSimAddOverscan.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimAddOverscan.c	(revision 14657)
@@ -29,5 +29,5 @@
 	psImageInit(variance, PS_SQR(readnoise));
 
-	ppSimAddNoise(signal, variance, config, readout->parent, rng);
+	ppSimAddNoise(signal, variance, readout->parent, config, rng);
 	psListAdd(readout->bias, PS_LIST_TAIL, signal);
 
Index: /trunk/ppSim/src/ppSimArguments.c
===================================================================
--- /trunk/ppSim/src/ppSimArguments.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimArguments.c	(revision 14657)
@@ -66,4 +66,5 @@
     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);
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1);
 
@@ -194,4 +195,8 @@
         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",
                          pa * M_PI / 180.0);
+
+
+	const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // Filter name
+	psMetadataAddStr(config->arguments, PS_LIST_TAIL, "PSF.MODEL", 0, "PSF model class", psfClass);
     }
 
Index: /trunk/ppSim/src/ppSimCreate.c
===================================================================
--- /trunk/ppSim/src/ppSimCreate.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimCreate.c	(revision 14657)
@@ -1,3 +1,9 @@
 # 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)
@@ -12,5 +18,4 @@
 
     pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE);
-    psFree(fpa);
     if (!file) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s.  "
@@ -20,4 +25,5 @@
     if (file->type != PM_FPA_FILE_IMAGE) {
         psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE);
+	psFree(fpa);
         psFree(file);
         return NULL;
@@ -25,13 +31,29 @@
     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");
-            return status;
+	    psFree(fpa);
+	    psFree(file);
+            return NULL;
         }
     }
+
+    // XXX only invoke this code for OBJECT types of images
+    // PPSIM.SOURCES carries the constructed, fake sources with their true parameters
+    pmFPAfile *simSources = pmFPAfileDefineOutput (config, file->fpa, "PPSIM.SOURCES");
+    if (!simSources) {
+        psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES");
+        return false;
+    }
+    simSources->save = true;
 
     pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes
@@ -44,4 +66,5 @@
             psFree(fpa);
             psFree(view);
+	    psFree(file);
             return NULL;
         }
@@ -55,4 +78,5 @@
                 psFree(fpa);
                 psFree(view);
+		psFree(file);
                 return NULL;
             }
@@ -66,4 +90,5 @@
                     psFree(fpa);
                     psFree(view);
+		    psFree(file);
                     return NULL;
                 }
@@ -71,4 +96,6 @@
         }
     }
+
+    psFree(fpa);
     psFree(view);
 
Index: /trunk/ppSim/src/ppSimInsertSource.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertSource.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimInsertSource.c	(revision 14657)
@@ -1,5 +1,5 @@
 # include "ppSim.h"
 
-bool ppSimInsertSources (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) {
+bool ppSimInsertSource (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) {
 
     bool mdok;
Index: /trunk/ppSim/src/ppSimInsertSources.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertSources.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimInsertSources.c	(revision 14657)
@@ -1,13 +1,22 @@
 # include "ppSim.h"
 
-bool ppSimInsertSources (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) {
+bool ppSimInsertSources (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {
 
     bool mdok;
 
     assert (stars);
+    assert (readout);
     
+    // XXX is this needed?
+    // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+
+    // XXX this is an estimate of the sky noise based on the inputs to the image simulation.
+    // XXX update this to allow the estimate based on the measured sky background
+    // XXX this is missing the gain.
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
 
-    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
@@ -23,4 +32,7 @@
     }
 
+    // Rough noise estimate, appropriate for entire cell (use for source radius?)
+    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+
     int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0");
     int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0");
@@ -35,25 +47,56 @@
     int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
 
-    // Rough noise estimate, appropriate for entire cell
-    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+    pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF");
+    assert (psf);
 
-    pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSF");
+    int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning);
+    int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning);
 
+    psArray *sources = psArrayAllocEmpty (stars->n);
+
+    // add sources to the readout image & weight
     for (long i = 0; i < stars->n; i++) {
 	ppSimStar *star = stars->data[i];
 
+	// star->x,y are in fpa coordinates
+
 	// Position on the cell and peak flux
-	float x = PPSIM_FPA_TO_CELL(star->x, x0Cell, xParityCell, binning, x0Chip, xParityChip);
-	float y = PPSIM_FPA_TO_CELL(star->y, y0Cell, yParityCell, binning, y0Chip, yParityChip);
+	float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
+	float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
 
-	// instantiate a model for the PSF at this location (normalized or not?)
-	// XXX 
-	pmModel *model = pmModelFromPSFforXY (psf, star->x, star->y, star->flux);
+	// Position on the cell and peak flux
+	float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
+	float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
+
+	if (xCell < 0) continue;
+	if (yCell < 0) continue;
+	if (xCell > readout->image->numCols) continue;
+	if (yCell > readout->image->numRows) continue;
+	// XXX need to apply col0, row0 if readout is a subarray
+
+	// XXX apply the expCorr to the star->flux before setting the model flux
+
+	// instantiate a model for the PSF at this location, set desired flux
+	pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
+	pmModelSetFlux (model, star->flux);
+
+	// XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
+	float radius = model->modelRadius (model->params, roughNoise);
+	radius = PS_MAX (radius, 1.0);
 
 	// construct a source, with model flux pixels set, based on the model
-	pmSource *source = pmSourceFromModel (model);
+	pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
 
-	ppSimInsertSource(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr);
+	// XXX add the sources to a source array
+
+	// insert the source flux in the image
+	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+	psArrayAdd (sources, 100,source);
     }
+
+    // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT"
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
+
+    // XXX many leaks in here, i think 
     return true;
 }
Index: /trunk/ppSim/src/ppSimLoop.c
===================================================================
--- /trunk/ppSim/src/ppSimLoop.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimLoop.c	(revision 14657)
@@ -1,5 +1,5 @@
 #include "ppSim.h"
 
-psExit ppSimLoop(pmConfig *config)
+bool ppSimLoop(pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(config, PS_EXIT_PROG_ERROR);
@@ -29,7 +29,10 @@
     pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA
 
+    // XXX if we include the psphot analysis, we will need to activate the correct pmFPAfiles
+    // at the correct times.  ppSim operates on PPSIM.OUTPUT and PPSIM.SOURCES
+
     // load any needed files (eg, input image, PSF)
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-	psError(PSPHOT_ERR_DATA, false, "failed IO for fpa in ppSim\n");
+	psError(PS_ERR_UNKNOWN, false, "failed IO for fpa in ppSim\n");
 	psFree(view);
 	return false;
@@ -43,12 +46,15 @@
 	// load any needed files (eg, input image, PSF)
 	if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-            psError(PSPHOT_ERR_DATA, false, "failed IO for chip %d in ppSim\n", view->chip);
+            psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in ppSim\n", view->chip);
 	    psFree (view);
 	    return false;
 	}
 
-	pmPSF *psf = NULL;
 	if (type == PPSIM_TYPE_OBJECT) {
-	    psf = ppSimSetPSF (config);
+	    if (!ppSimSetPSF (chip, config)) {
+		psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in ppSim\n", view->chip);
+		psFree (view);
+		return false;
+	    }
 	}
 
@@ -64,5 +70,5 @@
                 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "CELL.READDIR = 1 is the only supported mode.");
                 psFree(rng);
-                return PS_EXIT_CONFIG_ERROR;
+                return false;
             }
 
@@ -74,6 +80,6 @@
 
 		// TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately
-                psImage *signal = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels
-                psImage *variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
+                readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels
+                readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
 
                 psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects
@@ -82,24 +88,20 @@
                 }
 
-		psVector *biasRows = ppSimMakeBias (signal, variance, cell, config, rng);
+		psVector *biasRows = ppSimMakeBias (readout, config, rng);
 		if (type == PPSIM_TYPE_BIAS) goto done;
 		
-		ppSimMakeDark (signal, variance, config);
+		ppSimMakeDark (readout, config);
 		if (type == PPSIM_TYPE_DARK) goto done;
 
-		ppSimMakeSky (signal, variance, expCorr, type, config, fpa, chip, cell);
+		ppSimMakeSky (readout, expCorr, type, config);
 		if (type == PPSIM_TYPE_FLAT) goto done;
 
 		if (type == PPSIM_TYPE_OBJECT) {
-		    ppSimInsertStars (signal, variance, expCorr, stars, config, chip, cell);
+		    ppSimInsertSources (readout, expCorr, stars, config);
 		}
 
 	    done:
-		ppSimAddNoise(signal, variance, config, cell, rng);
-                ppSimSaturate(signal, config, cell);
-
-                readout->image = signal;
-                psFree(variance);
-                psFree(expCorr);
+		ppSimAddNoise(readout->image, readout->weight, cell, config, rng);
+                ppSimSaturate(readout, config);
 
 		ppSimAddOverscan (readout, config, biasCols, biasRows, rng);
@@ -123,5 +125,6 @@
                 psFree(rng);
                 psFree(view);
-                return PS_EXIT_SYS_ERROR;
+                // return PS_EXIT_SYS_ERROR;
+                return false;
             }
         }
@@ -135,5 +138,6 @@
             psFree(rng);
             psFree(view);
-            return PS_EXIT_SYS_ERROR;
+            // return PS_EXIT_SYS_ERROR;
+            return false;
         }
 
@@ -148,5 +152,6 @@
         psFree(rng);
         psFree(view);
-        return PS_EXIT_SYS_ERROR;
+        // return PS_EXIT_SYS_ERROR;
+        return false;
     }
 
@@ -154,4 +159,4 @@
     psFree(view);
 
-    return 0;
+    return true;
 }
Index: /trunk/ppSim/src/ppSimMakeBias.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeBias.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimMakeBias.c	(revision 14657)
@@ -1,7 +1,9 @@
 # include "ppSim.h"
 
-psVector *ppSimMakeBias (psImage *signal, psImage *variance, pmCell *cell, pmConfig *config, const psRandom *rng) {
+psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) {
 
     bool mdok;
+
+    pmCell *cell = readout->parent;
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
@@ -20,4 +22,7 @@
 	}
     }
+
+    psImage *signal = readout->image;
+    psImage *variance = readout->weight;
 
     int numRows = signal->numRows;
Index: /trunk/ppSim/src/ppSimMakeDark.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeDark.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimMakeDark.c	(revision 14657)
@@ -2,5 +2,8 @@
 
 // XXX add bounds to the inputs?
-bool ppSimMakeDark (psImage *signal, psImage *variance, pmConfig *config) {
+bool ppSimMakeDark (pmReadout *readout, pmConfig *config) {
+
+    psImage *signal = readout->image;
+    psImage *variance = readout->weight;
 
     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
Index: /trunk/ppSim/src/ppSimMakeSky.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeSky.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimMakeSky.c	(revision 14657)
@@ -2,5 +2,12 @@
 
 // XXX add bounds to the inputs?
-bool ppSimMakeSky (psImage *signal, psImage *variance, psImage *expCorr, ppSimType type, pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell) {
+bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config) {
+
+    psImage *signal = readout->image;
+    psImage *variance = readout->weight;
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent; 
+    pmFPA  *fpa  = chip->parent;
 
     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
Index: /trunk/ppSim/src/ppSimMakeSources.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeSources.c	(revision 14657)
+++ /trunk/ppSim/src/ppSimMakeSources.c	(revision 14657)
@@ -0,0 +1,75 @@
+# include "ppSim.h"
+
+bool ppSimMakeSources(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
+
+    assert (stars);
+
+    bool mdok;
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+
+    float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
+    float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude
+    float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes
+
+    float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
+    float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
+    float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
+
+    float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
+    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
+    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
+
+    if (starsDensity <= 0) return true;
+
+    // Size of FPA
+    psRegion *bounds = ppSimFPABounds (fpa);
+    // Size of focal plane (can't I get this from the config?)
+    int xSize = bounds->x1 - bounds->x0;
+    int ySize = bounds->y1 - bounds->y0;
+    psFree(bounds);
+
+    // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it
+    // shouldn't be too bad.
+    float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
+    if (!mdok) {
+	psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
+	psFree(bounds);
+	return false;
+    }
+
+    // Peak fluxes: faintest and brightest levels for random stars
+    float faint = 0.1 * 10.0 * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+    float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
+    if (bright < faint) {
+	psLogMsg("ppSim", PS_LOG_INFO,
+		 "Image noise is above brightest random star --- no random stars added.");
+    } else {
+
+	// Normalisation, set by the specified stellar density at the specified bright magnitude
+	float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
+	    powf(bright, starsLum);
+
+	// Total number of stars down to the faint flux end
+	long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
+
+	psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
+		 -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
+
+	long oldSize = stars->n;
+	psArrayRealloc (stars, stars->n + num);
+
+	for (long i = 0; i < num; i++) {
+	    ppSimStar *star = ppSimStarAlloc ();
+
+	    // make fpa center of distribution
+	    star->x    = psRandomUniform(rng) * xSize; // x position
+	    star->y    = psRandomUniform(rng) * ySize; // y position
+	    star->peak = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
+	    star->flux = star->peak * sqrt(2.0*M_PI)*seeing;
+	    stars->data[oldSize + i] = star;
+	}
+	stars->n = stars->nalloc;
+    }
+    return true;
+}
Index: /trunk/ppSim/src/ppSimSaturate.c
===================================================================
--- /trunk/ppSim/src/ppSimSaturate.c	(revision 14656)
+++ /trunk/ppSim/src/ppSimSaturate.c	(revision 14657)
@@ -2,10 +2,12 @@
 
 // Apply saturation limit to image
-bool ppSimSaturate(psImage *image, // Image to apply saturation
-		   const pmConfig *config,		       
-		   const pmCell *cell // Saturation level
+bool ppSimSaturate(pmReadout *readout, // Image to apply saturation
+		   const pmConfig *config // configuration data
     )
 {
     bool mdok;
+
+    psImage *image = readout->image;
+    pmCell *cell = readout->parent;
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
Index: /trunk/ppSim/src/ppSimSetPSF.c
===================================================================
--- /trunk/ppSim/src/ppSimSetPSF.c	(revision 14657)
+++ /trunk/ppSim/src/ppSimSetPSF.c	(revision 14657)
@@ -0,0 +1,78 @@
+# include "ppSim.h"
+static char *defaultModel = "PS_MODEL_GAUSS";
+
+bool ppSimSetPSF (pmChip *chip, pmConfig *config) {
+
+    bool status;
+    pmPSF *psf = NULL;
+    psPolynomial2D *param = NULL;
+
+    // the pmPSF IO functions stores the PSF on the chip->analysis
+    psf = psMetadataLookupPtr (&status, chip->analysis, "PSPHOT.PSF");
+    if (psf) {
+	return true;
+    }
+
+    // no supplied PSF, build one using supplied value for seeing
+    // XXX need to correct for the pixel scale
+    float seeing   = psMetadataLookupF32(&status, config->arguments, "SEEING"); // Seeing sigma (pix)
+
+    char *psfModelName = psMetadataLookupStr(&status, config->arguments, "PSF.MODEL"); // Seeing sigma (pix)
+    if (psfModelName == NULL) {
+	psfModelName = defaultModel;
+    }
+
+    pmModelType type = pmModelClassGetType (psfModelName);
+    if (type == -1) {
+	psError (PS_ERR_UNKNOWN, false, "invalid model name");
+        return false;
+    }
+
+    // generate a psf with no spatial variation
+    psf = pmPSFAlloc (type, true, NULL);
+    for (int i = 0; i < psf->params->n; i++) {
+	
+            if (i == PM_PAR_SKY)  continue;
+            if (i == PM_PAR_I0)   continue;
+            if (i == PM_PAR_XPOS) continue;
+            if (i == PM_PAR_YPOS) continue;
+
+	    // generate a 'constant' polynomial, unmasked
+            param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0);
+	    param->mask[0][0] = 0;
+
+            psf->params->data[i] = param;
+    }
+
+    psEllipseAxes axes;
+    psEllipsePol pol;
+    
+    // supply the semi-major axis
+    axes.major = seeing / 2.0;
+    axes.minor = seeing / 2.0;
+    axes.theta = 0.0;
+
+    pol = psEllipseAxesToPol (axes);
+
+    param = psf->params->data[PM_PAR_E0];
+    param->coeff[0][0] = pol.e0;
+
+    param = psf->params->data[PM_PAR_E1];
+    param->coeff[0][0] = pol.e1;
+
+    param = psf->params->data[PM_PAR_E2];
+    param->coeff[0][0] = pol.e2;
+
+    if (!strcasecmp (psfModelName, "PS_MODEL_QGAUSS")) {
+	param = psf->params->data[PM_PAR_7];
+	param->coeff[0][0] = 1.0;
+    }
+
+    if (!strcasecmp (psfModelName, "PS_MODEL_RGAUSS")) {
+	param = psf->params->data[PM_PAR_7];
+	param->coeff[0][0] = 1.0;
+    }
+
+    psMetadataAdd (chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,  "psphot psf", psf);
+    return true;
+}