Index: trunk/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- trunk/psphot/src/psphotExtendedSourceFits.c	(revision 35559)
+++ trunk/psphot/src/psphotExtendedSourceFits.c	(revision 35769)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -52,4 +53,6 @@
     psTimerStart ("psphot.extended");
 
+    psphotFitInitExtended();
+
     // find the currently selected readout
     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
@@ -70,4 +73,5 @@
     if (!sources->n) {
 	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+	psphotSersicModelClassCleanup();
 	return true;
     }
@@ -80,5 +84,5 @@
     // do not thread if we are trying to study the fitting process
     if (psTraceGetLevel ("psphot.psphotFitEXT") >= 6) {
-      nThreads = 0;
+	nThreads = 0;
     }
 
@@ -107,17 +111,93 @@
     }
 
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+	// default to the old method (chisqConvergence)
+	chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+	// default to the old method (chisqConvergence)
+	gainFactorMode = 0;
+    }
+
+    // Define source fitting parameters for extended source fits
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->saveCovariance = true;  // XXX make this a user option?
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    fitOptions->nIter          = fitIter;
+    fitOptions->minTol         = fitMinTol;
+    fitOptions->maxTol         = fitMaxTol;
+
+    fitOptions->gainFactorMode   = gainFactorMode;
+    fitOptions->chisqConvergence = chisqConvergence;
+
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
 
     // select the collection of desired models
-    psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
+    psMetadata *allModels = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
     if (!status) {
         psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
+	psFree (fitOptions);
+	psphotSersicModelClassCleanup();
         return true;
     }
+    if (allModels->list->n == 0) {
+        psWarning ("extended source model fits requested but no models are specified\n");
+	psFree (fitOptions);
+	psphotSersicModelClassCleanup();
+        return true;
+    }
+
+    psMetadata *models = NULL;
+    char *modelSelection = psMetadataLookupStr (&status, recipe, "EXTENDED_SOURCE_MODELS_SELECTION");
+    if (!modelSelection || !status) {
+	models = psMemIncrRefCounter(allModels);
+	goto got_models;
+    } 
+
+    if (!strcasecmp(modelSelection, "ALL")) {
+	models = psMemIncrRefCounter(allModels);
+	goto got_models;
+    }
+
+    if (!strcasecmp(modelSelection, "NONE")) {
+	psWarning ("extended source model fits requested but no models are selected (EXTENDED_SOURCE_MODELS_SELECTION = NONE)\n");
+	psFree (fitOptions);
+	psphotSersicModelClassCleanup();
+	return true;
+    }
+
+    psArray *selection = psStringSplitArray (modelSelection, ",", false);
+    if (selection->n == 0) {
+	psWarning ("extended source model fits requested but model selection string is empty (EXTENDED_SOURCE_MODELS_SELECTION)\n");
+	psFree (fitOptions);
+	psFree (selection);
+	psphotSersicModelClassCleanup();
+	return true;
+    }
+    models = psMetadataAlloc();
+    for (int i = 0; i < selection->n; i++) {
+	psMetadata *model = psMetadataLookupMetadata (&status, allModels, selection->data[i]);
+	if (!model) {
+	    psWarning ("extended source model selection string includes an invalid name %s\n", (char *) selection->data[i]);
+	    continue;
+	}
+	psMetadataAddMetadata (models, PS_LIST_TAIL, selection->data[i], PS_META_REPLACE, "extended model", model);
+    }
+    psFree (selection);
+
     if (models->list->n == 0) {
-        psWarning ("extended source model fits requested but no models are specified\n");
-        return true;
-    }
+	psWarning ("extended source model fits requested but no valid models in selection string (EXTENDED_SOURCE_MODELS_SELECTION = %s)\n", modelSelection);
+	psFree (fitOptions);
+	psFree (models);
+	psphotSersicModelClassCleanup();
+	return true;
+    }
+
+got_models:
 
     psphotInitRadiusEXT (recipe, readout);
@@ -128,34 +208,40 @@
     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
 
-      if (item->type != PS_DATA_METADATA) {
-        // XXX we could cull the bad entries or build a validated model folder
-        psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
-      }
-
-      psMetadata *model = (psMetadata *) item->data.md;
-
-      // check on the model type
-      char *modelName = psMetadataLookupStr (&status, model, "MODEL");
-      int modelType = pmModelClassGetType (modelName);
-      if (modelType < 0) {
-        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
-      }
-      psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
-
-      // check on the SNLIM, set a float value
-      char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
-      if (!status) {
-        psAbort("SNLIM not defined for extended source model %s\n", item->name);
-      }
-      float SNlim = atof (SNword);
-      psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
-
-      // check on the PSF-Convolution status
-      char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
-      if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
-        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
-      }
-      bool convolved = !strcasecmp (convolvedWord, "true");
-      psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
+	if (item->type != PS_DATA_METADATA) {
+	    // XXX we could cull the bad entries or build a validated model folder
+	    psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
+	}
+
+	psMetadata *model = (psMetadata *) item->data.md;
+
+	// check on the model type
+	char *modelName = psMetadataLookupStr (&status, model, "MODEL");
+	int modelType = pmModelClassGetType (modelName);
+	if (modelType < 0) {
+	    psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
+	}
+	psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
+
+	// check on the SNLIM, set a float value
+	char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
+	if (!status) {
+	    psAbort("SNLIM not defined for extended source model %s\n", item->name);
+	}
+	float SNlim = atof (SNword);
+	psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
+
+	// check on the PSF-Convolution status
+	char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
+	if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
+	    psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
+	}
+	bool convolved = !strcasecmp (convolvedWord, "true");
+	psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
+
+	if (convolved) {
+	    psLogMsg ("psphot", PS_LOG_INFO, "using convolved model class %s (%s)", modelName, item->name);
+	} else {
+	    psLogMsg ("psphot", PS_LOG_INFO, "using simple    model class %s (%s)", modelName, item->name);
+	}
     }
     psFree (iter);
@@ -195,11 +281,8 @@
             psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
             psArrayAdd(job->args, 1, iter);
-            psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            psArrayAdd(job->args, 1, fitOptions);
 
             PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32);
-
             PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
@@ -218,4 +301,7 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
 		psFree(AnalysisRegion);
+		psFree (fitOptions);
+		psFree (models);
+		psphotSersicModelClassCleanup();
                 return false;
             } 
@@ -224,20 +310,23 @@
 		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
 		psFree(AnalysisRegion);
+		psFree (fitOptions);
+		psFree (models);
+		psphotSersicModelClassCleanup();
 		return false;
 	    }
 	    psScalar *scalar = NULL;
-	    scalar = job->args->data[8];
+	    scalar = job->args->data[9];
 	    Next += scalar->data.S32;
-	    scalar = job->args->data[9];
+	    scalar = job->args->data[10];
 	    Nconvolve += scalar->data.S32;
-	    scalar = job->args->data[10];
+	    scalar = job->args->data[11];
 	    NconvolvePass += scalar->data.S32;
-	    scalar = job->args->data[11];
+	    scalar = job->args->data[12];
 	    Nplain += scalar->data.S32;
-	    scalar = job->args->data[12];
+	    scalar = job->args->data[13];
 	    NplainPass += scalar->data.S32;
-	    scalar = job->args->data[13];
+	    scalar = job->args->data[14];
 	    Nfaint += scalar->data.S32;
-	    scalar = job->args->data[14];
+	    scalar = job->args->data[15];
 	    Nfail += scalar->data.S32;
             psFree(job->args->data[3]); // iterator allocated above
@@ -250,4 +339,7 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
 	    psFree(AnalysisRegion);
+	    psFree (fitOptions);
+	    psFree (models);
+	    psphotSersicModelClassCleanup();
             return false;
         }
@@ -260,17 +352,17 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[8];
+                scalar = job->args->data[9];
                 Next += scalar->data.S32;
-                scalar = job->args->data[9];
+                scalar = job->args->data[10];
                 Nconvolve += scalar->data.S32;
-                scalar = job->args->data[10];
+                scalar = job->args->data[11];
                 NconvolvePass += scalar->data.S32;
-                scalar = job->args->data[11];
+                scalar = job->args->data[12];
                 Nplain += scalar->data.S32;
-                scalar = job->args->data[12];
+                scalar = job->args->data[13];
                 NplainPass += scalar->data.S32;
-                scalar = job->args->data[13];
+                scalar = job->args->data[14];
                 Nfaint += scalar->data.S32;
-                scalar = job->args->data[14];
+                scalar = job->args->data[15];
                 Nfail += scalar->data.S32;
                 psFree(job->args->data[3]); // metadata iterator allocated above
@@ -281,4 +373,6 @@
     psFree (cellGroups);
     psFree(AnalysisRegion);
+    psFree (fitOptions);
+    psFree (models);
 
     psphotSersicModelClassCleanup();
@@ -290,4 +384,7 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit, %d failed\n", Nfaint, Nfail);
+
+    psphotFitSummaryExtended();
+
     return true;
 }
@@ -308,24 +405,16 @@
 
     // arguments: readout, sources, models, region, psfSize, maskVal, markVal
-    pmReadout *readout      = job->args->data[0];
-    psArray *sources        = job->args->data[1];
-    psMetadata *models      = job->args->data[2];
+    pmReadout *readout       = job->args->data[0];
+    psArray *sources         = job->args->data[1];
+    psMetadata *models       = job->args->data[2];
     psMetadataIterator *iter = job->args->data[3];
-    psRegion *region        = job->args->data[4];
-    int psfSize             = PS_SCALAR_VALUE(job->args->data[5],S32);
-    int fitIter             = PS_SCALAR_VALUE(job->args->data[6],S32);
-    float fitMinTol         = PS_SCALAR_VALUE(job->args->data[7],F32);
-    float fitMaxTol         = PS_SCALAR_VALUE(job->args->data[8],F32);
-    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA);
-
-    // Define source fitting parameters for extended source fits
-    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode           = PM_SOURCE_FIT_EXT;
-    fitOptions->saveCovariance = true;  // XXX make this a user option?
-    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-    fitOptions->nIter          = fitIter;
-    fitOptions->minTol         = fitMinTol;
-    fitOptions->maxTol         = fitMaxTol;
+    psRegion *region         = job->args->data[4];
+    pmSourceFitOptions *fitOptions = job->args->data[5];
+
+    int psfSize             = PS_SCALAR_VALUE(job->args->data[6],S32);
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[8],PS_TYPE_IMAGE_MASK_DATA);
+
+    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
 
     // choose the sources of interest
@@ -348,4 +437,10 @@
         if (source->peak->y > region->y1) continue;
 
+
+	// XXX for a test, just do the obvious trail
+	// XXX if (source->peak->xf < 1100) continue;
+	// XXX if (source->peak->xf > 1400) continue;
+	// XXX if (source->peak->yf >  245) continue;
+
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -384,6 +479,4 @@
           assert (status);
 
-          // limit selection to some SN limit
-          // assert (source->peak); // how can a source not have a peak?
 	  // limit selection to some SN limit
 	  bool skipSource = false;
@@ -424,35 +517,49 @@
               }
           } else {
-              psFree (source->modelFlux);
-              source->modelFlux = NULL;
-              modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal);
-              if (!modelFit) {
-                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
-		  Nfail ++;
-                  continue;
-              }
-              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
-		       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
-              Nplain ++;
-              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
-                  NplainPass ++;
-		  source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
-              }
+	      bool doneFits = false;
+	      while (!doneFits) {
+		  psFree (source->modelFlux);
+		  source->modelFlux = NULL;
+		  modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
+		  if (!modelFit) {
+		      psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
+		      Nfail ++;
+		      doneFits = true;
+		      continue;
+		  }
+		  psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
+		  Nplain ++;
+		  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+		      NplainPass ++;
+		      source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
+		  }
+		  doneFits = true;
+
+		  if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+		      // calculate the object end points:
+		      float *PAR = modelFit->params->data.F32;
+		      float Xmin = PAR[PM_PAR_XPOS] - 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
+		      float Ymin = PAR[PM_PAR_YPOS] - 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
+		      float Xmax = PAR[PM_PAR_XPOS] + 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
+		      float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
+
+		      if (false && (source->peak->xf > 1100) && 
+			  (source->peak->xf < 1400) && 
+			  (source->peak->yf < 245)) {
+			  fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", 
+				   source->pixels->col0, source->pixels->row0, 
+				   source->pixels->col0 + source->pixels->numCols, source->pixels->row0 + source->pixels->numRows, 
+				   Xmin, Ymin, Xmax, Ymax);
+		      }
+		      if (PAR[PM_PAR_LENGTH] > 0.9*fitRadius) {
+			  doneFits = false;
+			  fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
+			  psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
+		      }
+		  }
+	      }
           }
-
-	  // XXX deprecate?
-	  // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
-	  // XXX extFitPars->Mxx = source->moments->Mxx;
-	  // XXX extFitPars->Mxy = source->moments->Mxy;
-	  // XXX extFitPars->Myy = source->moments->Myy;
-	  // XXX extFitPars->Mrf = source->moments->Mrf;
-	  // XXX extFitPars->Mrh = source->moments->Mrh;
-	  // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
-	  
-	  // save kron mag, but assign apMag & psfMag on output (not yet calculated)
-	  // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
-
-          // psArrayAdd (source->extFitPars, 4, extFitPars);
-	  // psFree (extFitPars);
+	  // XXX really need to do this in a cleaner way:
+	  if (!modelFit) continue;
 
           // test for fit quality / result
@@ -519,28 +626,54 @@
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     }
-    psFree (fitOptions);
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[9];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[10];
+    scalar->data.S32 = Nconvolve;
+
     scalar = job->args->data[11];
-    scalar->data.S32 = Next;
+    scalar->data.S32 = NconvolvePass;
 
     scalar = job->args->data[12];
-    scalar->data.S32 = Nconvolve;
+    scalar->data.S32 = Nplain;
 
     scalar = job->args->data[13];
-    scalar->data.S32 = NconvolvePass;
+    scalar->data.S32 = NplainPass;
 
     scalar = job->args->data[14];
-    scalar->data.S32 = Nplain;
+    scalar->data.S32 = Nfaint;
 
     scalar = job->args->data[15];
-    scalar->data.S32 = NplainPass;
-
-    scalar = job->args->data[16];
-    scalar->data.S32 = Nfaint;
-
-    scalar = job->args->data[17];
     scalar->data.S32 = Nfail;
 
     return true;
 }
+
+# define PAD_WINDOW 3.0
+bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius) {
+
+    psRegion newRegion;
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = newRadius;
+    *windowRadius = *fitRadius + PAD_WINDOW;
+
+    // check to see if new region is completely contained within old region
+    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
+    newRegion = psRegionForImage (readout->image, newRegion);
+
+    // redefine the pixels to match
+    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}