Index: trunk/psModules/src/camera/pmFPABin.c
===================================================================
--- trunk/psModules/src/camera/pmFPABin.c	(revision 34750)
+++ trunk/psModules/src/camera/pmFPABin.c	(revision 34800)
@@ -43,4 +43,8 @@
     }
 
+    int Nbits = (int) (ceil(log(maskVal)/log(2)) + 1);
+    int *bitcounter = malloc(sizeof(int) * Nbits);
+    int pxlcount;
+
     int xLast = numColsIn - 1, yLast = numRowsIn - 1; // Last index
     int yStart = psImageBinningGetFineY(binning, 0); // Starting input y for binning
@@ -55,6 +59,21 @@
             float sum = 0.0;            // Sum of pixels
             int numPix = 0;             // Number of pixels
+
+	    for (int j = 0; j < Nbits; j++) { // Reset bit counter
+	      bitcounter[j] = 0;
+	    }
+	    pxlcount = 0;
+	    
             for (int y = yStart; y < yStop; y++) {
                 for (int x = xStart; x < xStop; x++) {
+		  if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] != 0)) {
+		      for (int j = 0; j < Nbits; j++) {
+			psImageMaskType M = (psImageMaskType) pow(2,j);
+			if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & M) {
+			  bitcounter[j]++;
+			}
+		      }
+		    }
+		  
                     if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
                         continue;
@@ -65,7 +84,10 @@
                     sum += inImage->data.F32[y][x];
                     numPix++;
+
+
                 }
             }
-
+	    
+	    
 	    // Values to set
             float imageValue;
@@ -79,5 +101,11 @@
             }
             outImage->data.F32[yOut][xOut] = imageValue;
-            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
+	    //            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
+	    outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = 0;
+	    for (int j = 0; j < Nbits; j++) {
+	      if (bitcounter[j] > 0.5 * pxlcount) {
+		outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] |= (int) pow(2,j);
+	      }
+	    }
             xStart = xStop;
         }
Index: trunk/psModules/src/camera/pmFPAfileIO.c
===================================================================
--- trunk/psModules/src/camera/pmFPAfileIO.c	(revision 34750)
+++ trunk/psModules/src/camera/pmFPAfileIO.c	(revision 34800)
@@ -76,5 +76,4 @@
     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
         pmFPAfile *file = item->data.V;
-
         switch (place) {
           case PM_FPA_BEFORE:
Index: trunk/psModules/src/detrend/pmDark.c
===================================================================
--- trunk/psModules/src/detrend/pmDark.c	(revision 34750)
+++ trunk/psModules/src/detrend/pmDark.c	(revision 34800)
@@ -353,12 +353,35 @@
 
     // retrieve the required parameter vectors
-    psArray *values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
-    psAssert(values, "values not supplied");
+    psArray *in_values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
+    psAssert(in_values, "values not supplied");
     psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK");
     psAssert(roMask, "roMask not supplied");
-    psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
-    psAssert(orders, "orders not supplied");
-
-    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting
+    psVector *max_orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
+    psAssert(max_orders, "orders not supplied");
+
+    psArray *values_set = psArrayAlloc(max_orders->n);
+    psArray *poly_set = psArrayAlloc(max_orders->n);
+    psVector *logL = psVectorAlloc(max_orders->n,PS_TYPE_F64);
+
+    for (int i = 0; i < max_orders->n; i++) {
+      psVector *orders = psVectorAlloc(i+1,PS_TYPE_U8);
+      for (int j = 0; j < orders->n; j++) {
+	orders->data.U8[j] = max_orders->data.U8[j];
+      }
+      poly_set->data[i] =  psPolynomialMDAlloc(orders); // Polynomial for fitting
+      
+      psArray *values = psArrayAlloc(in_values->n);
+      
+      for (int j = 0; j < values->n; j++) {
+	psVector *these_values = psVectorAlloc(i+1,PS_TYPE_F32);
+	psVector *input_values = in_values->data[j];
+
+	for (int k = 0; k < orders->n; k++) {
+	  these_values->data.F32[k] = input_values->data.F32[k];
+	}
+	values->data[j] = these_values;
+      }
+      values_set->data[i] = values;
+    }
 
     // retrieve the norm vector, if supplied
@@ -383,5 +406,5 @@
     }
 
-    pmDarkVisualInit(values);
+    pmDarkVisualInit(values_set->data[max_orders->n - 1]);
 
     pmReadout *outReadout = output->readouts->data[0];
@@ -423,15 +446,53 @@
             }
 
-            if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
+	    int k_best = 0;
+	    for (int k = 0; k < max_orders->n; k++) {
+	      psPolynomialMD *poly = poly_set->data[k];
+	      psArray *values = values_set->data[k];
+	      
+	      if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
                 psErrorClear();         // Nothing we can do about it
                 psVectorInit(poly->coeff, NAN);
-            }
-
-            pmDarkVisualPixelFit(pixels, mask);
-            pmDarkVisualPixelModel(poly, values);
-
-            for (int k = 0; k < poly->coeff->n; k++) {
+	      }
+
+	      pmDarkVisualPixelFit(pixels, mask);
+	      pmDarkVisualPixelModel(poly, values);
+
+	      // Insert math here to choose optimum model.
+	      logL->data.F64[k] = 0.0;
+	      psPolynomialMD *polySig = poly_set->data[0];
+	      for (int m = 0; m < poly->deviations->n; m++) {
+		logL->data.F64[k] += pow(poly->deviations->data.F32[m] / polySig->stdevFit,2);
+/* 		if ((xOut == 20) && (yOut == 256)) { */
+/* 		  psTrace("psModules.detrend",3,"pmDarkCombine DEV: %d %d: input %d models: Norders: %d logL: %g value: %g\n", */
+/* 			  xOut,yOut,m,k,logL->data.F64[k],poly->deviations->data.F32[m]); */
+/* 		} */
+	      }
+	      if (k > 0) {
+		if ( ( logL->data.F64[k - 1] - logL->data.F64[k] ) > 1) { // Hard coded criterion for a ~5% limit with one degree of freedom
+		  k_best = k;
+		}
+	      }
+	      if ((xOut <= 600) && (yOut <= 600)) {
+		psTrace("psModules.detrend",3,"pmDarkCombine: %d %d: models: Norders: %d logL: %g BestOrders: %d\n",
+			xOut,yOut,k,logL->data.F64[k],k_best);
+	      }
+	    }
+	    if (k_best > 1) {
+	      k_best = 1;
+	    }
+/* 	    k_best = 1; */
+	    // Select the polynomial that seems best.
+	    psPolynomialMD *poly = poly_set->data[k_best];
+	      
+	    //            for (int k = 0; k < poly->coeff->n; k++) {
+	    for (int k = 0; k < max_orders->n + 1; k++) { // There is one more coefficient than is stored here.
                 pmReadout *ro = output->readouts->data[k]; // Readout of interest
-                ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+		if (k < poly->coeff->n) {
+		  ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+		}
+		else {
+		  ro->image->data.F32[yOut][xOut] = 0.0;
+		}
             }
             counts->data.U16[yOut][xOut] = poly->numFit;
Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 34750)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 34800)
@@ -1358,4 +1358,57 @@
 }
 
+bool pmStackSimpleMedianCombine(
+				pmReadout *combined,
+				psArray *input) {
+  int num = input->n;
+  //  int numCols, numRows;
+  int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+  int xSize, ySize;                   // Size of the output image
+
+  psArray *stack = psArrayAlloc(num); // Stack of readouts  
+  for (int i = 0; i < num; i++) {
+    //    pmStackData *data = input->data[i]; // Stack data for this input
+    pmReadout *ro = input->data[i]; // data->readout;  // Readout of interest
+    if (!ro) {
+      continue;
+    }
+    stack->data[i] = psMemIncrRefCounter(ro);
+  }    
+
+  if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+			      stack)) {
+    psError(psErrorCodeLast(), false, "Input stack is not valid.");
+    psFree(stack);
+    return false;
+  }
+
+  psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
+  psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+
+  for (int y = minInputRows; y < maxInputRows; y++) {
+    for (int x = minInputCols; x < maxInputCols; x++) {
+      for (int i = 0; i < input->n; i++) {
+	pmReadout *ro  = stack->data[i];
+	psImage *image = ro->image;
+	pixelData->data.F32[i] = image->data.F32[y][x];
+      }
+      if (!psVectorStats(stats,pixelData,NULL,NULL,0)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
+	psFree(stats);
+	psFree(pixelData);
+	psFree(stack);
+	return(false);
+      }
+      combined->image->data.F32[y][x] = stats->robustMedian;
+      combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+    }
+  }
+  
+  psFree(stats);
+  psFree(pixelData);
+  psFree(stack);
+  return (true);
+}
+
 /// Stack input images
 bool pmStackCombine(
Index: trunk/psModules/src/imcombine/pmStack.h
===================================================================
--- trunk/psModules/src/imcombine/pmStack.h	(revision 34750)
+++ trunk/psModules/src/imcombine/pmStack.h	(revision 34800)
@@ -41,4 +41,8 @@
                               float addVariance ///< Additional variance when rejecting
     );
+/// Stack input images simply
+bool pmStackSimpleMedianCombine(pmReadout *combined, ///< Combined readout (output)
+				psArray *input       ///< Input array of pmStackData
+				);
 
 /// Stack input images
Index: trunk/psModules/src/objects/Makefile.am
===================================================================
--- trunk/psModules/src/objects/Makefile.am	(revision 34750)
+++ trunk/psModules/src/objects/Makefile.am	(revision 34800)
@@ -111,4 +111,5 @@
 	pmSourceOutputs.h \
 	pmSourceIO.h \
+	pmSourceSatstar.h \ 
 	pmSourcePlots.h \
 	pmSourceVisual.h \
Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c
===================================================================
--- trunk/psModules/src/objects/pmFootprintCullPeaks.c	(revision 34750)
+++ trunk/psModules/src/objects/pmFootprintCullPeaks.c	(revision 34800)
@@ -91,5 +91,16 @@
 
 	// max flux is above threshold for brightest peak
-	pmPeak *maxPeak = fp->peaks->data[0];
+      pmPeak *maxPeak = NULL;
+      for (int i = 0; i < fp->peaks->n; i++) {
+	pmPeak *testPeak = fp->peaks->data[i];
+	float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
+	
+	if (isfinite(this_peak)) {
+	  maxPeak = fp->peaks->data[i];
+	  break;
+	}
+      }
+      psAssert(maxPeak,"maxPeak was not set in these peaks");
+      //      = fp->peaks->data[0];
 	float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;