diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 9131f69..4b659be 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -36,7 +36,7 @@ target_link_libraries(testPool ${tolink}) if (HDF5_FOUND) include_directories(${HDF5_INCLUDE_PATH}) - SET(tolink ${tolink} ${ZLIB}) + SET(tolink ${tolink} ${HDF5_CPP_LIBRARY} ${HDF5_LIBRARY} ${ZLIB}) add_executable(testReadFlash testReadFlash.cpp) target_link_libraries(testReadFlash ${tolink}) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 659d4cc..21f2691 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,6 +10,7 @@ SET(CosmoTool_SRCS growthFactor.cpp cosmopower.cpp cic.cpp + tf_fit.c ) IF (ENABLE_OPENMP) diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp index 1593db1..eb36b2c 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -108,6 +108,23 @@ static double powG(double y) return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1))); } + + +extern "C" { + void TFset_parameters(float omega0hh, float f_baryon, float Tcmb); + float TFfit_onek(float k, float *tf_baryon, float *tf_cdm); + + void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb, + int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm); + + float TFsound_horizon_fit(float omega0, float f_baryon, float hubble); + float TFk_peak(float omega0, float f_baryon, float hubble); + float TFnowiggles(float omega0, float f_baryon, float hubble, + float Tcmb, float k_hmpc); + float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc); +} + + double CosmoPower::powerEfstathiou(double k) { @@ -122,6 +139,12 @@ double CosmoPower::powerEfstathiou(double k) return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu)); } +void CosmoPower::updateHuWigglesOriginal() +{ + TFset_parameters( (OMEGA_C+OMEGA_B)*h*h, + OMEGA_B/(OMEGA_C+OMEGA_B), Theta_27*2.7); +} + void CosmoPower::updateHuWigglesConsts() { double f_b = OMEGA_B / OMEGA_0; @@ -343,6 +366,9 @@ double CosmoPower::power(double k) return (this->*eval)(k); } +double CosmoPower::powerHuWigglesOriginal(double k) +{ +} void CosmoPower::setFunction(CosmoFunction f) { @@ -355,6 +381,10 @@ void CosmoPower::setFunction(CosmoFunction f) updateHuWigglesConsts(); eval = &CosmoPower::powerHuWiggles; break; + case HU_WIGGLES_ORIGINAL: + updateHuWigglesOriginal(); + eval = &CosmoPower::powerHuWigglesOriginal; + break; case HU_BARYON: eval = &CosmoPower::powerHuBaryons; break; @@ -382,3 +412,4 @@ void CosmoPower::setNormalization(double A_K) { normPower = A_K;///power(0.002); } + diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp index b32c6a7..21ef45a 100644 --- a/src/cosmopower.hpp +++ b/src/cosmopower.hpp @@ -83,7 +83,8 @@ namespace CosmoTool { POWER_BARDEEN, POWER_SUGIYAMA, POWER_BDM, - POWER_TEST + POWER_TEST, + HU_WIGGLES_ORIGINAL }; CosmoPower(); @@ -95,6 +96,7 @@ namespace CosmoTool { void normalize(double k_min = -1, double k_max = -1); void setNormalization(double A_K); void updateHuWigglesConsts(); + void updateHuWigglesOriginal(); double eval_theta_theta(double k); double power(double k); @@ -111,7 +113,7 @@ namespace CosmoTool { double powerSugiyama(double k); double powerBDM(double k); double powerTest(double k); - + double powerHuWigglesOriginal(double k); }; }; diff --git a/src/loadSimu.hpp b/src/loadSimu.hpp index 6f848ce..2eb5ac8 100644 --- a/src/loadSimu.hpp +++ b/src/loadSimu.hpp @@ -36,6 +36,7 @@ knowledge of the CeCILL license and that you accept its terms. #ifndef __COSMOTOOLBOX_HPP #define __COSMOTOOLBOX_HPP +#include #include #include diff --git a/src/tf_fit.c b/src/tf_fit.c new file mode 100644 index 0000000..5a59649 --- /dev/null +++ b/src/tf_fit.c @@ -0,0 +1,337 @@ +/* The following routines implement all of the fitting formulae in +Eisenstein \& Hu (1997) */ + +/* There are two sets of routines here. The first set, + + TFfit_hmpc(), TFset_parameters(), and TFfit_onek(), + +calculate the transfer function for an arbitrary CDM+baryon universe using +the fitting formula in Section 3 of the paper. The second set, + + TFsound_horizon_fit(), TFk_peak(), TFnowiggles(), and TFzerobaryon(), + +calculate other quantities given in Section 4 of the paper. */ + +#include +#include +void TFset_parameters(float omega0hh, float f_baryon, float Tcmb); +float TFfit_onek(float k, float *tf_baryon, float *tf_cdm); + +void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb, + int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm); + +float TFsound_horizon_fit(float omega0, float f_baryon, float hubble); +float TFk_peak(float omega0, float f_baryon, float hubble); +float TFnowiggles(float omega0, float f_baryon, float hubble, + float Tcmb, float k_hmpc); +float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc); + +/* ------------------------ DRIVER ROUTINE --------------------------- */ +/* The following is an example of a driver routine you might use. */ +/* Basically, the driver routine needs to call TFset_parameters() to +set all the scalar parameters, and then call TFfit_onek() for each +wavenumber k you desire. */ + +/* While the routines use Mpc^-1 units internally, this driver has been +written to take an array of wavenumbers in units of h Mpc^-1. On the +other hand, if you want to use Mpc^-1 externally, you can do this by +altering the variables you pass to the driver: + omega0 -> omega0*hubble*hubble, hubble -> 1.0 */ + +/* INPUT: omega0 -- the matter density (baryons+CDM) in units of critical + f_baryon -- the ratio of baryon density to matter density + hubble -- the Hubble constant, in units of 100 km/s/Mpc + Tcmb -- the CMB temperature in Kelvin. T<=0 uses the COBE value 2.728. + numk -- the length of the following zero-offset array + k[] -- the array of wavevectors k[0..numk-1] */ + +/* INPUT/OUTPUT: There are three output arrays of transfer functions. +All are zero-offset and, if used, must have storage [0..numk-1] declared +in the calling program. However, if you substitute the NULL pointer for +one or more of the arrays, then that particular transfer function won't +be outputted. The transfer functions are: + + tf_full[] -- The full fitting formula, eq. (16), for the matter + transfer function. + tf_baryon[] -- The baryonic piece of the full fitting formula, eq. 21. + tf_cdm[] -- The CDM piece of the full fitting formula, eq. 17. */ + +/* Again, you can set these pointers to NULL in the function call if +you don't want a particular output. */ + +/* Various intermediate scalar quantities are stored in global variables, +so that you might more easily access them. However, this also means that +you would be better off not simply #include'ing this file in your programs, +but rather compiling it separately, calling only the driver, and using +extern declarations to access the intermediate quantities. */ + +/* ----------------------------- DRIVER ------------------------------- */ + +void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb, + int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm) +/* Remember: k[0..numk-1] is in units of h Mpc^-1. */ +{ + int j; + float tf_thisk, baryon_piece, cdm_piece; + + TFset_parameters(omega0*hubble*hubble, f_baryon, Tcmb); + + for (j=0;j 1 +and omega0 -> omega0*hubble^2. */ +{ + float omhh, sound_horizon_fit_mpc; + omhh = omega0*hubble*hubble; + sound_horizon_fit_mpc = + 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(omhh*f_baryon,0.75)); + return sound_horizon_fit_mpc*hubble; +} + +float TFk_peak(float omega0, float f_baryon, float hubble) +/* Input: omega0 -- CDM density, in units of critical density + f_baryon -- Baryon fraction, the ratio of baryon to CDM density. + hubble -- Hubble constant, in units of 100 km/s/Mpc +/* Output: The approximate location of the first baryonic peak, in h Mpc^-1 */ +/* Note: If you prefer to have the answer in units of Mpc^-1, use hubble -> 1 +and omega0 -> omega0*hubble^2. */ +{ + float omhh, k_peak_mpc; + omhh = omega0*hubble*hubble; + k_peak_mpc = 2.5*3.14159*(1+0.217*omhh)/TFsound_horizon_fit(omhh,f_baryon,1.0); + return k_peak_mpc/hubble; +} + +float TFnowiggles(float omega0, float f_baryon, float hubble, + float Tcmb, float k_hmpc) +/* Input: omega0 -- CDM density, in units of critical density + f_baryon -- Baryon fraction, the ratio of baryon to CDM density. + hubble -- Hubble constant, in units of 100 km/s/Mpc + Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of + COBE FIRAS value of 2.728 K + k_hmpc -- Wavenumber in units of (h Mpc^-1). */ +/* Output: The value of an approximate transfer function that captures the +non-oscillatory part of a partial baryon transfer function. In other words, +the baryon oscillations are left out, but the suppression of power below +the sound horizon is included. See equations (30) and (31). */ +/* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1 +and omega0 -> omega0*hubble^2. */ +{ + float k, omhh, theta_cmb, k_equality, q, xx, alpha_gamma, gamma_eff; + float q_eff, T_nowiggles_L0, T_nowiggles_C0; + + k = k_hmpc*hubble; /* Convert to Mpc^-1 */ + omhh = omega0*hubble*hubble; + if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */ + theta_cmb = Tcmb/2.7; + + k_equality = 0.0746*omhh/SQR(theta_cmb); + q = k/13.41/k_equality; + xx = k*TFsound_horizon_fit(omhh, f_baryon, 1.0); + + alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)* + SQR(f_baryon); + gamma_eff = omhh*(alpha_gamma+(1-alpha_gamma)/(1+POW4(0.43*xx))); + q_eff = q*omhh/gamma_eff; + + T_nowiggles_L0 = log(2.0*2.718282+1.8*q_eff); + T_nowiggles_C0 = 14.2 + 731.0/(1+62.5*q_eff); + return T_nowiggles_L0/(T_nowiggles_L0+T_nowiggles_C0*SQR(q_eff)); +} + +/* ======================= Zero Baryon Formula =========================== */ + +float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc) +/* Input: omega0 -- CDM density, in units of critical density + hubble -- Hubble constant, in units of 100 km/s/Mpc + Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of + COBE FIRAS value of 2.728 K + k_hmpc -- Wavenumber in units of (h Mpc^-1). */ +/* Output: The value of the transfer function for a zero-baryon universe. */ +/* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1 +and omega0 -> omega0*hubble^2. */ +{ + float k, omhh, theta_cmb, k_equality, q, T_0_L0, T_0_C0; + + k = k_hmpc*hubble; /* Convert to Mpc^-1 */ + omhh = omega0*hubble*hubble; + if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */ + theta_cmb = Tcmb/2.7; + + k_equality = 0.0746*omhh/SQR(theta_cmb); + q = k/13.41/k_equality; + + T_0_L0 = log(2.0*2.718282+1.8*q); + T_0_C0 = 14.2 + 731.0/(1+62.5*q); + return T_0_L0/(T_0_L0+T_0_C0*q*q); +}