Changeset 25033 for branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.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/galradius.c (modified) (6 diffs)

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);
   }