Index: /trunk/psphot/src/psphotSourceSize.c
===================================================================
--- /trunk/psphot/src/psphotSourceSize.c	(revision 13081)
+++ /trunk/psphot/src/psphotSourceSize.c	(revision 13081)
@@ -0,0 +1,78 @@
+# include "psphotInternal.h"
+
+// we need to call this function after sources have been fitted to the PSF
+// model and subtracted.  this function examines the 9 pixels in the 3x3 square containing the peak
+// and compares the observed flux to the model.  For the 4 pixel lines which contain the peak, we
+// calculate the value 
+
+/*
+ * 
+ *
+ *
+ *
+ *
+ */
+
+bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) {
+
+    // loop over all source
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+
+	// source must have been subtracted
+	if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
+
+	psF32 *resid = source->pixels->data.F32;
+	psF32 *weight = source->weight->data.F32;
+	psU8 *mask = source->maskObj->data.U8;
+
+	int xPeak = source->peak->xf - source->pixels->col0;
+	int yPeak = source->peak->yf - source->pixels->row0;
+
+	// XXX need to deal with edge peaks...
+	float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
+	float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0];
+	float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1];
+	float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];
+
+	float dcX = 2*weight[yPeak][xPeak] - weight[yPeak+0][xPeak-1] - weight[yPeak+0][xPeak+1];
+	float dcY = 2*weight[yPeak][xPeak] - weight[yPeak+1][xPeak+0] - weight[yPeak+1][xPeak+0];
+	float dcL = 2*weight[yPeak][xPeak] - weight[yPeak-1][xPeak-1] - weight[yPeak+1][xPeak+1];
+	float dcR = 2*weight[yPeak][xPeak] - weight[yPeak+1][xPeak-1] - weight[yPeak-1][xPeak+1];
+
+	// XXX plot this as a function of flux / magnitude?
+	// XXX overlay outliers on an image?
+	
+	// some possible classes:
+	// all 4 curvatures are highly negative : CR
+	// all 4 curvatures are highly positive : EXT
+
+	// at least 2 are significantly negative, none are significantly positive : CR
+	// at least 2 are significantly positive, none are significantly negative : EXT
+
+	// any are significantly negative, some may be significantly positive : CR
+	// any are significantly positive, none may be significantly positive : EXT
+	
+	/* Nn  Np  No
+	   4   0   0   CR_1
+	   3   1   0   CR_1
+	   3   0   1   CR_1
+	   2   2   0   CR_2
+	   2   1   1   CR_2
+	   2   0   2   CR_2
+	   1   3   0   CR_3
+	   1   2   1   CR_3
+	   1   1   2   CR_3
+	   1   0   3   CR_3
+	   0   4   0   EXT
+	   0   3   1   EXT
+	   0   2   2   EXT
+	   0   1   3   PSF
+	   0   0   4   PSF
+	*/
+
+	/* I can write the formal probability that the 4 measurements are consistent with a PSF
+	 * based on how large the error values are (Nsigma -> erf -> P(Nsigma)).
+	 * I should examine this value as a function of flux and also examine the distribution of
+	 * the probabilities for the PSF sources.
+	 */