Changeset 3378 in Sophya for trunk/Cosmo/SimLSS/pkspectrum.cc

Timestamp:

Nov 9, 2007, 7:04:12 PM (18 years ago)

Author:

cmv

Message:

modif sur remplissage NTuple/Minos, possibilite de changer les parametres O0,Om,Ol etc... , cmv 09/11/2007

File:

• trunk/Cosmo/SimLSS/pkspectrum.cc (modified) (12 diffs)

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -51 +51 @@
 TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon,int lp)
   : lp_(lp)
-  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
+  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h100_(h100) , tcmb_(tcmb)
   , nobaryon_(nobaryon) , nooscenv_(0), retpart_(ALL)
 {
@@ -60 +60 @@
 TransfertEisenstein::TransfertEisenstein(TransfertEisenstein& tf)
   : lp_(tf.lp_)
-  ,Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_)
+  ,Oc_(tf.Oc_) , Ob_(tf.Ob_) , h100_(tf.h100_) , tcmb_(tf.tcmb_)
   , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_), retpart_(tf.retpart_)
 {
@@ -67 +67 @@
 }
 
+TransfertEisenstein::~TransfertEisenstein(void)
+{
+}
+
 void TransfertEisenstein::zero_(void)
 {
@@ -78 +82 @@
   O0_ = Oc_ + Ob_;
   if(nobaryon_) {O0_ = Oc_; Ob_ = 0.;}
-  double H0 = 100. * h_, h2 = h_*h_;
-  if(lp_) cout<<"h100="<<h_<<" H0="<<H0<<") Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
+  double H0 = 100. * h100_, h2 = h100_*h100_;
+  if(lp_) cout<<"h100="<<h100_<<" H0="<<H0<<") Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
 
 
@@ -173 +177 @@
 }
 
-TransfertEisenstein::~TransfertEisenstein(void)
-{
+bool TransfertEisenstein::SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon)
+// Changement des valeurs des parametres (suivi de re-init eventuel)
+{
+ bool haschanged = false;
+
+ if(h100>0.) {h100_ = h100; haschanged = true;}
+ if(OmegaCDM0>0.) {Oc_ = OmegaCDM0; haschanged = true;}
+ if(OmegaBaryon0>0.) {Ob_ = OmegaBaryon0; haschanged = true;}
+ if(tcmb>0.) {tcmb_ = tcmb; haschanged = true;}
+ if(nobaryon!=nobaryon_) {nobaryon_ = nobaryon; haschanged = true;}
+
+ // et on recalcule les initialisations
+ if(haschanged) Init_();
+
+ return haschanged;
 }
 
@@ -202 +219 @@
 }
 
+void TransfertEisenstein::SetPrintLevel(int lp)
+{
+ lp_ = lp;
+}
+
+
 double TransfertEisenstein::T0tild(double k,double alphac,double betac)
 {
   // Formule 10 p 608
-  //double q = k*th2p7_*th2p7_/(O0_*h_*h_);
+  //double q = k*th2p7_*th2p7_/(O0_*h100_*h100_);
   double q = k/(13.41*keq_);
   // Formule 20 p 610
@@ -220 +243 @@
   //  OU Pour function lissee sans oscillation baryon
   if(nobaryon_  || nooscenv_ == 2) {
-    double gamma = O0_*h_;
+    double gamma = O0_*h100_;
     // Calcul de Gamma_eff, formule 30 page 612 (pour fct lissee)
     if( nobaryon_==false && nooscenv_ == 2 )
-      gamma = O0_*h_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*sfit_,4.))); // Gamma_eff
+      gamma = O0_*h100_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*sfit_,4.))); // Gamma_eff
     // Formule 28 page 612 : qui est est equivalent a:
-    //    q = k / h_ * th2p7_*th2p7_ / gamma;
+    //    q = k / h100_ * th2p7_*th2p7_ / gamma;
     // qui est est equivalent a:
     //    q = k / (13.41 * keq)                   pour Ob=0
@@ -231 +254 @@
     // Les resultats sont legerement differents a cause des valeurs approx.
     // des constantes numeriques: on prend comme W.Hu (tf_fit.c)
-    //double q = k / h_ * th2p7_*th2p7_ / gamma;  // Mpc^-1
-    double q = k/(13.41*keq_) * (O0_*h_/gamma);  // Mpc^-1
+    //double q = k / h100_ * th2p7_*th2p7_ / gamma;  // Mpc^-1
+    double q = k/(13.41*keq_) * (O0_*h100_/gamma);  // Mpc^-1
     // Formules 29 page 612
     double l0 = log(2.*M_E + 1.8*q);
@@ -284 +307 @@
 
 TransfertTabulate::TransfertTabulate(double h100,double OmegaCDM0,double OmegaBaryon0)
-: Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , kmin_(1.) , kmax_(-1.)
+: Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h100_(h100) , kmin_(1.) , kmax_(-1.)
 , interptyp_(0)
 {
@@ -290 +313 @@
 
 TransfertTabulate::TransfertTabulate(TransfertTabulate& tf)
-: Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , kmin_(tf.kmin_) , kmax_(tf.kmax_)
+: Oc_(tf.Oc_) , Ob_(tf.Ob_) , h100_(tf.h100_) , kmin_(tf.kmin_) , kmax_(tf.kmax_)
 , interptyp_(tf.interptyp_) , k_(tf.k_) , tf_(tf.tf_)
 {
@@ -324 +347 @@
    double k,tc,tb,tf;
    sscanf(line,"%lf %lf %lf",&k,&tc,&tb);
-   k *= h_;     // convert h^-1 Mpc  ->  Mpc
+   k *= h100_;     // convert h^-1 Mpc  ->  Mpc
    tf = (Oc_*tc+Ob_*tb)/(Oc_+Ob_);
    if(tf>tmax) tmax = tf;
@@ -382 +405 @@
 
  return D1z / norm_;
+}
+
+bool GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
+{
+ bool haschanged = false;
+
+ if(OmegaMatter0>0.) {O0_ = OmegaMatter0; haschanged = true;}
+ if(fabs(OmegaLambda0+12345.)>1e-6) {Ol_ = OmegaLambda0; haschanged = true;}
+
+ // et on recalcule les initialisations
+ if(haschanged) {
+   Ok_ = 1. - O0_ - Ol_;
+   norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
+   norm_ = (*this)(0.);
+ }
+
+ return haschanged;
 }