Index: trunk/psphot/src/fit_galaxies.c
===================================================================
--- trunk/psphot/src/fit_galaxies.c	(revision 4582)
+++ trunk/psphot/src/fit_galaxies.c	(revision 4630)
@@ -1,9 +1,7 @@
 # include "psphot.h"
-
-// fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
 
 bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
 { 
-    bool  status;
+    bool  status, goodfit;
     float x;
     float y;
@@ -13,28 +11,42 @@
     int   Niter = 0;
 
-    float MOMENT_R = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
-    // float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
-    float snFaint  = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
-    float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
-    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
-    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+    float  MOMENT_R    = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+    float  snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    float  OUTER       = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float  FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
+    float  FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+    char  *modelName   = psMetadataLookupPtr (&status, config, "GAL_MODEL");
 
     float FLUX_LIMIT  = FIT_NSIGMA * skyStats->sampleStdev;
 
-    psModelType   modelType   = psModelSetType ("PS_MODEL_SGAUSS");
+    psModelType   modelType   = psModelSetType (modelName);
     psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
-
-    psTraceSetLevel (".psModules.pmSourceMoments", 5);
 
     psTimerStart ("psphot");
     for (int i = 0; i < sources->n; i++) {
 	psSource *source = sources->data[i];
-	if (source->type != PS_SOURCE_GALAXY) continue;
-	if (source->moments->SN < snFaint) continue;
+
+	// sources which should not be fitted
+	// skip all valid stars
+	if (source->type == PS_SOURCE_PSFSTAR) continue;
+	if (source->type == PS_SOURCE_SATSTAR) continue;
+	if (source->type == PS_SOURCE_GOODSTAR) continue;
+	// skip all likely defects
+	if (source->type == PS_SOURCE_DEFECT) continue;
+	if (source->type == PS_SOURCE_SATURATED) continue;
+	// 
+	if (source->type == PS_SOURCE_FAINTSTAR) continue;
+	if (source->type == PS_SOURCE_POOR_FIT_PSF) continue;
+
+	// XXX when do we pick these up again?
+	if (source->moments->SN < snFaint) {
+	  source->type = PS_SOURCE_FAINT_GALAXY;  // better choice?
+	  continue;
+	}
 
 	// recalculate the source moments using the galaxy radius (larger)
 	status = pmSourceMoments_EAM (source, MOMENT_R);
 	if (!status) {
-	  fprintf (stderr, "invalid moments, skipping\n");
+	  source->type = PS_SOURCE_DROP_GALAXY;  // better choice?
 	  continue;
 	}
@@ -46,15 +58,9 @@
 	y = model->params->data.F32[3];
 
-	// need a better model guess and a better radius choice
-	// when radius is not fixed, we will need to check if new radius fits on image
-
 	// set the fit radius based on the object flux limit and the model
 	// FLUX_LIMIT should be set based on local sky model (not global median)
-	// model->radius = 25.0;
 	model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING;
-	if (isnan(model->radius)) {
-	  fprintf (stderr, "error in radius\n");
-	  continue;
-	}
+	if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
+
 	if (model->radius > OUTER) {
 	  // allocate image, noise, mask arrays for each peak (square of radius OUTER)
@@ -64,14 +70,26 @@
 	// fit as FLT, not PSF (skip poor fits)
 	psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
-	status = pmSourceFitModel (source, model, false);
+	status = pmSourceFitModel_EAM (source, model, false);
 	psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
-	if (!status || (model->params->data.F32[1] < 0)) {
+	if (!status) {
 	  // if the fit fails, we need to change the classification
 	  psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
-	  source->type = PS_SOURCE_OTHER;  // better choice?
+	  source->type = PS_SOURCE_FAIL_FIT_GAL;  // better choice?
+	  source->modelFLT = model;
 	  Nfail ++;
 	  continue;
 	}
 
+	goodfit = pmModelFitStatus (model);
+	if (!goodfit) {
+	  // if the fit fails, we need to change the classification
+	  psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
+	  source->type = PS_SOURCE_POOR_FIT_GAL;  // better choice?
+	  source->modelFLT = model;
+	  Nfail ++;
+	  continue;
+	}
+
+	source->type = PS_SOURCE_GALAXY;
 	source->modelFLT = model;
 	Niter += model[0].nIter;