Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 5802)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 5828)
@@ -1,3 +1,49 @@
 # include "psphot.h"
+
+psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) {
+
+    int ix = x - image->col0;
+    int iy = y - image->row0;
+
+    psF32 *Fm = &image->data.F32[iy - 1][ix];
+    psF32 *Fo = &image->data.F32[iy + 0][ix];
+    psF32 *Fp = &image->data.F32[iy + 1][ix];
+
+    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
+    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
+    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
+    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
+    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
+    poly->mask[2][2] = 1;
+    poly->mask[1][2] = 1;
+    poly->mask[2][1] = 1;
+
+    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
+
+    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
+    poly->coeff[0][1] = (Fyp - Fym)/6.0;
+    
+    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
+    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
+    
+    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
+    
+    return (poly);
+}
+
+psPlane psImageBicubeMin (psPolynomial2D *poly) {
+
+    psPlane min;
+
+    min.xErr = min.yErr = 0;
+
+    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
+
+    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
+    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
+    return (min);
+}
 
 bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
@@ -27,4 +73,15 @@
     index->n = 0;
 
+    bool UseAnalysisRegion = false;
+    psRegion AnalysisRegion;
+    char *region = psMetadataLookupPtr (&status, config, "ANALYSIS_REGION");
+    if (status) {
+	UseAnalysisRegion = true;
+	AnalysisRegion = psRegionFromString (region);
+	psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", 
+		  AnalysisRegion.x0, AnalysisRegion.y0, 
+		  AnalysisRegion.x1, AnalysisRegion.y1);
+    }
+
     for (int i = 0; i < sources->n; i++) {
 	pmSource *inSource = sources->data[i];
@@ -35,4 +92,11 @@
 	if (inSource->type == PM_SOURCE_DEFECT) continue; 
 	if (inSource->type == PM_SOURCE_SATURATED) continue;
+
+	if (UseAnalysisRegion) {
+	    if (inSource->moments->x < AnalysisRegion.x0) continue;
+	    if (inSource->moments->y < AnalysisRegion.y0) continue;
+	    if (inSource->moments->x > AnalysisRegion.x1) continue;
+	    if (inSource->moments->y > AnalysisRegion.y1) continue;
+	}
 
 	pmSource *otSource = pmSourceAlloc ();
@@ -51,6 +115,23 @@
 	    modelFLT->params->data.F32[2] = inSource->moments->x;
 	    modelFLT->params->data.F32[3] = inSource->moments->y;
-	}
-	// XXX EAM : add option to peak-up on peak (for non-sat objects)
+	} else {
+	    // peak-up on peak (for non-sat objects)
+
+	    int ix = inSource->peak->x;
+	    int iy = inSource->peak->y;
+
+	    psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
+	    psPlane min = psImageBicubeMin (bicube);
+
+	    psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", 
+		     modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy);
+	    
+	    // if min point is too deviant, keep the old value
+	    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
+		modelFLT->params->data.F32[2] = min.x + ix;
+		modelFLT->params->data.F32[3] = min.y + iy;
+	    }
+	    psFree (bicube);
+	}
 
 	// set PSF parameters for this model
@@ -149,9 +230,11 @@
 	// subtract object
 	pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
+	Fi->mode |= PM_SOURCE_SUBTRACTED;
+	Fi->mode |= PM_SOURCE_TEMPSUB;
     }
 
     // XXX EAM : need to free up many things here
 
-    psLogMsg ("psphot.emsemble", 4, "apply models: %f\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f\n", psTimerMark ("psphot"));
     return true;
 }