guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Imported code for flexible gadget units. Imported load ramses code. Always invoke metricTransform. Add the option to import box parameters from an already existing mock catalog

Browse source
This commit is contained in:
Guilhem Lavaux 2012-11-08 10:00:08 -05:00
parent 058b12bb3f
commit 9b8ce57a97
2 changed files with 149 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

166
c_tools/mock/generateMock.cpp
View file

@ -18,9 +18,67 @@ using namespace CosmoTool;
#define LIGHT_SPEED 299792.458
static double gadgetUnit=1e-3;
SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRedshift)
{
SimuData *d, *outd;
d = loadRamsesSimu(basename, baseid, -1, true, 0);
outd = new SimuData;
outd->NumPart = d->TotalNumPart;
outd->BoxSize = d->BoxSize;
outd->TotalNumPart = outd->NumPart;
outd->Hubble = d->Hubble;
outd->Omega_Lambda = d->Omega_Lambda;
outd->Omega_M = d->Omega_M;
outd->time = d->time;
for (int k = 0; k < 3; k++)
outd->Pos[k] = new float[outd->NumPart];
outd->Vel[2] = new float[outd->NumPart];
delete d;
int curCpu = 0;
cout << "loading cpu 0 " << endl;
while (d = loadRamsesSimu(basename, baseid, curCpu, true, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID))
{
for (int k = 0; k < 3; k++)
for (int i = 0; i < d->NumPart; i++)
{
assert(d->Id[i] >= 1);
assert(d->Id[i] <= outd->TotalNumPart);
outd->Pos[k][d->Id[i]-1] = d->Pos[k][i];
outd->Vel[2][d->Id[i]-1] = d->Vel[velAxis][i];
}
if (goRedshift)
for (int i = 0; i < d->NumPart; i++)
outd->Pos[velAxis][d->Id[i]-1] += d->Vel[velAxis][i]/100.;
delete d;
curCpu++;
cout << "loading cpu " << curCpu << endl;
}
return outd;
}
SimuData *myLoadGadget(const char *fname, int id, int flags)
{
return loadGadgetMulti(fname, id, flags);
SimuData *sim = loadGadgetMulti(fname, id, flags);
sim->BoxSize *= gadgetUnit*1000;
for (int j = 0; j < 3; j++)
{
if (sim->Pos[j] != 0) {
for (long i = 0; i < sim->NumPart; i++)
sim->Pos[j][i] *= gadgetUnit*1000;
}
}
return sim;
}
SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, SimuData *(*loadFunction)(const char *fname, int id, int flags))
@ -111,7 +169,7 @@ SimuData *doLoadMultidark(const char *multidarkname)
fscanf(fp, "%f\n", &outd->Omega_M);
fscanf(fp, "%f\n", &outd->Hubble);
fscanf(fp, "%f\n", &outd->time);
fscanf(fp, "%d\n", &outd->NumPart);
fscanf(fp, "%ld\n", &outd->NumPart);
outd->time = 1./(1.+outd->time); // convert to scale factor
outd->TotalNumPart = outd->NumPart;
@ -181,7 +239,7 @@ Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z
return buildFromVector(pairs);
}
void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact)
void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact, bool cosmo_flag)
{
int x0, x1, x2;
@ -200,7 +258,11 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
}
double z0 = 1/data->time - 1;
Interpolate z_vs_D = make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, 0., z0+8*data->BoxSize*100/LIGHT_SPEED); // Redshift 2*z0 should be sufficient ?
Interpolate z_vs_D =
make_cosmological_redshift(data->Omega_M, data->Omega_Lambda,
0., z0+8*data->BoxSize*100/LIGHT_SPEED);
// Redshift 2*z0 should be sufficient ? This is fragile.
//A proper solver is needed here.
double z_base = reshift ? z0 : 0;
TotalExpansion e_computer;
double baseComovingDistance;
@ -229,13 +291,15 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
// Distorted redshift
if (reduced_red == 0)
z = 0;
else {
z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
}
else if (cosmo_flag)
z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
else
z = reduced_red*LIGHT_SPEED/100.0;
expfact[i] = z / z_old;
// Add peculiar velocity
if (pecvel)
z += v/100;
// Add peculiar velocity
if (pecvel)
z += v/100;
}
catch(const InvalidRangeException& e) {
cout << "Trying to interpolate out of the tabulated range." << endl;
@ -404,6 +468,63 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info&
delete[] expansion_fac;
}
void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info)
{
NcFile f(args_info.inputParameter_arg);
NcVar *v;
int *particle_id;
double *expansion_fac;
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;
NcVar *v_id = f.get_var("particle_ids");
long *edges1;
double ranges[3][2];
double mul[3];
edges1 = v_id->edges();
assert(v_id->num_dims()==1);
boxed->NumPart = edges1[0];
delete[] edges1;
particle_id = new int[boxed->NumPart];
v_id->get(particle_id, boxed->NumPart);
ranges[0][0] = f.get_att("range_x_min")->as_double(0);
ranges[0][1] = f.get_att("range_x_max")->as_double(0);
ranges[1][0] = f.get_att("range_y_min")->as_double(0);
ranges[1][1] = f.get_att("range_y_max")->as_double(0);
ranges[2][0] = f.get_att("range_z_min")->as_double(0);
ranges[2][1] = f.get_att("range_z_max")->as_double(0);
for (int j = 0; j < 3; j++)
{
boxed->Pos[j] = new float[boxed->NumPart];
boxed->Vel[j] = 0;
mul[j] = 1.0/(ranges[j][1] - ranges[j][0]);
}
uint32_t k = 0;
for (uint32_t i = 0; i < boxed->NumPart; i++)
{
int id = particle_id[i];
for (int j = 0; j < 3; j++)
{
boxed->Pos[j][i] = (simu->Pos[j][id]-ranges[j][0])*mul[j];
}
}
delete[] particle_id;
}
int main(int argc, char **argv)
{
generateMock_info args_info;
@ -434,13 +555,15 @@ int main(int argc, char **argv)
generateMock_conf_print_version();
gadgetUnit=args_info.gadgetUnit_arg;
if (args_info.ramsesBase_given || args_info.ramsesId_given)
{
if (args_info.ramsesBase_given && args_info.ramsesId_given) {
//simu = doLoadRamses(args_info.ramsesBase_arg,
// args_info.ramsesId_arg,
// args_info.axis_arg, false);
}
simu = doLoadRamses(args_info.ramsesBase_arg,
args_info.ramsesId_arg,
args_info.axis_arg, false);
}
else
{
cerr << "Both ramsesBase and ramsesId are required to be able to load snapshots" << endl;
@ -490,14 +613,15 @@ int main(int argc, char **argv)
double *expfact;
if (args_info.cosmo_flag){
metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, args_info.peculiarVelocities_flag, expfact);
} else {
expfact = new double[simu->NumPart];
for (int j = 0; j < simu->NumPart; j++) expfact[j] = 1.0;
}
metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag,
args_info.peculiarVelocities_flag, expfact,
args_info.cosmo_flag);
if (args_info.inputParameter_given)
makeBoxFromParameter(simu, expfact, simuOut, args_info);
else
makeBox(simu, expfact, simuOut, args_info);
makeBox(simu, expfact, simuOut, args_info);
delete simu;
generateOutput(simuOut, args_info.axis_arg, args_info.output_arg);

4
c_tools/mock/generateMock.ggo
View file

@ -29,3 +29,7 @@ option "peculiarVelocities" - "Added peculiar velocities distortion" flag off
option "cosmo" - "Apply cosmological redshift" flag on
option "subsample" - "Subsample the input simulation by the specified amount" double optional
option "inputParameter" - "Input geometry (optional, warning!)" string optional
option "gadgetUnit" - "Unit of length in gadget file in Mpc/h" double optional default="0.001"