Index: trunk/ppSim/src/ppSimMakeSky.c
===================================================================
--- trunk/ppSim/src/ppSimMakeSky.c	(revision 17557)
+++ trunk/ppSim/src/ppSimMakeSky.c	(revision 18011)
@@ -1,5 +1,7 @@
 # include "ppSim.h"
 
-// XXX add bounds to the inputs?
+// this function sets the skyRate to a value for the night sky (SKY.RATE or SKY.MAGS) or for a
+// flat-field image (FLAT.RATE).  Include a shutter correction and a scattered light source
+
 bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config) {
 
@@ -15,14 +17,22 @@
     psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe
 
-    float expTime     = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time
-    float flatSigma   = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient
-    float flatRate    = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate
-    float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time
-    float skyRate     = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate
-    if (isnan(skyRate)) {
-	float zp      = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status);
-	float scale   = psMetadataLookupF32(&status, recipe, "SCALE");     assert (status);
-	float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS");  assert (status);
+    bool sky  = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux?
+    bool flat = psMetadataLookupBool(&status, recipe, "FLAT"); // Apply flat-field term?
+ 
+    float expTime      = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time
+
+    float flatSigma    = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient
+    float flatRate     = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate
+    float shutterTime  = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time
+    float scatterFrac  = psMetadataLookupF32(&status, recipe, "SCATTER.FRAC"); // scattered light fraction (max)
+    float skyRate      = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate
+    float skyMags      = psMetadataLookupF32(&status, recipe, "SKY.MAGS");  assert (status);
+    if (!isnan(skyMags)) {
+	float zp       = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status);
+	float scale    = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status);
 	skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
+    }
+    if (type == PPSIM_TYPE_FLAT) {
+      skyRate = flatRate;
     }
 
@@ -37,12 +47,14 @@
     int yParityCell   = psMetadataLookupS32(&status, cell->concepts, "CELL.YPARITY");
 
-    int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
+    int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y
 
     // Size of FPA
     psRegion *bounds = ppSimFPABounds (fpa);
+    int dXfpa = bounds->x1 - bounds->x0;
+    int dYfpa = bounds->y1 - bounds->y0;
 
     // Correct chip offsets so that boresight is in the middle of the FPA
-    x0Chip -= 0.5 * (bounds->x1 - bounds->x0);
-    y0Chip -= 0.5 * (bounds->y1 - bounds->y0);
+    x0Chip -= 0.5 * dXfpa;
+    y0Chip -= 0.5 * dYfpa;
 
     for (int y = 0; y < signal->numRows; y++) {
@@ -58,21 +70,30 @@
             float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;
 
-            // Gaussian flat-field over the FPA
-            float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) /
-                flatSigma / sqrtf(2.0 * M_PI);
+            // Gaussian flat-field over the FPA with flatValue = 1.0 at the field center
+            float flatValue = 1.0;
+	    if (flat) {
+		// we make the flat-field have a response of 1.0 at the field center (like a vignetting)
+		flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA)));
+	    }
 
-            if (type == PPSIM_TYPE_FLAT) {
-                float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field
-                signal->data.F32[y][x] += flatFlux;
-                variance->data.F32[y][x] += flatFlux;
-                continue;
-            }
+	    float scatterRate = 0.0;
 
-            expCorr->data.F32[y][x] = realExpTime / expTime;
+	    if (sky) {
+	      // add a scattered light term to the flat-field images (no
+	      if (type == PPSIM_TYPE_FLAT) {
+		float xF = 2.0*(xFPA / dXfpa) - 1.0;
+		scatterRate = scatterFrac * PS_SQR(xF);
+	      }
 
-            // Sky background
-            float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky
-            signal->data.F32[y][x] += skyFlux;
-            variance->data.F32[y][x] += skyFlux;
+	      // Sky background
+	      float skyFlux = (skyRate * (flatValue + scatterRate)) * realExpTime; // Flux from sky
+	      signal->data.F32[y][x] += skyFlux;
+	      variance->data.F32[y][x] += skyFlux;
+	    }
+
+	    // used later to modify the star and galaxy photometry 
+	    if (expCorr) {
+	      expCorr->data.F32[y][x] = realExpTime / expTime;
+	    }
 
             // TO DO: Add fringes