Better mock generation

mytools/generateMock.cpp

@@ -9,6 +9,7 @@
 #include <CosmoTool/fortran.hpp>
 #include "generateMock_conf.h"
 #include "gslIntegrate.hpp"
+#include <netcdfcpp.h>
 
 using namespace std;
 using namespace CosmoTool;
@@ -198,11 +199,93 @@ void generateOutput(SimuData *data, int axis,
   f.endCheckpoint();
 }
 
+void makeBox(SimuData *simu, SimuData *&boxed, generateMock_info& args_info)
+{
+  uint32_t goodParticles = 0;
+  double ranges[3][2] = {
+    { args_info.rangeX_min_arg, args_info.rangeX_max_arg },
+    { args_info.rangeY_min_arg, args_info.rangeY_max_arg },
+    { args_info.rangeZ_min_arg, args_info.rangeZ_max_arg }
+  };
+  double mul[3];
+  int *particle_id;
+
+  boxed = new SimuData;
+  boxed->Hubble = simu->Hubble;
+  boxed->Omega_M = simu->Omega_M;
+  boxed->Omega_Lambda = simu->Omega_Lambda;
+  boxed->time = simu->time;
+  boxed->BoxSize = simu->BoxSize;
+
+  for (uint32_t i = 0; i < simu->NumPart; i++)
+    {
+      bool acceptance = true;
+
+      for (int j = 0; j < 3; j++)
+	acceptance = 
+	  acceptance &&
+	  (simu->Pos[j][i] > ranges[j][0]) && 
+	  (simu->Pos[j][i] < ranges[j][1]);
+      
+      if (acceptance)
+	goodParticles++;
+    }
+  
+  for (int j = 0; j < 3; j++)
+    {
+      boxed->Pos[j] = new float[goodParticles];
+      boxed->Vel[j] = 0;
+      mul[j] = 1.0/(ranges[j][1] - ranges[j][0]);
+    }
+
+  boxed->NumPart = goodParticles;
+
+  particle_id = new int[goodParticles];
+
+  uint32_t k = 0;
+  for (uint32_t i = 0; i < simu->NumPart; i++)
+    {
+      bool acceptance = true;
+
+      for (int j = 0; j < 3; j++)
+	acceptance = 
+	  acceptance &&
+	  (simu->Pos[j][i] > ranges[j][0]) && 
+	  (simu->Pos[j][i] < ranges[j][1]);
+      
+      if (acceptance)
+	{
+	  for (int j = 0; j < 3; j++)
+	    {
+	      boxed->Pos[j][k] = (simu->Pos[j][i]-ranges[j][0])*mul[j];
+	      assert(boxed->Pos[j][k] > 0);
+	      assert(boxed->Pos[j][k] < 1);
+	    }
+	  particle_id[k] = i;
+	  k++;
+	}
+    }
+
+  NcFile f(args_info.outputParameter_arg);
+
+  f.add_att("range_x_min", ranges[0][0]);
+  f.add_att("range_x_max", ranges[0][1]);
+  f.add_att("range_y_min", ranges[1][0]);
+  f.add_att("range_y_max", ranges[1][1]);
+  f.add_att("range_z_min", ranges[2][0]);
+  f.add_att("range_z_max", ranges[2][1]);
+
+  NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
+  NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim);
+
+  v->put(particle_id, boxed->NumPart);
+}
+
 int main(int argc, char **argv)
 {
   generateMock_info args_info;
   generateMock_conf_params args_params;
-  SimuData *simu;
+  SimuData *simu, *simuOut;
  
   generateMock_conf_init(&args_info);
   generateMock_conf_params_init(&args_params);
@@ -239,7 +322,12 @@ int main(int argc, char **argv)
 
   metricTransform(simu, args_info.axis_arg);
 
-  generateOutput(simu, args_info.axis_arg, args_info.output_arg);
+  makeBox(simu, simuOut, args_info);
+  delete simu;
+
+  generateOutput(simuOut, args_info.axis_arg, args_info.output_arg);
+
+  delete simuOut;
   
   return 0;
 }

mytools/generateMock.ggo

@@ -10,4 +10,11 @@ option "ramsesId"   -   "Ramses snapshot id" int required
 option "axis"	    -	"Redshift axis (X=0, Y=1, Z=2)" int optional default="2"
 
 option "output"	    -   "Output filename for particles" string required
+option "outputParameter" - "Output geometry parameter file for postprocessing" string required
 
+option "rangeX_min"   -   "Minimum range in X for making the box" double required
+option "rangeX_max"   -   "Maximum range in X for making the box" double required
+option "rangeY_min"   -   "Minimum range in Y for making the box" double required
+option "rangeY_max"   -   "Maximum range in Y for making the box" double required
+option "rangeZ_min"   -   "Minimum range in Z for making the box (after distortion)" double required
+option "rangeZ_max"   -   "Maximum range in Z for making the box (after distortion)" double required

mytools/generateTestMock.cpp

@@ -0,0 +1,47 @@
+#include <iostream>
+#include <cstdlib>
+#include <CosmoTool/fortran.hpp>
+
+using namespace CosmoTool;
+using namespace std;
+
+#define LX 1.0
+#define LY 1.0
+#define LZ 1.0
+
+#define NUMPART (16*16*16)
+
+int main(int argc, char **argv)
+{
+  UnformattedWrite f("particles.bin");
+  
+  f.beginCheckpoint();
+  f.writeInt32(NUMPART);
+  f.endCheckpoint();
+
+  cout << "Writing X components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < NUMPART; i++)
+    {
+      f.writeReal32(drand48()*LX);
+    }
+  f.endCheckpoint();
+
+  cout << "Writing Y components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < NUMPART; i++)
+    {
+      f.writeReal32(drand48()*LY);
+    }
+  f.endCheckpoint();
+  
+  cout << "Writing Z components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < NUMPART; i++)
+    {
+      f.writeReal32(drand48()*LZ);
+    }
+  f.endCheckpoint();
+
+  return 0;
+}

voz.h

@@ -1,4 +1,4 @@
-#define MAXVERVER 2000
+#define MAXVERVER 4000
 #define NGUARD 42 /*Actually, the number of SPACES between guard points
 		    in each dim */
 

voz1b1.c

@@ -4,6 +4,8 @@
 #define DL for (d=0;d<3;d++)
 #define BF 1e30
 
+#define MAX(a,b) ( ((a) < (b)) ? (a) : (b) )
+
 int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs);
 int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol);
 int posread(char *posfile, float ***p, float fact);
@@ -35,7 +37,7 @@ int main(int argc, char *argv[]) {
     printf("Wrong number of arguments.\n");
     printf("arg1: position file\n");
     printf("arg2: border size\n");
-    printf("arg3: box size\n");
+    printf("arg3: boxsize\n");
     printf("arg4: suffix\n");
     printf("arg5: number of divisions\n");
     printf("arg6-8: b[0-2]\n\n");
@@ -78,6 +80,7 @@ int main(int argc, char *argv[]) {
   /* Boxsize should be the range in r, yielding a range 0-1 */
   np = posread(posfile,&r,1./boxsize);
   printf("%d particles\n",np);fflush(stdout);
+
   xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
   for (i=0; i<np;i++) {
     if (r[i][0]<xmin) xmin = r[i][0]; if (r[i][0]>xmax) xmax = r[i][0];

vozinit.c

@@ -141,8 +141,8 @@ int main(int argc, char *argv[]) {
   for (b[0]=0;b[0]<numdiv; b[0]++) {
    for (b[1] = 0; b[1] < numdiv; b[1]++) {
     for (b[2] = 0; b[2] < numdiv; b[2]++) {
-      fprintf(scr,"./voz1b1 %s %f %f %s %d %d %d %d &\n",
+      fprintf(scr,"./voz1b1 %s %f %f,%f,%f %s %d %d %d %d &\n",
-	     posfile,border,boxsize,suffix,numdiv,b[0],b[1],b[2]);
+	      posfile,border,boxsize,boxsize,boxsize,suffix,numdiv,b[0],b[1],b[2]);
       p++;
       if ((p == NUMCPU)) { fprintf(scr, "wait\n"); p = 0; }
     }