Index: trunk/ppBackground/src/ppBackgroundStack.h
===================================================================
--- trunk/ppBackground/src/ppBackgroundStack.h	(revision 36647)
+++ trunk/ppBackground/src/ppBackgroundStack.h	(revision 36649)
@@ -23,4 +23,5 @@
 
   psImageMap *modelMap;
+  psS32 model_iteration;
   // These are the full extent of the input data
   psF32 ra_min;
Index: trunk/ppBackground/src/ppBackgroundStackData.c
===================================================================
--- trunk/ppBackground/src/ppBackgroundStackData.c	(revision 36647)
+++ trunk/ppBackground/src/ppBackgroundStackData.c	(revision 36649)
@@ -50,4 +50,5 @@
 
     data->modelMap = NULL;
+    data->model_iteration = 0;
     data->ra_min   = 1e9;
     data->ra_max   = -1e9;
Index: trunk/ppBackground/src/ppBackgroundStackLoop.c
===================================================================
--- trunk/ppBackground/src/ppBackgroundStackLoop.c	(revision 36647)
+++ trunk/ppBackground/src/ppBackgroundStackLoop.c	(revision 36649)
@@ -9,4 +9,7 @@
 
 #include "ppBackgroundStack.h"
+
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+
 
 bool ppBackgroundStackLoop(ppBackgroundStackData *data // Run-time data
@@ -120,16 +123,8 @@
       if (tp->y > data->y_max) { data->y_max = tp->y; }
 
-/*       data->x_min -= data->ra_min; */
-/*       data->x_max -= data->ra_min; */
-/*       data->y_min -= data->dec_min; */
-/*       data->y_max -= data->dec_min; */
-      
-      
       psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
       psImageBinning *binning = psImageBinningAlloc();
       binning->nXruff = 13; // Number of samples
       binning->nYruff = 13; 
-      //    binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
-      //    binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
       binning->nXfine = ceil(data->x_max - data->x_min) + 1;
       binning->nYfine = ceil(data->y_max - data->y_min) + 1;
@@ -146,12 +141,4 @@
       P_PSIMAGE_SET_ROW0(sizeImage, data->y_min);
       data->modelMap = psImageMapAlloc(sizeImage,binning,stats);
-/*       psFree(sizeImage); */
-/*       data->modelMap = psImageMapNoImageAlloc( binning,stats); */
-/*       P_PSIMAGE_SET_COL0(data->modelMap->map, (data->x_min - binning->nXskip) / binning->nXbin); */
-/*       P_PSIMAGE_SET_ROW0(data->modelMap->map, (data->y_min - binning->nYskip) / binning->nYbin); */
-
-      // force col0/row0
-/*       data->modelMap->map->col0 = data->modelMap->binning->nXskip; */
-/*       data->modelMap->map->row0 = data->modelMap->binning->nYskip; */
       
       // PART 2:
@@ -173,5 +160,5 @@
 
       // This is where an iterative solution loop would likely start.
-      for (int iterator = 0; iterator < 2; iterator++) {
+      for (int iterator = 0; iterator < 4; iterator++) {
 	// Construct the offset information
 	printf("Model fit!\n");
@@ -310,4 +297,45 @@
 	    psFree(fp);
 	    psFree(tp);
+
+	    // Copy WCS (from ppStackUpdateHeader)
+	    pmHDU *inHDU = pmHDUFromCell(readout->parent);
+	    model->parent->hdu = pmHDUAlloc(NULL);
+	    corr->parent->hdu = pmHDUAlloc(NULL);
+	    pmHDU *modHDU= pmHDUFromCell(model->parent);
+	    pmHDU *corHDU= pmHDUFromCell(corr->parent);
+
+	    if (!modHDU || !inHDU) {
+	      psWarning("Unable to find HDU at FPA level to copy wcs!");
+	    }
+	    else {
+	      if (!pmAstromReadWCS(stack_model->fpa,model_cell->parent,inHDU->header,1.0)) {
+		psErrorClear();
+		psWarning("Unable to read WCS astrometry from input FPA!");
+	      }
+	      else {
+		if (!modHDU->header) {
+		  modHDU->header = psMetadataAlloc();
+		}
+		if (!pmAstromWriteWCS(modHDU->header, stack_model->fpa,model_cell->parent, WCS_TOLERANCE)) {
+		  psErrorClear();
+		  psWarning("Unable to read WCS astrometry from input FPA!");
+		}
+	      }
+	      if (!pmAstromReadWCS(stack_corr->fpa,corr_cell->parent,inHDU->header,1.0)) {
+		psErrorClear();
+		psWarning("Unable to read WCS astrometry from input FPA!");
+	      }
+	      else {
+		if (!corHDU->header) {
+		  corHDU->header = psMetadataAlloc();
+		}
+		if (!pmAstromWriteWCS(corHDU->header, stack_corr->fpa,corr_cell->parent, WCS_TOLERANCE)) {
+		  psErrorClear();
+		  psWarning("Unable to read WCS astrometry from input FPA!");
+		}
+	      }
+	    } // End WCS saving.
+
+	    
 	  } // Close readout
 	  printf("    I'm done with that readout\n");
@@ -339,4 +367,5 @@
       psFree(view);
       psFree(data->modelMap);
+      psFree(sizeImage);
     }
 		
Index: trunk/ppBackground/src/ppBackgroundStackMath.c
===================================================================
--- trunk/ppBackground/src/ppBackgroundStackMath.c	(revision 36647)
+++ trunk/ppBackground/src/ppBackgroundStackMath.c	(revision 36649)
@@ -55,5 +55,7 @@
 	  }
 	  psMetadata *chipData = psMetadataLookupPtr(NULL, expItem->data.md, workingChip);
+	  if (!chipData) { continue; }
 	  psImage *image      = psMetadataLookupPtr(NULL, chipData, "bkg image");
+	  if (!image) { continue; }
 	  psVectorAppend(tmp,image->data.F32[v][u]);
 	} // End loop over exposures
@@ -63,6 +65,17 @@
 
 	psFree(expIter);
+	psFree(tmp);
       } // End u
     } // End v
+    // Remove the median value from this data.  We just want the tilts, not the offsets
+    psStatsInit(stats);
+    psImageStats(stats,solution,NULL,0);
+    for (v = 0; v < solution->numRows; v++) {
+      for (u = 0; u < solution->numCols; u++) {
+	solution->data.F32[v][u] -= stats->robustMedian;
+      }
+    }
+    
+    
   } // End working chip scan
 
@@ -90,27 +103,37 @@
     psMetadataItem *chipItem;
     while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
+      //      const char *chipName = chipItem->name;
+      
       psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
       psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
       psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
-      
+      //      psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName);      
       psVector *obs  = psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
       psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
-      
+
+      int j = 0;
+      int used = 0;
       for (v = 0; v < image->numRows; v++) {
 	for (u = 0; u < image->numCols; u++) {
 	  if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)||
-	      (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
-	  psVectorAppend(obs,image->data.F32[v][u]);
+	      (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { j++; continue; }
+	  psVectorAppend(obs,image->data.F32[v][u]);// - camera->data.F32[v][u]);
 	  psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]));
-	}
-      }
-
-      psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
-      int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
-      if (!status) {
-	psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
-	psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
-      }
-      psFree(poly);
+	  j++;
+	  used++;
+	}
+      }
+
+      if (used > 0) {
+
+	psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
+	int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
+	printf("in model fit loop: %d %d %d %f %f\n",status,j,used,poly->coeff[0],poly->coeff[1]);
+	if (status && (poly->coeff[1] != 0.0)) {
+	  psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
+	  psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
+	}
+	psFree(poly);
+      }
     } // End OTA loop
     psFree(chipIter);
@@ -153,5 +176,5 @@
 	      (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
 
-	  model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u];
+	  model->data.F32[v][u] = scale * image->data.F32[v][u] + offset - camera->data.F32[v][u];
 	}
       }
@@ -168,8 +191,8 @@
 // This "averaging" is done using the psImageMapClipFit.
 bool ppBackgroundStackModelFit(ppBackgroundStackData *data) {
-  long j;
-  int u,v;
-
-  long used;
+  long j = 0;
+  int u,v;
+
+  long used = 0;
   psS16 N = psMetadataLookupS16(NULL, data->models, "N");
   psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
@@ -195,5 +218,12 @@
       psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
       psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
-      //      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
+#define DUMP_DATA 0
+#if DUMP_DATA
+      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
+
+      psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset");
+      psF32 scale  = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale");
+
+#endif
       for (v = 0; v < calib->numRows; v++) {
 	for (u = 0; u < calib->numCols; u++) {
@@ -214,12 +244,13 @@
 	  used++;
 	  j++;
-	  
-/* 	  printf("DATA %ld %ld %f %f %f %f\n", */
-/* 		 j,used, */
-/* 		 ra->data.F32[v][u], */
-/* 		 dec->data.F32[v][u], */
-/* 		 calib->data.F32[v][u], */
-/* 		 model->data.F32[v][u]); */
-		 
+#if DUMP_DATA
+	  printf("DATA %d %ld %ld %f %f %f %f %f %f\n",
+		 data->model_iteration,j,used,
+		 offset,scale,
+		 ra->data.F32[v][u],
+		 dec->data.F32[v][u],
+		 calib->data.F32[v][u],
+		 model->data.F32[v][u]);
+#endif 
 	}
       }
@@ -230,5 +261,5 @@
   bool fitStatus;
   bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E);
-
+  data->model_iteration++;
   psFree(expIter);
   psFree(X);