Index: /trunk/ppSim/src/ppSimInsertStars.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertStars.c	(revision 14815)
+++ /trunk/ppSim/src/ppSimInsertStars.c	(revision 14816)
@@ -1,3 +1,9 @@
 # include "ppSim.h"
+
+// Reset a pointer: free and set to NULL
+#define RESET(PTR) \
+    psFree(PTR); \
+    PTR = NULL;
+
 
 bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {
@@ -26,10 +32,10 @@
     float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
     if (isnan(readnoise)) {
-	psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
-	readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
-	if (!mdok) {
-	    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
-	    return false;
-	}
+        psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
+        readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
+        if (!mdok) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
+            return false;
+        }
     }
 
@@ -58,55 +64,71 @@
 
     // add sources to the readout image & weight
+    psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n);
     for (long i = 0; i < stars->n; i++) {
-	ppSimStar *star = stars->data[i];
+        ppSimStar *star = stars->data[i];
+        psTrace("ppSim", 10, "Inserting star at %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux);
 
-	// star->x,y are in fpa coordinates
+        // star->x,y are in fpa coordinates
 
-	// Position on the cell and peak flux
-	float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
-	float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
+        // Position on the cell and peak flux
+        float 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);
+        // 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 Note, the below does not put the edges of stars on the readout if they fall slightly off
+        // This will be visible as cut-off stars at amplifier boundaries (e.g., Megacam)
+        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
+        // Apply the expCorr to the star->flux before setting the model flux
+        float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell];
 
-	// instantiate a model for the PSF at this location, set desired flux
-	pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
-	pmModelSetFlux (model, star->flux);
+        // instantiate a model for the PSF at this location, set desired flux
+        pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
+        pmModelSetFlux (model, 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);
+        // 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);
+        // 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;
+        // 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
+        source->psfMag = -2.5*log10(flux);
+        source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
 
-	// insert the source flux in the image
-	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
-	psArrayAdd (sources, 100,source);
+        // 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
+
+        // Blow away the image parts of the source, which makes the memory explode
+        RESET(source->pixels);
+        RESET(source->weight);
+        RESET(source->maskObj);
+        RESET(source->maskView);
+        RESET(source->modelFlux);
+        RESET(source->psfFlux);
+        RESET(source->blends);
     }
 
     // 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);
+    psFree(sources);
 
-    // XXX many leaks in here, i think 
+    // XXX many leaks in here, i think
     return true;
 }
Index: /trunk/ppSim/src/ppSimMakeSky.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeSky.c	(revision 14815)
+++ /trunk/ppSim/src/ppSimMakeSky.c	(revision 14816)
@@ -8,5 +8,5 @@
 
     pmCell *cell = readout->parent;
-    pmChip *chip = cell->parent; 
+    pmChip *chip = cell->parent;
     pmFPA  *fpa  = chip->parent;
 
@@ -38,35 +38,35 @@
     for (int y = 0; y < signal->numRows; y++) {
 
-	float yFPA = PPSIM_CELL_TO_FPA(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
-	    (bounds->y1 - bounds->y0); // Relative y position in FPA
+        float yFPA = PPSIM_CELL_TO_FPA(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
+            (bounds->y1 - bounds->y0); // Relative y position in FPA
 
-	for (int x = 0; x < signal->numCols; x++) {
-	    float xFPA = PPSIM_CELL_TO_FPA(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 /
-		(bounds->x1 - bounds->x0); // Relative x position in FPA
+        for (int x = 0; x < signal->numCols; x++) {
+            float xFPA = PPSIM_CELL_TO_FPA(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 /
+                (bounds->x1 - bounds->x0); // Relative x position in FPA
 
-	    // Shutter: adjust exposure time
-	    float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;
+            // Shutter: adjust exposure time
+            float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;
 
-	    // Gaussian flat-field over the FPA
-	    float flatValue = exp(-0.5 / PS_SQR(flatSigma) *
-				  (PS_SQR(yFPA) + PS_SQR(xFPA))) / flatSigma / sqrt(M_PI);
+            // 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);
 
-	    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;
-	    }
+            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;
+            }
 
-	    expCorr->data.F32[y][x] = flatValue * realExpTime / expTime;
+            expCorr->data.F32[y][x] = realExpTime / expTime;
 
-	    // 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 * realExpTime; // Flux from sky
+            signal->data.F32[y][x] += skyFlux;
+            variance->data.F32[y][x] += skyFlux;
 
-	    // TO DO: Add fringes
+            // TO DO: Add fringes
 
-	}
+        }
     }
     psFree(bounds);
Index: /trunk/ppSim/src/ppSimMakeStars.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeStars.c	(revision 14815)
+++ /trunk/ppSim/src/ppSimMakeStars.c	(revision 14816)
@@ -1,3 +1,6 @@
 # include "ppSim.h"
+
+#define FAINT_FUDGE_FACTOR   1.0        // Fraction of the noise in a (boxcar) PSF for the faint limit
+
 
 bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
@@ -36,27 +39,30 @@
     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;
+        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;
+    // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars
+    float faint = FAINT_FUDGE_FACTOR * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime) *
+        seeing; // Faint limit is related to the noise in a (boxcar) PSF
+    float bright = powf(10.0, -0.4 * (starsMag - zp)) * expTime; // Bright limit is specified by user as mag
+    psTrace("ppSim", 6, "Faint limit: %f\n", faint);
+    psTrace("ppSim", 6, "Bright limit: %f\n", bright);
     if (bright < faint) {
-	psLogMsg("ppSim", PS_LOG_INFO,
-		 "Image noise is above brightest random star --- no random stars added.");
-	return true;
+        psLogMsg("ppSim", PS_LOG_INFO,
+                 "Image noise is above brightest random star --- no random stars added.");
+        return true;
     }
 
     // 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);
+        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);
+    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n",
+             num, -2.5 * log10(faint / expTime) + zp);
 
     long oldSize = stars->n;
@@ -64,16 +70,16 @@
 
     for (long i = 0; i < num; i++) {
-	ppSimStar *star = ppSimStarAlloc ();
+        ppSimStar *star = ppSimStarAlloc ();
 
-	// make fpa center of distribution
-	star->x    = psRandomUniform(rng) * xSize; // x position
-	star->y    = psRandomUniform(rng) * ySize; // y position
+        // make fpa center of distribution
+        star->x    = psRandomUniform(rng) * xSize; // x position
+        star->y    = psRandomUniform(rng) * ySize; // y position
 
-	star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
-	star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
+        star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
+        star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
 
-	stars->data[oldSize + i] = star;
+        stars->data[oldSize + i] = star;
     }
-    stars->n = stars->nalloc;
+    stars->n = oldSize + num;
 
     return true;