Index: trunk/psphot/src/psphotModelTestReadout.c
===================================================================
--- trunk/psphot/src/psphotModelTestReadout.c	(revision 32348)
+++ trunk/psphot/src/psphotModelTestReadout.c	(revision 35559)
@@ -5,4 +5,6 @@
 
     bool status;
+
+    psTimerStart ("modelTest");
 
     pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
@@ -41,4 +43,11 @@
     }
 
+    // load the psf model, if suppled.  FWHM_MAJ,FWHM_MIN,etc are determined and saved on
+    // readout->analysis. NOTE: this function currently only loads from PSPHOT.PSF.LOAD
+    if (!psphotLoadPSF (config, view, filerule)) { // ??? need to supply 2 ?
+        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
     float MIN_KRON_RADIUS = 5.0;
 
@@ -53,4 +62,8 @@
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
+
+    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
 
     // find the currently selected readout
@@ -84,5 +97,5 @@
 
     // define the source of interest
-   float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
+    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
@@ -99,9 +112,12 @@
     source->type = PM_SOURCE_TYPE_EXTENDED;
 
+    bool TEST_FIT_CONVOLVED = psMetadataLookupBool (&status, recipe, "TEST_FIT_CONVOLVED");
+
     // find the local sky
     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
     if (!status) psAbort("pmSourceLocalSky error");
 
-    { 
+    bool testMomentsRadius = psMetadataLookupBool (&status, config->arguments, "TEST_MOMENTS_RADIUS");
+    if (testMomentsRadius) { 
 	// XXX I want to test an iterative aperture for brighter sources
 	float radius = mRADIUS;
@@ -109,5 +125,5 @@
 
 	    // get the source moments
-	    status = pmSourceMoments (source, radius, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+	    status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal);
 	    if (!status) psAbort("psSourceMoments error");
 
@@ -121,5 +137,5 @@
 
     // get the source moments
-    status = pmSourceMoments (source, mRADIUS, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+    status = pmSourceMoments (source, mRADIUS, 0.25*mRADIUS, 0.0, MIN_KRON_RADIUS, maskVal);
     if (!status) psAbort("psSourceMoments error");
 
@@ -135,6 +151,14 @@
     fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
 
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    if (psf) {
+      // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
+      source->peak->rawFlux = source->moments->Peak;
+      float Io = source->moments->Peak;
+      source->modelPSF = pmModelFromPSFforXY(psf, xObj, yObj, Io);
+    }
+
     // get the initial model parameter guess
-    pmModel *model = pmSourceModelGuess (source, modelType);
+    pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal);
     source->modelEXT = model;
 
@@ -159,12 +183,15 @@
     fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
 
-    // list model input shape
-    psEllipseShape shape;
-    shape.sx  = 1.4 / model->params->data.F32[4];
-    shape.sy  = 1.4 / model->params->data.F32[5];
-    shape.sxy = model->params->data.F32[6];
-    axes = psEllipseShapeToAxes (shape, 20.0);
-
-    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+	fprintf (stderr, "guess: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]);
+    } else {
+	// list model input shape
+	psEllipseShape shape;
+	shape.sx  = 1.4 / model->params->data.F32[4];
+	shape.sy  = 1.4 / model->params->data.F32[5];
+	shape.sxy = model->params->data.F32[6];
+	axes = psEllipseShapeToAxes (shape, 20.0);
+	fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+    }
 
     fprintf (stderr, "input parameters: \n");
@@ -190,5 +217,24 @@
     fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
 
-    status = pmSourceFitModel (source, model, fitOptions, maskVal);
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+	fitOptions->mode = PM_SOURCE_FIT_NO_INDEX;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+	fitOptions->mode = PM_SOURCE_FIT_TRAIL;
+    }
+
+    if (TEST_FIT_CONVOLVED) {
+      pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+      if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+	psphotSersicModelClassGuessPCM (pcm, source);
+      } else {
+	pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
+      }
+      pmPCMupdate(pcm, source, fitOptions, model);
+      pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
+      psFree (pcm);
+    } else {
+      status = pmSourceFitModel (source, model, fitOptions, maskVal);
+    }
 
     // measure the source mags