Index: /branches/eam_branches/ipp-20110213/ppImage/src/Makefile.am
===================================================================
--- /branches/eam_branches/ipp-20110213/ppImage/src/Makefile.am	(revision 30679)
+++ /branches/eam_branches/ipp-20110213/ppImage/src/Makefile.am	(revision 30680)
@@ -2,5 +2,5 @@
 
 noinst_HEADERS = \
-	ppImage.h 
+	ppImage.h
 
 if HAVE_SVNVERSION
@@ -61,5 +61,7 @@
 	ppImageFileCheck.c \
 	ppImageVersion.c \
-	ppImageMemory.c
+	ppImageMemory.c \
+        ppImageAddNoise.c \
+        ppImageRandomGaussian.c
 
 CLEANFILES = *~
Index: /branches/eam_branches/ipp-20110213/ppImage/src/ppImage.h
===================================================================
--- /branches/eam_branches/ipp-20110213/ppImage/src/ppImage.h	(revision 30679)
+++ /branches/eam_branches/ipp-20110213/ppImage/src/ppImage.h	(revision 30680)
@@ -53,4 +53,5 @@
     bool doCrosstalkMeasure;            // measure crosstalk signal
     bool doCrosstalkCorrect;            // apply crosstalk correction
+    bool addNoise;                      // Add noise to degrade MD image to 3pi
 
     // output files requested
@@ -155,5 +156,4 @@
 bool ppImageDetrendBias(pmReadout *inputReadout, pmReadout *bias, pmReadout *dark, ppImageOptions *options);
 
-//bool ppImageDetrendNonLinear(pmReadout *input, ppImageOptions *options);
 bool ppImageDetrendNonLinear(pmReadout *input, pmFPAview *linearity, pmConfig *config);
 bool ppImageDetrendNonLinearLookup(pmReadout *input, psMetadataItem *dataItem);
@@ -291,3 +291,12 @@
 void ppImageMemoryDump(const char *description);
 
+
+//Functions needed to degrade MD exposures to 3pi exposures
+
+bool ppImageAddNoise(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmFPA *fpa) ;
+double ppImageRandomGaussian (const psRandom *rnd, double mean, double sigma);
+double ppImageRandomGaussianNorm (const psRandom *rnd);
+void ppImageRandomGaussianFree(void);
+
+
 #endif
Index: /branches/eam_branches/ipp-20110213/ppImage/src/ppImageAddNoise.c
===================================================================
--- /branches/eam_branches/ipp-20110213/ppImage/src/ppImageAddNoise.c	(revision 30680)
+++ /branches/eam_branches/ipp-20110213/ppImage/src/ppImageAddNoise.c	(revision 30680)
@@ -0,0 +1,98 @@
+# include "ppImage.h"
+
+/* This function degrades MD exposures to 3Pi exposures, by adding appropriate noise
+   and scaling the values in the image - Daniel Farrow.
+
+   EAM : Currently, this function is totally hard-wired to use values for GPC1 MD fields.
+   The sigmas should be included in a recipe, and probably calculated based on a target
+   exptime.
+
+ */
+
+bool ppImageAddNoise(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmFPA *fpa) 
+{
+
+  
+  // this step is optional.
+  if (!options->addNoise) {
+    return true;
+  }
+
+  // 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};
+  
+  // 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
+
+  // Something to choose the band, g,r,i,z,y = 0,1,2,3,4 respectively
+  char *filter = psMetadataLookupStr (NULL, fpa->concepts, "FPA.FILTERID");
+  
+  int band = 0;
+  if (!strcmp(filter, "g")) {
+    band = 1;
+  }
+  if (!strcmp(filter, "r")) {
+    band = 2;
+  }
+  if (!strcmp(filter, "i")) {
+    band = 3;
+  }
+  if (!strcmp(filter, "z")) {
+    band = 4;
+  }
+  if (!strcmp(filter, "y")) {
+    band = 5;
+  }
+  if (!band) {
+    psError(PS_ERR_UNKNOWN, true, "ppImageAddNoise doesn't recognise the filter %s, aborting", filter);
+    return false;
+  }
+
+  pmReadout *inReadout = pmFPAfileThisReadout(config->files, view, "PPIMAGE.INPUT");
+  
+  // find the currently selected readout
+  psImage *image = inReadout->image;
+  psImage *variance = inReadout->variance;
+  
+  // Warning, just in case this is left on accidently!
+  psWarning("addNoise is set; adding noise to the images.");
+
+  // Add in appropriate variance, and scale the image 
+  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));
+
+  fprintf (stderr, "mean sky: %f, scaling by %f, adding %f before re-scaling\n", MSKY_MN, rho, add_sigma);
+
+  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_sigma);
+      image->data.F32[iy][ix] /= rho;
+      variance->data.F32[iy][ix] += PS_SQR(add_sigma);
+      variance->data.F32[iy][ix] /= rho2;
+    }
+  }
+	  
+  // Update the metadata about exposure time
+  psMetadataAddF32(inReadout->parent->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "the modified exposure time", expTimes3Pi[band]);
+  psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "EXPTIME", PS_META_REPLACE, "the modified exposure time", expTimes3Pi[band]);
+  ppImageRandomGaussianFree();
+  psFree(rng);
+  psFree(stats);
+  
+  return true;
+}
Index: /branches/eam_branches/ipp-20110213/ppImage/src/ppImageLoop.c
===================================================================
--- /branches/eam_branches/ipp-20110213/ppImage/src/ppImageLoop.c	(revision 30679)
+++ /branches/eam_branches/ipp-20110213/ppImage/src/ppImageLoop.c	(revision 30680)
@@ -125,4 +125,9 @@
                     ESCAPE("Unable to measure CTE");
                 }
+
+		// optionally degrade a MD image to 3pi exposure times
+		if (!ppImageAddNoise(config, options, view, input->fpa)){
+                    ESCAPE("Unable to degrade MD image to 3pi");
+		}
             }
 
@@ -152,5 +157,5 @@
             ESCAPE("Unable to free detrend images");
         }
-    
+
         // Apply the fringe correction
         if (options->doFringe) {
Index: /branches/eam_branches/ipp-20110213/ppImage/src/ppImageOptions.c
===================================================================
--- /branches/eam_branches/ipp-20110213/ppImage/src/ppImageOptions.c	(revision 30679)
+++ /branches/eam_branches/ipp-20110213/ppImage/src/ppImageOptions.c	(revision 30680)
@@ -42,5 +42,6 @@
     options->applyParity     = false;   // Apply Cell parities
     options->doMaskStats     = false;   // Calculate mask fractions
-    
+    options->addNoise        = false;  //Degrade an MD image to a 3pi image
+
     // output files requested
     options->BaseFITS        = false;   // create output image
@@ -239,4 +240,5 @@
     options->doPatternCell = psMetadataLookupBool(NULL, recipe, "PATTERN.CELL");
     options->doMaskStats = psMetadataLookupBool(NULL, recipe, "MASK.STATS");
+    options->addNoise = psMetadataLookupBool(NULL, recipe, "ADDNOISE");
 
     options->doStats = false;
Index: /branches/eam_branches/ipp-20110213/ppImage/src/ppImageRandomGaussian.c
===================================================================
--- /branches/eam_branches/ipp-20110213/ppImage/src/ppImageRandomGaussian.c	(revision 30680)
+++ /branches/eam_branches/ipp-20110213/ppImage/src/ppImageRandomGaussian.c	(revision 30680)
@@ -0,0 +1,102 @@
+# include "ppImage.h"
+
+static int Ngaussint = 0;
+static double *gaussint = NULL;
+
+extern double drand48();
+
+double p_ppImageGaussian (double x, double mean, double sigma) {
+
+  double f;
+
+  f = exp (-0.5 * PS_SQR(x - mean) / PS_SQR(sigma)) / sqrt(2 * M_PI * PS_SQR(sigma));
+
+  return (f);
+
+}
+
+void ppImageRandomGaussianFree()
+{
+    psFree (gaussint);
+    return;
+}
+
+void ppImageRandomGaussianAlloc (int Nbin) {
+
+    gaussint = (double *) psAlloc(Nbin*sizeof(double));
+    return;
+}
+
+/* integrate a gaussian from -5 sigma to +5 sigma */
+void p_ppImageRandomGaussianInit (void) {
+
+  int i;
+  long A, B;
+  double val, x, dx, dx1, dx2, dx3, df;
+  double mean, sigma;
+
+  /* no need to generate this if it already exists */
+  if (gaussint) return;
+
+  A = time(NULL);
+  for (B = 0; A == time(NULL); B++);
+  srand48(B);
+
+  Ngaussint = 0x1000;
+  ppImageRandomGaussianAlloc (Ngaussint + 1);
+
+  val = 0;
+  dx = 1.0 / Ngaussint;
+  dx1 = dx / 3.0;
+  dx2 = 2.0*dx/3.0;
+  dx3 = dx;
+  mean = 0.0;
+  sigma = 1.0;
+
+  for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
+    df = (3.0*p_ppImageGaussian(x    , mean, sigma) +
+          9.0*p_ppImageGaussian(x+dx1, mean, sigma) +
+          9.0*p_ppImageGaussian(x+dx2, mean, sigma) +
+          3.0*p_ppImageGaussian(x+dx3, mean, sigma)) * (dx1/8.0);
+    val += df;
+    if (val > (i + 0.5) / (double) Ngaussint) {
+      gaussint[i] = x + dx / 2.0;
+      i++;
+    }
+  }
+}
+
+// XXX we are using drand48() rather than the random var supplied by rnd
+double ppImageRandomGaussian (const psRandom *rnd, double mean, double sigma) {
+
+  int i;
+  double y;
+
+  if (gaussint == NULL) {
+      p_ppImageRandomGaussianInit ();
+  }
+
+  y = drand48();
+  i = Ngaussint*y;
+  y = gaussint[i]*sigma + mean;
+
+  return (y);
+
+}
+
+// XXX we are using drand48() rather than the random var supplied by rnd
+double ppImageRandomGaussianNorm (const psRandom *rnd) {
+
+  int i;
+  double y;
+
+  if (gaussint == NULL) {
+      p_ppImageRandomGaussianInit ();
+  }
+
+  y = drand48();
+  i = Ngaussint*y;
+  y = gaussint[i];
+
+  return (y);
+}