Index: /branches/eam_branch_20070830/psLib/test/imageops/.cvsignore
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/.cvsignore	(revision 14862)
+++ /branches/eam_branch_20070830/psLib/test/imageops/.cvsignore	(revision 14863)
@@ -34,2 +34,5 @@
 tap_psImageStructManip
 tap_psImageMap
+
+tap_psImageMapFit
+tap_psImageMapFit2
Index: /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am	(revision 14862)
+++ /branches/eam_branch_20070830/psLib/test/imageops/Makefile.am	(revision 14863)
@@ -24,4 +24,5 @@
 	tap_psImageMap \
 	tap_psImageMapFit \
+	tap_psImageMapFit2 \
 	tap_psImageMaskOps
 
Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c	(revision 14862)
+++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit.c	(revision 14863)
@@ -6,4 +6,5 @@
 #include "pstap.h"
 
+// save function used to dump out test images while debugging algorithm
 int SaveImage (psMetadata *header, psImage *image, char *filename) {
 
@@ -14,22 +15,14 @@
 }
 
-# define TEST_4PT_0 0
-# define TEST_4PT_1 0
-# define TEST_4PT_2 0
-# define TEST_4PT_3 0
-
-# define TEST_6PT_0 0
-
-# define TEST_9PT_0 0
-# define TEST_9PT_1 0
+# define TEST_4PT_0 1
+# define TEST_4PT_1 1
+# define TEST_4PT_2 1
+# define TEST_4PT_3 1
+# define TEST_6PT_0 1
+# define TEST_9PT_0 1
+# define TEST_9PT_1 1
 # define TEST_9PT_2 1
-
-# define TEST1 0
-# define TEST2 0
-# define TEST3 0
-
-# define C00 +0.0
-# define C01 +0.0
-# define C10 +1.0
+# define TEST_9PT_3 1
+# define TEST_9PT_4 1
 
 int main (void)
@@ -41,6 +34,370 @@
 
     // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
+    # if (TEST_9PT_4)
+    {
+        psMemId id = psMemGetId();
+
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 6;
+	binning->nYfine = 6;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (50, PS_TYPE_F32);
+
+	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
+	// place the measurement points exactly on the ruff reference pixel centers
+	int n = 0;
+	for (float ix = 0.5; ix < 6.0; ix += 1.0) {
+	    for (float iy = 0.5; iy < 6.0; iy += 1.0) {
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		n++;
+	    }
+	}
+	x->n = n;
+	y->n = n;
+	f->n = n;
+
+	psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
+	for (int ix = 0; ix < 6; ix++) {
+	    for (int iy = 0; iy < 6; iy++) {
+		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
+
+	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
+	for (int ix = 0; ix < model->numCols; ix++) {
+	    for (int iy = 0; iy < model->numRows; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs");
+	    }
+	}
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
+    # if (TEST_9PT_3)
+    {
+        psMemId id = psMemGetId();
+
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 6;
+	binning->nYfine = 6;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (50, PS_TYPE_F32);
+
+	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
+	// place the measurement points exactly on the ruff reference pixel centers
+	int n = 0;
+	for (float ix = 1.0; ix < 6; ix += 2.0) {
+	    for (float iy = 1.0; iy < 6; iy += 2.0) {
+		y->data.F32[n] = iy;
+		x->data.F32[n] = ix;
+		if ((ix == 1.0) && (iy == 1.0)) {
+		    x->data.F32[n] = ix - 0.1;
+		}
+		if ((ix == 3.0) && (iy == 1.0)) {
+		    y->data.F32[n] = iy - 0.1;
+		}
+		if ((ix == 5.0) && (iy == 3.0)) {
+		    x->data.F32[n] = ix + 0.1;
+		}
+		if ((ix == 3.0) && (iy == 5.0)) {
+		    y->data.F32[n] = iy + 0.1;
+		}
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		n++;
+	    }
+	}
+	x->n = n;
+	y->n = n;
+	f->n = n;
+
+	psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
+	for (int ix = 0; ix < 6; ix++) {
+	    for (int iy = 0; iy < 6; iy++) {
+		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
+
+	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
+	for (int ix = 0; ix < model->numCols; ix++) {
+	    for (int iy = 0; iy < model->numRows; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs");
+	    }
+	}
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+	
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
+# if (TEST_9PT_2)
+    {
+        psMemId id = psMemGetId();
+
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 3;
+	binning->nYfine = 3;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (50, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (50, PS_TYPE_F32);
+
+	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
+	// place the measurement points exactly on the ruff reference pixel centers
+	int n = 0;
+	for (float ix = 0.0; ix < 3; ix += 1.0) {
+	    for (float iy = 0.0; iy < 3; iy += 1.0) {
+		x->data.F32[n] = ix + 0.5;
+		y->data.F32[n] = iy + 0.5;
+
+		if (ix == 0.0) {
+		    x->data.F32[n] -= 0.1;
+		}
+		if (iy == 0.0) {
+		    y->data.F32[n] -= 0.1;
+		}
+		if (ix == 2.0) {
+		    x->data.F32[n] += 0.1;
+		}
+		if (iy == 2.0) {
+		    y->data.F32[n] += 0.1;
+		}
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		n++;
+	    }
+	}
+	x->n = n;
+	y->n = n;
+	f->n = n;
+
+	psImage *field = psImageAlloc(3, 3, PS_TYPE_F32);
+	for (int ix = 0; ix < 3; ix++) {
+	    for (int iy = 0; iy < 3; iy++) {
+		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
+
+	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
+	for (int ix = 0; ix < model->numCols; ix++) {
+	    for (int iy = 0; iy < model->numRows; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], 1e-5, "model matches inputs");
+	    }
+	}
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+# endif
+
+    // still a simple test: 3x3 grid fitted to 9 points with simple slope and scale difference
+    # if (TEST_9PT_1)
+    {
+        psMemId id = psMemGetId();
+
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 6;
+	binning->nYfine = 6;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (9, PS_TYPE_F32);
+
+	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
+	// place the measurement points exactly on the ruff reference pixel centers
+	int n = 0;
+	for (int ix = 1; ix < 6; ix += 2) {
+	    for (int iy = 1; iy < 6; iy += 2) {
+		x->data.F32[n] = ix;
+		y->data.F32[n] = iy;
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		n++;
+	    }
+	}
+
+	psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
+	for (int ix = 0; ix < 6; ix++) {
+	    for (int iy = 0; iy < 6; iy++) {
+		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
+
+	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
+	for (int ix = 0; ix < model->numCols; ix++) {
+	    for (int iy = 0; iy < model->numRows; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // very simple test: 3x3 grid fitted to 9 points with simple slope
+    # if (TEST_9PT_0)
+    {
+	// function is defined over the range 0-1000, 0-1000
+        psMemId id = psMemGetId();
+
+	psImageBinning *binning = psImageBinningAlloc();
+	binning->nXfine = 3;
+	binning->nYfine = 3;
+	binning->nXruff = 3;
+	binning->nYruff = 3;
+
+	// generate a grid of test data points
+	psVector *x = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (9, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (9, PS_TYPE_F32);
+
+	int n = 0;
+	psImage *field = psImageAlloc(3, 3, PS_TYPE_F32);
+	for (int ix = 0; ix < 3; ix++) {
+	    for (int iy = 0; iy < 3; iy++) {
+		x->data.F32[n] = ix + 0.5;
+		y->data.F32[n] = iy + 0.5;
+		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+		field->data.F32[iy][ix] = f->data.F32[n];
+		n++;
+	    }
+	}
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
+
+	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
+	for (int ix = 0; ix < model->numCols; ix++) {
+	    for (int iy = 0; iy < model->numRows; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
+	    }
+	}
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    # endif
+
+    // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
     # if (TEST_6PT_0)
     {
+        psMemId id = psMemGetId();
+
 	psImageBinning *binning = psImageBinningAlloc();
 	binning->nXfine = 3;
@@ -94,9 +451,4 @@
 	// fit the data to the map
 	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
 
 	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
@@ -107,186 +459,18 @@
 	    }
 	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    # endif
-
-    // test for more points: 3x3 grid fitted to 9 points with simple slope and scale difference
-    # if (TEST_9PT_2)
-    {
-	psImageBinning *binning = psImageBinningAlloc();
-	binning->nXfine = 6;
-	binning->nYfine = 6;
-	binning->nXruff = 3;
-	binning->nYruff = 3;
-
-	// generate a grid of test data points
-	psVector *x = psVectorAlloc (50, PS_TYPE_F32);
-	psVector *y = psVectorAlloc (50, PS_TYPE_F32);
-	psVector *f = psVectorAlloc (50, PS_TYPE_F32);
-
-	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
-	// place the measurement points exactly on the ruff reference pixel centers
-	int n = 0;
-	for (float ix = 1.0; ix < 6; ix += 2.0) {
-	    for (float iy = 1.0; iy < 6; iy += 2.0) {
-		if ((ix == 3.0) && (iy == 1.0)) {
-		    x->data.F32[n] = ix + 0.0;
-		    y->data.F32[n] = iy + 0.1; // add in both points.  
-		    f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
-		    n++;
-		    x->data.F32[n] = ix + 0.0;
-		    y->data.F32[n] = iy - 0.1;
-		} else {
-		    x->data.F32[n] = ix;
-		    y->data.F32[n] = iy;
-		}
-		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
-		n++;
-	    }
-	}
-	x->n = n;
-	y->n = n;
-	f->n = n;
-
-	psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
-	for (int ix = 0; ix < 6; ix++) {
-	    for (int iy = 0; iy < 6; iy++) {
-		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
-	    }
-	}
-
-	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-
-	// scale defines both field and map image sizes (nXfine, nXruff)
-	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
-
-	// fit the data to the map
-	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
-
-	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
-	for (int ix = 0; ix < model->numCols; ix++) {
-	    for (int iy = 0; iy < model->numRows; iy++) {
-		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
-		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
-	    }
-	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    # endif
-
-    // still a simple test: 3x3 grid fitted to 9 points with simple slope and scale difference
-    # if (TEST_9PT_1)
-    {
-	psImageBinning *binning = psImageBinningAlloc();
-	binning->nXfine = 6;
-	binning->nYfine = 6;
-	binning->nXruff = 3;
-	binning->nYruff = 3;
-
-	// generate a grid of test data points
-	psVector *x = psVectorAlloc (9, PS_TYPE_F32);
-	psVector *y = psVectorAlloc (9, PS_TYPE_F32);
-	psVector *f = psVectorAlloc (9, PS_TYPE_F32);
-
-	// the underlying field is f = ix + iy, where ix,iy are fine pixel coordinates
-	// place the measurement points exactly on the ruff reference pixel centers
-	int n = 0;
-	for (int ix = 1; ix < 6; ix += 2) {
-	    for (int iy = 1; iy < 6; iy += 2) {
-		x->data.F32[n] = ix;
-		y->data.F32[n] = iy;
-		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
-		n++;
-	    }
-	}
-
-	psImage *field = psImageAlloc(6, 6, PS_TYPE_F32);
-	for (int ix = 0; ix < 6; ix++) {
-	    for (int iy = 0; iy < 6; iy++) {
-		field->data.F32[iy][ix] = (ix + 0.5) + (iy + 0.5);
-	    }
-	}
-
-	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-
-	// scale defines both field and map image sizes (nXfine, nXruff)
-	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
-
-	// fit the data to the map
-	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
-
-	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
-	for (int ix = 0; ix < model->numCols; ix++) {
-	    for (int iy = 0; iy < model->numRows; iy++) {
-		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
-		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
-	    }
-	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    # endif
-
-    // very simple test: 3x3 grid fitted to 9 points with simple slope
-    # if (TEST_9PT_0)
-    {
-	// function is defined over the range 0-1000, 0-1000
-	psImageBinning *binning = psImageBinningAlloc();
-	binning->nXfine = 3;
-	binning->nYfine = 3;
-	binning->nXruff = 3;
-	binning->nYruff = 3;
-
-	// generate a grid of test data points
-	psVector *x = psVectorAlloc (9, PS_TYPE_F32);
-	psVector *y = psVectorAlloc (9, PS_TYPE_F32);
-	psVector *f = psVectorAlloc (9, PS_TYPE_F32);
-
-	int n = 0;
-	psImage *field = psImageAlloc(3, 3, PS_TYPE_F32);
-	for (int ix = 0; ix < 3; ix++) {
-	    for (int iy = 0; iy < 3; iy++) {
-		x->data.F32[n] = ix + 0.5;
-		y->data.F32[n] = iy + 0.5;
-		f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
-		field->data.F32[iy][ix] = f->data.F32[n];
-		n++;
-	    }
-	}
-
-	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-
-	// scale defines both field and map image sizes (nXfine, nXruff)
-	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
-
-	// fit the data to the map
-	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
-
-	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
-	for (int ix = 0; ix < model->numCols; ix++) {
-	    for (int iy = 0; iy < model->numRows; iy++) {
-		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
-		is_float_tol (model->data.F32[iy][ix], field->data.F32[iy][ix], FLT_EPSILON, "model matches inputs");
-	    }
-	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
     # endif
@@ -298,4 +482,6 @@
     {
 	// function is defined over the range 0-1000, 0-1000
+        psMemId id = psMemGetId();
+
 	psImageBinning *binning = psImageBinningAlloc();
 	binning->nXfine = 4;
@@ -334,9 +520,4 @@
 	// fit the data to the map
 	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
 
 	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
@@ -347,6 +528,18 @@
 	    }
 	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
     # endif
@@ -358,4 +551,6 @@
     {
 	// function is defined over the range 0-1000, 0-1000
+        psMemId id = psMemGetId();
+
 	psImageBinning *binning = psImageBinningAlloc();
 	binning->nXfine = 2;
@@ -365,7 +560,7 @@
 
 	// generate a grid of test data points
-	psVector *x = psVectorAlloc (4, PS_TYPE_F32);
-	psVector *y = psVectorAlloc (4, PS_TYPE_F32);
-	psVector *f = psVectorAlloc (4, PS_TYPE_F32);
+	psVector *x = psVectorAlloc (5, PS_TYPE_F32);
+	psVector *y = psVectorAlloc (5, PS_TYPE_F32);
+	psVector *f = psVectorAlloc (5, PS_TYPE_F32);
 
 	int n = 0;
@@ -374,13 +569,19 @@
 	for (int ix = 0; ix < 2; ix++) {
 	    for (int iy = 0; iy < 2; iy++) {
-		// add in one extra point on the surface
+		# if (1)
 		if (!ix && !iy) {
-		    x->data.F32[n] = ix + 0.5;
-		    y->data.F32[n] = iy + 0.5 + 0.1;
+		    x->data.F32[n] = ix + 0.5 - 0.1;
+		    y->data.F32[n] = iy + 0.5 - 0.0;
+// 		    f->data.F32[n] = x->data.F32[n] + y->data.F32[n];
+// 		    n++;
+// 		    x->data.F32[n] = ix + 0.5;
+//		    y->data.F32[n] = iy + 0.5 - 0.1;
 		} else {
 		    x->data.F32[n] = ix + 0.5;
 		    y->data.F32[n] = iy + 0.5;
 		}
+		# endif
 		# if (0)
+		// offset points from centers
 		if (ix) {
 		    x->data.F32[n] = ix + 0.5 - 0.1;
@@ -399,17 +600,15 @@
 	    }
 	}
-
-	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-
-	// scale defines both field and map image sizes (nXfine, nXruff)
-	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
-
-	// fit the data to the map
-	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
+	x->n = n;
+	y->n = n;
+	f->n = n;
+
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+	// scale defines both field and map image sizes (nXfine, nXruff)
+	psImageMap *map = psImageMapAlloc (NULL, binning, stats);
+
+	// fit the data to the map
+	psImageMapFit (map, x, y, f, NULL);
 
 	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
@@ -420,6 +619,18 @@
 	    }
 	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
     # endif
@@ -429,4 +640,6 @@
     {
 	// function is defined over the range 0-1000, 0-1000
+        psMemId id = psMemGetId();
+
 	psImageBinning *binning = psImageBinningAlloc();
 	binning->nXfine = 4;
@@ -465,9 +678,4 @@
 	// fit the data to the map
 	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
 
 	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
@@ -478,6 +686,18 @@
 	    }
 	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
     # endif
@@ -487,4 +707,6 @@
     {
 	// function is defined over the range 0-1000, 0-1000
+        psMemId id = psMemGetId();
+
 	psImageBinning *binning = psImageBinningAlloc();
 	binning->nXfine = 2;
@@ -517,9 +739,4 @@
 	// fit the data to the map
 	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	SaveImage (NULL, map->map, "map.fits");
-	
-	SaveImage (NULL, field, "field.fits");
 
 	psImage *model = psImageAlloc(field->numCols, field->numRows, PS_TYPE_F32);
@@ -530,6 +747,18 @@
 	    }
 	}
-	SaveImage (NULL, model, "model.fits");
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+
+	// SaveImage (NULL, map->map, "map.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, model, "model.fits");
+
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
     # endif
Index: /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit2.c
===================================================================
--- /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit2.c	(revision 14862)
+++ /branches/eam_branch_20070830/psLib/test/imageops/tap_psImageMapFit2.c	(revision 14863)
@@ -1,6 +1,35 @@
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+
+#include "tap.h"
+#include "pstap.h"
+
+int SaveImage (psMetadata *header, psImage *image, char *filename) {
+
+    psFits *fits = psFitsOpen (filename, "w");
+    psFitsWriteImage (fits, NULL, image, 0, NULL);
+    psFitsClose (fits);
+    return (TRUE);
+}
+
+# define TEST1 0
+# define TEST2 1
+# define TEST3 0
+
+# define C00 +0.0
+# define C01 -2.0
+# define C10 +1.0
+
+int main (void)
+{
+
+    plan_tests(1);
 
     // make a model for a well-sampled field of a simple function (f = ax + by + c)
     # if (TEST1)
     {
+        psMemId id = psMemGetId();
+
 	// function is defined over the range 0-1000, 0-1000
 	psImageBinning *binning = psImageBinningAlloc();
@@ -33,12 +62,5 @@
 	// fit the data to the map
 	psImageMapFit (map, x, y, f, NULL);
-	psFree (binning);
-
-	// fprintf (stderr, "nGood: %d\n", map->nGood);
-	// fprintf (stderr, "nPoor: %d\n", map->nPoor);
-	// fprintf (stderr, "nBad:  %d\n", map->nBad);
-	
-	SaveImage (NULL, map->map, "map.fits");
-	
+
 	psImage *field = psImageAlloc(1000, 1000, PS_TYPE_F32);
 	for (int ix = 0; ix < 1000; ix++) {
@@ -47,5 +69,4 @@
 	    }
 	}
-	SaveImage (NULL, field, "field.fits");
 
 	// measure difference between model (map) and data (field)
@@ -54,9 +75,9 @@
 		int xo = psImageBinningGetFineX(map->binning, ix + 0.5);
 		int yo = psImageBinningGetFineY(map->binning, iy + 0.5);
-		float df = field->data.F32[yo][xo] - map->map->data.F32[iy][ix];
-		fprintf (stderr, "%d %d -> %d %d  :  %f = %f - %f\n", ix, iy, xo, yo, df, field->data.F32[yo][xo], map->map->data.F32[iy][ix]);
-	    }
-	}
-
+		is_float_tol (field->data.F32[yo][xo], map->map->data.F32[iy][ix], 1e-3, "model matches data");
+	    }
+	}
+
+	// XXX test on the model or don't bother
 	psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32);
 	for (int ix = 0; ix < 1000; ix++) {
@@ -65,5 +86,18 @@
 	    }
 	}
-	SaveImage (NULL, model, "model.fits");
+
+	// SaveImage (NULL, model, "model.fits");
+	// SaveImage (NULL, field, "field.fits");
+	// SaveImage (NULL, map->map, "map.fits");
+	
+	psFree (model);
+	psFree (binning);
+	psFree (map);
+	psFree (stats);
+	psFree (field);
+	psFree (x);
+	psFree (y);
+	psFree (f);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
     # endif
@@ -90,7 +124,7 @@
 	for (int ix = 0; ix < 1000; ix += 20) {
 	    for (int iy = 0; iy < 1000; iy += 20) {
-		x->data.F32[x->n] = ix;
-		y->data.F32[y->n] = iy;
-		f->data.F32[f->n] = C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500);
+		x->data.F32[x->n] = ix + 0,5;
+		y->data.F32[y->n] = iy + 0.5;
+		f->data.F32[f->n] = PS_SQR(x->data.F32[x->n]) + PS_SQR(y->data.F32[x->n]);
 		psVectorExtend (x, 100, 1);
 		psVectorExtend (y, 100, 1);
@@ -107,11 +141,7 @@
 	// XXX this function needs to correct for the mean superpixel position 
 	// which is sampled...
-	psImageMapGenerate (map, x, y, f, 0.1);
+	psImageMapFit (map, x, y, f, NULL);
 	psFree (binning);
 
-	fprintf (stderr, "nGood: %d\n", map->nGood);
-	fprintf (stderr, "nPoor: %d\n", map->nPoor);
-	fprintf (stderr, "nBad:  %d\n", map->nBad);
-	
 	SaveImage (NULL, map->map, "map.fits");
 	
@@ -119,5 +149,5 @@
 	for (int ix = 0; ix < 1000; ix++) {
 	    for (int iy = 0; iy < 1000; iy++) {
-		field->data.F32[iy][ix] = C10*PS_SQR(ix-500) + C01*PS_SQR(iy-500);
+		field->data.F32[iy][ix] = PS_SQR(ix+0.5) + PS_SQR(iy+0.5);
 	    }
 	}
@@ -133,4 +163,13 @@
 	    }
 	}
+
+	// XXX test on the model or don't bother
+	psImage *model = psImageAlloc(1000, 1000, PS_TYPE_F32);
+	for (int ix = 0; ix < 1000; ix++) {
+	    for (int iy = 0; iy < 1000; iy++) {
+		model->data.F32[iy][ix] = psImageUnbinPixel_V2 (ix + 0.5, iy + 0.5, map->map, map->binning);
+	    }
+	}
+	SaveImage (NULL, model, "model.fits");
     }
     # endif
@@ -215,2 +254,4 @@
     // choose a model scale for a poorly-sampled field of a non-polynomial function (f = (a(x-xo)^2 + b(y-yo)^2)^-1)
     
+    return exit_status();
+}