Changeset 36375 for trunk/psphot/src/psphotPetrosianStats.c

Timestamp:

Dec 10, 2013, 2:55:11 PM (13 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20130904

Location:

trunk

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotPetrosianStats.c (modified) (7 diffs)

trunk/psphot/src

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff

trunk/psphot/src/psphotPetrosianStats.c

@@ -6 +6 @@
 // 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 +21 @@
     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 +55 @@
     float dFsum2 = 0.0;
 
-    float nSigma = 3.0;
+    float nSigma = 2.0;
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
@@ -117 +119 @@
                 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 +162 @@
             fprintf (stderr, "nan pet radius\n");
         }
+        source->mode2 |= PM_SOURCE_MODE2_PETRO_INSIG_RATIO;
     }
 
@@ -180 +195 @@
     bool found50 = false;
     bool found90 = false;
+
     // XXX use bisection to do this faster:
     for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
@@ -235 +251 @@
 }
 
+// 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);

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff