#include "sopnamsp.h" #include "machdefs.h" #include #include #include #include #include #include #include "pexceptions.h" #include "constcosmo.h" #include "geneutils.h" #include "pkspectrum.h" /////////////////////////////////////////////////////////// //******************** InitialSpectrum ******************// /////////////////////////////////////////////////////////// InitialSpectrum::InitialSpectrum(double n,double a) : n_(n), A_(a) { } InitialSpectrum::InitialSpectrum(InitialSpectrum& pkinf) : n_(pkinf.n_), A_(pkinf.A_) { } InitialSpectrum::~InitialSpectrum(void) { } void InitialSpectrum::SetNorm(double a) { A_ = a; } void InitialSpectrum::SetSlope(double n) { n_ = n; } /////////////////////////////////////////////////////////// //****************** TransfertEisenstein ****************// /////////////////////////////////////////////////////////// // From Eisenstein & Hu ApJ 496:605-614 1998 April 1 TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon) : Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb) , nobaryon_(nobaryon) , nooscenv_(0) { if(nobaryon_ == false && Ob_>0.) Init_With_Baryons(); else Init_Without_Baryon(); } TransfertEisenstein::TransfertEisenstein(TransfertEisenstein& tf) : Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_) , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_) { if(nobaryon_ == false && Ob_>0.) Init_With_Baryons(); else Init_Without_Baryon(); } void TransfertEisenstein::Init_With_Baryons(void) { int lp = 1; nobaryon_ = false; O0_ = Oc_ + Ob_; if(lp) cout<<"h100="<0.038&&k<0.07) cout<<"k="<2) { cout<<"VarianceSpectrum::SetFilter: Error bad value for type of filter"<0 a l'ordre O(x^6) // DL(x) = 1-x^2*(1-x^2/2) // pour x<0.01 |DL(x)-G(X)^2|<2.0e-13 if(typfilter_ == 1) if(x<0.01) return 1.-x2*(1.-x2/2.); else return exp(-x2); // Filtre top-hat T(x) = 3*(sin(x)-x*cos(x))/x^3 // --- Gestion de la pseudo-divergence pour x->0 // on prend le DL de T(x)^2 pour x->0 a l'ordre O(x^7) // DL(x) = 1-x^2/5*(1-3*x^2/35*(1-4*x^2/81)) // pour x<0.1 |DL(x)-T(X)^2|<2.5e-13 double f2=0.; if(x<0.1) { f2 = 1.-x2/5.*(1.-3.*x2/35.*(1.-4.*x2/81.)); } else { f2 = 3.*(sin(x)-x*cos(x))/(x2*x); f2 *= f2; } return f2; } double VarianceSpectrum::Variance(double R,double kmin,double kmax) // Compute variance of spectrum pk_ by integration // Input: // R = taille du filter top-hat ou gaussien // kmin,kmax = bornes en k de l'integrale pour calculer la variance // Return: // valeur de la variance (sigma^2) // Remarque: // la meilleure approximation du filtre top-hat (R) est un filtre gaussien avec (Rg=R/sqrt(5)) // la variance renvoyee est la variance de la masse { if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) { cout<<"VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax"<=kmax"); } R_ = R; double lkmin = log10(kmin), lkmax = log10(kmax); double var = IntegrateFuncLog(*this,lkmin,lkmax,dperc_,dlogkinc_,dlogkmax_,glorder_); return var; } //------------------------------------ double VarianceSpectrum::FindMaximum(double R,double kmin,double kmax,double eps) // Retourne le maximum de la fonction a integrer // La recherche a lieu entre [kmin,kmax] par pas logarithmiques // Input: // R : taille du filter top-hat ou gaussien // kmin,kmax : intervalle de recherche // eps : precision requise sur les valeurs // Return: // position (en k) du maximum { if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) { cout<<"VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<=kmax || eps<=0"); } R_ = R; int n = 10; // toujours >2 double lkmin = log10(kmin), lkmax = log10(kmax), dlk = (lkmax-lkmin)/n; double lkfind=lkmin, pkfind=-1.; while(1) { for(int i=0; i<=n; i++) { double lk = lkmin + i*dlk; double v = (*this)(pow(10.,lk)); if(vlkmin) lkmin = lkfind-dlk; if(lkfind+dlk= high" (bit0 =1) // rc |= 2 "f(kmax) >= high" (bit1 =1) // rc |= 4 "f(k) < high pour tout k" (bit2 =1) { if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax || eps<=0.) { cout<<"VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<=kmax || eps<=0"); } R_ = R; int n = 10; // toujours >2 int rc = 0; double lkmin,lkmax,dlk,lkfind; // --- Find kmin lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n; while(1) { lkfind = lkmin; for(int i=0;i<=n;i++) { if( (*this)(pow(10,lkfind)) >= high ) break; //cmvbug lkfind = lkmin + i*dlk; } //cout<<"VarianceSpectrum::FindLimits[kmin]: lkfind="<