Index: /branches/eam_branch_20080324/psphot/src/psphotSourceSize.c
===================================================================
--- /branches/eam_branch_20080324/psphot/src/psphotSourceSize.c	(revision 17238)
+++ /branches/eam_branch_20080324/psphot/src/psphotSourceSize.c	(revision 17239)
@@ -1,4 +1,6 @@
 # include "psphotInternal.h"
 # include <gsl/gsl_sf_gamma.h>
+
+float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -14,5 +16,8 @@
     psTimerStart ("psphot");
 
-    float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT");
+    float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
+    assert (status);
+
+    float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
     assert (status);
 
@@ -25,5 +30,5 @@
 
 	source->crNsigma  = -1.0;
-	source->extNsigma = -1.0;
+	source->extNsigma = 0.0;
 
 	// source must have been subtracted
@@ -31,7 +36,7 @@
 	if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
 
-	psF32 **resid = source->pixels->data.F32;
+	psF32 **resid  = source->pixels->data.F32;
 	psF32 **weight = source->weight->data.F32;
-	psU8 **mask = source->maskObj->data.U8;
+	psU8 **mask    = source->maskObj->data.U8;
 
 	int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
@@ -56,15 +61,10 @@
 
 	// measure the flux at the 1 sigma contour
-	psVector *contour = psphotModelContour (source->pixels, source->maskObj, source->modelPSF, 1.0);
-
-	// XXX for now, just skip any masked pixels
-	extSum = 0.0;
-	bool keep = true;
-	for (int iy = -1; (iy <= +1) && keep; iy++) {
-	    for (int ix = -1; (ix <= +1) && keep; ix++) {
-		if (mask[yPeak+iy][xPeak+ix]) { keep &= false; }
-	    }
-	}
-	if (!keep) continue;
+	// XXX prevent a source from being both CR and EXT?
+	source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0);
+
+	if (source->extNsigma > EXT_NSIGMA_LIMIT) {
+	  source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
+	}
 
 	// XXX need to deal with edge peaks... and mask
@@ -126,5 +126,4 @@
 	}
 	source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
-	source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0;
 	// NOTE: abs needed to make the Nsigma value positive 
 
@@ -135,5 +134,5 @@
 
 	if (source->crNsigma > CR_NSIGMA_LIMIT) {
-	  source->mode |= PM_SOURCE_MODE_CRLIMIT;
+	  source->mode |= PM_SOURCE_MODE_CR_LIMIT;
 	}
     }
@@ -183,13 +182,10 @@
 
 
-// given the PSF ellipse parameters, navigate around the 1sigma contour
-// XXX return the Nsigma total deviation?  
-// XXX return just the sum around the contour?  
-// this is measure on the residual image - should we ignore negative deviations?
-psVector *psphotModelContour (psImage *image, psImage *mask, pmModel *model, float Ro) {
+// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
+// deviation in sigmas.  This is measure on the residual image - should we ignore negative
+// deviations?
+float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) {
 
     psF32 *PAR = model->params->data.F32;
-    psVector *contour = psVectorAllocEmpty (50, PS_TYPE_F32);
-    int nPts = 0;
 
     // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY;
@@ -201,10 +197,13 @@
     // solve this for x2:
     float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0);
-    if (Q < 0.0) return contour; // ellipse is imaginary
+    if (Q < 0.0) return NAN; // ellipse is imaginary
 
     int xMax = sqrt(Q);
     int xMin = -1.0*xMax;
 
-    for (int x = MIN; x <= MAX; x++) {
+    int nPts = 0;
+    float nSigma = 0.0;
+
+    for (int x = xMin; x <= xMax; x++) {
 	float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]);
 	float B = x * PAR[PM_PAR_SXY];
@@ -217,16 +216,16 @@
 	float yM = (-B - sqrt (T)) / (2.0 * A);
 
-	int xPix  = x  + PAR[PM_PAR_X0] - source->pixels->col0 + 0.5;
-	int yPixM = yM + PAR[PM_PAR_Y0] - source->pixels->row0 + 0.5;
-	int yPixP = yP + PAR[PM_PAR_Y0] - source->pixels->row0 + 0.5;
+	int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
+	int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
+	int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
 
 	if (xPix < 0) continue;
-	if (xPix >= source->pixels->numCols) continue;
-
-	if ((yPixM >= 0) && (yPixM < source->pixels->numRows)) {
-	    if (!mask || !mask->data.U8[xPix][yPixM]) {
-		contour->data.F32[nPts] = image->data.F32[xPix][yPixM];
-		psVectorExtend (contour, 100, 1);
-		nPts++;
+	if (xPix >= image->numCols) continue;
+
+	if ((yPixM >= 0) && (yPixM < image->numRows)) {
+	    if (!mask || !mask->data.U8[yPixM][xPix]) {
+		float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]);
+		nSigma += dSigma;
+		nPts ++;
 	    }
 	}
@@ -234,13 +233,13 @@
 	if (yPixM == yPixP) continue;
 
-	if ((yPixP >= 0) && (yPixP < source->pixels->numRows)) {
-	    if (!mask || !mask->data.U8[xPix][yPixP]) {
-		contour->data.F32[nPts] = image->data.F32[xPix][yPixP];
-		psVectorExtend (contour, 100, 1);
-		nPts++;
+	if ((yPixP >= 0) && (yPixP < image->numRows)) {
+	    if (!mask || !mask->data.U8[yPixP][xPix]) {
+		float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
+		nSigma += dSigma;
+		nPts ++;
 	    }
 	}
     }	
-
-    return contour;
+    nSigma /= nPts;
+    return nSigma;
 }