bug in uns_rectify, please compare with snaprect

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@111 f264a43e-d52d-4b82-913b-c2bd5215a18a

lib/projects/nemodep/crectify.cc

@@ -58,7 +58,7 @@ CRectify::CRectify()
 bool CRectify::rectify(const int _nbody,const float _time,
                        float * _pos, float * _vel,float * _mass, float * _rho,
                        const bool _use_rho, const bool _rectify,
-                       std::string _cod_file,std::string _rect_file, const float _threshold)
+                       std::string _cod_file,std::string _rect_file, const float _radius,const float _dmin, const float _dmax)
 {
   bool ok = true;
   nbody= _nbody;
@@ -68,7 +68,9 @@ bool CRectify::rectify(const int _nbody,const float _time,
   rho  = _rho;
   time = _time;
 
-  threshold = _threshold;
+  radius = _radius;
+  dmin = _dmin;
+  dmax = _dmax;
   cod_file  = _cod_file;
   rect_file = _rect_file;
 
@@ -194,7 +196,8 @@ void CRectify::processRho()
   int ii=0;
   for (std::vector<CDataIndex>::iterator it=rho_di.begin(); it != rho_di.end(); it++) {
     float logrho = log((*it).data);
-    if (((logrho - minrho)*100./(maxrho-minrho))>=threshold) { // we keep particles
+    if (((logrho - minrho)*100./(maxrho-minrho))>=dmin &&
+        ((logrho - minrho)*100./(maxrho-minrho))<=dmax) { // we keep particles
       int sindex=(*it).index;
       vpos.push_back(pos[sindex*3+0]);
       vpos.push_back(pos[sindex*3+1]);
@@ -228,36 +231,51 @@ void CRectify::findMoment()
 {
   CLRM(w_qpole);
   float w_sum=0.0;
+  float distance=0.0;
 
   if ( use_rho) { // data are stored in Vectors
     for (unsigned int i=0; i<vmass.size(); i++) {
+      ABSV(distance,&vpos[i*3]);
+      if (distance <=radius) {
         w_sum += vrho[i] * vmass[i];
       }
+    }
 
     float tmpv[3], pos_b[3];
     matrix tmpm;
     for (unsigned int i=0; i<vmass.size(); i++) {
-      float w_b = vrho[i] * vmass[i];
       vectutils::subv(pos_b, &vpos[i*3], codf); // shift to cod
+      ABSV(distance,pos_b);
+      if (distance <=radius) {
+        float w_b = vrho[i] * vmass[i];
+        //vectutils::subv(pos_b, &vpos[i*3], codf); // shift to cod
         MULVS(tmpv, pos_b, w_b);
         OUTVP(tmpm, tmpv, pos_b);
         ADDM(w_qpole, w_qpole, tmpm);
       }
+    }
   } else {  // data are stored in arrays
     for ( int i=0; i<nbody; i++) {
+      ABSV(distance,&pos[i*3]);
+      if (distance <=radius) {
         w_sum += mass[i];
       }
+    }
 
     float tmpv[3], pos_b[3];
     matrix tmpm;
     for ( int i=0; i<nbody; i++) {
-      float w_b = mass[i];
       vectutils::subv(pos_b, &pos[i*3], codf); // shift to cod
+      ABSV(distance,pos_b);
+      if (distance <=radius) {
+        float w_b = mass[i];
+        //vectutils::subv(pos_b, &pos[i*3], codf); // shift to cod
         MULVS(tmpv, pos_b, w_b);
         OUTVP(tmpm, tmpv, pos_b);
         ADDM(w_qpole, w_qpole, tmpm);
       }
     }
+  }
   DIVMS(w_qpole, w_qpole, w_sum);
   std::cerr << " find moment, done\n";
 }

lib/projects/nemodep/crectify.h

@@ -74,7 +74,7 @@ public:
   bool rectify(const int nbody,const float time,
                float * pos, float * vel,float * mass, float * weight,
                const bool _use_rho, const bool rectify=false,
-               std::string cod_file="",std::string rect_file="", const float threshold=0.0);
+               std::string cod_file="",std::string rect_file="", const float radius=50.0, const float dmin=0.0, const float dmax=100.0);
 
   void process();
   //------------------------------------------------------------------------------
@@ -96,7 +96,7 @@ public:
 private:
   jclut::CDensity * density;
   float time, * pos, * vel, * mass, * rho;
-  float threshold;
+  float radius,dmin, dmax;
   float oldframe[3][3];
   std::string cod_file, rect_file;
   bool rect, use_rho;

py/mains/g3cpu_parse.py

@@ -131,7 +131,7 @@ def parseCPU(file,dt,tf,tl):
 
     time=np.array([], dtype=np.float32)
     cpu =np.array([], dtype=np.float32)
-    
+    print "cpu.size = ", cpu.size
     first=True
     tlast=0.0
     t=0.0
@@ -160,10 +160,17 @@ def parseCPU(file,dt,tf,tl):
                 cputotlast=cputot
                 #total=a.readline()
                 #print "diff =",t, diff_t,diff_cpu
-                time=np.append(time,t)       # save time in numpy array
-                cpu=np.append(cpu,diff_cpu)  # save cputime in numpy array
+
+                idx=time.size
+                time.resize(idx+1)
+                time[idx]=t
+                cpu.resize(idx+1)
+                cpu[idx]=diff_cpu
+                #time=np.append(time,t)       # save time in numpy array
+                #cpu=np.append(cpu,diff_cpu)  # save cputime in numpy array
 
     a.close()
+    print "cpu.size = ", cpu.size
     return time,cpu,ncores
 
 def printHelp(prog,file,dt,tf,tl,out):

src/nemodep/uns_rectify.cc

@@ -35,7 +35,9 @@ const char * defv[] = {
   "rectf=\n            rectified info file",
   "rho=t\n             use rho as weight factor, otherwise mass",
   "codfile=\n          use cod from file to recenter",
-  "threshold=0.0\n    percentage of density kept >= threshold",
+  "radius=50.0\n       cutting radius",
+  "dmin=0.0\n          keep particles if %log(density) >= dmin",
+  "dmax=100.0\n        keep particles if %log(density) <= dmax",
   "times=all\n		     selected time",
   "first=f\n           add a trailing numbering to the first output file",
   "verbose=f\n         verbose on/off",
@@ -63,7 +65,9 @@ int main(int argc, char ** argv )
   bool        use_rho  (getbparam((char *) "rho"     ));
   std::string select_c (getparam ((char *) "select"  ));
   std::string select_t (getparam ((char *) "times"    ));
-  float       threshold(getdparam ((char *) "threshold"));
+  float       radius(getdparam ((char *) "radius"));
+  float       dmin(getdparam ((char *) "dmin"));
+  float       dmax(getdparam ((char *) "dmax"));
   bool        first    (getbparam((char *) "first"     ));
   bool        verbose = getbparam((char *) "verbose"  );
 
@@ -120,7 +124,7 @@ int main(int argc, char ** argv )
 
       std::cerr << "Snapshot ="<< uns->snapshot->getFileName() << "\n";
 
-      crectify->rectify(nbody,time,pos,vel,mass,rho,use_rho,true,codf,rectf, threshold);
+      crectify->rectify(nbody,time,pos,vel,mass,rho,use_rho,true,codf,rectf, radius,dmin,dmax);
 
       if (outname.length()!=0) {
         stringstream number;