Changeset 25033

Timestamp:

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

Author:

eugene

Message:

more tests of robust petrosian analysis

Location:

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

Files:

• Makefile (modified) (1 diff)
• galradbins.c (added)
• galradius.c (modified) (6 diffs)
• init.c (modified) (2 diffs)
• petrosian.c (modified) (3 diffs)

branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile

@@ -38 +38 @@
 $(SRC)/mksersic.$(ARCH).o          \
 $(SRC)/galradius.$(ARCH).o         \
+$(SRC)/galradbins.$(ARCH).o        \
 $(SRC)/galsectors.$(ARCH).o        \
 $(SRC)/galprofiles.$(ARCH).o       \

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

@@ -11 +11 @@
   opihi_flt *flux, *radius, *values;
   Buffer *buf;
-  Vector *fvec, *rvec, *tvec;
+  Vector *fvec, *rvec, *tvec, *Fvec, *Rvec;
   char name[128];
 
@@ -22 +22 @@
   if ((rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((fvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if ((Rvec = SelectVector ("rg", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fvec = SelectVector ("fg", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
   Fmin = atof(argv[3]);
   Fmax = atof(argv[4]);
@@ -50 +54 @@
     Rbin = 1;
   } else {
-    Rbin = MAX(1, Rsum / Rnpt);
+    Rbin = MAX(1, 0.5*(Rsum / Rnpt));
   }
 
@@ -56 +60 @@
   if (Rbin <= 2) {
     Nout   = fvec[0].Nelements;
+    ResetVector (Rvec, OPIHI_FLT, Nout);
+    ResetVector (Fvec, OPIHI_FLT, Nout);
+    memcpy (Fvec[0].elements.Flt, fvec[0].elements.Flt, Fvec[0].Nelements*sizeof(opihi_flt));
+    memcpy (Rvec[0].elements.Flt, rvec[0].elements.Flt, Rvec[0].Nelements*sizeof(opihi_flt));
     flux   = fvec[0].elements.Flt;
     radius = rvec[0].elements.Flt;
-    // XXX handle the necessary free below
   } else {
     // rebin the vectors by Rbin values:
     Nout = fvec[0].Nelements / Rbin + 1;
-    ALLOCATE (flux,   opihi_flt, Nout);
-    ALLOCATE (radius, opihi_flt, Nout);
     ALLOCATE (values, opihi_flt, fvec[0].Nelements);
   
+    ResetVector (Rvec, OPIHI_FLT, Nout);
+    ResetVector (Fvec, OPIHI_FLT, Nout);
+    flux   = Fvec[0].elements.Flt;
+    radius = Rvec[0].elements.Flt;
+
     for (i = 0; i < Nout; i++) {
       radius[i] = (i + 0.5)*Rbin;
@@ -73 +83 @@
       N = 0;
       for (j = 0; j < fvec[0].Nelements; j++) {
-        if (rvec[0].elements.Flt[i] < Rmin) continue;
-        if (rvec[0].elements.Flt[i] > Rmax) continue;
-        values[N] = fvec[0].elements.Flt[i];
+        if (rvec[0].elements.Flt[j] < Rmin) continue;
+        if (rvec[0].elements.Flt[j] > Rmax) continue;
+        values[N] = fvec[0].elements.Flt[j];
         N++;
       }
@@ -87 +97 @@
       }
     }
+    free (values);
   }
 

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

@@ -25 +25 @@
 int galprofiles             PROTO((int, char **));
 int galradius               PROTO((int, char **));
+int galradbins              PROTO((int, char **));
 int elliprofile             PROTO((int, char **));
 int petrosian               PROTO((int, char **));
@@ -71 +72 @@
   {1, "galprofiles", galprofiles,  "generate radial vectors with interpolation along paths"},
   {1, "galradius",   galradius,    "generate radial vectors with interpolation along paths"},
+  {1, "galradbins",  galradbins,   "generate radial vectors with interpolation along paths"},
   {1, "elliprofile", elliprofile,  "generate radial vectors with interpolation along paths"},
   {1, "petrosian",   petrosian,    "petrosian parameters given radial bins"},

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