Index: trunk/psModules/src/objects/Makefile.am
===================================================================
--- trunk/psModules/src/objects/Makefile.am	(revision 32632)
+++ trunk/psModules/src/objects/Makefile.am	(revision 32633)
@@ -45,4 +45,5 @@
 	pmSourceIO_CMF_PS1_V2.c \
 	pmSourceIO_CMF_PS1_V3.c \
+	pmSourceIO_CMF_PS1_V4.c \
 	pmSourceIO_CMF_PS1_SV1.c \
 	pmSourceIO_CMF_PS1_DV1.c \
@@ -130,5 +131,5 @@
 
 # pmSourceID_CMF_* functions use a common framework
-BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c
+BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c pmSourceIO_CMF_PS1_V4.c
 
 pmSourceIO_CMF_PS1_V1.c : pmSourceIO_CMF.c.in mksource.pl
@@ -141,3 +142,6 @@
 	mksource.pl pmSourceIO_CMF.c.in PS1_V3 pmSourceIO_CMF_PS1_V3.c
 
+pmSourceIO_CMF_PS1_V4.c : pmSourceIO_CMF.c.in mksource.pl
+	mksource.pl pmSourceIO_CMF.c.in PS1_V4 pmSourceIO_CMF_PS1_V4.c
+
 # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in
Index: trunk/psModules/src/objects/mksource.pl
===================================================================
--- trunk/psModules/src/objects/mksource.pl	(revision 32632)
+++ trunk/psModules/src/objects/mksource.pl	(revision 32633)
@@ -16,5 +16,6 @@
 %cmfmodes = ("PS1_V1", 1,
 	     "PS1_V2", 2,
-	     "PS1_V3", 3);
+	     "PS1_V3", 3,
+	     "PS1_V4", 4);
 
 print "1: $cmfmodes{1}\n";
Index: trunk/psModules/src/objects/pmSource.c
===================================================================
--- trunk/psModules/src/objects/pmSource.c	(revision 32632)
+++ trunk/psModules/src/objects/pmSource.c	(revision 32633)
@@ -137,4 +137,8 @@
     source->apFlux    	     = NAN;
     source->apFluxErr 	     = NAN; 
+
+    source->skyRadius  	     = NAN;
+    source->skyFlux    	     = NAN;
+    source->skySlope   	     = NAN;
 
     source->pixWeightNotBad  = NAN;
Index: trunk/psModules/src/objects/pmSource.h
===================================================================
--- trunk/psModules/src/objects/pmSource.h	(revision 32632)
+++ trunk/psModules/src/objects/pmSource.h	(revision 32633)
@@ -36,6 +36,8 @@
     PM_SOURCE_TMPF_MOMENTS_MEASURED  = 0x0010,
     PM_SOURCE_TMPF_CANDIDATE_PSFSTAR = 0x0020,
-    PM_SOURCE_TMPF_STACK_KEEP        = 0x0040,
-    PM_SOURCE_TMPF_STACK_SKIP        = 0x0080,
+    PM_SOURCE_TMPF_RADIAL_KEEP       = 0x0040,
+    PM_SOURCE_TMPF_RADIAL_SKIP       = 0x0080,
+    PM_SOURCE_TMPF_PETRO_KEEP        = 0x0100,
+    PM_SOURCE_TMPF_PETRO_SKIP        = 0x0200,
 } pmSourceTmpF;
 
@@ -95,4 +97,8 @@
     float apFluxErr;                    ///< apFluxErr corresponding to psfMag or extMag (depending on type)
 
+    float skyRadius;			///< radius at which profile hits local sky (or goes flat)
+    float skyFlux;			///< mean flux per pixel in aperture at which profile hits local sky (or goes flat)
+    float skySlope;			///< mean flux slope at which profile hits local sky (or goes flat)
+
     float pixWeightNotBad;              ///< PSF-weighted coverage of unmasked (not BAD) pixels
     float pixWeightNotPoor;             ///< PSF-weighted coverage of unmasked (not POOR) pixels
Index: trunk/psModules/src/objects/pmSourceIO.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO.c	(revision 32632)
+++ trunk/psModules/src/objects/pmSourceIO.c	(revision 32633)
@@ -569,4 +569,5 @@
 	    PM_SOURCES_WRITE("PS1_V2",    CMF_PS1_V2);
 	    PM_SOURCES_WRITE("PS1_V3",    CMF_PS1_V3);
+	    PM_SOURCES_WRITE("PS1_V4",    CMF_PS1_V4);
 	    PM_SOURCES_WRITE("PS1_SV1",   CMF_PS1_SV1);
 	    PM_SOURCES_WRITE("PS1_DV1",   CMF_PS1_DV1);
@@ -1025,4 +1026,7 @@
                 sources = pmSourcesRead_CMF_PS1_V3 (file->fits, hdu->header);
             }
+            if (!strcmp (exttype, "PS1_V4")) {
+                sources = pmSourcesRead_CMF_PS1_V4 (file->fits, hdu->header);
+            }
             if (!strcmp (exttype, "PS1_SV1")) {
                 sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header);
Index: trunk/psModules/src/objects/pmSourceIO.h
===================================================================
--- trunk/psModules/src/objects/pmSourceIO.h	(revision 32632)
+++ trunk/psModules/src/objects/pmSourceIO.h	(revision 32633)
@@ -62,4 +62,9 @@
 bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 
+bool pmSourcesWrite_CMF_PS1_V4(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_V4_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_V4_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_V4_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
+
 bool pmSourcesWrite_CMF_PS1_SV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_CMF_PS1_SV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
@@ -86,4 +91,5 @@
 psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header);
+psArray *pmSourcesRead_CMF_PS1_V4 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_CMF.c.in	(revision 32632)
+++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in	(revision 32633)
@@ -127,5 +127,5 @@
 
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
         @<PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
@@ -138,5 +138,5 @@
         @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
 
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
@@ -150,5 +150,5 @@
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
@@ -158,26 +158,23 @@
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
 
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
-
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
-
-	// XXX do not keep this long term, just a TEST:
-        // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF",    PS_DATA_F32, "Kron Flux",                                  moments.KronPSF);
-        // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF_SIG",PS_DATA_F32, "Kron Flux",                                  moments.KronPSFErr);
-	// Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d
-        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX",   PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.KronCore);
-	// psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
+
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
+
+        @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
+        @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
+        @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
 
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
-        @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
+        @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
 
         // XXX not sure how to get this : need to load Nimages with weight?
@@ -286,9 +283,9 @@
         @ALL@     source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
-        @=PS1_V3@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
+        @>PS1_V2@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
 
         // XXX use these to determine PAR[PM_PAR_I0] if they exist?
-        @=PS1_V3@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
-        @=PS1_V3@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
+        @>PS1_V2@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
+        @>PS1_V2@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
@@ -311,5 +308,5 @@
 
         @ALL@     source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
-        @=PS1_V3@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
+        @>PS1_V2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
         @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
         @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
@@ -329,27 +326,31 @@
         @ALL@     source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
 
-        @=PS1_V3@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
-        @=PS1_V3@ source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
-        @=PS1_V3@ source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
-        @=PS1_V3@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
-
-        @=PS1_V3@ source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
-        @=PS1_V3@ source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
+        @>PS1_V2@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
+        @>PS1_V2@ source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
+        @>PS1_V2@ source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
+        @>PS1_V2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
+
+        @>PS1_V2@ source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
+        @>PS1_V2@ source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
+
+        @>PS1_V3@ source->skyRadius            = psMetadataLookupF32 (&status, row, "SKY_LIMIT_RAD");
+        @>PS1_V3@ source->skyFlux              = psMetadataLookupF32 (&status, row, "SKY_LIMIT_FLUX");
+        @>PS1_V3@ source->skySlope             = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE");
 
 	// XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
 	// we are storing enough information so the output will be consistent with the input
-        @=PS1_V3@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
-        @=PS1_V3@ source->moments->Mxxy = 0.0;
-        @=PS1_V3@ source->moments->Mxyy = 0.0;
-        @=PS1_V3@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
-
-        @=PS1_V3@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
-        @=PS1_V3@ source->moments->Mxxxy = 0.0;
-        @=PS1_V3@ source->moments->Mxxyy = 0.0;
-        @=PS1_V3@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
-        @=PS1_V3@ source->moments->Myyyy = 0.0;
+        @>PS1_V2@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
+        @>PS1_V2@ source->moments->Mxxy = 0.0;
+        @>PS1_V2@ source->moments->Mxyy = 0.0;
+        @>PS1_V2@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
+
+        @>PS1_V2@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
+        @>PS1_V2@ source->moments->Mxxxy = 0.0;
+        @>PS1_V2@ source->moments->Mxxyy = 0.0;
+        @>PS1_V2@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
+        @>PS1_V2@ source->moments->Myyyy = 0.0;
 
         @ALL@     source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
-        @=PS1_V3@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
+        @>PS1_V2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
         assert (status);
 
Index: trunk/psphot/src/Makefile.am
===================================================================
--- trunk/psphot/src/Makefile.am	(revision 32632)
+++ trunk/psphot/src/Makefile.am	(revision 32633)
@@ -185,4 +185,5 @@
 	psphotKronMasked.c             \
 	psphotKronIterate.c            \
+	psphotRadialProfileWings.c     \
 	psphotDeblendSatstars.c	       \
 	psphotMosaicSubimage.c	       \
@@ -201,5 +202,4 @@
 	psphotRadialBins.c	       \
 	psphotRadialApertures.c	       \
-	psphotRadialAperturesByObject.c \
 	psphotPetrosian.c	       \
         psphotPetrosianRadialBins.c    \
Index: trunk/psphot/src/psphot.h
===================================================================
--- trunk/psphot/src/psphot.h	(revision 32632)
+++ trunk/psphot/src/psphot.h	(revision 32633)
@@ -428,5 +428,6 @@
 bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry);
 bool psphotRadialApertures_Threaded (psThreadJob *job);
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVer);
+// bool psphotRadialApertureSource (pmSource *source, int entry);
 
 bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
@@ -466,4 +467,8 @@
 bool psphotKronIterate_Threaded (psThreadJob *job);
 
+bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources);
+bool psphotRadialProfileWings_Threaded (psThreadJob *job);
+
 bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
 
Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c
===================================================================
--- trunk/psphot/src/psphotExtendedSourceAnalysis.c	(revision 32632)
+++ trunk/psphot/src/psphotExtendedSourceAnalysis.c	(revision 32633)
@@ -209,6 +209,6 @@
 	// if we have checked the source validity on the basis of the object set, then 
 	// we either skip these tests below or we skip the source completely
-	if (source->tmpFlags & PM_SOURCE_TMPF_STACK_SKIP) continue;
-	if (source->tmpFlags & PM_SOURCE_TMPF_STACK_KEEP) goto keepSource;
+	if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+	if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
 	// skip PSF-like and non-astronomical objects
Index: trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- trunk/psphot/src/psphotFitSourcesLinear.c	(revision 32632)
+++ trunk/psphot/src/psphotFitSourcesLinear.c	(revision 32633)
@@ -134,4 +134,7 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+
+	// psAssert (source->peak, "source without peak??");
+	// psAssert (source->peak->footprint, "peak without footprint??");
 
         // turn this bit off and turn it on again if we pass this test
Index: trunk/psphot/src/psphotGuessModels.c
===================================================================
--- trunk/psphot/src/psphotGuessModels.c	(revision 32632)
+++ trunk/psphot/src/psphotGuessModels.c	(revision 32633)
@@ -170,4 +170,6 @@
         if (!source->peak) continue;
 
+	// psAssert (source->peak->footprint, "peak without footprint??");
+
         nSrc ++;
 
Index: trunk/psphot/src/psphotKronIterate.c
===================================================================
--- trunk/psphot/src/psphotKronIterate.c	(revision 32632)
+++ trunk/psphot/src/psphotKronIterate.c	(revision 32633)
@@ -151,5 +151,5 @@
             }
 # else
-	    if (!psphotExtendedSourceFits_Threaded(job)) {
+	    if (!psphotKronIterate_Threaded(job)) {
 		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
 		psFree(AnalysisRegion);
@@ -192,5 +192,6 @@
     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
 
-    for (int j = 0; j < 5; j++) {
+    // XXX TEST : set iteration to 1
+    for (int j = 0; j < 1; j++) {
 	for (int i = 0; i < sources->n; i++) {
 
@@ -209,13 +210,17 @@
 
 	    // use S/N to control max window size
-	    float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
+	    // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
 
 	    // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux)
 	    // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2)
 	    // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0)
-	    float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
+	    // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
 
 	    // iterate to the window radius
-	    float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
+	    // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
+
+	    // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
+	    float maxWindow = source->skyRadius;
+	    float windowRadius = PS_MAX(RADIUS, maxWindow);
 
 	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
Index: trunk/psphot/src/psphotRadialApertures.c
===================================================================
--- trunk/psphot/src/psphotRadialApertures.c	(revision 32632)
+++ trunk/psphot/src/psphotRadialApertures.c	(revision 32633)
@@ -26,6 +26,12 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
 	if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
@@ -35,4 +41,11 @@
     return true;
 }
+
+// these values are used by all threads repeatedly (and are not modified)
+static psVector *aperRadii = NULL;
+static psVector *aperRadii2 = NULL;
+static float outerRadius = NAN;
+static float SN_LIM = NAN;
+static psImageMaskType maskVal = 0;
 
 // aperture-like measurements for extended sources
@@ -75,4 +88,25 @@
         nThreads = 0;
     }
+
+    // aperRadii stores the upper bounds of the annuli
+    // XXX keep the same name here as for the petrosian / elliptical apertures?
+    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
+    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+
+    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
+
+    // save the R^2 values as well for quicker comparison
+    aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
+    for (int i = 0; i < aperRadii->n; i++) {
+	aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // S/N limit to perform full non-linear fits
+    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     // source analysis is done in S/N order (brightest first)
@@ -117,12 +151,9 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, AnalysisRegion);
-            psArrayAdd(job->args, 1, recipe);
-
             PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
-
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
 
-// set this to 0 to run without threading
+	    // set this to 0 to run without threading
 # if (1)	    
             if (!psThreadJobAddPending(job)) {
@@ -138,5 +169,5 @@
 	    }
 	    psScalar *scalar = NULL;
-	    scalar = job->args->data[6];
+	    scalar = job->args->data[5];
 	    Nradial += scalar->data.S32;
 	    psFree(job);
@@ -158,5 +189,5 @@
             } else {
 		psScalar *scalar = NULL;
-		scalar = job->args->data[6];
+		scalar = job->args->data[5];
 		Nradial += scalar->data.S32;
             }
@@ -166,4 +197,5 @@
     psFree (cellGroups);
     psFree(AnalysisRegion);
+    psFree (aperRadii2);
 
     psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
@@ -173,5 +205,4 @@
 bool psphotRadialApertures_Threaded (psThreadJob *job) {
 
-    bool status;
     int Nradial = 0;
 
@@ -180,23 +211,11 @@
     psArray *sources        = job->args->data[1];
     psRegion *region        = job->args->data[2];
-    psMetadata *recipe      = job->args->data[3];
-
-    int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
-    int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
-
-    // radMax stores the upper bounds of the annuli
-    // XXX keep the same name here as for the petrosian / elliptical apertures?
-    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
-    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
-
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    assert (maskVal);
-
-    // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
-
-    float outerRadius = radMax->data.F32[radMax->n - 1];
+    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
+
+    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
 
     // choose the sources of interest
@@ -204,4 +223,9 @@
 
 	pmSource *source = sources->data[i];
+
+	// if we have checked the source validity on the basis of the object set, then 
+	// we either skip these tests below or we skip the source completely
+	if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
+	if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
 
 	// skip PSF-like and non-astronomical objects
@@ -222,4 +246,6 @@
 	if (source->peak->x > region->x1) continue;
 	if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
 	// allocate pmSourceExtendedParameters, if not already defined
@@ -240,17 +266,7 @@
 	}
 
-	// we need to change the view for the radial aperture analysis, but we want to recover exactly
-	// the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
-	psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
-	psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
-	psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
-	psRegion oldRegion    = source->region;
-
 	Nradial ++;
 
-	// force source image to be a bit larger...
-	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
-
-	if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
+	if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
 	    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 	} else {
@@ -258,20 +274,19 @@
 	}
 
-	pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
-	psFree(source->maskObj);   source->maskObj   = oldMaskObj;
-	psFree(source->modelFlux); source->modelFlux = oldModelFlux;
-	psFree(source->psfImage);  source->psfImage  = oldPSFimage;
-	
 	// re-subtract the object, leave local sky
 	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
-    psScalar *scalar = job->args->data[6];
+    psScalar *scalar = job->args->data[5];
     scalar->data.S32 = Nradial;
+
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
 
     return true;
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
-
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
+					    
     // if we are a child source, save the results to the parent source radial aperture array
     psArray *radialAperSet = source->radialAper;
@@ -285,59 +300,68 @@
     radialAperSet->data[entry] = radialAper;
 
-    // storage for the derived pixel values
-    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
+    // find the largest aperture of interest (use only apertures with inner radii <=
+    // source->skyRadius)
+    int lastAp = aperRadii->n;
+    for (int i = 0; i < aperRadii->n; i++) {
+	if (aperRadii->data.F32[i] < source->skyRadius) continue;
+	lastAp = i + 1;
+	break;
+    }
 
     // outer-most radius for initial truncation
-    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax  = aperRadii->data.F32[lastAp - 1];
     float Rmax2 = PS_SQR(Rmax);
 
-    // store the R^2 values for the apertures
-    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
-    for (int i = 0; i < aperRadii->n; i++) {
-	aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
-    }
+    // in this function, the operatins are relative to the full image (readout->image, etc)
 
     float xCM = NAN, yCM = NAN;
     if (pmSourcePositionUseMoments(source)) {
-	xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-	yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+	xCM = source->moments->Mx; // index coord of peak in readout
+	yCM = source->moments->My; // index coord of peak in readout
     } else {
-	xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
-	yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
-    }
+	xCM = source->peak->xf; // index coord of peak in readout
+	yCM = source->peak->yf; // index coord of peak in readout
+    }
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    pixRadius2->n = 0;
+    pixFlux->n = 0;
+    pixVar->n = 0;
 
     // one pass through the pixels to select the valid pixels and calculate R^2
-    for (int iy = 0; iy < source->pixels->numRows; iy++) {
-
-	float yDiff = iy - yCM;
-	if (fabs(yDiff) > Rmax) continue;
-
-	float *vPix = source->pixels->data.F32[iy];
-	float *vWgt = source->variance->data.F32[iy];
-	psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
-
-	for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
-
-	    if (vMsk) {
-		if (*vMsk & maskVal) {
-		    vMsk++;
-		    continue;
-		}
-		vMsk++;
-	    }
-	    if (isnan(*vPix)) continue;
-
-	    float xDiff = ix - xCM;
-	    if (fabs(xDiff) > Rmax) continue;
+    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
+
+	float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
+	int yPix = (int) yDiff;
+
+	if (yPix < 0) continue;
+	if (yPix > Ny - 1) continue;
+	if (fabs(iy) > Rmax) continue;
+
+	float *vPix = readout->image->data.F32[yPix];
+	float *vWgt = readout->variance->data.F32[yPix];
+	psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
+
+	for (int ix = -Rmax; ix < Rmax + 1; ix++) {
+
+	    float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
+	    int xPix = (int) xDiff;
+	    
+	    if (xPix < 0) continue;
+	    if (xPix > Nx - 1) continue;
+	    if (fabs(ix) > Rmax) continue;
+	    
+	    if (vMsk[xPix] & maskVal) continue;
+	    if (isnan(vPix[xPix])) continue;
 
 	    // radius is just a function of (xDiff, yDiff)
-	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    float r2  = PS_SQR(ix) + PS_SQR(iy);
 	    if (r2 > Rmax2) continue;
 
 	    psVectorAppend(pixRadius2, r2);
-	    psVectorAppend(pixFlux, *vPix);
-	    psVectorAppend(pixVar, *vWgt);
+	    psVectorAppend(pixFlux, vPix[xPix]);
+	    psVectorAppend(pixVar, vWgt[xPix]);
 	}
     }
@@ -348,8 +372,15 @@
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
+    // init the apertures of interest to 0.0, the rest go to NAN
     psVectorInit (flux,    0.0);
     psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
+    for (int i = lastAp; i < flux->n; i++) {
+	flux->data.F32[i] = NAN;
+	fluxStd->data.F32[i] = NAN;
+	fluxErr->data.F32[i] = NAN;
+	fill->data.F32[i] = NAN;
+    }
 
     float *rPix2 = pixRadius2->data.F32;
@@ -358,6 +389,8 @@
 	int j = 0;
 	float *aRad2 = aperRadii2->data.F32;
-	for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
-	for (; j < aperRadii2->n; j++, aRad2++) {
+	for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
+
+	// XXX I can speed this up by only saving this single aperture
+	for (; j < lastAp; j++, aRad2++) {
 	    flux->data.F32[j]    += pixFlux->data.F32[i];
 	    fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
@@ -371,7 +404,7 @@
        2) the fractional fill factor (count of valid pixels / effective area of the aperture
        3) the error on the flux within that aperture
-     */
-
-    for (int i = 0; i < flux->n; i++) {
+    */
+
+    for (int i = 0; i < lastAp; i++) {
 	// calculate the total flux for bin 'nOut'
 	float Area = M_PI*aperRadii2->data.F32[i];
@@ -383,5 +416,5 @@
 	// XXX report the total flux or the mask-corrected flux?
 	// flux->data.F32[i]    = SBmean * Area;
-	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
 
 	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
@@ -393,16 +426,21 @@
     }
     
+# if (1)
     radialAper->flux = flux;
     radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
-
-    psFree (aperRadii2);
-    psFree (pixRadius2);
-    psFree (pixFlux);
-    psFree (pixVar);
+# else
+    // XXX TEST
+    psFree(flux);
+    psFree(fluxStd);
+    psFree(fluxErr);
+    psFree(fill);
+# endif
 
     return true;
 }
+
+/*** below is a test to use a sort to speed this up, not very successfully ***/
 
 static int nCalls = 0;
@@ -495,18 +533,18 @@
 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
 # define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_VECT(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B];	\
-    radius->data.F32[B] = tmp; \
-    tmp = pixFlux->data.F32[A]; \
-    pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
-    pixFlux->data.F32[B] = tmp; \
-    tmp = pixVar->data.F32[A]; \
-    pixVar->data.F32[A] = pixVar->data.F32[B]; \
-    pixVar->data.F32[B] = tmp; \
-  } \
-}
+# define SWAP_VECT(TYPE,A,B) {					\
+	float tmp;						\
+	if (A != B) {						\
+	    tmp = radius->data.F32[A];				\
+	    radius->data.F32[A] = radius->data.F32[B];		\
+	    radius->data.F32[B] = tmp;				\
+	    tmp = pixFlux->data.F32[A];				\
+	    pixFlux->data.F32[A] = pixFlux->data.F32[B];	\
+	    pixFlux->data.F32[B] = tmp;				\
+	    tmp = pixVar->data.F32[A];				\
+	    pixVar->data.F32[A] = pixVar->data.F32[B];		\
+	    pixVar->data.F32[B] = tmp;				\
+	}							\
+    }
 
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {
Index: trunk/psphot/src/psphotRadialProfileWings.c
===================================================================
--- trunk/psphot/src/psphotRadialProfileWings.c	(revision 32633)
+++ trunk/psphot/src/psphotRadialProfileWings.c	(revision 32633)
@@ -0,0 +1,386 @@
+# include "psphotInternal.h"
+
+// measure the radial profile outside of the core.  the goal is to find the point at which we
+// reach sky + X sigma
+
+bool psphotRadialProfileWingsSource (pmSource *source, pmReadout *readout, psImageMaskType maskVal);
+bool psphotRadialProfileFluxAtRadius (psVector *flux, psVector *fluxVar, pmSource *source, pmReadout *readout, float Radius, float dRadius, psImageMaskType maskVal);
+
+bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    // return true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Radial Profile Wings ---");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        if (!psphotRadialProfileWingsReadout (config, recipe, view, readout, sources)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// these are set before we fork off to threads and used by all threads as constant values
+static float MAX_RADIUS = NAN;
+static float MIN_RADIUS = NAN;
+static float SKY_STDEV  = NAN;
+// static FILE *file = NULL;
+
+bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources) {
+
+    bool status = false;
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping radial profile wings");
+        return true;
+    }
+
+    psTimerStart ("psphot.wings");
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // XXX is this a good recipe value to use for MAX RADIUS??
+    MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
+    if (!status) {
+        MAX_RADIUS = 50.0;
+    }
+
+    MIN_RADIUS = 0.25*psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
+    if (!status) {
+        MIN_RADIUS = 0.25*psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    }
+
+    // SKY_STDEV is the sigma of the sky model (ie, smoothed on large scales)
+    SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_SIG");
+    if (!status) {
+	SKY_STDEV = 1.0; // a crude default value (why would this not exist?)
+    }
+
+    // bit-masks to test for good/bad pixels
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortByFlux);
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        return true;
+    }
+
+    // file = fopen ("radii.dat", "w");
+
+    // threaded measurement of the source magnitudes
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_RADIAL_PROFILE_WINGS");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,       PS_TYPE_IMAGE_MASK);
+
+// set this to 0 to run without threading
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return false;
+            }
+# else
+	    if (!psphotRadialProfileWings_Threaded(job)) {
+		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		return false;
+	    }
+	    psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+
+    // fclose (file);
+
+    psLogMsg ("psphot.wings", PS_LOG_WARN, "measure radial profile wings : %f sec for %ld objects\n", psTimerMark ("psphot.wings"), sources->n);
+    return true;
+}
+
+bool psphotRadialProfileWings_Threaded (psThreadJob *job) {
+
+    pmReadout *readout              = job->args->data[0];
+    psArray *sources                = job->args->data[1];
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
+
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+	if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+	// allocate space for moments
+	if (!source->moments) continue;
+
+	// replace object in image
+	bool reSubtract = false;
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	    reSubtract = true;
+	}
+
+	// re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+	// XXX don't measure on subraster images:
+	// pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+	// psAssert (source->pixels, "WTF?");
+
+	// this function populates moments->Mrf,KronFlux,KronFluxErr
+	psphotRadialProfileWingsSource (source, readout, maskVal);
+
+	// if we subtracted it above, re-subtract the object, leave local sky
+	if (reSubtract) {
+	    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+	}
+    }
+    return true;
+}
+
+# define TEST_X 3158
+# define TEST_Y 3096
+
+// XXX use integer radius values?  the rings assume integer values, right? or do they?
+bool psphotRadialProfileWingsSource (pmSource *source, pmReadout *readout, psImageMaskType maskVal) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+
+    // psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    // XXX unclear if I should run this analysis on both 1st and 2nd pass for 1st pass objects,
+    // or just use the 1st pass value.  I think I should just use the 1st pass value, so skip
+    // any that have already been assigned
+    if (isfinite(source->skyRadius)) return true;
+
+    // radii will be MIN_RADIUS to MAX_RADIUS in NN log steps:
+    float NSTEP = 10.0;
+    float MIN_DR = 2;
+    float NSIGMA = 1.0;
+    float THRESHOLD = 2.0*SKY_STDEV;
+    float alpha = pow ((MAX_RADIUS / MIN_RADIUS), 1.0/NSTEP) - 1.0;
+    float dRmax = MAX_RADIUS * alpha / (1.0 + alpha); // approximate size of last annulus, to get a rough size for vector allocation
+
+    int iter = 0;
+
+    psVector *flux = psVectorAllocEmpty(7*MAX_RADIUS*dRmax, PS_TYPE_F32);
+    psVector *fluxVar = psVectorAllocEmpty(7*MAX_RADIUS*dRmax, PS_TYPE_F32);
+
+    // should I just use sample median here?
+    psStats *fluxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    psStats *varStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+
+    float lastFlux = NAN; 
+    float lastRadius = NAN; 
+    float lastSlope = NAN; 
+    bool limit = false;
+    float limitRadius = NAN;
+    float limitFlux = NAN;
+    float limitSlope = NAN;
+
+    // note: radius is the inner radius of the annulus; outer radius = inner radius + dR
+    for (float radius = MIN_RADIUS; !limit && (radius < MAX_RADIUS); iter ++) {
+
+	float dR = (int)(radius * alpha);
+	if (dR < MIN_DR) dR = MIN_DR;
+	float outerRadius = radius + dR;
+	float meanRadius = (2.0/3.0) * (outerRadius*outerRadius*outerRadius - radius*radius*radius) / (PS_SQR(outerRadius) - PS_SQR(radius));
+
+	// extract a vector of the pixel values (signal, variance) at this radius + dR
+	psphotRadialProfileFluxAtRadius (flux, fluxVar, source, readout, radius, dR, maskVal);
+
+	psStatsInit (fluxStats);
+	psStatsInit (varStats);
+
+	psVectorStats (fluxStats, flux, NULL, NULL, 0);
+	psVectorStats (varStats, fluxVar, NULL, NULL, 0);
+	
+	// get the mean / median flux statistic and significance here
+	float meanFlux = fluxStats->robustMedian;
+	float meanFluxError = sqrt(varStats->sampleMean / fluxVar->n);
+	// this is a bit crude on the flux error: the error. technically, it should be
+	// sqrt(\sum(variance)) / Npts.  I am using the identity \sum(variance) =
+	// \mean(variance) * Npts and cancelling the Npts term inside and out of the sqrt()
+	
+	float slope = NAN;
+	if (isfinite(lastFlux)) {
+	    slope = (meanFlux - lastFlux) / (meanRadius - lastRadius);
+	}
+
+	// fprintf (stderr, "%f %f : %f : %f %f  :  %f\n", source->peak->xf, source->peak->yf, radius, meanFlux, meanFluxError, slope);
+
+	if (!limit) {
+	    limit |= (meanFlux - NSIGMA * meanFluxError < THRESHOLD); // dropped to sky level
+	    limit |= isfinite(slope) && (fabs(slope) < 3.0); // SB no longer changing.
+	    limitRadius = meanRadius;
+	    limitFlux = meanFlux;
+	    limitSlope = slope;
+	}
+
+	// completion criteria:
+	// 1) flux - NSIGMA * dflux <= sky
+	// 2) flux rising?
+	// 3) flux flat?
+
+	lastFlux = meanFlux;
+	lastRadius = meanRadius;
+	lastSlope = slope;
+
+	// reset the flux & fluxVar vector length to zero for re-use above:
+	flux->n = 0;
+	fluxVar->n = 0;
+	radius += dR;
+    }
+
+    if (!limit) {
+	limitRadius = lastRadius;
+	limitFlux = lastFlux;
+	limitSlope = lastSlope;
+    }
+    // fprintf (file, "%f %f : %f %f : %f\n", source->peak->xf, source->peak->yf, limitRadius, limitFlux, limitSlope);
+
+    psFree (flux);
+    psFree (fluxVar);
+    psFree (fluxStats);
+    psFree (varStats);
+
+    source->skyRadius = limitRadius;
+    source->skyFlux   = limitFlux;
+    source->skySlope  = limitSlope;
+
+    // save the max radius (and anything else?)
+    // source->moments->Mrf = Mrf;
+
+    return true;
+}
+
+// inline this?  macro this?
+bool psphotRadialProfileGetFlux (psVector *flux, psVector *fluxVar, pmReadout *readout, int xc, int yc, psImageMaskType maskVal) {
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    if (xc < 0) return false;
+    if (xc >= Nx) return false;
+    if (yc < 0) return false;
+    if (yc >= Ny) return false;
+    
+    if (readout->mask && (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yc][xc] & maskVal)) return false;
+    psVectorAppend (flux, readout->image->data.F32[yc][xc]);
+    psVectorAppend (fluxVar, readout->variance->data.F32[yc][xc]);
+    return true;
+}
+
+// select the pixels at the given radius and extract the flux and variance into the given vectors
+// XXX should radius be 'int'?
+bool psphotRadialProfileFluxAtRadius (psVector *flux, psVector *fluxVar, pmSource *source, pmReadout *readout, float Radius, float dRadius, psImageMaskType maskVal) {
+
+    psAssert (flux, "must allocate output flux vector");
+    psAssert (fluxVar, "must allocate output fluxVar vector");
+
+    // the peak position is less accurate but less subject to extreme deviations
+    float dX = source->moments->Mx - source->peak->xf;
+    float dY = source->moments->My - source->peak->yf;
+    float dR = hypot(dX, dY);
+    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
+    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
+
+    for (int radius = Radius; radius < Radius + dRadius; radius ++) {
+
+	int x = 0;
+	int y = radius;
+	int d = 5 - 4*radius;
+
+	while (x <= y) {
+	    psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo + x), (Yo + y), maskVal);
+	    psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo + x), (Yo - y), maskVal);
+	    psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo + y), (Yo + x), maskVal);
+	    psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo - y), (Yo + x), maskVal);
+	    
+	    if (x > 0) {
+		psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo - x), (Yo + y), maskVal);
+		psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo - x), (Yo - y), maskVal);
+		psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo - y), (Yo - x), maskVal);
+		psphotRadialProfileGetFlux(flux, fluxVar, readout, (Xo + y), (Yo - x), maskVal);
+	    }
+	    
+	    if (d < 0) {
+		d = d + 8*x + 4;
+	    } else {
+		d = d + 8*(x-y) + 8;
+		y--;
+	    }
+	    x++;
+	}
+    }
+    return true;
+}
Index: trunk/psphot/src/psphotReadout.c
===================================================================
--- trunk/psphot/src/psphotReadout.c	(revision 32632)
+++ trunk/psphot/src/psphotReadout.c	(revision 32633)
@@ -190,4 +190,7 @@
     psphotDumpChisqs (config, view, filerule);
 
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, filerule);
+
     // re-measure the kron mags with models subtracted.  this pass uses a circular window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     
@@ -305,4 +308,7 @@
 pass1finish:
 
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, filerule);
+
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, filerule);
Index: trunk/psphot/src/psphotSetThreads.c
===================================================================
--- trunk/psphot/src/psphotSetThreads.c	(revision 32632)
+++ trunk/psphot/src/psphotSetThreads.c	(revision 32633)
@@ -50,6 +50,11 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 7);
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 6);
     task->function = &psphotRadialApertures_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_PROFILE_WINGS", 3);
+    task->function = &psphotRadialProfileWings_Threaded;
     psThreadTaskAdd(task);
     psFree(task);
Index: trunk/psphot/src/psphotSourceFits.c
===================================================================
--- trunk/psphot/src/psphotSourceFits.c	(revision 32632)
+++ trunk/psphot/src/psphotSourceFits.c	(revision 32633)
@@ -357,4 +357,6 @@
     pmSource *newSrc = pmSourceCopy (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
+    newSrc->peak->footprint = source->peak->footprint; // just a reference; the peak does not own the footprint
+    psArrayAdd(newSrc->peak->footprint->peaks, 1, newSrc->peak); // the footprint owns the peak
 
     // build cached models and subtract
Index: trunk/psphot/src/psphotSourceMatch.c
===================================================================
--- trunk/psphot/src/psphotSourceMatch.c	(revision 32632)
+++ trunk/psphot/src/psphotSourceMatch.c	(revision 32633)
@@ -67,6 +67,9 @@
     objects = psArraySort (objects, pmPhotObjSortByX); 
  
-    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
-    psVectorInit (found, 0);
+    psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundSrc, 0);
+
+    psVector *foundObj = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundObj, 0);
 
     // match sources to existing objects
@@ -93,4 +96,11 @@
         if (dx > +1.02*RADIUS) NEXT2; 
  
+	/* this block will match a given detection to the closest object within range of that detection.
+	   XXX note that this matches ALL detections within range of the single object to that same object 
+	   this is bad, but I cannot just go in linear order (ie, mark off each object as they are
+	   used).  I should make a list of all Nobj * Ndet pairs in range and choose the matches
+	   based on their separations.  UGH
+	*/
+    
         // we are within match range, look for matches: 
 	int Jmin = -1;
@@ -98,4 +108,6 @@
         for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
  
+	    // skip objects that are already assigned:
+	    if (foundObj->data.U8[J]) continue;
 	    obj = objects->data[J]; 
 	    
@@ -117,8 +129,9 @@
 	}
 	obj = objects->data[Jmin]; 
+	foundObj->data.U8[Jmin] = 1;
 
 	// add to object
 	pmPhotObjAddSource (obj, src);
-	found->data.U8[i] = 1;
+	foundSrc->data.U8[i] = 1;
         i++; 
     } 
@@ -128,5 +141,5 @@
     for (i = 0; i < sources->n; i++) {
 
-	if (found->data.U8[i]) continue;
+	if (foundSrc->data.U8[i]) continue;
 
         pmSource *src = sources->data[i]; 
@@ -139,5 +152,6 @@
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
 
-    psFree (found);
+    psFree (foundSrc);
+    psFree (foundObj);
     return true;
 } 
@@ -234,10 +248,12 @@
 	    
 	    // assign to a footprint on this readout->image
-	    peak->footprint = pmFootprintCopyData(footprint, readout->image);
-
-	    // the peak does not claim ownership of the footprint (it does not free it). save a copy of this 
-	    // footprint on detections->footprints so we can free it later
-	    psArrayAdd(detections->footprints, 100, peak->footprint); 
-	    psFree (peak->footprint);
+	    if (footprint) {
+		peak->footprint = pmFootprintCopyData(footprint, readout->image);
+
+		// the peak does not claim ownership of the footprint (it does not free it). save a copy of this 
+		// footprint on detections->footprints so we can free it later
+		psArrayAdd(detections->footprints, 100, peak->footprint); 
+		psFree (peak->footprint);
+	    }
 	    
 	    // create a new source
@@ -268,4 +284,5 @@
         pmPhotObj *obj = objects->data[i]; 
 	nSources += obj->sources->n;
+	psAssert (obj->sources->n == nImages, "failed to match sources?");
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
Index: trunk/psphot/src/psphotSourceSize.c
===================================================================
--- trunk/psphot/src/psphotSourceSize.c	(revision 32632)
+++ trunk/psphot/src/psphotSourceSize.c	(revision 32633)
@@ -156,5 +156,5 @@
 
     // XXX fix this (was source->n  - first)
-    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
+    psLogMsg ("psphot.size", PS_LOG_WARN, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
 
     psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
Index: trunk/psphot/src/psphotSourceStats.c
===================================================================
--- trunk/psphot/src/psphotSourceStats.c	(revision 32632)
+++ trunk/psphot/src/psphotSourceStats.c	(revision 32633)
@@ -124,4 +124,6 @@
         // add the peak
         source->peak = psMemIncrRefCounter(peak);
+
+	// psAssert (source->peak->footprint, "peak without footprint??");
 
         // allocate space for moments
Index: trunk/psphot/src/psphotStackObjects.c
===================================================================
--- trunk/psphot/src/psphotStackObjects.c	(revision 32632)
+++ trunk/psphot/src/psphotStackObjects.c	(revision 32633)
@@ -70,5 +70,6 @@
 
     // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_PETRO  = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
@@ -81,5 +82,6 @@
 	// we check each source for an object and keep the object if any source is valid
 
-	bool keepObject = false;
+	bool keepObjectRadial = false;
+	bool keepObjectPetro = false;
 	for (int j = 0; j < object->sources->n; j++) {
 
@@ -94,21 +96,4 @@
 	    if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
 	    
-	// optionally allow non-extended objects to get petrosians as well
-	if (!doPetroStars) {
-	    if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-	    if (source->type == PM_SOURCE_TYPE_STAR) continue;
-	}
-
-	    // limit selection to some SN limit
-	    // assert (source->peak); // how can a source not have a peak?
-	    // limit selection to some SN limit
-	    bool skipSource = false;
-	    if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-		skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
-	    } else {
-		skipSource = (sqrt(source->peak->detValue) < SN_LIM);
-	    }
-	    if (skipSource) continue;
-
 	    // limit selection by analysis region (this automatically apply
 	    if (source->peak->x < AnalysisRegion.x0) continue;
@@ -117,5 +102,23 @@
 	    if (source->peak->y > AnalysisRegion.y1) continue;
 	    
-	    keepObject = true;
+	    // SN limit tests for RADIAL APERTURES:
+	    bool skipSourceRadial = false;
+	    if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+		skipSourceRadial = (source->moments->KronFlux < SN_LIM_RADIAL * source->moments->KronFluxErr);
+	    } else {
+		skipSourceRadial = (sqrt(source->peak->detValue) < SN_LIM_RADIAL);
+	    }
+	    if (!skipSourceRadial) keepObjectRadial = true;
+
+	    // SN limit tests for PETRO
+	    bool skipSourcePetro = false;
+	    if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+		skipSourcePetro = (source->moments->KronFlux < SN_LIM_PETRO * source->moments->KronFluxErr);
+	    } else {
+		skipSourcePetro = doPetroStars ? (sqrt(source->peak->detValue) < SN_LIM_PETRO) : true;
+	    }
+	    if (!skipSourcePetro) keepObjectPetro = true;
+
+	    keepObjectPetro = true;
 	}
 
@@ -128,10 +131,18 @@
 	    // avoid the single-detection tests
 
-	    if (keepObject) {
-		source->tmpFlags |=  PM_SOURCE_TMPF_STACK_KEEP;
-		source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_SKIP;
+	    if (keepObjectPetro) {
+		source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_KEEP;
+		source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_SKIP;
 	    } else {
-		source->tmpFlags |=  PM_SOURCE_TMPF_STACK_SKIP;
-		source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_KEEP;
+		source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_SKIP;
+		source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_KEEP;
+	    }	    
+
+	    if (keepObjectRadial) {
+		source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_KEEP;
+		source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP;
+	    } else {
+		source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_SKIP;
+		source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP;
 	    }	    
 	}
Index: trunk/psphot/src/psphotStackReadout.c
===================================================================
--- trunk/psphot/src/psphotStackReadout.c	(revision 32632)
+++ trunk/psphot/src/psphotStackReadout.c	(revision 32633)
@@ -43,4 +43,6 @@
 
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
+
+    psArray *objects = NULL; // used below after 'pass1finish' label
 
     // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
@@ -169,4 +171,7 @@
     psphotStackVisualFilerule(config, view, STACK_SRC);
 
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
     // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
@@ -184,10 +189,10 @@
     psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
+    // if we only do one pass, skip to extended source analysis
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
-
-    // if we only do one pass, skip to extended source analysis
-    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -238,10 +243,23 @@
 	// XXX check on free of sources...
 	psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
-
-	// NOTE: apply to ALL sources
-	psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
     }
 
 pass1finish:
+
+    // generate the objects (objects unify the sources from the different images) NOTE: could
+    // this just match the detections for the chisq image, and not bother measuring the source
+    // stats in that case...?
+    objects = psphotMatchSources (config, view, STACK_SRC);
+    psMemDump("matchsources");
+
+    psphotStackObjectsUnifyPosition (objects);
+
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+
+    // NOTE: apply to ALL sources
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
+
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, STACK_SRC);
 
     // re-measure the kron mags with models subtracted
@@ -254,14 +272,4 @@
 
     psMemDump("psfstats");
-
-    // generate the objects (objects unify the sources from the different images)
-    // XXX this could just match the detections for the chisq image, and not bother measuring the
-    // source stats in that case...
-    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
-    psMemDump("matchsources");
-
-    psphotStackObjectsUnifyPosition (objects);
-
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
@@ -281,7 +289,5 @@
 
     // measure circular, radial apertures (objects sorted by S/N)
-    // XXX can we just use psphotRadialApertures
-    // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none')
-    // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF); 
+    // this forces photometry on the undetected sources from other images
     psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     psMemDump("extmeas");
@@ -323,5 +329,7 @@
     psphotMagnitudes(config, view, STACK_SRC);
 
-    if (0 && !psphotEfficiency(config, view, STACK_DET)) {
+    // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the 
+    // PSF of the SRC image.  can we fake it by generating the PSF for DET as well (up above)? 
+    if (false && !psphotEfficiency(config, view, STACK_DET)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();