Changeset 17335

Timestamp:

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

Author:

eugene

Message:

fix algorithms

Location:

branches/eam_branch_20080324/psphot/src

Files:

• psphotAnnuli.c (modified) (2 diffs)
• psphotIsophotal.c (modified) (4 diffs)
• psphotKron.c (modified) (1 diff)
• psphotPetrosian.c (modified) (5 diffs)
• psphotRadialProfile.c (modified) (1 diff)

branches/eam_branch_20080324/psphot/src/psphotAnnuli.c

@@ -1 @@
-# include "psphot.h"
+# include "psphotInternal.h"
 
 bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
@@ -55 +55 @@
   source->extpars->annuli->fluxVar = fluxSquare;
 
+  psFree (pixelCount);
+
   return true;
 }

branches/eam_branch_20080324/psphot/src/psphotIsophotal.c

@@ -1 @@
-# include "psphot.h"
+# include "psphotInternal.h"
+
+static float ISOPHOT_FLUX = NAN;
 
 bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
@@ -14 +16 @@
 
   // flux at which to measure isophotal parameters
-  // XXX cache this?
-  float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
+  if (!isfinite(ISOPHOT_FLUX)) {
+    // XXX ISOPHOTAL_FLUX should be specified in mags, need the zero point to get counts/sec
+    ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
+    assert (status);
+  }
 
   // find the first bin below the flux level and the last above the level
-  // XXX can this be done faster with disection?
+  // XXX can this be done faster with bisection?
   // XXX do I need to worry about crazy outliers?  
   // XXX should i be smoothing or fitting the curve?
@@ -27 +32 @@
     if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i;
   }
-
+  // if we don't go out far enough, we have a problem...
+  if (lastAbove == flux->n - 1) {
+    psTrace ("psphot", 5, "did not go out far enough to reach isophotal magnitude");
+    // XXX raise a flag ?
+    return false;
+  }
+  if (firstBelow < 0) {
+    psTrace ("psphot", 5, "did not go out far enough to bound isophotal magnitude: error unmeasured");
+    // XXX raise a flag ?
+    lastAbove = firstBelow;
+    return false;
+  }
+  
   // need to examine pixels in this vicinity
   float isophotalFluxFirst = 0;
@@ -56 +73 @@
   source->extpars->isophot->radErr = isophotalRadErr;
 
+  psTrace ("psphot", 5, "Isophot flux:%f +/- %f @ %f +/- %f for %f, %f\n",
+           source->extpars->isophot->mag, source->extpars->isophot->magErr,
+           source->extpars->isophot->rad, source->extpars->isophot->radErr,
+           source->peak->xf, source->peak->yf);
+
   return true;
 

branches/eam_branch_20080324/psphot/src/psphotKron.c

@@ -1 @@
-# include "psphot.h"
+# include "psphotInternal.h"
 
 bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal) {

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;
 

branches/eam_branch_20080324/psphot/src/psphotRadialProfile.c

@@ -41 +41 @@
 
     int n = 0;
-    float Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
-    float Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+    
+    float Xo = 0.0;
+    float Yo = 0.0;
+
+    if (source->modelEXT) {
+      Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
+      Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+    } else {
+      Xo = source->peak->xf - source->pixels->col0;
+      Yo = source->peak->yf - source->pixels->row0;
+    }
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
         for (int ix = 0; ix < source->pixels->numCols; ix++) {