Changeset 25033 for branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c

Timestamp:

Aug 10, 2009, 3:43:50 AM (17 years ago)

Author:

eugene

Message:

more tests of robust petrosian analysis

File:

• branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c (modified) (3 diffs)

branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c

@@ -5 +5 @@
   int i, above;
   float *in;
-  double Fsum, R_90, rad_90, flux_90;
-  Vector *rvec, *fvec, *Rvec, *Fvec;
+  double Fsum, Asum, Area, R_90, rad_90, flux_90;
+  Vector *rvec, *avec, *fvec, *Rvec, *Fvec, *Svec, *Avec;
 
-  if (argc != 2) {
-    gprint (GP_ERR, "USAGE: petrosian (radius) (flux) (P_ratio) (P_flux)\n");
+  if (argc != 8) {
+    gprint (GP_ERR, "USAGE: petrosian (radius) (area) (flux) (P_ratio) (P_flux) (P_sb) (P_area)\n");
     return (FALSE);
   }
@@ -15 +15 @@
   /* select input vectors */
   if ((rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((fvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((Rvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((Fvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((avec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((fvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Rvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fvec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Svec = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Avec = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
   // difficult work goes into cleaning the input galaxy profile:
@@ -34 +37 @@
   ResetVector (Fvec, OPIHI_FLT, fvec[0].Nelements);
   ResetVector (Rvec, OPIHI_FLT, fvec[0].Nelements);
+  ResetVector (Svec, OPIHI_FLT, fvec[0].Nelements);
+  ResetVector (Avec, OPIHI_FLT, fvec[0].Nelements);
 
   above = TRUE;
   Fsum = 0.0;
+  Asum = 0.0;
+  Area = avec[0].elements.Flt[0];
   for (i = 0; i < fvec[0].Nelements; i++) {
-    Fsum += fvec[0].elements.Flt[i] * M_PI * SQ(rvec[0].elements.Flt[i]);
-    Rvec[0].elements.Flt[i] = Fsum / fvec[0].elements.Flt[i];
+    // for nan bins, we keep the area for use with the next valid bin
+    if (isnan(fvec[0].elements.Flt[i])) {
+      Area += avec[0].elements.Flt[i];
+      continue;
+    } 
+    Fsum += fvec[0].elements.Flt[i] * Area;
+    Asum += Area;
+    if (i+1 < fvec[0].Nelements) {
+      Area = avec[0].elements.Flt[i+1];
+    }
+
+    Rvec[0].elements.Flt[i] = Asum * fvec[0].elements.Flt[i] / Fsum;
     Fvec[0].elements.Flt[i] = Fsum;
+    Svec[0].elements.Flt[i] = Fsum / Asum;
+    Avec[0].elements.Flt[i] = Asum;
 
     // anytime we transition below the petrosian ratio R_90, calculate the radius and flux