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

Store the expansion factor of particles in params.nc

Browse source
This commit is contained in:
Guilhem Lavaux 2011-03-03 10:06:28 -06:00
parent cfa81cf3f0
commit 4cade55a44
1 changed files with 17 additions and 4 deletions
Download patch file Download diff file Expand all files Collapse all files

21
mytools/generateMock.cpp
View file

@ -100,7 +100,7 @@ Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z
return buildFromVector(pairs); return buildFromVector(pairs);
} }
void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel) void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact)
{ {
int x0, x1, x2; int x0, x1, x2;
@ -125,6 +125,8 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel)
TotalExpansion e_computer; TotalExpansion e_computer;
double baseComovingDistance; double baseComovingDistance;
expfact = new double[data->NumPart];
cout << "Using base redshift z=" << z0 << endl; cout << "Using base redshift z=" << z0 << endl;
e_computer.Omega_M = data->Omega_M; e_computer.Omega_M = data->Omega_M;
@ -138,11 +140,13 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel)
float& y = data->Pos[x1][i]; float& y = data->Pos[x1][i];
float& z = data->Pos[x2][i]; float& z = data->Pos[x2][i];
float& v = data->Vel[2][i]; float& v = data->Vel[2][i];
float z_old = z;
double reduced_red = (z + baseComovingDistance)*100./LIGHT_SPEED; double reduced_red = (z + baseComovingDistance)*100./LIGHT_SPEED;
// Distorted redshift // Distorted redshift
z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.; z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
expfact[i] = z / z_old;
// Add peculiar velocity // Add peculiar velocity
if (pecvel) if (pecvel)
z += v/100; z += v/100;
@ -201,7 +205,7 @@ void generateOutput(SimuData *data, int axis,
f.endCheckpoint(); f.endCheckpoint();
} }
void makeBox(SimuData *simu, SimuData *&boxed, generateMock_info& args_info) void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info)
{ {
uint32_t goodParticles = 0; uint32_t goodParticles = 0;
double ranges[3][2] = { double ranges[3][2] = {
@ -243,6 +247,7 @@ void makeBox(SimuData *simu, SimuData *&boxed, generateMock_info& args_info)
boxed->NumPart = goodParticles; boxed->NumPart = goodParticles;
particle_id = new int[goodParticles]; particle_id = new int[goodParticles];
double *expansion_fac = new double[goodParticles];
uint32_t k = 0; uint32_t k = 0;
for (uint32_t i = 0; i < simu->NumPart; i++) for (uint32_t i = 0; i < simu->NumPart; i++)
@ -264,6 +269,7 @@ void makeBox(SimuData *simu, SimuData *&boxed, generateMock_info& args_info)
assert(boxed->Pos[j][k] < 1); assert(boxed->Pos[j][k] < 1);
} }
particle_id[k] = i; particle_id[k] = i;
expansion_fac[k] = efac[i];
k++; k++;
} }
} }
@ -279,8 +285,13 @@ void makeBox(SimuData *simu, SimuData *&boxed, generateMock_info& args_info)
NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart); NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim); NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim);
NcVar *v2 = f.add_var("expansion", ncDouble, NumPart_dim);
v->put(particle_id, boxed->NumPart); v->put(particle_id, boxed->NumPart);
v2->put(expansion_fac, boxed->NumPart);
delete[] particle_id;
delete[] expansion_fac;
} }
int main(int argc, char **argv) int main(int argc, char **argv)
@ -322,9 +333,11 @@ int main(int argc, char **argv)
cout << "Omega_M = " << simu->Omega_M << endl; cout << "Omega_M = " << simu->Omega_M << endl;
cout << "Omega_Lambda = " << simu->Omega_Lambda << endl; cout << "Omega_Lambda = " << simu->Omega_Lambda << endl;
metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, args_info.peculiarVelocities_flag); double *expfact;
makeBox(simu, simuOut, args_info); metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag, args_info.peculiarVelocities_flag, expfact);
makeBox(simu, expfact, simuOut, args_info);
delete simu; delete simu;
generateOutput(simuOut, args_info.axis_arg, args_info.output_arg); generateOutput(simuOut, args_info.axis_arg, args_info.output_arg);