Index: trunk/psModules/src/detrend/pmDark.c
===================================================================
--- trunk/psModules/src/detrend/pmDark.c	(revision 34085)
+++ 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;