Changeset 17335 for branches/eam_branch_20080324/psphot/src/psphotPetrosian.c

Timestamp:

Apr 6, 2008, 10:11:39 AM (18 years ago)

Author:

eugene

Message:

fix algorithms

File:

• branches/eam_branch_20080324/psphot/src/psphotPetrosian.c (modified) (5 diffs)

branches/eam_branch_20080324/psphot/src/psphotPetrosian.c

@@ -1 @@
-# include "psphot.h"
+# include "psphotInternal.h"
+
+static float PETROSIAN_R0 = NAN;
+static float PETROSIAN_RF = NAN;
 
 bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
@@ -14 +17 @@
 
   // flux at which to measure isophotal parameters
-  // XXX cache this?
-  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
-  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
+  if (!isfinite(PETROSIAN_R0)) {
+    PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
+    PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
+    assert (status);
+  }
 
   // first find flux at R0
-  int firstBelow = -1;
-  int lastAbove = -1;
+  int firstAbove = -1;
+  int lastBelow = -1;
   for (int i = 0; i < radius->n; i++) {
-    if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
-    if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i;
+    if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i;
+    if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i;
+  }
+  // if we don't go out far enough, we have a problem...
+  if (lastBelow == radius->n - 1) {
+    psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius...");
+    // XXX skip object? raise a flag ?
+    return false;
+  }
+  if (firstAbove < 0) {
+    psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius");
+    // XXX raise a flag ?
+    return false;
   }
 
@@ -29 +45 @@
   float fluxR0 = 0.0;
   int fluxRn = 0;
-  for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) {
+  for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) {
     fluxR0 += flux->data.F32[i];
     fluxRn ++;
@@ -42 +58 @@
   // XXX do I need to worry about crazy outliers?  
   // XXX should i be smoothing or fitting the curve?
-  firstBelow = -1;
-  lastAbove = -1;
+  int firstBelow = -1;
+  int lastAbove = -1;
   for (int i = 0; i < flux->n; i++) {
     if (flux->data.F32[i] > fluxRP) lastAbove = i;
     if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
+  }
+  // if we don't go out far enough, we have a problem...
+  if (lastAbove == radius->n - 1) {
+    psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius...");
+    // XXX skip object? raise a flag ?
+    return false;
+  }
+  if (firstBelow < 0) {
+    psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius");
+    // XXX raise a flag ?
+    return false;
   }
 
@@ -78 +105 @@
   source->extpars->petrosian->radErr = radErr;
 
+  psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n",
+           source->extpars->petrosian->mag, source->extpars->petrosian->magErr,
+           source->extpars->petrosian->rad, source->extpars->petrosian->radErr,
+           source->peak->xf, source->peak->yf);
+
   return true;