Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 4573)
+++ /trunk/psphot/Makefile	(revision 4574)
@@ -39,4 +39,5 @@
 $(SRC)/subtract_galaxies.$(ARCH).o \
 $(SRC)/LocalSky.$(ARCH).o \
+$(SRC)/basic_classes.$(ARCH).o \
 $(SRC)/find_defects.$(ARCH).o
 
Index: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4573)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4574)
@@ -65,5 +65,10 @@
 	status = pmSourceFitModel (source, model, true);
 	psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
-	if (!status) continue;
+	if (!status || (model->params->data.F32[1] < 0)) {
+	  // if the fit fails, we need to change the classification
+	  psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
+	  source->type = PS_SOURCE_OTHER;  // better choice?
+	  continue;
+	}
 
 	source->modelPSF = model;
@@ -71,5 +76,5 @@
 	Nfit ++;
 
-	mark_psf_source (source, shapeNsigma);
+	mark_psf_source (source, shapeNsigma, SATURATE);
 	if (subtract_psf_source (source)) {
 	  Nsub ++;
Index: /trunk/psphot/src/basic_classes.c
===================================================================
--- /trunk/psphot/src/basic_classes.c	(revision 4573)
+++ /trunk/psphot/src/basic_classes.c	(revision 4574)
@@ -2,4 +2,6 @@
 
 bool basic_classes (psArray *sources, psMetadata *config) {
+
+    pmPSFClump psfClump;
 
     // group into STAR, COSMIC, GALAXY, SATURATED
Index: /trunk/psphot/src/fit_galaxies.c
===================================================================
--- /trunk/psphot/src/fit_galaxies.c	(revision 4573)
+++ /trunk/psphot/src/fit_galaxies.c	(revision 4574)
@@ -3,13 +3,27 @@
 // fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
 
-bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources) 
+bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
 { 
     bool  status;
+    float x;
+    float y;
+    float sky;
     int   Nfit = 0;
+    int   Nfail = 0;
     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");
-    psModelType modelType = psModelSetType ("PS_MODEL_PGAUSS");
+    float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
+    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+
+    float FLUX_LIMIT  = FIT_NSIGMA * skyStats->sampleStdev;
+
+    psModelType   modelType   = psModelSetType ("PS_MODEL_RGAUSS");
+    psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
+
+    psTraceSetLevel (".psModules.pmSourceMoments", 5);
 
     psTimerStart ("psphot");
@@ -19,22 +33,53 @@
 	if (source->moments->SN < snFaint) continue;
 
+	// recalculate the source moments using the galaxy radius (larger)
+	status = pmSourceMoments (source, MOMENT_R);
+
+	// use the source moments, etc to guess basic model parameters
+	psModel  *model  = pmSourceModelGuess (source, modelType); 
+
+	x = model->params->data.F32[2];
+	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
-# if (0)
+
+	// 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 (model->radius > OUTER) {
 	  // allocate image, noise, mask arrays for each peak (square of radius OUTER)
-	  pmSourceDefinePixels (source, imdata, x, y, OUTER);
+	  pmSourceDefinePixels (source, imdata, x, y, model->radius);
 	}
-# endif
-	psModel  *model  = pmSourceModelGuess (source, modelType); 
 
 	// fit as FLT, not PSF (skip poor fits)
-	if (!pmSourceFitModel (source, model, false)) continue;
+	psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
+	status = pmSourceFitModel (source, model, false);
+	psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
+	if (!status || (model->params->data.F32[1] < 0)) {
+	  // 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?
+	  Nfail ++;
+	  continue;
+	}
+
 	source->modelFLT = model;
-
 	Niter += model[0].nIter;
 	Nfit++;
+
+	// subtract object, leave local sky
+	sky = model->params->data.F32[0];
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
     }
-    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);
     return (true);
 }
Index: /trunk/psphot/src/mark_psf_source.c
===================================================================
--- /trunk/psphot/src/mark_psf_source.c	(revision 4573)
+++ /trunk/psphot/src/mark_psf_source.c	(revision 4574)
@@ -15,5 +15,5 @@
 // PS_SOURCE_BRIGHTSTAR 
 # define MIN_DS 0.01
-bool mark_psf_source (psSource *source, float shapeNsigma)
+bool mark_psf_source (psSource *source, float shapeNsigma, float SATURATE)
 { 
     float dSX, dSY, SX, SY, SN;
@@ -22,6 +22,17 @@
     if (source->modelPSF == NULL) return (false);
 
-    // SATSTAR objects are fitted with PSF regardless
-    if (source->type == PS_SOURCE_SATSTAR) return (true);
+    // if object has fitted peak above saturation, label as SATSTAR
+    // remember: fit does not use saturated pixels (masked)
+    if (source->modelPSF->params->data.F32[1] >= SATURATE) {
+	if (source->type == PS_SOURCE_PSFSTAR) {
+	    psLogMsg ("psphot", 3, "PSFSTAR marked saturated\n");
+	}
+	source->type = PS_SOURCE_SATSTAR;
+	return (true);
+    } 
+    if (source->type == PS_SOURCE_SATSTAR) {
+	psLogMsg ("psphot", 4, "SATSTAR marked bright (fitted peak below saturation)\n");
+	source->type = PS_SOURCE_BRIGHTSTAR;
+    }
 
     SN  = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1];
@@ -35,25 +46,25 @@
     // sigma = 1 / SX
     // dsx = 1 / (SX * SN)
-    // dsx_o = hypot (1/SX*SN + MIN_DSX)
+    // dsx_o = hypot (1/(SX*SN), MIN_DSX)
 
     // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
     if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
-      if (source->type == PS_SOURCE_PSFSTAR) return (true);
-      source->type = PS_SOURCE_BRIGHTSTAR;
-      return (true);
+	if (source->type == PS_SOURCE_PSFSTAR) return (true);
+	source->type = PS_SOURCE_BRIGHTSTAR;
+	return (true);
     }
     
     if (source->type == PS_SOURCE_PSFSTAR) {
-      psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
+	psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
     }
 
     if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
-      source->type = PS_SOURCE_GALAXY;
-      return (false);
+	source->type = PS_SOURCE_GALAXY;
+	return (false);
     }
     // replace DEFECT with COSMIC?
     if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
-      source->type = PS_SOURCE_DEFECT;
-      return (false);
+	source->type = PS_SOURCE_DEFECT;
+	return (false);
     }
     psTrace (".psphot.mark_psf_source", 2, "unexpected result: unmarked object\n");
Index: /trunk/psphot/src/psphot-utils.c
===================================================================
--- /trunk/psphot/src/psphot-utils.c	(revision 4573)
+++ /trunk/psphot/src/psphot-utils.c	(revision 4574)
@@ -96,4 +96,7 @@
 	model = (psModel  *) source->modelPSF;
 	if (model == NULL) continue;
+	if (source->type == PS_SOURCE_GALAXY) continue;
+	if (source->type == PS_SOURCE_DEFECT) continue;
+	if (source->type == PS_SOURCE_SATURATED) continue;
 	params = model->params;
 	dparams = model->dparams;
@@ -107,5 +110,7 @@
 	    fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
 	}
-	fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
+	fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", 
+		 source[0].type, log10(model[0].chisq), source[0].moments->SN, 
+		 model[0].radius, model[0].nDOF, model[0].nIter);
     }
     fclose (f);
@@ -134,4 +139,12 @@
 	model = (psModel  *) source->modelFLT;
 	if (model == NULL) continue;
+	if (source->type == PS_SOURCE_PSFSTAR) continue;
+	if (source->type == PS_SOURCE_SATSTAR) continue;
+	if (source->type == PS_SOURCE_BRIGHTSTAR) continue;
+	if (source->type == PS_SOURCE_FAINTSTAR) continue;
+	if (source->type == PS_SOURCE_OTHER) continue;
+	if (source->type == PS_SOURCE_DEFECT) continue;
+	if (source->type == PS_SOURCE_SATURATED) continue;
+
 	params = model->params;
 	dparams = model->dparams;
@@ -287,2 +300,84 @@
 }
 
+/* this code is in pmObjects.h / psModules:psEllipse */
+# if (0) 
+typedef struct {
+  double major;
+  double minor;
+  double theta;
+} EllipseAxes;
+
+typedef struct {
+  double x2;
+  double y2;
+  double xy;
+} EllipseMoments;
+
+typedef struct {
+  double sx;
+  double sy;
+  double sxy;
+} EllipseShape;
+
+EllipseAxes EllipseMomentsToAxes (EllipseMoments moments) {
+
+  EllipseAxes axes;
+
+  double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.sy));
+
+  axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f);
+  axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f)
+  axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2;
+  // theta in radians
+
+  return (axes);
+}
+
+EllipseShape EllipseMomentsToShape (EllipseMoments moments) {
+
+  EllipseShape shape;
+
+  double f = sqrt(0.25*PS_SQR(axes.x2 - axes.y2) + PS_SQR(axes.sxy));
+
+  shape.sx = moments.y2 / f;
+  shape.sy = moments.x2 / f;
+  shape.sxy = -2.0*moments.xy / f;
+  // XXX - these are probably wrong
+
+  return (shape);
+}
+
+EllipseMoments EllipseShapeToMoments (EllipseShape shape) {
+
+  EllipseMoments moments;
+
+  return (moments);
+}
+
+EllipseAxes EllipseShapeToAxes (EllipseShape shape) {
+
+  EllipseAxes axes;
+
+
+  return (axes);
+}
+
+EllipseShape EllipseAxesToShape (EllipseAxes axes) {
+
+  EllipseShape shape;
+
+  shape.sx = PS_SQR(cos(axes.theta) / axes.major) + PS_SQR(sin(axes.theta) / axes.minor);
+  shape.sy = PS_SQR(sin(axes.theta) / axes.major) + PS_SQR(cos(axes.theta) / axes.minor);
+  shape.sxy = cos(axes.theta)*sin(axes.theta)*(1.0/PS_SQR(axes.major) - 1.0/PS_SQR(axes.minor));
+
+  return (shape);
+}
+
+EllipseMoments EllipseAxesToMoments (EllipseAxes axes) {
+
+  EllipseMoments moments;
+
+  return (moments);
+}
+
+# endif
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 4573)
+++ /trunk/psphot/src/psphot.c	(revision 4574)
@@ -41,5 +41,5 @@
 
     // fit extended objects with galaxy models
-    // fit_galaxies (imdata, config, sources);
+    fit_galaxies (imdata, config, sources, sky);
 
     // subtract_galaxies (sources, config);
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4573)
+++ /trunk/psphot/src/psphot.h	(revision 4574)
@@ -46,9 +46,10 @@
 bool mark_psf_sources (psArray *sources, psMetadata *config);
 bool subtract_psf_sources (psArray *sources);
-bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
+bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
 bool subtract_galaxies (psArray *sources, psMetadata *config);
 bool subtract_psf_source (psSource *source);
-bool mark_psf_source (psSource *source, float shapeNsigma);
+bool mark_psf_source (psSource *source, float shapeNsigma, float SATURATE);
 bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
+bool basic_classes (psArray *sources, psMetadata *config);
 
 // psf utilities
Index: /trunk/psphot/src/source_moments.c
===================================================================
--- /trunk/psphot/src/source_moments.c	(revision 4573)
+++ /trunk/psphot/src/source_moments.c	(revision 4574)
@@ -5,5 +5,4 @@
     bool     status  = false;
     psArray *sources = NULL;
-    pmPSFClump psfClump;
 
     psTimerStart ("psphot");
Index: /trunk/psphot/src/subtract_psf_source.c
===================================================================
--- /trunk/psphot/src/subtract_psf_source.c	(revision 4573)
+++ /trunk/psphot/src/subtract_psf_source.c	(revision 4574)
@@ -9,6 +9,8 @@
 
   // non-stellar sources are ignored
+  if (source->type == PS_SOURCE_OTHER) return (false);
   if (source->type == PS_SOURCE_GALAXY) return (false);
   if (source->type == PS_SOURCE_DEFECT) return (false);
+  if (source->type == PS_SOURCE_SATURATED) return (false);
 
   psModel *model = source->modelPSF;
@@ -18,5 +20,5 @@
   }	    
 
-  psImage  *pixels = source->pixels;
+  psImage *pixels = source->pixels;
 
   // subtract object, leave local sky