Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c	(revision 34408)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c	(revision 34409)
@@ -1,4 +1,6 @@
 # include "psphotInternal.h"
 bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
+bool psphotBlendSetSourceSatstar (psImage *image, pmSource *source);
+float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -385,5 +387,6 @@
     if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
 	// find the radius at which we hit something like 0.1*SATURATION
-	return false;
+	psphotBlendSetSourceSatstar (IDimage, source);
+	return true;
     }
 
@@ -431,6 +434,6 @@
     float rMxx = 0.5 / PS_SQR(shape.sx);
     float rMyy = 0.5 / PS_SQR(shape.sy);
-    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
-                                    // implementation. XXX: Is this correct?
+    float Sxy = shape.sxy;    // factor of -1 is included to match the previous window function
+    // implementation. XXX: Is this correct?
 
     int ID = source->id;
@@ -443,4 +446,5 @@
 
 	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+	    if (z > 2.311) continue;
 
 	    float f = Io*exp(-z);
@@ -453,2 +457,71 @@
     return true;
 }
+
+bool psphotBlendSetSourceSatstar (psImage *IDimage, pmSource *source) {
+
+    float Xo = source->satstar->Xo;
+    float Yo = source->satstar->Yo;
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+
+    float lPeak = logFmodel->data.F32[0];
+    float fPeak = pow(10.0, lPeak);
+    float threshold = 0.1*fPeak;
+    float logThresh = log10(threshold);
+
+    float radius = NAN;
+
+    // what is the radius for a specific peak fraction?
+    for (int i = 0; i < logFmodel->n; i++) {
+	float logF = logFmodel->data.F32[i];
+	float flux = pow(10.0, logF);
+	if (flux > threshold) continue;
+	if (i == 0) continue;
+	
+	float logF0 = logFmodel->data.F32[i - 1];
+	float logR0 = logRmodel->data.F32[i - 1];
+
+	float logF1 = logFmodel->data.F32[i];
+	float logR1 = logRmodel->data.F32[i];
+
+	float logR = InterpolateValues (logF0, logR0, logF1, logR1, logThresh);
+	radius = pow(10.0, logR);
+	break;
+    }
+
+    if (!isfinite(radius)) {
+	for (int i = logFmodel->n - 1; i >= 0; i--) {
+	    if (!isfinite(logFmodel->data.F32[i])) continue;
+	    float logR = logRmodel->data.F32[i];
+	    radius = pow(10.0, logR);
+	    break;
+	}
+    }
+
+    if (!isfinite(radius)) return false;
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    // region to mask
+    int minX = PS_MIN(PS_MAX(Xo - radius, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + radius, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - radius, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + radius, 0), Ny - 1);
+
+    int ID = source->id;
+
+    for (int iy = minY; iy < maxY; iy++) {
+	for (int ix = minX; ix < maxX; ix++) {
+
+	    float dX = (ix - Xo);
+	    float dY = (iy - Yo);
+	    float R = hypot (dX, dY) ;
+	    if (R > radius) continue;
+
+            IDimage->data.S32[iy][ix] = ID;
+	}
+    }
+    return true;
+}
+