From 793e285c6fd8401963480ed103b4ee89ff9707f2 Mon Sep 17 00:00:00 2001 From: Your Name Date: Tue, 22 Feb 2011 14:35:03 -0500 Subject: [PATCH] Generate distorted particle distribution --- mytools/generateMock.cpp | 155 ++++++++++++++++++++++++++++++++++++++- mytools/generateMock.ggo | 5 +- mytools/gslIntegrate.hpp | 33 +++++++++ 3 files changed, 189 insertions(+), 4 deletions(-) create mode 100644 mytools/gslIntegrate.hpp diff --git a/mytools/generateMock.cpp b/mytools/generateMock.cpp index d7331f3..9c62796 100644 --- a/mytools/generateMock.cpp +++ b/mytools/generateMock.cpp @@ -1,13 +1,20 @@ +#include #include #include #include +#include #include #include +#include +#include #include "generateMock_conf.h" +#include "gslIntegrate.hpp" using namespace std; using namespace CosmoTool; +#define LIGHT_SPEED 299792.458 + SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRedshift) { SimuData *d, *outd; @@ -19,6 +26,8 @@ SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRed outd->BoxSize = d->BoxSize; outd->TotalNumPart = outd->NumPart; outd->Hubble = d->Hubble; + outd->Omega_Lambda = d->Omega_Lambda; + for (int k = 0; k < 3; k++) outd->Pos[k] = new float[outd->NumPart]; outd->Vel[2] = new float[outd->NumPart]; @@ -49,6 +58,144 @@ SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRed return outd; } +static double cubic(double a) +{ + return a*a*a; +} + +struct TotalExpansion +{ + double Omega_M, Omega_L; + + double operator()(double z) + { + return 1/sqrt(Omega_M*cubic(1+z) + Omega_L); + } +}; + +Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z1) +{ + TotalExpansion e_computer; + double D_tilde, Q, Qprime; + InterpolatePairs pairs; + + e_computer.Omega_M = OM; + e_computer.Omega_L = OL; + + pairs.resize(100); + ofstream f("comoving_distance.txt"); + + for (int i = 0; i < 100; i++) + { + double z = z0 + (z1-z0)/100*i; + + pairs[i].second = z; + pairs[i].first = gslIntegrate(e_computer, 0, z, 1e-3); + f << z << " " << pairs[i].first << endl; + } + + return buildFromVector(pairs); +} + +void metricTransform(SimuData *data, int axis) +{ + int x0, x1, x2; + + switch (axis) { + case 0: + x0 = 1; x1 = 2; x2 = 0; + break; + case 1: + x0 = 0; x1 = 2; x2 = 1; + break; + case 2: + x0 = 0; x1 = 1; x2 = 2; + break; + default: + abort(); + } + + Interpolate z_vs_D = make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, 0., 2.0); // Redshift 2 should be sufficient ? + + double z0 = 1/data->time - 1; + TotalExpansion e_computer; + double baseComovingDistance; + + cout << "Using base redshift z=" << z0 << endl; + + e_computer.Omega_M = data->Omega_M; + e_computer.Omega_L = data->Omega_Lambda; + baseComovingDistance = LIGHT_SPEED/100.* gslIntegrate(e_computer, 0, z0, 1e-3); + cout << "Comoving distance = " << baseComovingDistance << " Mpc/h" << endl; + + for (uint32_t i = 0; i < data->NumPart; i++) + { + float& x = data->Pos[x0][i]; + float& y = data->Pos[x1][i]; + float& z = data->Pos[x2][i]; + float& v = data->Vel[2][i]; + + double reduced_red = (z + baseComovingDistance)*100./LIGHT_SPEED; + + // Distorted redshift + z = z_vs_D.compute(reduced_red)*LIGHT_SPEED/100.; + // Add peculiar velocity + z += v; + } +} + +void generateOutput(SimuData *data, int axis, + const std::string& fname) +{ + UnformattedWrite f(fname); + + cout << "Generating output particles to " << fname << endl; + + int x0, x1, x2; + + switch (axis) { + case 0: + x0 = 1; x1 = 2; x2 = 0; + break; + case 1: + x0 = 0; x1 = 2; x2 = 1; + break; + case 2: + x0 = 0; x1 = 1; x2 = 2; + break; + default: + abort(); + } + + f.beginCheckpoint(); + f.writeInt32(data->NumPart); + f.endCheckpoint(); + + cout << "Writing X components..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->Pos[x0][i]); + } + f.endCheckpoint(); + + cout << "Writing Y components..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->Pos[x1][i]); + } + f.endCheckpoint(); + + cout << "Writing Z components..." << endl; + f.beginCheckpoint(); + for (uint32_t i = 0; i < data->NumPart; i++) + { + f.writeReal32(data->Pos[x2][i]); + } + f.endCheckpoint(); +} + int main(int argc, char **argv) { generateMock_info args_info; @@ -81,10 +228,16 @@ int main(int argc, char **argv) simu = doLoadRamses(args_info.ramsesBase_arg, args_info.ramsesId_arg, - args_info.axis_arg, true); + args_info.axis_arg, false); cout << "Hubble = " << simu->Hubble << endl; cout << "Boxsize = " << simu->BoxSize << endl; + cout << "Omega_M = " << simu->Omega_M << endl; + cout << "Omega_Lambda = " << simu->Omega_Lambda << endl; + + metricTransform(simu, args_info.axis_arg); + + generateOutput(simu, args_info.axis_arg, args_info.output_arg); return 0; } diff --git a/mytools/generateMock.ggo b/mytools/generateMock.ggo index 8b1ecc1..27b5d76 100644 --- a/mytools/generateMock.ggo +++ b/mytools/generateMock.ggo @@ -8,7 +8,6 @@ option "ramsesBase" - "Base directory for ramses" string required option "ramsesId" - "Ramses snapshot id" int required option "axis" - "Redshift axis (X=0, Y=1, Z=2)" int optional default="2" -option "bufferSize" - "Size of the buffer regions in unit of X-Y size" double optional default="0.1" -option "divisions" - "number of divisions" int optional default="2" -option "suffix" - "suffix for the run" string optional default="run" \ No newline at end of file +option "output" - "Output filename for particles" string required + diff --git a/mytools/gslIntegrate.hpp b/mytools/gslIntegrate.hpp new file mode 100644 index 0000000..f3f2dd3 --- /dev/null +++ b/mytools/gslIntegrate.hpp @@ -0,0 +1,33 @@ +#ifndef __MYGSL_INTEGRATE_HPP +#define __MYGSL_INTEGRATE_HPP + +#include + +template +double gslSpecialFunction(double x, void *param) +{ + FunT *f = (FunT *)param; + + return (*f)(x); +} + +template +double gslIntegrate(FunT& v, double a, double b, double prec, int NPTS = 1024) +{ + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NPTS); + gsl_function f; + double result; + double abserr; + + f.function = &gslSpecialFunction; + f.params = &v; + + gsl_integration_qag(&f, a, b, 0, prec, NPTS, GSL_INTEG_GAUSS61, + w, &result, &abserr); + + gsl_integration_workspace_free(w); + + return result; +} + +#endif