Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c
===================================================================
--- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c	(revision 25351)
+++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c	(revision 25352)
@@ -245,4 +245,6 @@
         source->modelPSF->radiusFit = options->radius;
 
+	// XXXX use a different radius for the aperture magnitude than for the PSF fit?
+
         // set object mask to define valid pixels
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
@@ -422,5 +424,23 @@
               psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
 
-
+    float dSys = psVectorSystematicError (psfTry->metric, psfTry->metricErr, 0.1);
+    fprintf (stderr, "systematic error: %f\n", dSys);
+
+    int n = 0;
+    psVector *bright = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32);
+    psVector *brightErr = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32);
+    for (int i = 0; i < psfTry->metric->n; i++) {
+	if (!isfinite(psfTry->metric->data.F32[i])) continue;
+	if (!isfinite(psfTry->metricErr->data.F32[i])) continue;
+	if (psfTry->metricErr->data.F32[i] <= 0.0) continue;
+	if (psfTry->metricErr->data.F32[i] > 0.005) continue;
+	bright->data.F32[n] = psfTry->metric->data.F32[i];
+	brightErr->data.F32[n] = psfTry->metricErr->data.F32[i];
+	n++;
+    }
+    bright->n = brightErr->n = n;
+
+    float dSysBright = psVectorSystematicError (bright, brightErr, 0.1);
+    fprintf (stderr, "bright systematic error: %f\n", dSysBright);
 
     // XXX test dump of fitted model (dump when tracing?)
@@ -582,5 +602,4 @@
 }
 
-
 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
 
@@ -596,4 +615,6 @@
     // parameters, with the constraint that the minor axis must be greater than a minimum
     // threshold.
+
+    // XXX re-read the sextractor manual on handling 'infinitely thin' sources...
 
     // convert the measured source shape paramters to polarization terms
@@ -1069,2 +1090,93 @@
     return true;
 }
+
+float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) {
+
+    psAssert(residuals, "residuals cannot be NULL");
+    psAssert(errors, "errors cannot be NULL");
+    psAssert(residuals->n == errors->n, "residuals and errors must be the same length");
+
+    // given a vector of residuals and their formal errors, calculated the necessary systematic
+    // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction
+    // highest chi-square contributors (allowed outliers)
+
+    psVector *mask  = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK);
+    psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32);
+
+    // calculate the chisq vector:
+    int Ngood = 0;
+    for (int i = 0; i < residuals->n; i++) {
+	chisq->data.F32[i] = PS_MAX_F32;
+	if (!isfinite(residuals->data.F32[i])) continue;
+	if (!isfinite(errors->data.F32[i])) continue;
+	if (errors->data.F32[i] <= 0.0) continue;
+	chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]);
+	Ngood ++;
+    }
+
+    psVector *index = psVectorSortIndex(NULL, chisq);
+
+    // toss out the clipFraction highest chisq values
+    for (int i = 0; i < residuals->n; i++) {
+	int n = index->data.S32[i];
+	if (i < (1.0 - clipFraction)*Ngood) {
+	    mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0;
+	} else {
+	    mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1;
+	}
+    }
+
+    // Ndof ~= Ngood
+    // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof
+    // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0
+    
+    // use Newton-Raphson to solve for S2:
+
+    // use median sigma to calculate the initial guess for S2:
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, errors, NULL, mask, 1);
+    float errorMedian = stats->sampleMedian;
+    
+    float nPts = 0.0;
+    float res2mean = 0.0;
+    float ChiSq = 0.0;
+    for (int i = 0; i < residuals->n; i++) {
+	int n = index->data.S32[i];
+	if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
+	res2mean += PS_SQR(residuals->data.F32[n]);
+	ChiSq += PS_SQR(residuals->data.F32[n] / errors->data.F32[n]);
+	nPts += 1.0;
+    }
+    res2mean /= nPts;
+    ChiSq /= nPts;
+    
+    float S2guess = res2mean - PS_SQR(errorMedian);
+
+    psLogMsg ("psModules", 3, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n", 
+	      ChiSq, residuals->n, Ngood, nPts, S2guess);
+
+    for (int iter = 0; iter < 10; iter++) {
+
+	ChiSq = 0.0;
+	float dRdS = 0.0;
+	for (int i = 0; i < residuals->n; i++) {
+	    int n = index->data.S32[i];
+	    if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
+	    float error2 = PS_SQR(errors->data.F32[n]) + S2guess;
+	    ChiSq += PS_SQR(residuals->data.F32[n]) / error2;
+	    dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2);
+	}
+	ChiSq /= nPts;
+	dRdS /= nPts;
+
+	// Note the sign on dS: dRdS above is -1 * dR/dS formally
+	float dS = (ChiSq - 1.0) / dRdS;
+	S2guess += dS;
+
+	psLogMsg ("psModules", 3, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess);
+    }
+
+    // free local allocations
+    return (sqrt(S2guess));
+}
+