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