Index: /trunk/psModules/test/objects/Makefile.am
===================================================================
--- /trunk/psModules/test/objects/Makefile.am	(revision 10078)
+++ /trunk/psModules/test/objects/Makefile.am	(revision 10079)
@@ -17,5 +17,6 @@
 check_PROGRAMS = \
 	tap_pmSourcePhotometry \
-	tap_pmGrowthCurve
+	tap_pmGrowthCurve \
+	tap_pmSourceFitModel
 
 test: check
Index: /trunk/psModules/test/objects/tap_pmSourceFitModel.c
===================================================================
--- /trunk/psModules/test/objects/tap_pmSourceFitModel.c	(revision 10078)
+++ /trunk/psModules/test/objects/tap_pmSourceFitModel.c	(revision 10079)
@@ -7,12 +7,112 @@
 #include "pstap.h"
 
+bool fitModelFlux (psRandom *seed, float flux, float radius);
+
 int main (void)
 {
-
     pmModelGroupInit ();
+    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
 
     plan_tests(240);
 
-    diag("measure the photometry bias as a function of magnitude for a sky offset");
+    // build a gauss-deviate vector (mean = 0.0, sigma = 1.0)
+    psRandom *seed = psRandomAlloc (PS_RANDOM_TAUS, 0);
 
+    # if (0)
+
+        fitModelFlux (seed, 10000.0);
+    fitModelFlux (seed, 3000.0);
+    fitModelFlux (seed, 1000.0);
+    fitModelFlux (seed, 300.0);
+    fitModelFlux (seed, 100.0);
+    fitModelFlux (seed, 30.0);
+    # endif
+
+    fitModelFlux (seed, 30.0, 4.0);
+    fitModelFlux (seed, 30.0, 5.0);
+    fitModelFlux (seed, 30.0, 7.0);
+    fitModelFlux (seed, 30.0, 10.0);
+    fitModelFlux (seed, 30.0, 15.0);
+    fitModelFlux (seed, 30.0, 20.0);
+    fitModelFlux (seed, 30.0, 25.0);
+    fitModelFlux (seed, 30.0, 35.0);
+
+    return exit_status();
 }
+
+bool fitModelFlux (psRandom *seed, float flux, float radius)
+{
+
+    psMemId id = psMemGetId();
+    diag("test model fit");
+
+    psVector *rnd = psVectorAlloc (1000, PS_TYPE_F32);
+    for (int i = 0; i < rnd->n; i++) {
+        rnd->data.F32[i] = psRandomGaussian (seed);
+    }
+
+    // construct a model
+    pmSource *source = pmSourceAlloc ();
+    source->moments = pmMomentsAlloc ();
+
+    pmModelGroupInit ();
+    pmModelType type = pmModelSetType ("PS_MODEL_GAUSS");
+    source->modelEXT = pmModelAlloc (type);
+
+    source->modelEXT->params->data.F32[0] = 0;
+    source->modelEXT->params->data.F32[1] = flux;
+    source->modelEXT->params->data.F32[2] = 50;
+    source->modelEXT->params->data.F32[3] = 50;
+    source->modelEXT->params->data.F32[4] = 2.0*sqrt(2.0);
+    source->modelEXT->params->data.F32[5] = 2.0*sqrt(2.0);
+    source->modelEXT->params->data.F32[6] = 0;
+
+    source->pixels = psImageAlloc (100, 100, PS_TYPE_F32);
+    source->weight = psImageAlloc (100, 100, PS_TYPE_F32);
+    source->mask   = psImageAlloc (100, 100, PS_TYPE_U8);
+    psImageInit (source->pixels, 0.0);
+    psImageInit (source->weight, 0.0);
+    psImageInit (source->mask, 0);
+
+    // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
+    pmModelAdd (source->pixels, source->mask, source->modelEXT, false, false);
+    int npix = 0;
+    for (int j = 0; j < source->pixels->numRows; j++) {
+        for (int i = 0; i < source->pixels->numCols; i++) {
+            float flux = source->pixels->data.F32[j][i];
+            float var = flux + 100 + PS_SQR(5);
+            source->pixels->data.F32[j][i] += rnd->data.F32[npix]*sqrt(var);
+            source->weight->data.F32[j][i] = var;
+            npix ++;
+            if (npix == rnd->n)
+                npix = 0;
+        }
+    }
+
+    psFits *fits = psFitsOpen ("test.fits", "w");
+    psFitsWriteImage (fits, NULL, source->pixels, 0, NULL);
+    psFitsClose (fits);
+
+    // save the original model, modify params
+    pmModel *guess = pmModelCopy (source->modelEXT);
+    guess->params->data.F32[1] *= 0.9;
+    guess->params->data.F32[2] += 1.0;
+    guess->params->data.F32[3] -= 1.0;
+    guess->params->data.F32[4] *= 0.9;
+    guess->params->data.F32[5] *= 0.9;
+
+    pmSourceFitModelInit (15, 0.1, 100, true);
+
+    psImageKeepCircle (source->mask, 50, 50, radius, "OR", PM_MASK_MARK);
+    pmSourceFitModel (source, guess, PM_SOURCE_FIT_EXT);
+    psImageKeepCircle (source->mask, 50, 50, radius, "AND", PS_NOT_U8(PM_MASK_MARK));
+
+    ok_float_tol (guess->params->data.F32[1], source->modelEXT->params->data.F32[1], 100, "Io fitted: %f vs %f", guess->params->data.F32[1], source->modelEXT->params->data.F32[1]);
+    ok_float_tol (guess->params->data.F32[2], source->modelEXT->params->data.F32[2], 0.01, "Xo fitted: %f vs %f", guess->params->data.F32[2], source->modelEXT->params->data.F32[2]);
+    ok_float_tol (guess->params->data.F32[3], source->modelEXT->params->data.F32[3], 0.01, "Yo fitted: %f vs %f", guess->params->data.F32[3], source->modelEXT->params->data.F32[3]);
+    ok_float_tol (guess->params->data.F32[4], source->modelEXT->params->data.F32[4], 0.01, "Sx fitted: %f vs %f", guess->params->data.F32[4], source->modelEXT->params->data.F32[4]);
+    ok_float_tol (guess->params->data.F32[5], source->modelEXT->params->data.F32[5], 0.01, "Sy fitted: %f vs %f", guess->params->data.F32[5], source->modelEXT->params->data.F32[5]);
+
+    ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    return true;
+}