Index: trunk/psphot/src/psphotPetrosianStats.c
===================================================================
--- trunk/psphot/src/psphotPetrosianStats.c	(revision 34086)
+++ trunk/psphot/src/psphotPetrosianStats.c	(revision 36375)
@@ -6,4 +6,5 @@
 // generate the Petrosian radius and flux from the mean surface brightness (r_i)
 
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X);
 float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
 float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
@@ -20,4 +21,5 @@
     if (!profile->binSB) {
 	psLogMsg ("psphot", PS_LOG_DETAIL, "no petrosian profile, skipping source %f, %f", source->peak->xf, source->peak->yf);
+	source->mode2 |= PM_SOURCE_MODE2_PETRO_NO_PROFILE;
 	return true;
     }
@@ -53,5 +55,5 @@
     float dFsum2 = 0.0;
 
-    float nSigma = 3.0;
+    float nSigma = 2.0;
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
@@ -117,7 +119,19 @@
 		petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
 		petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
-	    } else {
-		petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
-		petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
+		source->mode2 |= PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN;
+	    } else {
+	      // petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+	      if (nOut > 1) {
+		petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-2], &refRadius->data.F32[nOut-2],   PETROSIAN_RATIO);
+	      } else {
+		petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-3], &refRadius->data.F32[nOut-3], PETROSIAN_RATIO);
+	      }
+# if (PS_TRACE_ON)
+	      float petRadiusLinear = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+	      if (fabs(petRadius - petRadiusLinear) > fabs(refRadius->data.F32[nOut] - refRadius->data.F32[nOut-1])) { 
+		fprintf (stderr, "big difference : %f vs %f\n", petRadius, petRadiusLinear); 
+	      }
+#endif
+	      petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
 	    }
 	    above = false;
@@ -148,4 +162,5 @@
 	    fprintf (stderr, "nan pet radius\n");
 	}
+	source->mode2 |= PM_SOURCE_MODE2_PETRO_INSIG_RATIO;
     }
 
@@ -180,4 +195,5 @@
     bool found50 = false;
     bool found90 = false;
+
     // XXX use bisection to do this faster:
     for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
@@ -235,4 +251,23 @@
 }
 
+// Lagrange's form of the interpolating polynomial...
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X) {
+
+  float dx01 = Xin[0] - Xin[1];
+  float dx02 = Xin[0] - Xin[2];
+  float dx12 = Xin[1] - Xin[2];
+
+  float dx0  = X - Xin[0];
+  float dx1  = X - Xin[1];
+  float dx2  = X - Xin[2];
+
+  float y0 = Yin[0]*dx1*dx2/(dx01*dx02);
+  float y1 = Yin[1]*dx0*dx2/(dx01*dx12); // need - sign
+  float y2 = Yin[2]*dx0*dx1/(dx02*dx12); 
+
+  float Y = y0 - y1 + y2;
+  return Y;
+}
+
 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
     float dydx = (Y1 - Y0) / (X1 - X0);