Index: /trunk/ippconfig/recipes/psphot.config
===================================================================
--- /trunk/ippconfig/recipes/psphot.config	(revision 34225)
+++ /trunk/ippconfig/recipes/psphot.config	(revision 34226)
@@ -233,4 +233,11 @@
 @RADIAL.ANNULAR.BINS.UPPER           F32      0.96 2.72 4.16 7.04 12.0 18.56 29.76 45.68 72.80 112.80 176.88 # perl parser fails without a comment ticket #1476 still?
 
+# In very dense fields it doesn't make sense to measure out to large radial apertures.
+# If the number of sources is above RADIAL.CUT.RADIUS.NUM.SOURCES.LIMIT 
+# bins with     RADIAL.ANNULAR.BINS.UPPER > RADIAL.SOURCES.OVER.LIMIT.RADIUS
+# will be removed from the recipe. (Note this is currently only used in psphotStack).
+RADIAL.NUM.SOURCES.LIMIT            S32         50000
+RAIDAL.SOURCES.OVER.LIMIT.RADIUS    F32         20.0    # pixels
+
 # Extended source fit parameters
 EXTENDED_TEST                       METADATA
@@ -354,4 +361,5 @@
 @PSPHOT.STACK.TARGET.PSF.FWHM       F32   6.0 8.0 # FWHM of target PSF (if NOT AUTO sized; pixels)
 PSPHOT.STACK.USE.RAW                BOOL  T
+PSPHOT.STACK.SPLIT.LINEAR.FIT       BOOL  T    # require positive flux for pass 3 sources then fit matched sources 
 
 RADIAL_APERTURES                    BOOL  F    # calculate flux in circular radial apertures?
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 34225)
+++ /trunk/psphot/src/psphot.h	(revision 34226)
@@ -264,4 +264,5 @@
 bool            psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS);
 bool            psphotEllipticalContour (pmSource *source);
+bool            psphotLimitRadialApertures(psMetadata *recipe, long nSources);
 
 // psphotVisual functions
Index: /trunk/psphot/src/psphotRadialBins.c
===================================================================
--- /trunk/psphot/src/psphotRadialBins.c	(revision 34225)
+++ /trunk/psphot/src/psphotRadialBins.c	(revision 34226)
@@ -202,4 +202,69 @@
 }
 
+// If the number of sources is large, edit the recipe to remove radial bins beyond
+// a value specified in the recipe
+bool psphotLimitRadialApertures(psMetadata *recipe, long nSources) {
+
+    bool status = false;
+    long sourceLimit = psMetadataLookupS32 (&status, recipe, "RADIAL.NUM.SOURCES.LIMIT");
+    if (!status) {
+        sourceLimit = 50000;
+    }
+    if (nSources < sourceLimit) {
+        return true;
+    }
+    psF32 maxRadius = psMetadataLookupF32 (&status, recipe, "RADIAL.SOURCES.OVER.LIMIT.RADIUS");
+    if (!status) {
+        maxRadius = 20.0;
+    }
+    psLogMsg ("psphot", PS_LOG_INFO, "Number of objects: %ld is greater than limit %ld. Limiting radial annular bins to %.3f pixels\n", 
+            nSources, sourceLimit, maxRadius);
+
+    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    if (!radMin || !radMin->n) {
+	psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)");
+	return false;
+    }
+    if (!radMax || !radMax->n) {
+	psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)");
+	return false;
+    }
+    if (radMax->n != radMin->n) {
+	psError (PSPHOT_ERR_CONFIG, true, "length of radMin %ld and radMax %ld not equal)", radMin->n, radMax->n);
+	return false;
+    }
+    int i = 0;
+    psF32 lastRadius = 0;
+    for (; i < radMax->n; i++) {
+        if (radMax->data.F32[i] > maxRadius) {
+            break;
+        }
+        lastRadius = radMax->data.F32[i];
+    }
+    if (i == radMax->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "Radius of all bins is within the limit. lastRadius: %.3f\n", lastRadius);
+        return true;
+    }
+    psLogMsg ("psphot", PS_LOG_INFO, "  radius limit exceeded at bin %d of %ld. New max radius is %.3f\n",
+            i, radMax->n, lastRadius);
+
+    psVector *radMinNew = psVectorRealloc(radMin, i);
+    if (radMinNew != radMin) {
+        psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.LOWER", PS_META_REPLACE, "", radMinNew);
+        // XXX: I don't need this so I?
+        // psFree(radMinNew);
+    }
+
+    psVector *radMaxNew = psVectorRealloc(radMax, i);
+    if (radMaxNew != radMax) {
+        psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.UPPER", PS_META_REPLACE, "", radMaxNew);
+        // XXX: I don't need to free this do I?
+        // psFree(radMaxNew);
+    }
+
+    return true;
+}
+
 // the area-weighted mean radius is given by:
 
Index: /trunk/psphot/src/psphotStackReadout.c
===================================================================
--- /trunk/psphot/src/psphotStackReadout.c	(revision 34225)
+++ /trunk/psphot/src/psphotStackReadout.c	(revision 34226)
@@ -286,11 +286,16 @@
     }
 
+    // gcc doesn't like the label to refer to a declaration so we declare this here
+    bool splitLinearFit;
+
 pass1finish:
 
-#ifdef notyet
-    // Split the fit of detected sources from matched sources. But not yet.
-    // NOTE: apply to ALL sources. Only include sources with postitive flux in the fit
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
-#endif
+    splitLinearFit = psMetadataLookupBool(NULL, recipe, "PSPHOT.STACK.SPLIT.LINEAR.FIT");
+    if (splitLinearFit) {
+        psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
+        // Fit the detected sources separately from matched that wea are about to create.
+        // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
+        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
+    }
 
     // generate the objects (objects unify the sources from the different images) NOTE: could
@@ -300,4 +305,8 @@
     psMemDump("matchsources");
 
+    // check the source density. If it too high change the number of radial bins
+    // in the recipe.
+    psphotLimitRadialApertures(recipe, objects->n);
+
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
@@ -308,11 +317,12 @@
     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
-#ifdef notyet
-    // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
-    // not be included in this fit.
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
-#else
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
-#endif
+    if (splitLinearFit) {
+        // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
+        // not be included in this fit.
+        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
+    } else {
+        // Fit all sources together
+        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
+    }
 
     // measure the radial profiles to the sky (only measures new objects)