Changeset 25105 for branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c

Timestamp:

Aug 18, 2009, 6:24:53 AM (17 years ago)

Author:

eugene

Message:

updates for petrosian analysis study

Location:

branches/eam_branches/20090715/psphot/src

Files:

• . (modified) (1 prop)
• psphotPetrosianRadialBins.c (modified) (2 diffs)

branches/eam_branches/20090715/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

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

@@ -2 +2 @@
 
 // convert the flux vs elliptical radius to annular bins
-bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian) {
 
-    psVector *radius = petrosian->radiusElli;
-    psVector *flux = petrosian->fluxElli;
+// we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels
+// themselves.  this function does not attempt to measure the radial profiles accurately
+// for radii that are smaller than a minimum (currently 1.0 pixels).  
 
-    sort (radius, flux);
+// for small radii, we are measuring the mean surface brightness in non-overlapping radial
+// bins.  for large radii (r > 2 pixels), we are measuring the surface brightness for a
+// radius range \alpha r_i < i < \beta r_i, but performing this measurement for radii more
+// finely spaced than r_{i+1} = r_i * \beta / \alpha.  for the integration, we need to
+// track the non-overlapping radius values.
+
+# define PETROSIAN_ALPHA 0.8
+# define PETROSIAN_BETA 1.25
+
+bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
+
+    psVector *radius = petrosian->radiusElliptical;
+    psVector *flux = petrosian->fluxElliptical;
+
+    // sort incoming vectors by radius
+    psphotPetrosianSortPair (radius, flux);
+
+    int nMax = radiusMax;
 
     // radBin stores the centers of the radial bins, 
     // radMin, radMax store the bounds
-    psVector *radBin = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *radMin = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *radMax = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *radMin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radMax  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radAlp  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radBet  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+
+    psVector *fluxBin = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radBin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *area    = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+
+    psVectorInit (fluxBin, 0.0);
+    psVectorInit (radBin, 0.0);
 
     // generate radial bin bounds
     radMin->data.F32[0] = 0.0;
     radMax->data.F32[0] = 1.0;
+    radAlp->data.F32[0] = 0.0;
+    radBet->data.F32[0] = 1.0;
     
     radMin->data.F32[1] = 1.0;
     radMax->data.F32[1] = 1.5;
+    radAlp->data.F32[1] = 1.0;
+    radBet->data.F32[1] = 1.5;
     
     radMin->data.F32[2] = 1.5;
     radMax->data.F32[2] = 2.0;
+    radAlp->data.F32[2] = 1.5;
+    radBet->data.F32[2] = 2.0;
     
-    for (int i = 3; i < rmax; i++) {
-        radMin->data.F32[i] = (i - 1);
-        radMax->data.F32[i] = i;
+    int nPts = 3;
+    for (int i = 3; i < radiusMax; i++) {
+        radMin->data.F32[nPts] = (i - 1);
+        radMax->data.F32[nPts] = i;
+        nPts++;
     }
+    radMin->n = radMax->n = radAlp->n = radBet->n = nPts;
 
     // generate radial area-weighted mean radius & non-overlapping areas
-    float Rprev = 0.0;
     for (int i = 0; i < radMin->n; i++) {
         float rMin = radMin->data.F32[i];
@@ -42 +75 @@
         float rMax3 = rMax2*rMax;
 
-        float rBin = (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0;
+        float rBin = 2.0 * (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0;
         
+        // XXX calculate area-weighted radius rather than asserting?
         radBin->data.F32[i] = rBin;
+        area->data.F32[i] = M_PI * (rMax2 - rMin2);
 
-        area->data.F32[i] = 
+        if (i > 2) {
+            radAlp->data.F32[i] = rBin*PETROSIAN_ALPHA;
+            radBet->data.F32[i] = rBin*PETROSIAN_BETA;
+        }
     }
+
+    // storage vector for stats
+    psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+
+    // integrate flux, radius for each of these bins.  since flux is sorted by radius, 
+    // we can do this fairly quickly
+
+    bool done = false;
+    int nOut = 0;
+    float Rmin = radAlp->data.F32[nOut];
+    float Rmax = radBet->data.F32[nOut];
+    float Rnxt = radAlp->data.F32[nOut+1];  // minimum radius for next range
+    int iNext = 0;
+    for (int i = 0; !done && (i < radius->n); i++) {
+        if (radius->data.F32[i] < Rnxt) {
+          iNext = i;
+        }
+        if (radius->data.F32[i] > Rmax) {
+            // calculate the value for the nOut bin
+            psVectorStats (stats, values, NULL, NULL, 0);
+            fluxBin->data.F32[nOut] = stats->sampleMedian;
+            nOut ++;
+            if (nOut >= nMax) break;
+            Rmin = radAlp->data.F32[nOut];
+            Rmax = radBet->data.F32[nOut];
+            Rnxt = (nOut < nMax - 1) ? radAlp->data.F32[nOut+1] : Rmax;  // minimum radius for next range
+            values->n = 0;
+            psStatsInit(stats);
+            i = iNext;
+        }
+        if (radius->data.F32[i] < Rmin) {
+            continue;
+        }
+        psVectorAppend (values, flux->data.F32[i]);
+    }
+    fluxBin->n = radBin->n = area->n = nOut;
+    // XXX I think this misses the last radial bin -- do we care?
+
+    // save the vectors
+    petrosian->radialBins = radBin;
+    petrosian->area = area;
+    petrosian->binnedFlux = fluxBin;
+
+    psphotPetrosianVisualProfileRadii (radius, flux, radBin, fluxBin, 0.0);
+
+    psFree(radMin);
+    psFree(radMax);
+    psFree(radAlp);
+    psFree(radBet);
+    psFree(values);
+    psFree(stats);
+
+    return true;
 }
 

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy