Changeset 32316

Timestamp:

Sep 5, 2011, 8:49:42 AM (15 years ago)

Author:

eugene

Message:

deal with the cases of small radii (may be an error in the PetRadiusErr)

File:

• branches/eam_branches/ipp-20110710/psphot/src/psphotPetrosianStats.c (modified) (4 diffs)

branches/eam_branches/ipp-20110710/psphot/src/psphotPetrosianStats.c

@@ -143 +143 @@
     // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!)
     if (!anyPetro) {
-        // 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);
-            petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]);
-
-        } else {
-            int n0 = lowestSignificantRadius-1;
-            int n1 = lowestSignificantRadius;
-            petRadius    = InterpolateValues     (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO);
-            petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]);
+        petRadius = refRadius->data.F32[lowestSignificantRadius];
+        petRadiusErr = NAN;
+        if (!isfinite(petRadius)) {
+            fprintf (stderr, "nan pet radius\n");
         }
     }
@@ -170 +163 @@
                 petArea    = InterpolateValues     (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius);
                 petApix    = InterpolateValues     (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius);
+                if (!isfinite(petFlux)) {
+                    fprintf (stderr, "nan pet flux\n");
+                }
                 break;
             }
@@ -188 +184 @@
         if (!found50 && (fluxSum->data.F32[i] > flux50)) {
             if (i == 0) {
-                psWarning ("does this case make any sense? (fluxSum[0] %f > flux50 %f)", fluxSum->data.F32[i], flux50);
-                continue;
+                R50    = InterpolateValues     (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux50);
+                R50err = InterpolateValuesErrX (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux50, sqrt(fluxSumErr2->data.F32[i]), sqrt(fluxSumErr2->data.F32[i+1]));
+                found50 = true;
             } else {
                 R50    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
@@ -198 +195 @@
         if (!found90 && (fluxSum->data.F32[i] > flux90)) {
             if (i == 0) {
-                psWarning ("does this case make any sense? (fluxSum[0] %f > flux90 %f)", fluxSum->data.F32[i], flux90);
-                continue;
+                R90    = InterpolateValues     (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux90);
+                R90err = InterpolateValuesErrX (fluxSum->data.F32[i], refRadius->data.F32[i], fluxSum->data.F32[i+1], refRadius->data.F32[i+1], flux90, sqrt(fluxSumErr2->data.F32[i]), sqrt(fluxSumErr2->data.F32[i+1]));
+                found90 = true;
             } else {
                 R90    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);