diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 17e01c4..765e1b3 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -55,3 +55,7 @@ if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES}) set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) + +add_executable(test_cosmopower test_cosmopower.cpp) +target_link_libraries(test_cosmopower ${tolink}) + diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp index e82c25e..19aec97 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -191,15 +191,15 @@ static double gslPowSpecNorm(double k, void *params) return c->integrandNormalize(k); } -double CosmoPower::integrandNormalize(double k) +double CosmoPower::integrandNormalize(double x) { + double k = (1-x)/x; double f = tophatFilter(k*8.0/h); - return power(k)*k*k*f*f; + return power(k)*k*k*f*f/(x*x); } void CosmoPower::normalize() { - int Nsteps = 30000; double normVal = 0; double abserr; gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); @@ -210,8 +210,15 @@ void CosmoPower::normalize() normPower = 1; - gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); + ofstream ff("PP_k.txt"); + for (int i = 0; i < 100; i++) + { + double k = pow(10.0, 4.0*i/100.-2); + ff << k << " " << power(k) << endl; + } + // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); + gsl_integration_qag(&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr); gsl_integration_workspace_free(w); normVal /= (2*M_PI*M_PI); @@ -226,6 +233,16 @@ void CosmoPower::updateCosmology() beta = pow(OMEGA_0, 5./9); OmegaEff = OMEGA_0 * h * h; Gamma0 = OMEGA_0 * h * h; + omega_B = OMEGA_B * h * h; + omega_C = OMEGA_C * h * h; +} + +void CosmoPower::updatePhysicalCosmology() +{ + OMEGA_B = omega_B / (h*h); + OMEGA_C = omega_C / (h*h); + OMEGA_0 = Gamma0 / (h*h); + beta = pow(OMEGA_0, 5./9); } double CosmoPower::eval_theta_theta(double k) @@ -282,3 +299,8 @@ void CosmoPower::setFunction(CosmoFunction f) abort(); } } + +void CosmoPower::setNormalization(double A_K) +{ + normPower = A_K/power(0.002); +} diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp index b700fc0..1edf5b6 100644 --- a/src/cosmopower.hpp +++ b/src/cosmopower.hpp @@ -17,6 +17,8 @@ namespace CosmoTool { double SIGMA8; double OMEGA_B; double OMEGA_C; + double omega_B; + double omega_C; double Theta_27; // DERIVED VARIABLES @@ -44,7 +46,9 @@ namespace CosmoTool { void setFunction(CosmoFunction f); void updateCosmology(); + void updatePhysicalCosmology(); void normalize(); + void setNormalization(double A_K); double eval_theta_theta(double k); diff --git a/src/fourier/details/euclidian_spectrum_1d.hpp b/src/fourier/details/euclidian_spectrum_1d.hpp index 6e2c21e..fbc4f7f 100644 --- a/src/fourier/details/euclidian_spectrum_1d.hpp +++ b/src/fourier/details/euclidian_spectrum_1d.hpp @@ -1,6 +1,7 @@ #ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D #define __DETAILS_EUCLIDIAN_SPECTRUM_1D +#include #include diff --git a/src/h5_readFlash.cpp b/src/h5_readFlash.cpp index 0c78f58..d8f6475 100644 --- a/src/h5_readFlash.cpp +++ b/src/h5_readFlash.cpp @@ -185,7 +185,7 @@ void h5_read_flash3_particles (H5File* file, float *vel1, float *vel2, float *vel3, - int *id) + long *id) { herr_t status; @@ -341,7 +341,7 @@ void h5_read_flash3_particles (H5File* file, if (id) { for(p=0; p < (pcount); p++) { - id[p+poffset] = (int) *(partBuffer+iptag-1+p*numProps); + id[p+poffset] = (long) *(partBuffer+iptag-1+p*numProps); } } if (pos1 && pos2 && pos3) { diff --git a/src/h5_readFlash.hpp b/src/h5_readFlash.hpp index 5187232..71392bc 100644 --- a/src/h5_readFlash.hpp +++ b/src/h5_readFlash.hpp @@ -29,7 +29,7 @@ void h5_read_flash3_particles (H5File* file, float *vel1, float *vel2, float *vel3, - int *id); + long *id); void h5_read_flash3_header_info(H5File* file, double* time, /* simulation time */ diff --git a/src/loadFlash.cpp b/src/loadFlash.cpp index ada498d..f00f098 100644 --- a/src/loadFlash.cpp +++ b/src/loadFlash.cpp @@ -71,7 +71,7 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags) } } if (loadflags & NEED_GADGET_ID) { - data->Id = new int[data->NumPart]; + data->Id = new long[data->NumPart]; if (data->Id == 0) { delete data; return 0; diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index ee30d80..ced42b6 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -195,7 +195,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, try { f->beginCheckpoint(); - data->Id = new int[data->NumPart]; + data->Id = new long[data->NumPart]; if (data->Id == 0) { delete f; diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp index 5cd2e1e..e05f426 100644 --- a/src/loadRamses.cpp +++ b/src/loadRamses.cpp @@ -384,7 +384,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI } if (flags & NEED_GADGET_ID) { - data->Id = new int[nPar]; + data->Id = new long[nPar]; } for (int k = 0; k < 3; k++) diff --git a/src/loadSimu.hpp b/src/loadSimu.hpp index 7a72259..e93fd64 100644 --- a/src/loadSimu.hpp +++ b/src/loadSimu.hpp @@ -38,7 +38,7 @@ namespace CosmoTool long NumPart; long TotalNumPart; - int *Id; + long *Id; float *Pos[3]; float *Vel[3]; int *type;