diff --git a/mytools/generateMock.cpp b/mytools/generateMock.cpp index 7e86645..4426ea9 100644 --- a/mytools/generateMock.cpp +++ b/mytools/generateMock.cpp @@ -9,6 +9,7 @@ #include #include "generateMock_conf.h" #include "gslIntegrate.hpp" +#include 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; } diff --git a/mytools/generateMock.ggo b/mytools/generateMock.ggo index 27b5d76..67d8e8f 100644 --- a/mytools/generateMock.ggo +++ b/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 \ No newline at end of file diff --git a/mytools/generateTestMock.cpp b/mytools/generateTestMock.cpp new file mode 100644 index 0000000..9853c36 --- /dev/null +++ b/mytools/generateTestMock.cpp @@ -0,0 +1,47 @@ +#include +#include +#include + +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; +} diff --git a/voz.h b/voz.h index 3a59560..99a20bc 100644 --- a/voz.h +++ b/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 */ diff --git a/voz1b1.c b/voz1b1.c index 1d97121..ebc7f0b 100644 --- a/voz1b1.c +++ b/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; ixmax) xmax = r[i][0]; diff --git a/vozinit.c b/vozinit.c index f1609c3..8b3db6c 100644 --- a/vozinit.c +++ b/vozinit.c @@ -141,8 +141,8 @@ int main(int argc, char *argv[]) { for (b[0]=0;b[0]