Changeset 28017 for branches/czw_branch/20100427/psphot/src/psphotPetrosianStats.c

Timestamp:

May 18, 2010, 5:06:01 PM (16 years ago)

Author:

watersc1

Message:

Almost working version of the mask-stats/software revision code.

Almost.

Location:

branches/czw_branch/20100427

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotPetrosianStats.c (modified) (4 diffs)

branches/czw_branch/20100427/psphot/src/psphotPetrosianStats.c

@@ -10 +10 @@
 bool psphotPetrosianStats (pmSource *source) {
 
-    pmSourceRadialProfile *profile = source->extpars->profile;
-
-    float petRadius = NAN;
-    float petFlux = NAN;
+    psAssert (source, "missing source");
+    psAssert (source->extpars, "missing extpars");
+    psAssert (source->extpars->petProfile, "missing petProfile");
+
+    pmSourceRadialProfile *profile = source->extpars->petProfile;
 
     psVector *binSB      = profile->binSB;
@@ -28 +29 @@
     psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
 
+    float petRadius = NAN;
+    float petFlux = NAN;
+
     bool anyPetro = false;
     bool manyPetro = false;
@@ -38 +42 @@
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
+
+    // find the Petrosian Radius and Petrosian Flux
 
     int nOut = 0;
@@ -142 +148 @@
     }
 
-    if (!source->extpars->petrosian_80) {
-        source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc ();
-    }
-    pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
+    // now measure the radii R90 and R50 where flux = 0.9 (or 0.5) * petFlux;
+    float flux90 = 0.9 * petFlux;
+    float flux50 = 0.5 * petFlux;
+    float R50 = NAN;
+    float R90 = NAN;
+    bool found50 = false;
+    bool found90 = false;
+    // XXX use bisection to do this faster:
+    for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
+        if (!found50 && (fluxSum->data.F32[i] > flux50)) {
+            if (i == 0) {
+                psWarning ("does this case make any sense? (fluxSum[0] > flux50)");
+                continue;
+            } else {
+                R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
+                found50 = true;
+            }
+        }
+        if (!found90 && (fluxSum->data.F32[i] > flux90)) {
+            if (i == 0) {
+                psWarning ("does this case make any sense? (fluxSum[0] > flux90)");
+                continue;
+            } else {
+                R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
+                found90 = true;
+            }
+        }
+    }
+
 
     // XXX save flags (anyPetro, manyPetro)
-    petrosian->radius = petRadius;
-    petrosian->flux   = petFlux;
-
-    // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
+    source->extpars->petrosianRadius = petRadius;
+    source->extpars->petrosianFlux   = petFlux;
+    source->extpars->petrosianR50    = R50;
+    source->extpars->petrosianR90    = R90;
+    
+    // XXX add the errors
+    source->extpars->petrosianRadiusErr = NAN;
+    source->extpars->petrosianFluxErr   = NAN;
+    source->extpars->petrosianR50Err    = NAN;
+    source->extpars->petrosianR90Err    = NAN;
+
+    fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
+    psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
 
     psFree(fluxSum);