Changeset 25275 for branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.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/psphotPetrosianRadialBins.c (modified) (4 diffs)

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

@@ -13 +13 @@
 // track the non-overlapping radius values.
 
-bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
+bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise) {
+
+    // XXX for testing, let's just set this to a value
+    
+    float skyModelErrorSQ = PS_SQR(skynoise);
+
+# if (0)
+    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
+    pmReadout *backModel = psphotSelectBackground (config, view);
+    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
+# endif
 
     psVector *radius = petrosian->radiusElliptical;
@@ -104 +114 @@
     psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    // psStats *stats = psStatsAlloc(PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4);
 
     // integrate flux, radius for each of these bins.  since flux is sorted by radius, 
@@ -123 +134 @@
             // binSB->data.F32[nOut] = stats->sampleMedian;
             binSB->data.F32[nOut] = stats->robustMedian;
-            binSBstdev->data.F32[nOut] = stats->robustStdev / sqrt(values->n);
+            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->robustStdev) / values->n + skyModelErrorSQ);
+            // binSB->data.F32[nOut] = stats->fittedMean;
+            // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ);
+
+            // error in the SB is the stdev per bin / sqrt (number of pixels) 
+            // added in quadrature to a fraction of the local sky (not the 
+            // residual flux, but the sky from the sky model)
 
             if (1) {
@@ -153 +170 @@
     petrosian->binSBstdev = binSBstdev;
 
-    psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, 0.0);
+    psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0);
 
     psFree(radMin);