Index: /branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceMoments.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceMoments.c	(revision 32215)
+++ /branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceMoments.c	(revision 32216)
@@ -80,8 +80,4 @@
     }
 
-    if (source->moments->nPixels != 0) {
-	fprintf (stderr, "remeasure moments: %f,%f\n", source->peak->xf, source->peak->yf);
-    }
-
     float Sum = 0.0;
     float Var = 0.0;
@@ -120,6 +116,6 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
-    float RFW = 0.0;
-    float RHW = 0.0;
+    float RFa = 0.0;
+    float RSa = 0.0;
 
     float RF = 0.0;
@@ -154,4 +150,5 @@
     float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
+    // calculate the higher-order moments using Xo,Yo
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
@@ -205,12 +202,5 @@
 	    Sum += pDiff;
 
-	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
 	    float r = sqrt(r2);
-	    float rf = r * fDiff;
-	    float rh = sqrt(r) * fDiff;
-	    float rs = fDiff;
-
-	    float rfw = r * pDiff;
-	    float rhw = sqrt(r) * pDiff;
 
 	    float x = xDiff * pDiff;
@@ -232,11 +222,4 @@
 	    float yyyy = yDiff * yyy / r2;
 
-	    RF  += rf;
-	    RH  += rh;
-	    RS  += rs;
-
-	    RFW  += rfw;
-	    RHW  += rhw;
-
 	    XX  += xx;
 	    XY  += xy;
@@ -253,10 +236,22 @@
 	    XYYY  += xyyy;
 	    YYYY  += yyyy;
+
+	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
+	    // XXX float r = sqrt(r2);
+	    // XXX float rf = r * fDiff;
+	    // XXX float rh = sqrt(r) * fDiff;
+	    // XXX float rs = fDiff;
+	    // XXX 
+	    // XXX float rfw = r * pDiff;
+	    // XXX float rhw = sqrt(r) * pDiff;
+	    // XXX 
+	    // XXX RF  += rf;
+	    // XXX RH  += rh;
+	    // XXX RS  += rs;
+	    // XXX 
+	    // XXX RFW  += rfw;
+	    // XXX RHW  += rhw;
 	}
     }
-
-    source->moments->Mrf = RF/RS;
-    source->moments->Mrh = RH/RS;
-
     source->moments->Mxx = XX/Sum;
     source->moments->Mxy = XY/Sum;
@@ -273,4 +268,86 @@
     source->moments->Mxyyy = XYYY/Sum;
     source->moments->Myyyy = YYYY/Sum;
+
+# define TEST_X1 167
+# define TEST_Y1 299
+# define TEST_X2 180
+# define TEST_Y2 300
+    if ((fabs(Xo - TEST_X1) < 3) && (fabs(Yo - TEST_Y1) < 3)) {
+	fprintf (stderr, "test obj 1\n");
+    }
+    if ((fabs(Xo - TEST_X2) < 3) && (fabs(Yo - TEST_Y2) < 3)) {
+	fprintf (stderr, "test obj 2\n");
+    }
+
+    float **vPix = source->pixels->data.F32;
+    float **vWgt = source->variance->data.F32;
+    psImageMaskType  **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    // calculate the 1st radial moment (for kron flux) -- symmetrical averaging
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+	float yDiff = row - yCM;
+	if (fabs(yDiff) > radius) continue;
+
+	// coordinate of mirror pixel
+	int yFlip = yCM - yDiff;
+	if (yFlip < 0) continue;
+	if (yFlip >= source->pixels->numRows) continue;
+
+	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+	    // check mask and value for this pixel
+	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
+	    if (isnan(vPix[row][col])) continue;
+
+	    float xDiff = col - xCM;
+	    if (fabs(xDiff) > radius) continue;
+
+	    // coordinate of mirror pixel
+	    int xFlip = xCM - xDiff;
+	    if (xFlip < 0) continue;
+	    if (xFlip >= source->pixels->numCols) continue;
+
+	    // check mask and value for mirror pixel
+	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+	    if (isnan(vPix[yFlip][xFlip])) continue;
+
+	    // radius is just a function of (xDiff, yDiff)
+	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    if (r2 > R2) continue;
+
+	    float fDiff1 = vPix[row][col] - sky;
+	    float fDiff2 = vPix[yFlip][xFlip] - sky;
+	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+
+	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
+	    float r = sqrt(r2);
+	    float rf = r * pDiff;
+	    float rh = sqrt(r) * pDiff;
+	    float rs = 0.5 * (fDiff1 + fDiff2);
+
+	    float rfa = r * fDiff1;
+	    float rsa = fDiff1;
+
+	    RF  += rf;
+	    RH  += rh;
+	    RS  += rs;
+
+	    RFa  += rfa;
+	    RSa  += rsa;
+	}
+    }
+
+    source->moments->Mrf = RF/RS;
+    source->moments->Mrh = RH/RS;
+
+    float R1 = RFa / RSa;
+    if ((fabs(Xo - TEST_X1) < 3) && (fabs(Yo - TEST_Y1) < 3)) {
+	fprintf (stderr, "R1: %f vs %f\n", R1, source->moments->Mrf);
+    }
+    if ((fabs(Xo - TEST_X2) < 3) && (fabs(Yo - TEST_Y2) < 3)) {
+	fprintf (stderr, "R2: %f vs %f\n", R1, source->moments->Mrf);
+    }
+
+    // fprintf (stderr, "Rad: %f vs %f\n", R1, source->moments->Mrf);
 
     // if Mrf (first radial moment) is very small, we are getting into low-significance
@@ -294,14 +371,155 @@
     float SumOuter = 0.0;
 
+    // calculate the Kron flux, and related fluxes (symmetrical averaging)
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
+	
 	float yDiff = row - yCM;
 	if (fabs(yDiff) > radKouter) continue;
+	
+	// coordinate of mirror pixel
+	int yFlip = yCM - yDiff;
+	if (yFlip < 0) continue;
+	if (yFlip >= source->pixels->numRows) continue;
+	
+	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+	    // check mask and value for this pixel
+	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
+	    if (isnan(vPix[row][col])) continue;
+	    
+	    float xDiff = col - xCM;
+	    if (fabs(xDiff) > radKouter) continue;
+	    
+	    // coordinate of mirror pixel
+	    int xFlip = xCM - xDiff;
+	    if (xFlip < 0) continue;
+	    if (xFlip >= source->pixels->numCols) continue;
+	    
+	    // check mask and value for mirror pixel
+	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+	    if (isnan(vPix[yFlip][xFlip])) continue;
+	    
+	    // radKron is just a function of (xDiff, yDiff)
+	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+
+	    float fDiff1 = vPix[row][col] - sky;
+	    float fDiff2 = vPix[yFlip][xFlip] - sky;
+	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+	    // float pDiff = vPix[row][col] - sky;
+	    float wDiff = vWgt[row][col];
+				    
+	    // skip pixels below specified significance level.  this is allowed, but should be
+	    // avoided -- the over-weights the wings of bright stars compared to those of faint
+	    // stars.
+	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
+	    
+# define WEIGHTED 0
+# if (WEIGHTED)
+	    float z = r2 * rsigma2 / 4.0;
+	    assert (z >= 0.0);
+	    float weight  = exp(-z);
+# else
+	    float weight  = 1.0;
+# endif
+
+	    float r  = sqrt(r2);
+	    if (r < radKron) {
+		Sum += pDiff*weight;
+		Var += wDiff*weight;
+		nKronPix ++;
+		// if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
+	    }
+
+	    // use sigma (fixed by psf) not a radKron based value
+	    if (r < sigma) {
+		SumCore += pDiff;
+		VarCore += wDiff;
+		nCorePix ++;
+	    }
+
+	    if ((r > radKinner) && (r < radKron)) {
+		SumInner += pDiff;
+		nInner ++;
+	    }
+	    if ((r > radKron)  && (r < radKouter)) {
+		SumOuter += pDiff;
+		nOuter ++;
+	    }
+	}
+    }
+    // *** should I rescale these fluxes by pi R^2 / nNpix?
+    // XXX source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
+    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
+    // XXX source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
+    // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
+    // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
+
+    source->moments->KronCore    = SumCore;
+    source->moments->KronCoreErr = sqrt(VarCore);
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    source->moments->KronFinner = SumInner;
+    source->moments->KronFouter = SumOuter;
+
+    // XXX not sure I should save this here...
+    source->moments->KronFluxPSF    = source->moments->KronFlux;
+    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+    source->moments->KronRadiusPSF  = source->moments->Mrf;
+
+    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
+	     source->moments->Mrf,   source->moments->KronFlux, 
+	     source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
+	     source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
+	     source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
+
+    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
+	     source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
+
+    return(true);
+}
+
+bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+
+    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
+    // Sum = SUM (z - sky)
+    // X1  = SUM (x - xc)*(z - sky)
+    // .. etc
+
+    float sky = 0.0;
+
+    float peakPixel = -PS_MAX_F32;
+    psS32 numPixels = 0;
+    float Sum = 0.0;
+    float Var = 0.0;
+    float X1 = 0.0;
+    float Y1 = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
+
+    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
+    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
+
+    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
+    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
+    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
+    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
+    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
+
+    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
+    // not depend on the fractional pixel location of the source.  However, the aperture
+    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
+    // position of the expected centroid
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+	float yDiff = row + 0.5 - yPeak;
+	if (fabs(yDiff) > radius) continue;
 
 	float *vPix = source->pixels->data.F32[row];
 	float *vWgt = source->variance->data.F32[row];
 
-	psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-	// psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
+	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+	// psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
 	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -315,134 +533,4 @@
 	    if (isnan(*vPix)) continue;
 
-	    float xDiff = col - xCM;
-	    if (fabs(xDiff) > radKouter) continue;
-
-	    // radKron is just a function of (xDiff, yDiff)
-	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-
-	    float pDiff = *vPix - sky;
-	    float wDiff = *vWgt;
-
-	    // skip pixels below specified significance level.  this is allowed, but should be
-	    // avoided -- the over-weights the wings of bright stars compared to those of faint
-	    // stars.
-	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-
-# define WEIGHTED 1
-# if (WEIGHTED)
-	    float z = r2 * rsigma2 / 4.0;
-	    assert (z >= 0.0);
-	    float weight  = exp(-z);
-# endif
-
-	    float r  = sqrt(r2);
-	    if (r < radKron) {
-# if (WEIGHTED)
-		Sum += pDiff*weight;
-		Var += wDiff*weight;
-# else
-		Sum += pDiff;
-		Var += wDiff;
-# endif
-		nKronPix ++;
-		// if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
-	    }
-
-	    // use sigma (fixed by psf) not a radKron based value
-	    if (r < sigma) {
-		SumCore += pDiff;
-		VarCore += wDiff;
-		nCorePix ++;
-	    }
-
-	    if ((r > radKinner) && (r < radKron)) {
-		SumInner += pDiff;
-		nInner ++;
-	    }
-	    if ((r > radKron)  && (r < radKouter)) {
-		SumOuter += pDiff;
-		nOuter ++;
-	    }
-	}
-    }
-    // *** should I rescale these fluxes by pi R^2 / nNpix?
-    source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
-    source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
-    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
-    source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
-
-    // XXX not sure I should save this here...
-    source->moments->KronFluxPSF    = source->moments->KronFlux;
-    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
-    source->moments->KronRadiusPSF  = source->moments->Mrf;
-
-    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
-	     source->moments->Mrf,   source->moments->KronFlux, 
-	     source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
-	     source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
-	     source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
-
-    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
-	     source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
-
-    return(true);
-}
-
-bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
-
-    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // .. etc
-
-    float sky = 0.0;
-
-    float peakPixel = -PS_MAX_F32;
-    psS32 numPixels = 0;
-    float Sum = 0.0;
-    float Var = 0.0;
-    float X1 = 0.0;
-    float Y1 = 0.0;
-    float R2 = PS_SQR(radius);
-    float minSN2 = PS_SQR(minSN);
-    float rsigma2 = 0.5 / PS_SQR(sigma);
-
-    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
-    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
-
-    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
-    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
-    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
-    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
-
-    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
-    // not depend on the fractional pixel location of the source.  However, the aperture
-    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
-    // position of the expected centroid
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-	float yDiff = row + 0.5 - yPeak;
-	if (fabs(yDiff) > radius) continue;
-
-	float *vPix = source->pixels->data.F32[row];
-	float *vWgt = source->variance->data.F32[row];
-
-	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-	// psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
-
-	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
-	    if (vMsk) {
-		if (*vMsk & maskVal) {
-		    vMsk++;
-		    continue;
-		}
-		vMsk++;
-	    }
-	    if (isnan(*vPix)) continue;
-
 	    float xDiff = col + 0.5 - xPeak;
 	    if (fabs(xDiff) > radius) continue;
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c	(revision 32215)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c	(revision 32216)
@@ -10,4 +10,6 @@
 {
     bool status = true;
+
+    // return true;
 
     // select the appropriate recipe information
@@ -127,93 +129,63 @@
     psphotSaveImage (NULL, kronWindow, "kron.window.v0.fits");
 
-    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
-	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+    for (int j = 0; j < 5; j++) {
+	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
+	    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+		pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	    }
+
+	    // iterate to the window radius
+	    float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+
+	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+	    pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+	    // clear the window function for this source based on the moments
+	    psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
+	    // psphotKronWindowSetSource (source, kronWindow, false, false);
+	    // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
+
+	    // 165, 539;
+	    if ((fabs(source->peak->xf - 165) < 3) && (fabs(source->peak->yf - 539) < 3)) {
+		fprintf (stderr, "test obj\n");
+	    }
+
+	    // this function populates moments->Mrf,KronFlux,KronFluxErr
+	    psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
+	    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+	    // set a window function for each source based on the moments
+	    psphotKronWindowSetSource (source, kronWindow, true, true);
+	    // psphotKronWindowSetSource (source, kronWindow, false, true);
+
+	    // test source fluxes
+	    pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+	    float kmag = -2.5*log10(source->moments->KronFlux);
+# define TEST_X1 167
+# define TEST_Y1 299
+# define TEST_X2 180
+# define TEST_Y2 300
+	    if ((fabs(source->peak->xf - TEST_X1) < 3) && (fabs(source->peak->yf - TEST_Y1) < 3)) {
+		fprintf (stderr, "R1: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
+	    }
+	    if ((fabs(source->peak->xf - TEST_X2) < 3) && (fabs(source->peak->yf - TEST_Y2) < 3)) {
+		fprintf (stderr, "R2: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
+	    }
+
+	    // re-subtract the object, leave local sky
+	    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 	}
-
-	// iterate to the window radius
-	float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
-
-	// re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-	pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
-
-	// clear the window function for this source based on the moments
-	psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
-	// psphotKronWindowSetSource (source, kronWindow, false, false);
-	// psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
-
-	// this function populates moments->Mrf,KronFlux,KronFluxErr
-	psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
-	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-	// set a window function for each source based on the moments
-	psphotKronWindowSetSource (source, kronWindow, true, true);
-	// psphotKronWindowSetSource (source, kronWindow, false, true);
-
-	// test source fluxes
-        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
-	float kmag = -2.5*log10(source->moments->KronFlux);
-	if (source->psfMag - kmag > 0.25) {
-	    // fprintf (stderr, "continue\n");
-	}
-
-	// re-subtract the object, leave local sky
-	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    }
-
-    psphotSaveImage (NULL, kronWindow, "kron.window.v1.fits");
-
-    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
-	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-	}
-
-	// iterate to the window radius
-	float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
-
-	// re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-	pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
-
-	// clear the window function for this source based on the moments
-	psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
-	// psphotKronWindowSetSource (source, kronWindow, false, false);
-	// psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
-
-	// this function populates moments->Mrf,KronFlux,KronFluxErr
-	psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
-	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-	// set a window function for each source based on the moments
-	psphotKronWindowSetSource (source, kronWindow, true, true);
-	// psphotKronWindowSetSource (source, kronWindow, false, true);
-
-	// test source fluxes
-        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
-	float kmag = -2.5*log10(source->moments->KronFlux);
-	if (source->psfMag - kmag > 0.25) {
-	    // fprintf (stderr, "continue\n");
-	}
-
-	// re-subtract the object, leave local sky
-	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    }
-
-    psphotSaveImage (NULL, kronWindow, "kron.window.v2.fits");
+	char name[64];
+	sprintf (name, "kron.window.v%d.fits", j+1);
+	psphotSaveImage (NULL, kronWindow, name);
+    }
     psFree (kronWindow);
 
@@ -231,4 +203,5 @@
 
     psF32 R2 = PS_SQR(radius);
+    float rsigma2 = 0.5 / PS_SQR(radius/2.0);
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -261,4 +234,10 @@
     int Ywo = source->pixels->row0;
 
+    psF32 **vPix = source->pixels->data.F32;
+    psF32 **vWin = kronWindow->data.F32;
+    psF32 **vWgt = source->variance->data.F32;
+    
+    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
@@ -266,21 +245,25 @@
 	if (fabs(yDiff) > radius) continue;
 
-	psF32 *vPix = source->pixels->data.F32[row];
-	psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
-
-	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-
-	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWin++) {
-	    if (vMsk) {
-		if (*vMsk & maskVal) {
-		    vMsk++;
-		    continue;
-		}
-		vMsk++;
-	    }
-	    if (isnan(*vPix)) continue;
+	// coordinate of mirror pixel
+	int yFlip = yCM - yDiff;
+	if (yFlip < 0) continue;
+	if (yFlip >= source->pixels->numRows) continue;
+
+	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+	    // check mask and value for this pixel
+	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
+	    if (isnan(vPix[row][col])) continue;
 
 	    psF32 xDiff = col - xCM;
 	    if (fabs(xDiff) > radius) continue;
+
+	    // coordinate of mirror pixel
+	    int xFlip = xCM - xDiff;
+	    if (xFlip < 0) continue;
+	    if (xFlip >= source->pixels->numCols) continue;
+
+	    // check mask and value for mirror pixel
+	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+	    if (isnan(vPix[yFlip][xFlip])) continue;
 
 	    // radius is just a function of (xDiff, yDiff)
@@ -289,11 +272,21 @@
 
 	    // flux * window
-	    float weight  = *vWin;
+	    float z = r2 * rsigma2;
+	    assert (z >= 0.0);
+
 	    // float weight  = 1.0;
-	    psF32 pDiff = *vPix * weight;
+	    float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
+	    float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
+	    // float weight1  = vWin[row+Ywo][col+Xwo];
+	    // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
+
+	    float fDiff1 = vPix[row][col]*weight1;
+	    float fDiff2 = vPix[yFlip][xFlip]*weight2;
+
+	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
 
 	    // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
 	    psF32 rf = pDiff * sqrt(r2);
-	    psF32 rs = pDiff;
+	    psF32 rs = 0.5 * (fDiff1 + fDiff2);
 
 	    RF  += rf;
@@ -322,22 +315,25 @@
 	if (fabs(yDiff) > radKron) continue;
 
-	psF32 *vPix = source->pixels->data.F32[row];
-	psF32 *vWgt = source->variance->data.F32[row];
-	psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
-
-	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-
-	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++, vWin++) {
-	    if (vMsk) {
-		if (*vMsk & maskVal) {
-		    vMsk++;
-		    continue;
-		}
-		vMsk++;
-	    }
-	    if (isnan(*vPix)) continue;
+	// coordinate of mirror pixel
+	int yFlip = yCM - yDiff;
+	if (yFlip < 0) continue;
+	if (yFlip >= source->pixels->numRows) continue;
+
+	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+	    // check mask and value for this pixel
+	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
+	    if (isnan(vPix[row][col])) continue;
 
 	    psF32 xDiff = col - xCM;
 	    if (fabs(xDiff) > radKron) continue;
+
+	    // coordinate of mirror pixel
+	    int xFlip = xCM - xDiff;
+	    if (xFlip < 0) continue;
+	    if (xFlip >= source->pixels->numCols) continue;
+
+	    // check mask and value for mirror pixel
+	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+	    if (isnan(vPix[yFlip][xFlip])) continue;
 
 	    // radKron is just a function of (xDiff, yDiff)
@@ -345,12 +341,22 @@
 	    if (r2 > radKron2) continue;
 
-	    float weight  = *vWin;
+	    // float z = r2 * rsigma2;
+	    // assert (z >= 0.0);
+
 	    // float weight  = 1.0;
-	    psF32 pDiff = *vPix * weight;
-	    psF32 wDiff = *vWgt * weight;
+	    // float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
+	    // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
+	    float weight1  = vWin[row+Ywo][col+Xwo];
+	    float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
+
+	    float fDiff1 = vPix[row][col]*weight1;
+	    float fDiff2 = vPix[yFlip][xFlip]*weight2;
+
+	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+	    psF32 wDiff = vWgt[row][col] * weight1;
 
 	    Sum += pDiff;
 	    Var += wDiff;
-	    Win += weight;
+	    Win += weight1;
 	    nKronPix ++;
 	}
@@ -387,7 +393,9 @@
     float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
 
+    // float kratio = source->moments->KronFinner / source->moments->KronFlux;
+
+    float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
     // float scale = useKronRadius ? 2.0 * source->moments->Mrf / Mmajor : 2.0;
-    float scale = 3.0 * source->moments->Mrf / Mmajor;
-    // float scale = 2.0;
+    // float scale = (kratio > 0.4) ? 9.0 * source->moments->Mrf / Mmajor : 3.0 * source->moments->Mrf / Mmajor;
 
     float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2