Index: /trunk/psphot/doc/notes.20100715.txt
===================================================================
--- /trunk/psphot/doc/notes.20100715.txt	(revision 29547)
+++ /trunk/psphot/doc/notes.20100715.txt	(revision 29548)
@@ -1,2 +1,18 @@
+
+2010.08.24
+
+  Remaining work to be done:
+
+  * test and turn on CR masking
+  * double-check source size analysis:
+    * are we doing the right thing with sources flagged as bad in 'rough'?
+    * are we doing the right thing for SAT, CR, EXT, PSF after SourceSize?
+  * convert CR_SIGMA and EXT_SIGMA to probabilities
+  * example results of model fits vs fake inputs
+  * what is the right solution for the stack PSF photometry?
+    * Nigel sees 4-6%
+    * we see 1.9% at the worst
+    * non-poisson errors are clearly better -- check on the chisq values
+    * why is QGAUSS so good for the one example (better for poisson, if trendy)
 
 2010.08.12
Index: /trunk/psphot/doc/notes.20101013.txt
===================================================================
--- /trunk/psphot/doc/notes.20101013.txt	(revision 29548)
+++ /trunk/psphot/doc/notes.20101013.txt	(revision 29548)
@@ -0,0 +1,80 @@
+
+20101019
+
+  We can identify CRs in the (M_minor, KronMag) plane as a fairly
+  tight clump (M_minor ~< 1.0, KronMag ~< -8.0).  The best bounds of
+  the clump vary somewhat from image to image.  I've added a dynamic
+  CR limit function which looks for the peak in M_minor and in KronMag
+  in those regions, then finds the valley floor.  It seems to work
+  fairly well.  Outstanding issues:
+
+  * only adjust the limits once
+  * record limits in the headers
+  * record CR counts in the headers
+  * mask the CRs.
+
+20101017
+  
+  pmSourceMoments:
+    * FWHM(src) vs FWHM(win) & 1st pass centroid accuracy (no noise)
+      * 4.0 : 6.0 -> 30% 
+      * 4.0 : 8.0 -> 20%
+      * 4.0 : 15.0 -> 6.7%
+   * FWHM(src) vs Tophat(radius)
+      * 4.0 : 5.0 -> ~30%
+      * 4.0 : 10.0 -> 0.2%
+      * 4.0 : 8.0 -> 3%
+      * 4.0 : 6.0 -> 18%
+      * 4.0 : 7.0 -> 9%
+
+  pmSourceMomentsGetCentroid(source, 8.0, 0.0, 0.0, 0);
+  pmSourceMomentsGetCentroid(source, 10.0, 8.0/2.35, 0.0, 0);
+
+  Io = 10000 : (S/N ~ 400)
+
+dX stats: -0.000708 +/- 0.004586
+dY stats: -0.000296 +/- 0.004678
+
+  Io = 1000 : (S/N ~ 100)
+
+dX stats: -0.000131 +/- 0.020428
+dY stats: 0.001038 +/- 0.020556
+
+  Io = 100 : (S/N ~ 10)
+
+dX stats: -0.025900 +/- 0.166111
+dY stats: 0.025336 +/- 0.169681
+
+raising bits : these functions may raise:
+
+psphotSourceStats : BELOW_MOMENTS_SN | SKY_FAILURE | SKYVAR_FAILURE | BIG_RADIUS | MOMENTS_FAILURE
+
+psphotBasicDeblend : BLEND 
+
+psphotDeblendSatstars : BLEND 
+
+psphotRoughClass : SATSTAR | SATURATED | DEFECT | PSFSTAR
+
+psphotChoosePSF : can remove PSFSTAR from failed candidates
+
+psphotGuessModels : BADPSF
+
+psphotMergeSources : PSFSTAR | EXTERNAL 
+  * note : this filters out 'bad' sources -- maybe too restrictive?
+  
+psphotFitSourcesLinearStack : LINEAR_FIT
+
+psphotSourceSize : SIZE_SKIPPED | DEFECT | EXT_LIMIT | CR_LIMIT 
+
+psphotApResid : AP_MAGS
+
+psphotRadialAperturesByObject RADIAL_FLUX | EXTENDED_STATS
+
+psphotExtendedSourceFits : EXTENDED_FIT
+
+****
+
+flags from psastro:
+
+ OFF_CHIP | ON_GHOST | ON_GLINT | ON_SPIKE
+
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 29547)
+++ /trunk/psphot/src/psphot.h	(revision 29548)
@@ -167,5 +167,5 @@
 
 // used by psphotFindDetections
-psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal);
+psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
 psArray        *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
 bool            psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
@@ -258,10 +258,12 @@
 bool            psphotVisualShowFlags (psArray *sources);
 bool            psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
+bool            psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources);
 bool            psphotVisualShowResidualImage (pmReadout *readout);
-bool            psphotVisualPlotApResid (psArray *sources, float mean, float error);
+bool            psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag);
 bool            psphotVisualPlotChisq (psArray *sources);
 bool            psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources);
 bool            psphotVisualShowPetrosians (psArray *sources);
 bool            psphotVisualEraseOverlays (int channel, char *overlay);
+bool            psphotVisualClose(void);
 
 bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
Index: /trunk/psphot/src/psphotApResid.c
===================================================================
--- /trunk/psphot/src/psphotApResid.c	(revision 29547)
+++ /trunk/psphot/src/psphotApResid.c	(revision 29548)
@@ -339,5 +339,5 @@
     psFree (dMag);
 
-    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
+    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid, true);
 
     return true;
Index: /trunk/psphot/src/psphotCleanup.c
===================================================================
--- /trunk/psphot/src/psphotCleanup.c	(revision 29547)
+++ /trunk/psphot/src/psphotCleanup.c	(revision 29548)
@@ -14,4 +14,5 @@
     psFree (config);
 
+    psphotVisualClose();
     psMemCheckCorruption (stderr, true);
     pmModelClassCleanup ();
@@ -19,4 +20,5 @@
     pmConceptsDone ();
     pmConfigDone ();
+    pmVisualCleanup ();
     psLibFinalize();
     // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
Index: /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c
===================================================================
--- /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c	(revision 29547)
+++ /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c	(revision 29548)
@@ -38,7 +38,5 @@
     // which extended source analyses should we perform?
     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
-    bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
-    bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
 
     // number of images used to define sources
@@ -117,13 +115,16 @@
 
 	    // if we request any of these measurements, we require the radial profile
-	    if (doPetrosian || doIsophotal || doAnnuli || doKron) {
+	    if (doPetrosian || doAnnuli) {
 		if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
 		    // all measurements below require the radial profile; skip them all
 		    // re-subtract the object, leave local sky
-		    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+		    psTrace ("psphot", 5, "FAILED radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 		    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 		    continue;
+		} else {
+		    psTrace ("psphot", 5, "measured radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+		    Nannuli ++;
+		    source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
 		}
-		source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
 	    }
 
@@ -138,5 +139,4 @@
 		}
 	    }
-
 
 	    // re-subtract the object, leave local sky
Index: /trunk/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotExtendedSourceFits.c	(revision 29547)
+++ /trunk/psphot/src/psphotExtendedSourceFits.c	(revision 29548)
@@ -85,4 +85,6 @@
         return true;
     }
+
+    psphotInitRadiusEXT (recipe, readout);
 
     // validate the model entries
@@ -167,9 +169,27 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
 
-            if (!psThreadJobAddPending(job)) {
+            if (false && !psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
 		psFree(AnalysisRegion);
                 return false;
-            }
+            } else {
+		if (!psphotExtendedSourceFits_Threaded(job)) {
+		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		    psFree(AnalysisRegion);
+		    return false;
+		}
+                psScalar *scalar = NULL;
+                scalar = job->args->data[7];
+                Next += scalar->data.S32;
+                scalar = job->args->data[8];
+                Nconvolve += scalar->data.S32;
+                scalar = job->args->data[9];
+                NconvolvePass += scalar->data.S32;
+                scalar = job->args->data[10];
+                Nplain += scalar->data.S32;
+                scalar = job->args->data[11];
+                NplainPass += scalar->data.S32;
+		psFree(job);
+	    }
         }
 
@@ -272,5 +292,5 @@
 	    fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
 	    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-	    // XXX raise an error flag of some kind
+	    // XXX raise an error of some kind
 	    continue;
 	}
@@ -296,5 +316,6 @@
 	    if (!modelFluxStart) {
 		fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
-		// XXX raise an error flag of some kind?
+		pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+		// XXX raise an error of some kind?
 		continue;
 	    }
Index: /trunk/psphot/src/psphotFindDetections.c
===================================================================
--- /trunk/psphot/src/psphotFindDetections.c	(revision 29547)
+++ /trunk/psphot/src/psphotFindDetections.c	(revision 29548)
@@ -85,5 +85,5 @@
 
     // generate the smoothed significance image
-    psImage *significance = psphotSignificanceImage (readout, recipe, pass, maskVal);
+    psImage *significance = psphotSignificanceImage (readout, recipe, maskVal);
 
     // display the significance image
Index: /trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 29547)
+++ /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 29548)
@@ -104,5 +104,5 @@
     float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
     if (!status) {
-        MIN_VALID_FLUX = 1e-8;
+        MIN_VALID_FLUX = 0.0;
     }
     float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");
Index: /trunk/psphot/src/psphotFitSourcesLinearStack.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinearStack.c	(revision 29547)
+++ /trunk/psphot/src/psphotFitSourcesLinearStack.c	(revision 29548)
@@ -35,4 +35,13 @@
     bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
+
+    float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
+    if (!status) {
+        MIN_VALID_FLUX = 0.0;
+    }
+    float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");
+    if (!status) {
+        MAX_VALID_FLUX = 1e+8;
+    }
 
     // XXX store a local static array of covar factors for each of the images (by image ID)
@@ -123,7 +132,7 @@
 
     psSparseConstraint constraint;
-    constraint.paramMin   = 0.0;
-    constraint.paramMax   = 1e8;
-    constraint.paramDelta = 1e8;
+    constraint.paramMin   = MIN_VALID_FLUX;
+    constraint.paramMax   = MAX_VALID_FLUX;
+    constraint.paramDelta = 1e7;
 
     // solve for normalization terms (need include local sky?)
Index: /trunk/psphot/src/psphotGuessModels.c
===================================================================
--- /trunk/psphot/src/psphotGuessModels.c	(revision 29547)
+++ /trunk/psphot/src/psphotGuessModels.c	(revision 29548)
@@ -176,9 +176,14 @@
 
         // We have two options to get a guess for the object position: the position from the
-        // peak and the position from the moments.  Use the peak position if (a) there are no
-        // moments and (b) the sources is not saturated
-
+        // peak and the position from the moments.  Use the peak position if there are no
+        // moments
+
+# if (0)
+        bool useMoments = true;
+# else
         bool useMoments = false;
         useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR);  // we only want to try if SATSTAR is set, but..
+# endif
+
         useMoments = (useMoments && source->moments);          // can't if there are no moments
         useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
Index: /trunk/psphot/src/psphotImageLoop.c
===================================================================
--- /trunk/psphot/src/psphotImageLoop.c	(revision 29547)
+++ /trunk/psphot/src/psphotImageLoop.c	(revision 29548)
@@ -30,4 +30,10 @@
     // files associated with the science image
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
 
     // for psphot, we force data to be read at the chip level
@@ -95,5 +101,5 @@
                 if (readout->mask) {
                     psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
-                    if (!pmReadoutMaskNonfinite(readout, maskSat)) {
+                    if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) {
                         psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
                         psFree(view);
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 29547)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 29548)
@@ -226,5 +226,5 @@
     // measure the source mags
     pmSourcePhotometryModel (&fitMag, model);
-    pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->maskObj, maskVal);
+    pmSourcePhotometryAper  (&obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
     fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
 
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 29547)
+++ /trunk/psphot/src/psphotOutput.c	(revision 29548)
@@ -270,4 +270,9 @@
     psMetadataItemSupplement (&status, header, analysis, "NDET_CR");
 
+    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID");
+    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID.SYSERR");
+    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.CR.MAX.SIZE");
+    psMetadataItemSupplement (&status, header, analysis, "PSPHOT.CR.MAX.MAG");
+
     // sky background model statistics
     psMetadataItemSupplement (&status, header, analysis, "MSKY_MN");
Index: /trunk/psphot/src/psphotSignificanceImage.c
===================================================================
--- /trunk/psphot/src/psphotSignificanceImage.c	(revision 29547)
+++ /trunk/psphot/src/psphotSignificanceImage.c	(revision 29548)
@@ -4,5 +4,5 @@
 // (S/N)^2.  If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the
 // smoothing kernel.
-psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal) {
+psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) {
 
     float SIGMA_SMTH, NSIGMA_SMTH;
@@ -67,4 +67,5 @@
     // XXX change these to recipe value checks
     if (psTraceGetLevel("psphot") > 5) {
+	static int pass = 0;
         char name[64];
         sprintf (name, "imsmooth.v%d.fits", pass);
@@ -72,4 +73,5 @@
         sprintf (name, "wtsmooth.v%d.fits", pass);
         psphotSaveImage(NULL, smooth_wt, name);
+	pass ++;
     }
 
@@ -108,6 +110,8 @@
     if (psTraceGetLevel("psphot") > 5) {
         char name[64];
+	static int pass = 0;
         sprintf (name, "snsmooth.v%d.fits", pass);
         psphotSaveImage (NULL, smooth_im, name);
+	pass ++;
     }
 
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 29547)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 29548)
@@ -551,4 +551,5 @@
 
     	if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
+	    psFree (pcm);
 	    model->flags |= PM_MODEL_STATUS_BADARGS;
 	    return model;
Index: /trunk/psphot/src/psphotSourceMatch.c
===================================================================
--- /trunk/psphot/src/psphotSourceMatch.c	(revision 29547)
+++ /trunk/psphot/src/psphotSourceMatch.c	(revision 29548)
@@ -228,4 +228,5 @@
 	    pmSource *source = pmSourceAlloc();
 	    source->imageID = index;
+	    source->mode2 |= PM_SOURCE_MODE2_MATCHED; // source is generated based on another image
 
 	    // add the peak
Index: /trunk/psphot/src/psphotSourceSize.c
===================================================================
--- /trunk/psphot/src/psphotSourceSize.c	(revision 29547)
+++ /trunk/psphot/src/psphotSourceSize.c	(revision 29548)
@@ -1,6 +1,4 @@
 # include "psphotInternal.h"
 # include <gsl/gsl_sf_gamma.h>
-
-# define KRON 1
 
 typedef struct {
@@ -16,15 +14,23 @@
     int grow;
     int xtest, ytest;
-    bool apply; // apply CR mask?
+    bool applyCRmask; // apply CR mask?
+    bool dynamicLimitsCR; // apply CR mask?
+    float sizeLimitCR;
+    float magLimitCR;
 } psphotSourceSizeOptions;
 
 // local functions:
-bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
+bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe);
+bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe);
 bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
-bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
+bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
 bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
 bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
 int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
+float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction);
+
+// save test output images?
+# define DUMPPICS 0
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -111,5 +117,5 @@
     assert (status);
 
-    // XXX recipe name is not great
+    // location of a single test source
     options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
     options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
@@ -128,5 +134,5 @@
     }
 
-    options.apply = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
+    options.applyCRmask = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
     if (!status) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined.");
@@ -134,21 +140,16 @@
     }
 
-    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
-    // of an object.  We need to model this distribution for the PSF stars before we can test
-    // the significance for a specific object
-    // XXX move this to the code that generates the PSF?
-    // XXX store the results on pmPSF?
-
-    // XXX this should only be done on the first pass (ie, if we have newSources or allSources?)
-    if (getPSFsize) {
-        psphotSourceSizePSF (&options, sources, psf);
-    }
-
-    // classify the sources based on ApResid and Moments (extended sources)
+    // determine the distribution of (PSF_mag - KRON_mag) for the PSF sources (saved on readout->analysis)
+    psphotSourceSizePSF (&options, readout, sources, psf, recipe);
+
+    // adjust the user-supplied limits based on the distribution of CRs in the (Mminor, mKron) space
+    psphotDynamicLimitsCR(&options, readout, sources, psf, recipe);
+
+    // classify the sources based on ApResid and Moments
     // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
     psphotSourceClass(readout, sources, recipe, psf, &options);
 
-    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
-    psphotSourceSizeCR (readout, sources, &options);
+    // attempt to mask the candidate CRs; flag if CR nature is confirmed
+    psphotSourceSelectCR(readout, sources, &options);
 
     // XXX fix this (was source->n  - first)
@@ -156,6 +157,7 @@
 
     psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
+    psphotVisualPlotSourceSizeAlt (recipe, readout->analysis, sources);
     psphotVisualShowSourceSize (readout, sources);
-    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
+    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr, false);
     psphotVisualShowSatStars (recipe, psf, sources);
 
@@ -163,9 +165,33 @@
 }
 
+# define SAVE_PSF_OPTIONS(RO,OPT)					\
+    /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */ \
+    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID", PS_META_REPLACE, "locus of PSF stars in PSF_MAG - KRON_MAG", OPT->ApResid); \
+    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.PSF.APRESID.SYSERR",  PS_META_REPLACE, "systematic error of PSF_MAG - KRON_MAG",  OPT->ApSysErr);
+
 // model the apmifit distribution for the psf stars:
-bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
+bool psphotSourceSizePSF (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) {
+
+    // We are using the value PSF_MAG - KRON_MAG as a measure of the extendedness of an object.
+    // We need to model this distribution for the PSF stars before we can test the significance
+    // for a specific object.
+
+    // NOTE: we require that objects have had moments measured, and we also require psfMags to
+    // be calculated.  but, we do not require aperture mags or any other photometry values that
+    // require pixel analysis.
+
+    bool status1 = false;
+    bool status2 = false;
+    options->ApResid = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.PSF.APRESID");
+    options->ApSysErr = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.PSF.APRESID.SYSERR");
+    psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of PSF ApResid values");
+
+    // if they are saved on readout->analysis, we have already calculated these values
+    if (status1 && status2) {
+	return true;
+    }
 
     // select stats from the psf stars
-    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *ApOff = psVectorAllocEmpty (100, PS_TYPE_F32);
     psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
 
@@ -173,6 +199,6 @@
     psImageMaskType maskVal = options->maskVal | options->markVal;
 
-    // XXX  why PHOT_WEIGHT??
-    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
+    // with PSFONLY, we do not need modify the pixels
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
 
     int num = 0;                        // Number of sources measured
@@ -182,61 +208,182 @@
         num++;
 
-        // replace object in image
-        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
-        }
-
-        // clear the mask bit and set the circular mask pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
-
-        // XXX can we test if psfMag is set and calculate only if needed?
         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
 
-        // clear the mask bit
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
-
-        // re-subtract the object, leave local sky
-        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
-
-	// XXX test: switch to kron flux
-# if (KRON)
-        float apMag = -2.5*log10(source->moments->KronFlux);
-# else
-        float apMag = -2.5*log10(source->moments->Sum);
-# endif
-        float dMag = source->psfMag - apMag;
-
-        psVectorAppend (Ap, dMag);
+        float kMag = -2.5*log10(source->moments->KronFlux);
+        float dMag = source->psfMag - kMag;
+
+        psVectorAppend (ApOff, dMag);
         psVectorAppend (ApErr, source->errMag);
     }
     if (num == 0) {
-        // Not raising an error, because errors aren't being checked elsewhere in this function
-        psFree(Ap);
+	// if we cannot determine the PSF distribution, call all objects PSFs...
+	options->ApResid = NAN;
+	options->ApSysErr = NAN;
+        psFree(ApOff);
         psFree(ApErr);
+	SAVE_PSF_OPTIONS(readout, options);
         return false;
     }
 
     // model the distribution as a mean or median value and a systematic error from that value:
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-    psVectorStats (stats, Ap, NULL, NULL, 0);
-
-    psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
-    for (int i = 0; i < Ap->n; i++) {
-        dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
-    }
-
-    options->ApResid = stats->robustMedian;
-    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
-    // XXX this is quite arbitrary...
-    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
+    // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);    
+    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN);    
+    psVectorStats (stats, ApOff, NULL, NULL, 0);
+
+    psVector *dAp = psVectorAlloc (ApOff->n, PS_TYPE_F32);
+    for (int i = 0; i < ApOff->n; i++) {
+        dAp->data.F32[i] = ApOff->data.F32[i] - stats->clippedMean;
+    }
+
+    options->ApResid = stats->clippedMean;
+    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.1);
+
+    // this is quite arbitrary... a large value means fewer things classified as extended.
+    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.05;
     psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
 
-    psFree (Ap);
+    psFree (ApOff);
     psFree (ApErr);
     psFree (stats);
     psFree (dAp);
 
+    SAVE_PSF_OPTIONS(readout, options);
     return true;
+}
+
+# define SAVE_CR_OPTIONS(RO,OPT)					\
+    /* save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass */  \
+    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.SIZE", PS_META_REPLACE, "dynamically-set minor axis size limit for cosmic rays", OPT->sizeLimitCR); \
+    psMetadataAddF32(RO->analysis, PS_LIST_TAIL, "PSPHOT.CR.MAX.MAG",  PS_META_REPLACE, "dynamically-set kron magnitude limit for cosmic rays",  OPT->magLimitCR);
+
+// model the size and magnitude distribution of the Cosmic Rays
+// ** CRs are reliably flagged by a combination on Mminor < X && mag (or flux) > Y
+bool psphotDynamicLimitsCR (psphotSourceSizeOptions *options, pmReadout *readout, psArray *sources, pmPSF *psf, psMetadata *recipe) {
+
+    /* attempt to describe the CR sources:
+       - input parameters are sizeLimit, magLimit
+       - first try to refine the sizeLimit:
+       -- select objects which meet the magLimit and exceed the sizeLimit by a factor of 1.5
+       -- generate the histogram
+       -- look for a peak in the histogram
+       -- look for the min between valley and upper limit
+       -- look for first bin within 5% of the valley floor after peak (new sizeLimit)
+       
+       - next try to refine the magLimit
+       -- select objects which meet the sizeLimit and go fainter than the magLimit by 1.0 mag
+       -- generate the histogram
+       -- look for a peak in the histogram
+       -- look for the min between valley and upper limit
+       -- look for first bin within 5% of the valley floor after peak (new magLimit)
+    */
+
+    bool status  = false;
+    bool status1 = false;
+    bool status2 = false;
+    options->sizeLimitCR = psMetadataLookupF32 (&status1, readout->analysis, "PSPHOT.CR.MAX.SIZE");
+    if (!status1) {
+	options->sizeLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.SIZE");
+	if (!status) {
+	    options->sizeLimitCR = 1.0;
+	}
+    } 
+    options->magLimitCR = psMetadataLookupF32 (&status2, readout->analysis, "PSPHOT.CR.MAX.MAG");
+    if (!status2) {
+	options->magLimitCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.MAX.MAG");
+	if (!status) {
+	    options->magLimitCR = -8.0;
+	}
+    }
+    psAssert ((status1 && status2) || (!status1 && !status2), "inconsistent record of dynamic CR limits");
+
+    // if they are saved on readout->analysis, we have already calculated these values
+    if (status1 && status2) {
+	return true;
+    }
+
+    // if we do not want to dynamically set these, save the user-supplied values and exit
+    options->dynamicLimitsCR = psMetadataLookupBool (&status, recipe, "PSPHOT.CR.AUTOSCALE");
+    if (!status) {
+	options->dynamicLimitsCR = true;
+    }
+    if (!options->dynamicLimitsCR) {
+	SAVE_CR_OPTIONS(readout, options); // macro defined above
+	return true;
+    }
+
+    psVector *minor = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *mKron = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *mask  = NULL;
+
+    psImageMaskType markVal = options->markVal;
+    psImageMaskType maskVal = options->maskVal | options->markVal;
+
+    // with PSFONLY, we do not need modify the pixels
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
+
+    // generate vectors for all the objects of possible interest (so we can just access those
+    // vectors in the next sections)
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        // XXX can we test if psfMag is set and calculate only if needed?
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
+
+        // convert to Mmaj, Mmin:
+        psF32 Mxx = source->moments->Mxx;
+        psF32 Myy = source->moments->Myy;
+        psF32 Mxy = source->moments->Mxy;
+
+        float KronMag = -2.5*log10(source->moments->KronFlux);
+
+	float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
+
+	if (Mminor > options->sizeLimitCR * 1.5) continue;
+	if (KronMag > options->magLimitCR + 2.5) continue;
+
+        psVectorAppend (mKron, KronMag);
+        psVectorAppend (minor, Mminor);
+    }
+
+    // if too few objects meet the criterion, give up..
+    if (mKron->n < 50) goto escape;
+
+    // set distinct masks for (minor > sizeLimit) or (mKron > magLimit)
+    mask = psVectorAlloc(mKron->n, PS_TYPE_U8);
+    psVectorInit(mask, 0);
+    for (int i = 0; i < mKron->n; i++) {
+	if (mKron->data.F32[i] > options->magLimitCR) {
+	    mask->data.U8[i] |= 0x01;
+	}
+	if (minor->data.F32[i] > options->sizeLimitCR) {
+	    mask->data.U8[i] |= 0x02;
+	}
+    }
+	
+    float delta1 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 0.5 / minor->n));
+    float newSizeLimit = psphotSourceSizeFindThreshold(minor, mask, 0x01, 0.0, options->sizeLimitCR * 1.5, delta1, options->sizeLimitCR, 0.05);
+    if (isfinite(newSizeLimit)) {
+	options->sizeLimitCR = newSizeLimit;
+    }
+
+    float delta2 = PS_MAX(0.02, PS_MIN(0.2, 5.0 * 2.0 / minor->n));
+    float newMagLimit = psphotSourceSizeFindThreshold(mKron, mask, 0x02, -15.0, options->magLimitCR + 2.5, delta2, options->magLimitCR, 0.05);
+    if (isfinite(newMagLimit)) {
+	options->magLimitCR = newMagLimit;
+    }
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "CR limits : %f mag | %f pix^2\n", options->magLimitCR, options->sizeLimitCR);
+
+    psFree (mKron);
+    psFree (minor);
+    psFree (mask);
+
+    // save these on readout->analysis (a) for output to the header and (b) to prevent re-calculating in another pass
+    SAVE_CR_OPTIONS(readout, options); // macro defined above
+    return true;
+
+ escape:
+    SAVE_CR_OPTIONS(readout, options); // macro defined above
+    return false;
 }
 
@@ -248,5 +395,12 @@
     char regionName[64];
 
-    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
+    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nskip");
+
+    if (!psphotSourceClassRegion (NULL, &psfClump, sources, recipe, psf, options)) {
+	psLogMsg ("psphot", 4, "Failed to determine source classification for full image\n");
+    } else {
+	psLogMsg ("psphot", 4, "source classification for full image\n");
+    }
+    return true;
 
     int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
@@ -274,4 +428,5 @@
             continue;
         }
+	psLogMsg ("psphot", 4, "source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
         // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     }
@@ -280,5 +435,4 @@
 }
 
-# define SIZE_SN_LIM 10
 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
 
@@ -290,9 +444,10 @@
     int Npsf  = 0;
     int Ncr   = 0;
-    int Nmiss = 0;
     int Nskip = 0;
 
     pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
-    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
+
+    // request the pixWeight values as well as the magnitudes
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; 
 
     psImageMaskType markVal = options->markVal;
@@ -305,8 +460,10 @@
         // psfClumps are found for image subregions:
         // skip sources not in this region
-        if (source->peak->x <  region->x0) continue;
-        if (source->peak->x >= region->x1) continue;
-        if (source->peak->y <  region->y0) continue;
-        if (source->peak->y >= region->y1) continue;
+	if (region) {
+	    if (source->peak->x <  region->x0) continue;
+	    if (source->peak->x >= region->x1) continue;
+	    if (source->peak->y <  region->y0) continue;
+	    if (source->peak->y >= region->y1) continue;
+	}
 
         // skip source if it was already measured
@@ -320,8 +477,9 @@
             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
             psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
-            continue;
-        }
-
-        // we are basically classifying by moments
+            Nskip ++;
+            continue;
+        }
+
+        // we are classifying by moments and PSF_MAG - KRON_MAG
         psAssert (source->moments, "why is this source missing moments?");
         if (source->mode & noMoments) {
@@ -333,6 +491,8 @@
         psF32 Mxx = source->moments->Mxx;
         psF32 Myy = source->moments->Myy;
-
-        // replace object in image
+        psF32 Mxy = source->moments->Mxy;
+	float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
+
+        // replace object in image to measure mag & psfWeights
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
             pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
@@ -343,5 +503,4 @@
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
 
-        // XXX can we test if psfMag is set and calculate only if needed?
         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
 
@@ -352,35 +511,16 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
 
-# if (KRON)
-        float apMag = -2.5*log10(source->moments->KronFlux);
-# else
-        float apMag = -2.5*log10(source->moments->Sum);
-# endif
-        float dMag = source->psfMag - apMag;
-
-        // set nSigma to include both systematic and poisson error terms
-        // XXX the 'poisson error' contribution for size is probably wrong...
-	// XXX add in a hard floor on the Ap Sys Err (to be a bit generous)
-        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025));
-        float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
-        float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
-
-	// fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f, nSigmaMxx: %f, nSigmaMyy: %f\n", 
-	// source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My, 
-	// source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY);
-
-        // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs
-        float sizeBias = 1.0;
-        if (source->pixWeightNotBad < 0.9) {
-            sizeBias = 3.0;
-        }
-        if (source->pixWeightNotPoor < 0.9) {
-            sizeBias = 3.0;
-        }
-
-        float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
-        float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
-
-        // include MAG, MXX, and MYY?
+        float kMag = -2.5*log10(source->moments->KronFlux);
+        float dMag = source->psfMag - kMag;
+
+	// sources without a valid magnitudes cannot have their size measured
+	if (!isfinite(kMag) || !isfinite(source->psfMag)) {
+            Nskip ++;
+            continue;
+	}
+
+        // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
+	// floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
+        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02));
         source->extNsigma = nSigmaMAG;
 
@@ -389,16 +529,4 @@
         // * CR & defect should have a faintess limit (min S/N)
         // * SAT stars should not be faint, but defects may?
-
-        // Anything within this region is a probably PSF-like object. Saturated stars may land
-        // in this region, but are detected elsewhere on the basis of their peak value.
-        bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
-        if (isPSF) {
-          psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
-                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-          source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
-          Npsf ++;
-          continue;
-        }
 
         // Defects may not always match CRs from peak curvature analysis
@@ -407,33 +535,45 @@
         // XXX only accept brightish detections as CRs
         // (nSigmaMAG < -options->nSigmaApResid) ||
-        bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
+
+        // saturated star (too many saturated pixels or peak above saturation limit).  These
+        // may also be saturated galaxies, or just large saturated regions.
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g SAT\n",
+                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
+            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+            Nsat ++;
+            continue;
+        }
+
+        // any sources missing a large fraction should just be treated as PSFs
+        if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
+            psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
+                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,kMag,dMag,nSigmaMAG,
+                    options->nSigmaApResid,options->nSigmaMoments);
+	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+	    Npsf ++;
+	    continue;
+        }
+
+	// CRs are flagged by a combination on Mminor < options->sizeLimitCR && kmag < options->magLimitCR
+	// NOTE: we only flag the CRs here; when we mask them we verify their CR nature (otherwise -> PSF)
+        // bool isCR = (kMag < options->magLimitCR) && (Mminor < options->sizeLimitCR);
+	// XXX skip if we have already marked it??
+        bool isCR = (source->moments->SN > 7.0) && (Mminor < options->sizeLimitCR);
         if (isCR) {
-            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
-                    source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g CR\n",
+                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
             source->mode |= PM_SOURCE_MODE_DEFECT;
             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
+	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
             Ncr ++;
             continue;
         }
 
-        // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
-        // just large saturated regions.
-        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
-                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
-            Nsat ++;
-            continue;
-        }
-
-        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
-        bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
+        // Likely extended source (PSF_MAG - KRON_MAG is larger than limit)
+        bool isEXT = (nSigmaMAG > options->nSigmaApResid);
         if (isEXT) {
-          psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
-                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-
+            psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g EXT\n",
+                    source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
@@ -441,15 +581,13 @@
             continue;
         }
-        psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
-                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-                options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-
-        // sources that reach here are probably too faint for a reasonable source size measurement
-        // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
-        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
-        Nmiss ++;
-    }
-
-    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
+
+        // Everything else should just be treated as a PSF
+	psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g PSF\n",
+		source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
+	source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+	Npsf ++;
+    }
+
+    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nskip);
 
     return true;
@@ -458,5 +596,5 @@
 // given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
 // if enough of the object is detected as 'sharp', consider the object a cosmic ray
-bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
+bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
 
     psTimerStart ("psphot.cr");
@@ -466,10 +604,4 @@
         pmSource *source = sources->data[i];
 
-        // skip source if it was already measured
-        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
-            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
-            continue;
-        }
-
         // only check candidates marked above
         if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
@@ -478,7 +610,6 @@
         }
 
-        // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
-        // XXX this may be degenerate with the above test
-        if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
+	// once we are here, remove the temporary flag (this allows us to try again if we do not mark it now)
+	source->tmpFlags &= ~PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
 
         // Integer position of peak
@@ -504,5 +635,5 @@
         // XXX this is running slowly and is too agressive, but it more-or-less works
         psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
-        if (options->apply) {
+        if (options->applyCRmask) {
             psphotMaskCosmicRay(readout, source, options->crMask);
         } else {
@@ -531,12 +662,11 @@
 
     // XXX test : save the mask image
-    if (0) {
-        psphotSaveImage (NULL, readout->mask,   "mask.fits");
-    }
+# if (DUMPPICS)
+    psphotSaveImage (NULL, readout->mask,   "crmask.fits");
+# endif
 
     return true;
 }
 
-# define DUMPPICS 0
 # define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
 # define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
@@ -721,5 +851,4 @@
     if (nCRpix > 1) {
         source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     }
     // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
@@ -1216,2 +1345,64 @@
 }
 
+float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction) {
+
+    // make a histogram of the sources with mKron < magLimit
+    int nValue = (int)((maxValue - minValue) / delta) + 1;
+
+    psVector *histogram = psVectorAlloc(nValue, PS_TYPE_S32);
+    psVectorInit (histogram, 0);
+
+    for (int i = 0; i < value->n; i++) {
+	if (mask->data.U8[i] & maskValue) continue;
+	int bin = (value->data.F32[i] - minValue) / delta;
+	if (bin < 0.0) continue;
+	if (bin > histogram->n - 1) continue;
+	histogram->data.S32[bin] ++;
+    }
+
+    // find the peak of this histogram, but stop search at guess
+    int last = (guess - minValue) / delta;
+    int nPeak = 0;
+    int vPeak = histogram->data.S32[0];
+    for (int i = 1; (i < last) && (i < histogram->n); i++) {
+	if (histogram->data.S32[i] < vPeak) continue;
+	nPeak = i;
+	vPeak = histogram->data.S32[i];
+    }
+
+    // start at the peak and find the valley between here and the end of the histogram
+    int nValley = nPeak;
+    int vValley = histogram->data.S32[nPeak];
+    for (int i = nPeak + 1; i < histogram->n; i++) {
+	if (histogram->data.S32[i] > vValley) continue;
+	nValley = i;
+	vValley = histogram->data.S32[i];
+    }
+
+    psLogMsg ("psphot", PS_LOG_MINUTIA, "CR limits threshold : peak %d @ %f, valley %d @ %f (%f sigma)\n", 
+	      vPeak, delta * nPeak + minValue, vValley, delta * nValley + minValue, vPeak / PS_MAX(1.0, sqrt(vValley)));
+
+    if (nValley == nPeak) {
+	psFree (histogram);
+	return NAN;
+    }
+
+    if (vPeak < 3.0*sqrt(vValley)) {
+	psFree (histogram);
+	return NAN;
+    }
+
+    /// search for the first bin after the peak with f < vValley + 0.05*(vPeak - vValley)
+    float vLimit = vValley + fraction*(vPeak - vValley);
+    int nLimit = nValley;
+    for (int i = nPeak; i < nValley; i++) {
+	if (histogram->data.S32[i] > vLimit) continue;
+	nLimit = i;
+	break;
+    }
+    float result = nLimit * delta + minValue;
+
+    psFree (histogram);
+    
+    return result;
+}
Index: /trunk/psphot/src/psphotSourceStats.c
===================================================================
--- /trunk/psphot/src/psphotSourceStats.c	(revision 29547)
+++ /trunk/psphot/src/psphotSourceStats.c	(revision 29548)
@@ -45,4 +45,12 @@
     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     psAssert (readout, "missing readout?");
+
+    if (psTraceGetLevel("psphot") > 5) {
+	static int pass = 0;
+        char name[64];
+        sprintf (name, "srstats.v%d.fits", pass);
+        psphotSaveImage(NULL, readout->image, name);
+	pass ++;
+    }
 
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
@@ -418,5 +426,5 @@
 	    
 	    float apMag = NAN;
-	    pmSourcePhotometryAper (&apMag, NULL, source->pixels, source->maskObj, maskVal);
+	    pmSourcePhotometryAper (&apMag, NULL, NULL, NULL, source->pixels, source->variance, source->maskObj, maskVal);
 	    fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux));
 
Index: /trunk/psphot/src/psphotStackMatchPSFs.c
===================================================================
--- /trunk/psphot/src/psphotStackMatchPSFs.c	(revision 29547)
+++ /trunk/psphot/src/psphotStackMatchPSFs.c	(revision 29548)
@@ -85,4 +85,5 @@
 
     // set NAN pixels to 'SAT'
+    // XXX replace this is pmReadoutMaskInvalid?
     psImageMaskType maskVal = pmConfigMaskGet("SAT", config);
     if (!pmReadoutMaskNonfinite(readoutSrc, maskVal)) {
@@ -106,5 +107,5 @@
     rescaleData(readoutOut, config, options, index);
 
-    dumpImage(readoutOut, readoutSrc, index, "convolved");
+    // dumpImage(readoutOut, readoutSrc, index, "convolved");
 
     return true;
Index: /trunk/psphot/src/psphotStackMatchPSFsUtils.c
===================================================================
--- /trunk/psphot/src/psphotStackMatchPSFsUtils.c	(revision 29547)
+++ /trunk/psphot/src/psphotStackMatchPSFsUtils.c	(revision 29548)
@@ -361,5 +361,9 @@
 	// Scale the input parameters
 	widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
-	if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, options->inputSeeing->data.F32[index], options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
+
+	// we need to register the FWHM values for use downstream 
+	pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing);
+
+	if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, scaleRef, scaleMin, scaleMax)) {
 	    psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
 	    goto escape;
Index: /trunk/psphot/src/psphotStackParseCamera.c
===================================================================
--- /trunk/psphot/src/psphotStackParseCamera.c	(revision 29547)
+++ /trunk/psphot/src/psphotStackParseCamera.c	(revision 29548)
@@ -91,5 +91,5 @@
 
 	if (!rawInputFile && !cnvInputFile) {
-	    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s (%d) lacks IMAGE of type STR", item->name, i);
+	    psError(PSPHOT_ERR_CONFIG, true, "Component %s (%d) lacks both RAW:IMAGE and CNV:IMAGE of type STR", item->name, i);
 	    return false;
 	}
Index: /trunk/psphot/src/psphotStackReadout.c
===================================================================
--- /trunk/psphot/src/psphotStackReadout.c	(revision 29547)
+++ /trunk/psphot/src/psphotStackReadout.c	(revision 29548)
@@ -154,7 +154,12 @@
 
     psphotStackObjectsUnifyPosition (objects);
+
+    // measure circular, radial apertures (objects sorted by S/N)
     psphotRadialAperturesByObject (config, objects, view, STACK_OUT); 
 
+    // measure elliptical apertures, petrosians (objects sorted by S/N)
     psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
+
+    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
     psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
 
@@ -162,5 +167,5 @@
     psphotMagnitudes(config, view, STACK_OUT);
 
-    if (!psphotEfficiency(config, view, STACK_OUT)) {
+    if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();
Index: /trunk/psphot/src/psphotVisual.c
===================================================================
--- /trunk/psphot/src/psphotVisual.c	(revision 29547)
+++ /trunk/psphot/src/psphotVisual.c	(revision 29548)
@@ -25,4 +25,13 @@
 static int kapa2 = -1;
 static int kapa3 = -1;
+
+/** destroy windows at the end of a run*/
+bool psphotVisualClose(void)
+{
+    if(kapa1 != -1) KiiClose(kapa1);
+    if(kapa2 != -1) KiiClose(kapa2);
+    if(kapa3 != -1) KiiClose(kapa3);
+    return true;
+}
 
 int psphotKapaChannel (int channel) {
@@ -462,5 +471,5 @@
     if (myKapa == -1) return false;
 
-    KapaClearPlots (myKapa);
+    KapaClearSections (myKapa);
     KapaInitGraph (&graphdata);
     KapaSetFont (myKapa, "courier", 14);
@@ -1389,5 +1398,5 @@
     section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
 
-    KapaClearPlots (myKapa);
+    KapaClearSections (myKapa);
     // first section : mag vs CR nSigma
     section.dx = 1.0;
@@ -1621,5 +1630,5 @@
     if (myKapa == -1) return false;
 
-    KapaClearPlots (myKapa);
+    KapaClearSections (myKapa);
     KapaInitGraph (&graphdata);
     KapaSetFont (myKapa, "courier", 14);
@@ -2125,4 +2134,271 @@
 }
 
+bool PlotSourceSizeAltSetVectors(psVector *m, psVector *k, psVector *v1, psVector *v2, psVector *v3, pmSource *source) {
+
+    float Mxx = source->moments->Mxx;
+    float Myy = source->moments->Myy;
+    float Mxy = source->moments->Mxy;
+    float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
+    float KronMag = -2.5*log10(source->moments->KronFlux);
+    
+    psVectorAppend(m,  source->psfMag);
+    psVectorAppend(k,  KronMag);
+    psVectorAppend(v1, Mminor);
+    psVectorAppend(v2, source->psfMag - KronMag);
+    psVectorAppend(v3, source->extNsigma);
+    return true;
+}
+
+bool PlotSourceSizeAltAddPoints(Graphdata *graphdata, int myKapa, psVector *x, psVector *y, char *colorname, int ptype, float size) {
+
+    graphdata->color = KapaColorByName (colorname);
+    graphdata->ptype = ptype;
+    graphdata->size = size;
+    graphdata->style = 2;
+    KapaPrepPlot   (myKapa, x->n, graphdata);
+    KapaPlotVector (myKapa, x->n, x->data.F32, "x");
+    KapaPlotVector (myKapa, x->n, y->data.F32, "y");
+    return true;
+}
+
+bool psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
+
+    Graphdata graphdata;
+    KapaSection section;
+
+    if (!pmVisualTestLevel("psphot.size", 2)) return true;
+
+    int myKapa = psphotKapaChannel (2);
+    if (myKapa == -1) return false;
+
+    KapaClearSections (myKapa);
+    KapaInitGraph (&graphdata);
+    KapaSetFont (myKapa, "courier", 14);
+
+    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
+
+    // storage vectors for data to be plotted
+    psVector *SATm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *SATk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *SAT1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *SAT2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *SAT3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *PSFm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *PSFk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *PSF1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *PSF2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *PSF3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *EXTm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *EXTk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *EXT1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *EXT2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *EXT3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *DEFm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *DEFk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *DEF1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *DEF2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *DEF3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *BADm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *BADk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *BAD1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *BAD2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *BAD3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    psVector *CRm = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *CRk = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *CR1 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *CR2 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *CR3 = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+
+    // int notPSF = PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT;
+    int nSkip = 0;
+
+    // construct the vectors
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+	if (source->moments == NULL) {
+	    nSkip ++;
+	    continue;
+	}
+
+	bool found = false;
+
+	// only plot the measured sources...
+	if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) {
+	    nSkip ++;
+	    continue;
+	}
+
+        // any sources missing a large fraction should just be treated as PSFs
+        if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
+	    PlotSourceSizeAltSetVectors(BADm, BADk, BAD1, BAD2, BAD3, source);
+        }
+	if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
+	    PlotSourceSizeAltSetVectors(CRm, CRk, CR1, CR2, CR3, source);
+	    found = true;
+	}
+	if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+	    PlotSourceSizeAltSetVectors(SATm, SATk, SAT1, SAT2, SAT3, source);
+	    found = true;
+	}
+	if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+	    PlotSourceSizeAltSetVectors(EXTm, EXTk, EXT1, EXT2, EXT3, source);
+	    found = true;
+	}
+	if (source->mode & PM_SOURCE_MODE_DEFECT) {
+	    PlotSourceSizeAltSetVectors(DEFm, DEFk, DEF1, DEF2, DEF3, source);
+	    found = true;
+	}
+	if (!found) {
+	    PlotSourceSizeAltSetVectors(PSFm, PSFk, PSF1, PSF2, PSF3, source);
+	}
+	// if (!(source->mode & notPSF) && !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
+	//     PlotSourceSizeAltSetVectors(PSFm, PSFk, PSF1, PSF2, PSF3, source);
+	// }
+    }
+    // three sections: kronMag vs Mminor, psfMag vs psfMag - KronMag, psfMag vs extNsigma
+
+    // --- first section: kronMag vs Mminor ---
+    section.dx = 1.00;
+    section.dy = 0.33;
+    section.x  = 0.00;
+    section.y  = 0.00;
+    section.name = psStringCopy ("Mminor");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = -17.1;
+    graphdata.xmax =  -6.9;
+    graphdata.ymin = -0.5;
+    graphdata.ymax = +7.1;
+    KapaSetLimits (myKapa, &graphdata);
+
+    graphdata.padXm = NAN;
+    graphdata.padYm = 5.0;
+    graphdata.padXp = NAN;
+    graphdata.padYp = NAN;
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "kron mag", KAPA_LABEL_XM);
+    KapaSendLabel (myKapa, "M_minor| (pixels^2|)", KAPA_LABEL_YM);
+
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFk, PSF1, "black", 0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTk, EXT1, "blue",  0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa,  CRk,  CR1, "red",   0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFk, DEF1, "red",   7, 1.0);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATk, SAT1, "blue",  7, 1.2);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADk, BAD1, "green", 7, 1.4);
+
+    // --- second section: dMag ----
+    section.dx = 1.00;
+    section.dy = 0.33;
+    section.x  = 0.00;
+    section.y  = 0.33;
+    section.name = psStringCopy ("dMag");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = -17.1;
+    graphdata.xmax =  -6.9;
+    graphdata.ymin = -0.75;
+    graphdata.ymax = +1.50;
+    KapaSetLimits (myKapa, &graphdata);
+
+    graphdata.padXm = NAN;
+    graphdata.padYm = 5.0;
+    graphdata.padXp = 0.0;
+    graphdata.padYp = NAN;
+    strcpy (graphdata.labels, "0200");
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "dMag", KAPA_LABEL_YM);
+
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFm, PSF2, "black", 0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTm, EXT2, "blue",  0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa,  CRm,  CR2, "red",   0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFm, DEF2, "red",   7, 1.0);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATm, SAT2, "blue",  7, 1.2);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADm, BAD2, "green", 7, 1.4);
+
+    // --- third section: nSigma ---
+    section.dx = 1.00;
+    section.dy = 0.33;
+    section.x  = 0.00;
+    section.y  = 0.66;
+    section.name = psStringCopy ("nSigma");
+    KapaSetSection (myKapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = -17.1;
+    graphdata.xmax = -6.9;
+    graphdata.ymin = -10.1;
+    graphdata.ymax = +10.1;
+    KapaSetLimits (myKapa, &graphdata);
+
+    graphdata.padXm = 0.0;
+    graphdata.padYm = 5.0;
+    graphdata.padXp = NAN;
+    graphdata.padYp = NAN;
+    strcpy (graphdata.labels, "0210");
+    KapaBox (myKapa, &graphdata);
+    KapaSendLabel (myKapa, "psf msg", KAPA_LABEL_XP);
+    KapaSendLabel (myKapa, "EXT nSigma", KAPA_LABEL_YM);
+
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, PSFm, PSF3, "black", 0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, EXTm, EXT3, "blue",  0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa,  CRm,  CR3, "red",   0, 0.5);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, DEFm, DEF3, "red",   7, 1.0);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, SATm, SAT3, "blue",  7, 1.0);
+    PlotSourceSizeAltAddPoints(&graphdata, myKapa, BADm, BAD3, "green", 7, 1.4);
+
+    fprintf (stderr, "PSF: %ld, EXT: %ld, CR: %ld, DEF: %ld, SAT: %ld, TOTAL: %ld, nSkip: %d\n", 
+	     PSFm->n, EXTm->n, CRm->n, DEFm->n, SATm->n, 
+	     PSFm->n+ EXTm->n+ CRm->n+ DEFm->n+ SATm->n, nSkip);
+
+    psFree (SATk);
+    psFree (SATm);
+    psFree (SAT1);
+    psFree (SAT2);
+    psFree (SAT3);
+
+    psFree (EXTk);
+    psFree (EXTm);
+    psFree (EXT1);
+    psFree (EXT2);
+    psFree (EXT3);
+
+    psFree (PSFk);
+    psFree (PSFm);
+    psFree (PSF1);
+    psFree (PSF2);
+    psFree (PSF3);
+
+    psFree (DEFk);
+    psFree (DEFm);
+    psFree (DEF1);
+    psFree (DEF2);
+    psFree (DEF3);
+
+    psFree (BADk);
+    psFree (BADm);
+    psFree (BAD1);
+    psFree (BAD2);
+    psFree (BAD3);
+
+    psFree (CRk);
+    psFree (CRm);
+    psFree (CR1);
+    psFree (CR2);
+    psFree (CR3);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
 bool psphotVisualShowResidualImage (pmReadout *readout) {
 
@@ -2138,5 +2414,5 @@
 }
 
-bool psphotVisualPlotApResid (psArray *sources, float mean, float error) {
+bool psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag) {
 
     Graphdata graphdata;
@@ -2148,5 +2424,5 @@
     if (myKapa == -1) return false;
 
-    KapaClearPlots (myKapa);
+    KapaClearSections (myKapa);
     KapaInitGraph (&graphdata);
 
@@ -2169,6 +2445,14 @@
 	if (!isfinite (source->psfMag)) continue;
 
+	float dMag;
+	if (useApMag) {
+	    dMag = source->apMag - source->psfMag;
+	} else {
+	    float kMag = -2.5*log10(source->moments->KronFlux);
+	    dMag = source->psfMag - kMag;
+	}
+
 	x->data.F32[n] = source->psfMag;
-	y->data.F32[n] = source->apMag - source->psfMag;
+	y->data.F32[n] = dMag;
 	dy->data.F32[n] = source->errMag;
 	graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
@@ -2266,5 +2550,5 @@
     if (myKapa == -1) return false;
 
-    KapaClearPlots (myKapa);
+    KapaClearSections (myKapa);
     KapaInitGraph (&graphdata);
 
Index: /trunk/psphot/test/notes.txt
===================================================================
--- /trunk/psphot/test/notes.txt	(revision 29548)
+++ /trunk/psphot/test/notes.txt	(revision 29548)
@@ -0,0 +1,5 @@
+
+psphot tests
+
+This is a long-overdue set of psphot unit tests (or are these system tests?) 
+
Index: /trunk/psphot/test/tap.pro
===================================================================
--- /trunk/psphot/test/tap.pro	(revision 29548)
+++ /trunk/psphot/test/tap.pro	(revision 29548)
@@ -0,0 +1,52 @@
+# -*-sh-*-
+
+macro tapOK
+  if ($0 != 3) 
+    echo "USAGE: tapOK (condition) (message)"
+    break
+  end
+
+  if ($1)
+    echo "ok : $2"
+    $TAP_LAST = 1
+  else
+    echo "not ok : $2"
+    $TAP_NFAIL ++
+    $TAP_LAST = 0
+    if ($TAP_BREAK)
+     break
+    end
+  end
+  $TAP_NDONE ++
+end
+
+macro tapPLAN
+  if ($0 != 2) 
+    echo "USAGE: tapPLAN (Ntests)"
+    break
+  end
+
+  $TAP_NTEST = $1
+  $TAP_NFAIL = 0
+  $TAP_NDONE = 0
+  if (not($?TAP_BREAK)) set TAP_BREAK = 0
+end
+
+macro tapDONE
+  if ($0 != 1) 
+    echo "USAGE: tapDONE"
+    break
+  end
+
+  if ($TAP_NDONE != $TAP_NTEST) 
+    echo "planned tests ($TAP_NTEST) not equal to done tests ($TAP_NDONE)"
+  end
+
+  if ($TAP_NFAIL) 
+    echo "failed $TAP_NFAIL of $TAP_NDONE"
+  end
+
+  if ($TAP_NFAIL == 0) 
+    echo "passed $TAP_NDONE tests"
+  end
+end
Index: /trunk/psphot/test/tap_psphot_basic.pro
===================================================================
--- /trunk/psphot/test/tap_psphot_basic.pro	(revision 29548)
+++ /trunk/psphot/test/tap_psphot_basic.pro	(revision 29548)
@@ -0,0 +1,211 @@
+#!/usr/bin/env mana 
+# -*-sh-*-
+
+# macros to support tap tests in mana/dvo
+input tap.pro
+
+# create a fake image (ppSim, all stars are fake)
+macro make.fake
+
+  local options
+
+  # basic options for these image (r, 100sec)
+  $options = -type OBJECT -filter r -exptime 100.0
+  $options = $options -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 -seeing 1.0
+  $options = $options -D PSF.MODEL PS_MODEL_GAUSS
+  
+  $options = $options -Df STARS.DENSITY 20.0
+
+  $options = $options -camera SIMTEST -recipe PPSIM STACKTEST.MAKE
+  exec ppSim $options testimage
+end
+
+macro run.psphot
+  exec psphot -file testimage.fits testimage.phot
+end
+
+macro test.init
+  data testimage.phot.cmf
+
+  # a bit hackish : cannot append to non-existent scalar
+  $fields = Chip.psf
+  for i 0 $pairs:n
+    list word -split $pairs:$i
+    $fields = $fields $word:0
+  end
+
+  read -fits $fields
+
+  data testimage.dat
+  read Xraw 1 Yraw 2 Fraw 3 Galaxy 4 Mraw 5 dMraw 6 SXXraw 7 SYYraw 8 SXYraw 9 S7raw 10
+  set dFraw = Fraw * dMraw
+  set dXraw = 4.0 * dMraw / 2.35
+  set dYraw = 4.0 * dMraw / 2.35
+
+  match2d X_PSF Y_PSF Xraw Yraw 2.0 -index1 index1 -index2 index2
+end
+
+macro test.plots
+  if ($0 != 2)
+    echo "USAGE: test.plots (outroot)"
+    break
+  end
+
+  clear -s 
+
+  resize 1200 1175
+
+  $nx = 0
+  $ny = 0
+  $np = 0
+  for i 0 $pairs:n
+    section a$i {0.5*$nx} {0.8 - 0.2*$ny} 0.5 0.2
+    show.pair $i
+    $ny ++
+    if ($ny == 5) 
+      $nx ++
+      $ny = 0
+    end
+    if ($nx == 2)
+      png -name $1.$np.png
+      clear -s
+      $np ++
+      $nx = 0
+    end
+  end  
+  png -name $1.$np.png
+end
+
+macro show.pair
+  if ($0 != 2)
+    echo "USAGE: show.pair (N)"
+    break
+  end
+
+  list word -split $pairs:$1
+  if ($word:n != 7)
+    echo "invalid pair $1"
+    break
+  end
+
+  $Nr = $word:3
+
+  reindex v1 = $word:0 using index1
+  if ("$word:6" == "V") 
+    reindex v2 = $word:1 using index2
+  end
+  if ("$word:6" == "S") 
+    set v2 = $word:1 + zero(index1)
+  end
+  reindex rv = $word:2 using index$Nr
+
+  set delta = v1 - v2
+  if (("$word:4" == "def") || ("$word:5" == "def"))
+    lim rv delta; box; plot rv delta
+  else
+    lim rv $word:4 $word:5; box; plot rv delta
+  end
+  label -y $word:0 -x $word:3
+end
+
+macro show.dpair
+  if ($0 != 2)
+    echo "USAGE: show.dpair (N)"
+    break
+  end
+
+  list word -split $dpairs:$1
+  if ($word:n != 7)
+    echo "invalid dpair $1"
+    break
+  end
+
+  $Nr = $word:4
+
+  reindex v1 = $word:0 using index1
+  reindex dv = $word:1 using index1
+  reindex v2 = $word:2 using index2
+  reindex rv = $word:3 using index$Nr
+
+  set delta = (v1 - v2) / dv
+  if (("$word:5" == "def") || ("$word:6" == "def"))
+    lim rv delta; box; plot rv delta
+  else
+    lim rv $word:5 $word:6; box; plot rv delta
+  end
+  label -y $word:0 -x $word:3
+end
+
+macro show.pair.old
+  set delta = $1 - $2
+  lim mag2 delta; clear; box; plot mag2 delta
+end
+
+# list the entries that should be compared and the vector to use as a range
+list dpairs
+  X_PSF      	    X_PSF_SIG 	      Xraw Mraw 2 -10.0 10.0
+  Y_PSF 	    Y_PSF_SIG 	      Yraw Mraw 2 -10.0 10.0
+  PSF_INST_MAG 	    PSF_INST_MAG_SIG  Mraw Mraw 2 -10.0 10.0
+  PSF_INST_FLUX	    PSF_INST_FLUX_SIG Fraw Mraw 2 -10.0 10.0
+  AP_MAG  	    PSF_INST_MAG_SIG  Mraw Mraw 2 -10.0 10.0
+  AP_MAG_RAW	    PSF_INST_MAG_SIG  Mraw Mraw 2 -10.0 10.0
+end
+
+# list the entries that should be compared and the vector to use as a range
+list pairs
+  X_PSF      	     Xraw Mraw 2 -1.01 1.01 V
+  Y_PSF 	     Yraw Mraw 2 -1.01 1.01 V
+  X_PSF_SIG 	    dXraw Mraw 2 -1.01 1.01 V
+  Y_PSF_SIG 	    dYraw Mraw 2 -1.01 1.01 V
+  PSF_INST_MAG 	     Mraw Mraw 2 -1.01 1.01 V
+  PSF_INST_MAG_SIG    0.0 Mraw 2 -0.01 1.01 S
+  PSF_INST_FLUX	     Fraw Mraw 2   def  def V
+  PSF_INST_FLUX_SIG dFraw Mraw 2   def  def V
+  AP_MAG  	     Mraw Mraw 2 -1.01 1.01 V
+  AP_MAG_RAW	     Mraw Mraw 2 -1.01 1.01 V
+  AP_MAG_RADIUS      0.0  Mraw 2 -0.01 20.1 S
+  SKY     	     0.0  Mraw 2   def  def S
+  SKY_SIGMA	     0.0  Mraw 2   def  def S
+  PSF_CHISQ	     1.0  Mraw 2   def  def S
+  CR_NSIGMA	     0.0  Mraw 2   def  def S
+  EXT_NSIGMA	     0.0  Mraw 2 -5.01 5.01 S
+  PSF_MAJOR	     0.0  Mraw 2 -0.01 5.01 S
+  PSF_MINOR	     0.0  Mraw 2 -0.01 5.01 S
+  PSF_THETA	     0.0  Mraw 2 -1.61 1.61 S
+  PSF_QF  	     0.0  Mraw 2 -0.10 1.10 S
+  PSF_QF_PERFECT     0.0  Mraw 2 -0.10 1.10 S
+  PSF_NDOF	     0.0  Mraw 2  def  def  S
+  PSF_NPIX	     0.0  Mraw 2  def  def  S
+  MOMENTS_XX	     0.0  Mraw 2 -0.01 3.01 S
+  MOMENTS_XY	     0.0  Mraw 2 -3.01 3.01 S
+  MOMENTS_YY	     0.0  Mraw 2 -0.01 3.01 S
+  MOMENTS_M3C	     0.0  Mraw 2 -3.01 3.01 S
+  MOMENTS_M3S	     0.0  Mraw 2 -3.01 3.01 S
+  MOMENTS_M4C	     0.0  Mraw 2 -2.01 2.01 S
+  MOMENTS_M4S	     0.0  Mraw 2 -2.01 2.01 S
+  MOMENTS_R1	     0.0  Mraw 2 -5.01 5.01 S
+  MOMENTS_RH	     0.0  Mraw 2 -5.01 5.01 S
+  KRON_FLUX	     0.0  Mraw 2  def  def  S
+  KRON_FLUX_ERR	     0.0  Mraw 2  def  def  S
+  KRON_FLUX_INNER    0.0  Mraw 2  def  def  S
+  KRON_FLUX_OUTER    0.0  Mraw 2  def  def  S
+#   CAL_PSF_MAG	     none Mraw 2 -1.01 1.01 V
+#   CAL_PSF_MAG_SIG  none Mraw 2 -1.01 1.01 V
+#   RA_PSF  	     none Mraw 2 -1.01 1.01 V
+#   DEC_PSF 	     none Mraw 2 -1.01 1.01 V
+#   PEAK_FLUX_AS_MAG none Mraw 2 -1.01 1.01 V
+#  FLAGS   	     0.0  Mraw 2 -1.01 1.01 S
+#  FLAGS2             0.0  Mraw 2 -1.01 1.01 S
+end
+
+if ($SCRIPT) 
+  make.fake
+  exit 0
+end
+
+
+  # write tmp.dat X_PSF Y_PSF X_PSF_SIG Y_PSF_SIG PSF_INST_MAG PSF_INST_MAG_SIG
+  # exec gcompare -m tmp.dat 1 2 testimage.dat 1 2 2.0 > tmp.match.dat
+  # 
+  # data tmp.match.dat
+  # read x1 1 y1 2 dx1 3 dy1 4 psfmag1 5 dpsfmag1 6 x2 7 y2 8 flux 9 galaxy 10 mag2 11 dmag2 12 sxx 13 syy 14 sxy 15 s7 16