Merge branch 'master' of bitbucket.org:cosmicvoids/void_identification

c_tools/mock/generateMock.cpp

@@ -125,6 +125,7 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
       float z_old = z;
 
       double reduced_red = (z + baseComovingDistance)*100./LIGHT_SPEED;
+      double reduced_base = (reshift ? (baseComovingDistance*100./LIGHT_SPEED) : 0);
       try
         {
 
@@ -134,7 +135,7 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
           else if (cosmo_flag)
 	    z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
           else
-            z = reduced_red*LIGHT_SPEED/100.0;
+            z = (reduced_red-reduced_base)*LIGHT_SPEED/100.0;
 
           if (expfact)
             expfact[i] = z / z_old; 

c_tools/mock/loaders/sdf_loader.cpp

@@ -69,7 +69,6 @@ public:
     return num_splitting;
   }
 
-
   int64_t getStart(int id)
   {
     return sdf_header->TotalNumPart * int64_t(id) / num_splitting;
@@ -87,7 +86,7 @@ public:
 #define INFINITY std::numeric_limits<float>::max()
     float shift = 0.5*d->BoxSize;
     rescale_position /= d->time;
-    float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
+//    float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
 
     if (d->Pos[0] != 0)
       {
@@ -96,10 +95,10 @@ public:
             for (int64_t i = 0; i < d->NumPart; i++)
               {
                 d->Pos[k][i] = (d->Pos[k][i]*rescale_position + shift);
-                min_pos[k] = std::min(min_pos[k], d->Pos[k][i]);
+//                min_pos[k] = std::min(min_pos[k], d->Pos[k][i]);
-                max_pos[k] = std::max(max_pos[k], d->Pos[k][i]);
+//                max_pos[k] = std::max(max_pos[k], d->Pos[k][i]);
               }
-            cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
+//            cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
           }
       }
     if (d->Vel[0] != 0)
@@ -115,6 +114,20 @@ public:
   }
 
 
+  void enforceBoxSize(SimuData *d)
+  {
+    for (int k = 0; k < 3; k++)
+      {
+        for (int64_t i = 0; i < d->NumPart; i++)
+         {
+           while (d->Pos[k][i]<0)
+             d->Pos[k][i] += d->BoxSize;
+           while (d->Pos[k][i]>=d->BoxSize)
+             d->Pos[k][i] -= d->BoxSize;
+         }
+       }
+  }
+
   SimuData *loadFile(int id) {
     SimuData *d;
     int64_t starts[7];
@@ -177,6 +190,7 @@ public:
     if (load_flags & (NEED_POSITION | NEED_VELOCITY))
       rescaleParticles(d, d->Hubble*1e-5, one_kpc/one_Gyr);
 
+    enforceBoxSize(d);
     applyTransformations(d);
     basicPreprocessing(d, preproc);
 

external/cosmotool/src/loadGadget.cpp

@@ -268,10 +268,10 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
 
 
 
-void CosmoTool::writeGadget(const char *fname, SimuData *data, int GadgetFormat)
+void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int GadgetFormat)
 {
   UnformattedWrite *f;
-  int npart[6];
+  int npart[6], npartTotal[6];
   float mass[6];
 
   if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0)
@@ -283,12 +283,12 @@ void CosmoTool::writeGadget(const char *fname, SimuData *data, int GadgetFormat)
 
   for (int i = 0; i < 6; i++)
     {
-      npart[i] = 0;
+      npart[i] = npartTotal[i] = 0;
       mass[i] = 0;
     }
   mass[1] = 1.0;
-
   npart[1] = data->NumPart;
+  npartTotal[1] = data->TotalNumPart;
   
   f->beginCheckpoint();
   for (int i = 0; i < 6; i++)
@@ -302,7 +302,7 @@ void CosmoTool::writeGadget(const char *fname, SimuData *data, int GadgetFormat)
   f->writeInt32(0);
 
   for (int i = 0; i < 6; i++)
-    f->writeInt32(npart[i]);
+    f->writeInt32(npartTotal[i]);
   f->writeInt32(0);
   f->writeInt32(1);
   f->writeReal64(data->BoxSize);

external/cosmotool/src/loadGadget.hpp

@@ -36,6 +36,8 @@ knowledge of the CeCILL license and that you accept its terms.
 #ifndef __COSMO_LOAD_GADGET_HPP
 #define __COSMO_LOAD_GADGET_HPP
 
+#include <string>
+
 #include "load_data.hpp"
 #include "loadSimu.hpp"
 
@@ -45,7 +47,7 @@ namespace CosmoTool {
 			    int GadgetFormat = 1, SimuFilter filter = 0);
  
   // Only single snapshot supported
-  void writeGadget(const char *fname, SimuData *data, int GadgetFormat = 1);
+  void writeGadget(const std::string& fname, SimuData *data, int GadgetFormat = 1);
  
 };
 

external/cosmotool/src/mykdtree.hpp

@@ -136,7 +136,13 @@ namespace CosmoTool {
     void setPeriodic(bool on, CoordType replicate)
      {
        periodic = on;
-       this->replicate = replicate;
+       std::fill(this->replicate, this->replicate+N, replicate);
+     }
+
+    void setPeriodic(bool on, const coords& replicate)
+     {
+       periodic = on;
+       std::copy(replicate, replicate+N, this->replicate);
      }
 
     uint32_t getIntersection(const coords& x, CoordType r, 
@@ -193,7 +199,7 @@ namespace CosmoTool {
     Cell *base_cell;
 
     bool periodic;
-    CoordType replicate;
+    coords replicate;
 
     Node *buildTree(Cell **cell0,
 		    uint32_t NumCells,

external/cosmotool/src/replicateGenerator.hpp

@@ -46,18 +46,29 @@ namespace CosmoTool
   {
   public:
     typedef Coord Coords[N];
-    Coord replicate;
+    Coords replicate;
 
     ReplicateGenerator(const Coords& x, Coord shift)
     {
       face = 0; 
-      replicate = shift;
+      std::fill(replicate, replicate+N, shift);
       numFaces = spower<N,long>(3);
       std::copy(x, x+N, x_base);
       if (!next())
 	abort();
     }
 
+    ReplicateGenerator(const Coords& x, Coords& shift)
+    {
+      face = 0;
+      std::copy(shift, shift+N, replicate);
+      numFaces = spower<N,long>(3);
+      std::copy(x, x+N, x_base);
+      if (!next())
+        abort();
+    }
+
+
     bool next()
     {
       if (face == numFaces)
@@ -72,7 +83,7 @@ namespace CosmoTool
           int c_face;
           c_face = q_face % 3;
           q_face /= 3;
-          x_shifted[i] = x_base[i] + (c_face-1)*replicate;
+          x_shifted[i] = x_base[i] + (c_face-1)*replicate[i];
           no_move = no_move && (c_face == 1);
         }
       if (no_move)