Index: /tags/ipp-20101215/ppImage/src/ppImageAddNoise.c
===================================================================
--- /tags/ipp-20101215/ppImage/src/ppImageAddNoise.c	(revision 30640)
+++ /tags/ipp-20101215/ppImage/src/ppImageAddNoise.c	(revision 30641)
@@ -20,9 +20,9 @@
 
   // grizy variances to add to turn MD exposure -> 3pi, calculated from the DRM 
-  static float add_sigmas[] = {NAN, 13.39, 23.79, 75.25, 110.19, 128.84};
+  // static float add_sigmas[] = {NAN, 13.39, 23.79, 75.25, 110.19, 128.84};
   
   // Target Exposure times for 3pi in grizy
   static float expTimes3Pi[] = {NAN, 43.0, 40.0, 45.0, 30.0, 30.0}; 
-  float expTime   = psMetadataLookupF32(NULL, fpa->concepts, "FPA.EXPOSURE"); // Exposure time for image
+  float expTime = psMetadataLookupF32(NULL, fpa->concepts, "FPA.EXPOSURE"); // Exposure time for image
 
   // Something to choose the band, g,r,i,z,y = 0,1,2,3,4 respectively
@@ -60,15 +60,27 @@
 
   // Add in appropriate variance, and scale the image 
-  float expTimeScaling =  expTime/expTimes3Pi[band];
-  // XXX why was this not used? float expTimeScalingSqd = PS_SQR(expTimeScaling);
+  float rho  =  expTime/expTimes3Pi[band];
+  float rho2 = PS_SQR(rho);
 
   psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
 
+  psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+  stats->nSubsample = 10000;
+
+  psImageBackground(stats, NULL, inReadout->image, inReadout->mask, 0xffff, rng);
+  double MSKY_MN = stats->sampleMedian;
+
+  // set GAIN and RDNOISE to nominal values:
+  double GAIN = 1.0; // electrons / DN
+  double RDNOISE = 6.0; // electrons (== DN)
+
+  double add_sigma = sqrt((MSKY_MN/GAIN)*(rho - 1.0) + PS_SQR(RDNOISE)*(rho2 - 1.0));
+
   for (int iy = 0; iy < image->numRows; iy++){
     for (int ix = 0; ix < image->numCols; ix++){
-      image->data.F32[iy][ix] += ppImageRandomGaussian(rng, 0.0, add_sigmas[band]);
-      image->data.F32[iy][ix] /= expTimeScaling;
-      variance->data.F32[iy][ix] += PS_SQR(add_sigmas[band]);
-      variance->data.F32[iy][ix] /= expTimeScaling;
+      image->data.F32[iy][ix] += ppImageRandomGaussian(rng, 0.0, add_sigma);
+      image->data.F32[iy][ix] /= rho;
+      variance->data.F32[iy][ix] += PS_SQR(add_sigma);
+      variance->data.F32[iy][ix] /= rho2;
     }
   }