load new ramses with families

src/cpart.cc

@@ -29,6 +29,7 @@ CPart::CPart(const std::string _indir, const bool _v)
   verbose=_v;
   indir = _indir;
   infile="";
+  exist_family=false;
   
   // keep filename untill last /
   int found=indir.find_last_of("/");
@@ -50,6 +51,15 @@ CPart::CPart(const std::string _indir, const bool _v)
     infile = indir + "/part_" + s_run_index + ".out00001";
     if (verbose)
       std::cerr << "infile =[" << infile <<"]\n";
+    // check if new ramses format with family
+    std::ifstream fi;
+    fi.open(std::string(indir+"/part_file_descriptor.txt").c_str());
+    if (fi.is_open()) {
+      exist_family=true;
+      fi.close();
+    } else {
+      exist_family=false;
+    }
   }
 }
 
@@ -125,6 +135,7 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
     part.readDataBlock((char *) tmp[6]);
     double * agetmp, * metaltmp;
     int * id=NULL;
+    char * family=NULL;
 
     if (req_bits&ID_BIT) {
       id = new int[npart];
@@ -136,6 +147,11 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
     if (nstar>0) { // || 1) { // there are stars
 
       part.skipBlock(); // skip level
+      if (exist_family) {
+        family = new char[npart];
+        part.readDataBlock((char *) family); // read family
+        part.skipBlock(); // skip tag
+      }
       agetmp = new double[npart];
       part.readDataBlock((char *) agetmp);
 
@@ -151,8 +167,23 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
 
       }
       for (int k=0; k<npart; k++) {
-        if ((agetmp[k]==0.&& (comp_bits&HALO_BIT))  ||  // it's DM
-            (agetmp[k]!=0.&& (comp_bits&STARS_BIT))) {  // its' stars
+        bool ok_stars=false, ok_dm=false;
+        if (exist_family) {
+          if (family[k]==2) { // stars
+            ok_stars=true;
+          }
+          if (family[k]==1) { // stars
+            ok_dm=true;
+          }
+        } else {
+          if (agetmp[k]!=0.) { // stars
+            ok_stars=true;
+          } else {            // dm
+            ok_dm=true;
+          }
+        }
+        if ((ok_dm && (comp_bits&HALO_BIT))  ||  // it's DM
+            (ok_stars&& (comp_bits&STARS_BIT))) {  // its' stars
           if ((tmp[0][k]>=xmin && tmp[0][k]<=xmax) &&
               (tmp[1][k]>=ymin && tmp[1][k]<=ymax) &&
               ((ndim<3)||(tmp[2][k]>=zmin && tmp[2][k]<=zmax))
@@ -184,16 +215,19 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
               particles->load_bits |= MASS_BIT;
               take=true;
             }
-            if (agetmp[k]!=0 && req_bits&AGE_BIT) { // stars only (age)
+            //if (agetmp[k]!=0 && req_bits&AGE_BIT) { // stars only (age)
+            if (ok_stars && req_bits&AGE_BIT) { // stars only (age)
               particles->age.push_back(agetmp[k]);
               particles->load_bits |= AGE_BIT;
             }
-            if (agetmp[k]!=0 && req_bits&METAL_BIT && found_metal) { // stars only (metallicity)
+            //if (agetmp[k]!=0 && req_bits&METAL_BIT && found_metal) { // stars only (metallicity)
+            if (ok_stars && req_bits&METAL_BIT && found_metal) { // stars only (metallicity)
 	      //std::cerr << "Stars => "<<metaltmp[k]<<"\n";
               particles->metal.push_back(metaltmp[k]);
               particles->load_bits |= METAL_BIT;
             }
-            if (agetmp[k]==0 && req_bits&METAL_BIT && found_metal) { // stars only (metallicity)
+            //if (agetmp[k]==0 && req_bits&METAL_BIT && found_metal) { // stars only (metallicity)
+            if (ok_dm && req_bits&METAL_BIT && found_metal) { // stars only (metallicity)
 	      //if ( tmp[6][k] == 0.0 ) 
 	      //std::cerr << "DM => "<<metaltmp[k]<<" mass=" << tmp[6][k] << "\n";
               //particles->metal.push_back(tmp[6][k]); // put mass for dark matter
@@ -205,11 +239,11 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
               particles->id.push_back(id[k]); // save real id for dm or stars
               particles->load_bits |= ID_BIT;
             }
-            if ( (take && agetmp[k]!=0) || (!req_bits && agetmp[k]!=0)) { // !req_bits for uns_info and siplay=f
+            if ( (take && ok_stars) || (!req_bits && ok_stars)) { // !req_bits for uns_info and siplay=f
                particles->indexes.push_back(4); // save star positions
                particles->nstars++;               
             }
-            if ((take && agetmp[k]==0) || (!req_bits && agetmp[k]==0)) {
+            if ((take && ok_dm) || (!req_bits && ok_dm)) {
                particles->indexes.push_back(1); // save DM positions
                particles->ndm++;
             }
@@ -223,6 +257,7 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
       }
       // garbage
       delete [] agetmp;
+      delete [] family;
       if (req_bits&ID_BIT) {
         delete [] id;
       }
@@ -232,7 +267,7 @@ template <class T> int CPart::loadData(uns::CParticles<T> * particles,
     }
     else {  // there are no stars
       if (comp_bits&HALO_BIT) { // DM sel
-        for (int k=0; k<npart; k++) {
+        for (int k=0; k<npart; k++) {        
           if ((tmp[0][k]>=xmin && tmp[0][k]<=xmax) &&
               (tmp[1][k]>=ymin && tmp[1][k]<=ymax) &&
               ((ndim<3)||(tmp[2][k]>=zmin && tmp[2][k]<=zmax))

src/cpart.h

@@ -60,7 +60,7 @@ private:
   int npart,nstar,ncpu,ndim,nbody,ndm;
   int ndm_box, nstar_box;
   std::string s_run_index;
-
+  bool exist_family;
   double xmin,xmax,ymin,ymax,zmin,zmax;
 
   CFortIO part;