Merge branch 'master' of bitbucket.org:glavaux/cosmotool

This commit is contained in:
Guilhem Lavaux 2015-04-29 10:18:15 +02:00
commit c50679c31d
5 changed files with 38 additions and 54 deletions

View File

@ -17,7 +17,8 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
POWER_BARDEEN "CosmoTool::CosmoPower::POWER_BARDEEN", POWER_BARDEEN "CosmoTool::CosmoPower::POWER_BARDEEN",
POWER_SUGIYAMA "CosmoTool::CosmoPower::POWER_SUGIYAMA", POWER_SUGIYAMA "CosmoTool::CosmoPower::POWER_SUGIYAMA",
POWER_BDM, POWER_BDM,
POWER_TEST POWER_TEST,
HU_WIGGLES_ORIGINAL "CosmoTool::CosmoPower::HU_WIGGLES_ORIGINAL"
cdef cppclass CosmoPower: cdef cppclass CosmoPower:
double n double n
@ -113,6 +114,8 @@ cdef class CosmologyPower:
f = POWER_BARDEEN f = POWER_BARDEEN
elif funcname=='SUGIYAMA': elif funcname=='SUGIYAMA':
f = POWER_SUGIYAMA f = POWER_SUGIYAMA
elif funcname=='HU_WIGGLES_ORIGINAL':
f = HU_WIGGLES_ORIGINAL
self.power.setFunction(f) self.power.setFunction(f)

View File

@ -9,7 +9,6 @@ SET(CosmoTool_SRCS
growthFactor.cpp growthFactor.cpp
cosmopower.cpp cosmopower.cpp
cic.cpp cic.cpp
tf_fit.c
) )
IF (Boost_FOUND) IF (Boost_FOUND)

View File

@ -41,6 +41,7 @@ knowledge of the CeCILL license and that you accept its terms.
#include <fstream> #include <fstream>
#include <gsl/gsl_integration.h> #include <gsl/gsl_integration.h>
#include "cosmopower.hpp" #include "cosmopower.hpp"
#include "tf_fit.hpp"
using namespace std; using namespace std;
using namespace CosmoTool; using namespace CosmoTool;
@ -69,9 +70,17 @@ CosmoPower::CosmoPower()
Theta_27 = 2.728/2.7; Theta_27 = 2.728/2.7;
ehu_params = 0;
updateCosmology(); updateCosmology();
} }
CosmoPower::~CosmoPower()
{
if (ehu_params)
delete ehu_params;
}
/* /*
* This is \hat{tophat} * This is \hat{tophat}
*/ */
@ -109,23 +118,6 @@ static double powG(double y)
return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1))); 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) double CosmoPower::powerEfstathiou(double k)
{ {
double a = 6.4/Gamma0; double a = 6.4/Gamma0;
@ -141,7 +133,10 @@ double CosmoPower::powerEfstathiou(double k)
void CosmoPower::updateHuWigglesOriginal() void CosmoPower::updateHuWigglesOriginal()
{ {
TFset_parameters( (OMEGA_C+OMEGA_B)*h*h, if (ehu_params == 0)
ehu_params = new TF_Transfer();
ehu_params->TFset_parameters( (OMEGA_C+OMEGA_B)*h*h,
OMEGA_B/(OMEGA_C+OMEGA_B), Theta_27*2.7); OMEGA_B/(OMEGA_C+OMEGA_B), Theta_27*2.7);
} }
@ -368,6 +363,13 @@ double CosmoPower::power(double k)
double CosmoPower::powerHuWigglesOriginal(double k) double CosmoPower::powerHuWigglesOriginal(double k)
{ {
float tfb, tfc;
ehu_params->TFfit_onek(k, &tfb, &tfc);
double T_k = OMEGA_B/OMEGA_0 * tfb + OMEGA_C/OMEGA_0 * tfc;
return normPower * pow(k,n) * T_k * T_k;
} }
void CosmoPower::setFunction(CosmoFunction f) void CosmoPower::setFunction(CosmoFunction f)

View File

@ -38,6 +38,8 @@ knowledge of the CeCILL license and that you accept its terms.
namespace CosmoTool { namespace CosmoTool {
struct TF_Transfer;
class CosmoPower class CosmoPower
{ {
public: public:
@ -74,6 +76,7 @@ namespace CosmoTool {
}; };
EHuParams ehu; EHuParams ehu;
TF_Transfer *ehu_params;
enum CosmoFunction { enum CosmoFunction {
POWER_EFSTATHIOU, POWER_EFSTATHIOU,
@ -88,6 +91,7 @@ namespace CosmoTool {
}; };
CosmoPower(); CosmoPower();
~CosmoPower();
void setFunction(CosmoFunction f); void setFunction(CosmoFunction f);

View File

@ -14,17 +14,6 @@ calculate other quantities given in Section 4 of the paper. */
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
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 --------------------------- */ /* ------------------------ DRIVER ROUTINE --------------------------- */
/* The following is an example of a driver routine you might use. */ /* The following is an example of a driver routine you might use. */
@ -65,26 +54,6 @@ 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 but rather compiling it separately, calling only the driver, and using
extern declarations to access the intermediate quantities. */ 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<numk;j++) {
tf_thisk = TFfit_onek(k[j]*hubble, &baryon_piece, &cdm_piece);
if (tf_full!=NULL) tf_full[j] = tf_thisk;
if (tf_baryon!=NULL) tf_baryon[j] = baryon_piece;
if (tf_cdm!=NULL) tf_cdm[j] = cdm_piece;
}
return;
}
/* ------------------------ FITTING FORMULAE ROUTINES ----------------- */ /* ------------------------ FITTING FORMULAE ROUTINES ----------------- */
/* There are two routines here. TFset_parameters() sets all the scalar /* There are two routines here. TFset_parameters() sets all the scalar
@ -97,6 +66,9 @@ global variables in case you wish to access them, e.g. by declaring
them as extern variables in your main program. */ them as extern variables in your main program. */
/* Note that all internal scales are in Mpc, without any Hubble constants! */ /* Note that all internal scales are in Mpc, without any Hubble constants! */
namespace CosmoTool {
struct TF_Transfer {
float omhh, /* Omega_matter*h^2 */ float omhh, /* Omega_matter*h^2 */
obhh, /* Omega_baryon*h^2 */ obhh, /* Omega_baryon*h^2 */
theta_cmb, /* Tcmb in units of 2.7 K */ theta_cmb, /* Tcmb in units of 2.7 K */
@ -117,11 +89,11 @@ float omhh, /* Omega_matter*h^2 */
alpha_gamma; /* Gamma suppression in approximate TF */ alpha_gamma; /* Gamma suppression in approximate TF */
/* Convenience from Numerical Recipes in C, 2nd edition */ /* Convenience from Numerical Recipes in C, 2nd edition */
static float sqrarg; float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg) #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static float cubearg; float cubearg;
#define CUBE(a) ((cubearg=(a)) == 0.0 ? 0.0 : cubearg*cubearg*cubearg) #define CUBE(a) ((cubearg=(a)) == 0.0 ? 0.0 : cubearg*cubearg*cubearg)
static float pow4arg; float pow4arg;
#define POW4(a) ((pow4arg=(a)) == 0.0 ? 0.0 : pow4arg*pow4arg*pow4arg*pow4arg) #define POW4(a) ((pow4arg=(a)) == 0.0 ? 0.0 : pow4arg*pow4arg*pow4arg*pow4arg)
/* Yes, I know the last one isn't optimal; it doesn't appear much */ /* Yes, I know the last one isn't optimal; it doesn't appear much */
@ -335,3 +307,7 @@ and omega0 -> omega0*hubble^2. */
T_0_C0 = 14.2 + 731.0/(1+62.5*q); T_0_C0 = 14.2 + 731.0/(1+62.5*q);
return T_0_L0/(T_0_L0+T_0_C0*q*q); return T_0_L0/(T_0_L0+T_0_C0*q*q);
} }
};
};