Changeset 25275 for branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c

Timestamp:

Sep 4, 2009, 2:18:54 PM (17 years ago)

Author:

eugene

Message:

updates to petrosian analysis

File:

• branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (modified) (8 diffs)

branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c

@@ -32 +32 @@
     float dFsum2 = 0.0;
 
+    float nSigma = 3.0;
+    int lowestSignificantRadius = 0;
+    float lowestSignificantRatio = 1.0;
+
     int nOut = 0;
     for (int i = 0; i < binSB->n; i++) {
@@ -48 +52 @@
             areaInner += areaSum->data.F32[nOut-1];
             fluxInner += fluxSum->data.F32[nOut-1];
-            fluxInnerErr2 = PS_SQR(fluxSumErr2->data.F32[nOut-1] * areaSum->data.F32[nOut-1]);
+            fluxInnerErr2 += fluxSumErr2->data.F32[nOut-1];
         }
 
@@ -64 +68 @@
 
         float meanSBerr = sqrt(fluxInnerErr2) / areaInner;
-        float ratioErr = ratio * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut]));
+        float ratioErr = fabs(ratio) * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut]));
 
         psVectorAppend(petRatioErr, ratioErr);
@@ -74 +78 @@
 
         if (1) {
-            fprintf (stderr, "%3d : %5.2f : %5.3f  %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", 
-                     i, refRadius->data.F32[nOut], binSB->data.F32[i], meanSB->data.F32[nOut], petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], fluxSum->data.F32[nOut], areaSum->data.F32[nOut]);
+            fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", 
+                     i, refRadius->data.F32[nOut], 
+                     binSB->data.F32[i], binSBstdev->data.F32[i], 
+                     meanSB->data.F32[nOut], meanSBerr, 
+                     petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], 
+                     fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
         }
     
         // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
         // we will keep and report the last (largest radius) value
-        if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) {
+        if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) {
             // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
             if (i == 0) { 
@@ -93 +101 @@
         }
     
+        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
+        // we will keep and report the last (largest radius) value
+        // find the last signficant measurement of the petrosian ratio
+        if (above && (petRatio->data.F32[nOut] < lowestSignificantRatio) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) {
+            lowestSignificantRadius = nOut;
+            lowestSignificantRatio = petRatio->data.F32[nOut];
+        }
+    
         // reset on transitions up, but do not re-calculate rad_90, flux_90
         if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) {
@@ -101 +117 @@
 
     if (!anyPetro) {
-        // set default radius:
+        // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
+        if (lowestSignificantRadius == 0) { 
+            // assume Fmax @ R = 0.0
+            petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
+        } else {
+            petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
+        }
     }
 
@@ -110 +132 @@
         if (refRadius->data.F32[i] > apRadius) {
             if (i == 0) {
-                psAbort ("does this case make any sense?");
+                psWarning ("does this case make any sense? (refRadius[0] > apRadius)");
+                continue;
             } else {
                 petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
@@ -119 +142 @@
 
     // save petRadius, petFlux
+    // XXX save flags (anyPetro, manyPetro)
     petrosian->petrosianRadius = petRadius;
     petrosian->petrosianFlux   = petFlux;