From fdd12c3d9c0eab0e966a6936204242397a2dfe72 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 16 Dec 2010 12:05:03 +0100 Subject: [PATCH] minor fixes (mainly constraints) --- constraints.cc | 193 +++++++++++++++++++++++++++++++++++++++-- log.cc | 2 + plugins/output_enzo.cc | 2 +- 3 files changed, 189 insertions(+), 8 deletions(-) diff --git a/constraints.cc b/constraints.cc index 951e316..47d8411 100644 --- a/constraints.cc +++ b/constraints.cc @@ -6,10 +6,172 @@ Copyright (C) 2010 Oliver Hahn + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + */ #include "constraints.hh" +double dsigma2_tophat( double k, void *params ); +double dsigma2_gauss( double k, void *params ); +double find_coll_z( const std::vector& z, const std::vector& sigma, double nu ); +void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, std::vector& z, std::vector& sigma ); +void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std::vector& z, std::vector& sigma ); + + + +double dsigma2_tophat( double k, void *pparams ) +{ + if( k<=0.0 ) + return 0.0; + + char **params = reinterpret_cast (pparams); + + transfer_function *ptf = reinterpret_cast(params[0]); + double x = k * *reinterpret_cast(params[1]); + double nspect = *reinterpret_cast(params[2]); + + double w = 3.0*(sin(x)-x*cos(x))/(x*x*x); + + double tfk = ptf->compute(k,total); + return k*k * w*w * pow(k,nspect) * tfk*tfk; +} + +double dsigma2_gauss( double k, void *pparams ) +{ + if( k<=0.0 ) + return 0.0; + + char **params = reinterpret_cast (pparams); + + transfer_function *ptf = reinterpret_cast (params[0]); + double x = k * *reinterpret_cast(params[1]); + double nspect = *reinterpret_cast(params[2]); + + double w = exp(-x*x*0.5); + + double tfk = ptf->compute(k,total); + return k*k * w*w * pow(k,nspect) * tfk*tfk; +} + +double find_coll_z( const std::vector& z, const std::vector& sigma, double nu ) +{ + double dcoll = 1.686/nu; + double zcoll = 0.0; + + gsl_interp_accel *acc = gsl_interp_accel_alloc (); + gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, z.size()); + + gsl_spline_init (spline, &sigma[0], &z[0], z.size() ); + + zcoll = gsl_spline_eval(spline, dcoll, acc ); + + gsl_spline_free (spline); + gsl_interp_accel_free (acc); + + return zcoll; +} + + +void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, std::vector& z, std::vector& sigma ) +{ + z.clear(); + sigma.clear(); + + cosmology cosm( cf ); + CosmoCalc ccalc( cosm, ptf ); + + double zmin = 0.0, zmax = 200.0; + int nz = 100; + for( int i=0; i ("cosmology","sigma_8"); + double nspec = cf.getValue("cosmology","nspec"); + + double sigma0 = 0.0; + { + double eight=8.0; + + char *params[3]; + params[0] = reinterpret_cast (ptf); + params[1] = reinterpret_cast (&eight); + params[2] = reinterpret_cast (&nspec); + + sigma0 = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast(params) )); + } + + for( int i=0; i (ptf); + params[1] = reinterpret_cast (&R); + params[2] = reinterpret_cast (&nspec); + + double sig = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast(params) )); + double Dz = ccalc.CalcGrowthFactor(1./(1.+z[i])); + sigma.push_back( sig*sigma8/sigma0*Dz/D0 ); + } +} + +void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std::vector& z, std::vector& sigma ) +{ + z.clear(); + sigma.clear(); + + cosmology cosm( cf ); + CosmoCalc ccalc( cosm, ptf ); + + double zmin = 0.0, zmax = 200.0; + int nz = 100; + for( int i=0; i ("cosmology","sigma_8"); + double nspec = cf.getValue("cosmology","nspec"); + + double sigma0 = 0.0; + { + double eight=8.0; + + char *params[3]; + params[0] = reinterpret_cast (ptf); + params[1] = reinterpret_cast (&eight); + params[2] = reinterpret_cast (&nspec); + + sigma0 = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast(params) )); + } + + for( int i=0; i (ptf); + params[1] = reinterpret_cast (&R); + params[2] = reinterpret_cast (&nspec); + + double sig = sqrt(4.0*M_PI*integrate( &dsigma2_gauss, 1e-4, 1e4, reinterpret_cast(params) )); + double Dz = ccalc.CalcGrowthFactor(1./(1.+z[i])); + + //std::cerr << z[i] << " " << sig << std::endl; + sigma.push_back( sig*sigma8/sigma0*Dz/D0 ); + } +} + constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) : pcf_( &cf ), ptf_( ptf ) @@ -26,11 +188,11 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) constr_level_ = std::max(constr_level_,levelmin_TF); - //double omegam = pcf_->getValue("cosmology","Omega_m"); - //double rhom = omegam*2.77519737e11; //... mean matter density + double omegam = pcf_->getValue("cosmology","Omega_m"); + double rhom = omegam*2.77519737e11; //... mean matter density //... use EdS density for estimation - double rhom = 2.77519737e11; + //double rhom = 2.77519737e11; std::map< std::string, constr_type> constr_type_map; constr_type_map.insert( std::pair("halo",halo) ); @@ -74,13 +236,30 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) else if( new_c.type == peak ) { //... peak type constraints take a scale and a peak height - sprintf(temp1,"constraint[%u].Rg",i); - new_c.Rg = cf.getValue( "constraints", temp1 ); + //sprintf(temp1,"constraint[%u].Rg",i); + //new_c.Rg = cf.getValue( "constraints", temp1 ); + //double mass = pow(new_c.Rg,3.0)*rhom*pow(2.*M_PI,1.5); + sprintf(temp1,"constraint[%u].mass",i); + double mass = cf.getValue( "constraints", temp1 ); + new_c.Rg = pow((mass/pow(2.*M_PI,1.5)/rhom),1./3.); + double Rtophat = pow(mass/4.0*3.0/M_PI/rhom,1./3.); sprintf(temp1,"constraint[%u].nu",i); - new_c.sigma = cf.getValue( "constraints", temp1 ); + double nu = cf.getValue( "constraints", temp1 ); - LOGINFO("Constraint %d : peak with Rg=%g h-1 Mpc and nu = %g",i,new_c.Rg,new_c.sigma); + std::vector z,sigma; + compute_sigma_tophat( cf, ptf, Rtophat, z, sigma ); + double zcoll = find_coll_z( z, sigma, nu ); + + //LOGINFO("Probable collapse redshift for constraint %d : z = %f @ M = %g", i, zcoll,mass ); + + compute_sigma_gauss( cf, ptf, new_c.Rg, z, sigma ); + new_c.sigma = nu*sigma.back(); + + //LOGINFO("Constraint %d : peak with Rg=%g h-1 Mpc and nu = %g",i,new_c.Rg,new_c.sigma); + LOGINFO("Constraint %3d : peak",i); + LOGINFO(" M = %g h-1 M_o, nu = %.2f sigma", mass, nu ); + LOGINFO(" estimated z_coll = %f, sigma = %f", zcoll, new_c.sigma ); } diff --git a/log.cc b/log.cc index 5a67d58..f0b019c 100644 --- a/log.cc +++ b/log.cc @@ -12,6 +12,8 @@ #include #include +std::string RemoveMultipleWhiteSpaces( std::string s ); + std::string MUSIC::log::outputFile_; std::ofstream MUSIC::log::outputStream_; diff --git a/plugins/output_enzo.cc b/plugins/output_enzo.cc index cfbde2a..9277cda 100644 --- a/plugins/output_enzo.cc +++ b/plugins/output_enzo.cc @@ -296,7 +296,7 @@ public: std::vector data; data.reserve( ng[0]*ng[1]*ng[2] ); - double dx = 1.0/(1<