Index: trunk/ppSim/src/ppSimInsertStars.c
===================================================================
--- trunk/ppSim/src/ppSimInsertStars.c	(revision 14463)
+++ trunk/ppSim/src/ppSimInsertStars.c	(revision 14668)
@@ -1,13 +1,24 @@
 # include "ppSim.h"
 
-bool ppSimInsertStars (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) {
+bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {
 
     bool mdok;
 
+    assert (readout);
     assert (stars);
-    
+
+    if (!stars->n) { return true; }
+
+    // XXX is this needed?
+    // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+
+    // XXX this is an estimate of the sky noise based on the inputs to the image simulation.
+    // XXX update this to allow the estimate based on the measured sky background
+    // XXX this is missing the gain.
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
 
-    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
@@ -23,4 +34,7 @@
     }
 
+    // Rough noise estimate, appropriate for entire cell (use for source radius?)
+    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+
     int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0");
     int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0");
@@ -35,15 +49,65 @@
     int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
 
-    // Rough noise estimate, appropriate for entire cell
-    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+    pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF");
+    assert (psf);
 
+    int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning);
+    int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning);
+
+    psArray *sources = psArrayAllocEmpty (stars->n);
+
+    // add sources to the readout image & weight
     for (long i = 0; i < stars->n; i++) {
 	ppSimStar *star = stars->data[i];
 
+	// star->x,y are in fpa coordinates
+
 	// Position on the cell and peak flux
-	float x = PPSIM_FPA_TO_CELL(star->x, x0Cell, xParityCell, binning, x0Chip, xParityChip);
-	float y = PPSIM_FPA_TO_CELL(star->y, y0Cell, yParityCell, binning, y0Chip, yParityChip);
-	ppSimInsertStar(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr);
+	float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
+	float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
+
+	// Position on the cell and peak flux
+	float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
+	float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
+
+	if (xCell < 0) continue;
+	if (yCell < 0) continue;
+	if (xCell > readout->image->numCols) continue;
+	if (yCell > readout->image->numRows) continue;
+	// XXX need to apply col0, row0 if readout is a subarray
+
+	// XXX apply the expCorr to the star->flux before setting the model flux
+
+	// instantiate a model for the PSF at this location, set desired flux
+	pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
+	pmModelSetFlux (model, star->flux);
+
+	// XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
+	float radius = model->modelRadius (model->params, roughNoise);
+	radius = PS_MAX (radius, 1.0);
+
+	// construct a source, with model flux pixels set, based on the model
+	pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
+
+	// XXX set the mag & err values (should this be done in pmSourceFromModel?)
+	// XXX i should be applying the gain and the correct effective area
+	psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+
+	source->psfMag = -2.5*log10(star->flux);
+	source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux;	
+	
+	// XXX add the sources to a source array
+
+	// insert the source flux in the image
+	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+	psArrayAdd (sources, 100,source);
     }
+
+    // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT"
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
+
+    // XXX many leaks in here, i think 
     return true;
 }
+