Changeset 3318 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Aug 28, 2007, 11:41:49 AM (18 years ago)

Author:

cmv

Message:

TransfertTabulate pour avoir les fct de transfert de CMBFast cmv 28/08/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvtransf.cc (modified) (5 diffs)
• cmvtstpk.cc (modified) (1 diff)
• pkspectrum.cc (modified) (2 diffs)
• pkspectrum.h (modified) (1 diff)

trunk/Cosmo/SimLSS/cmvtransf.cc

@@ -14 +14 @@
 
 void usage(void);
-void usage(void) {cout<<"cmvtuniv [k1,k2,npt] [h100,Om0,Ob0,tcmb]"<<endl;}
+void usage(void) {
+cout
+<<"cmvtransf [options]"<<endl
+<<"  -k k1,k2,npt : k range in Mpc^-1"<<endl
+<<"  -U h100,Om0,Ob0,tcmb : cosmology"<<endl
+<<"  -F filename : read also CMBFast file"<<endl
+<<endl;
+}
 
 int main(int narg,char *arg[])
@@ -20 +27 @@
  double k1 = 1e-6, k2 = 10.; int_4 npt = 100;
  double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;
+ string fcmbfast = "";
 
- if(narg>1) sscanf(arg[1],"%lf,%lf,%d",&k1,&k2,&npt);
+ char c;
+  while((c = getopt(narg,arg,"hk:U:F:")) != -1) {
+  switch (c) {
+  case 'k' :
+    sscanf(optarg,"%d,%lf,%lf",&npt,&k1,&k2);
+    break;
+  case 'U' :
+    sscanf(optarg,"%lf,%lf,%lf,%lf",&h100,&Om0,&Ob0,&tcmb);
+    break;
+  case 'F' :
+    fcmbfast = optarg;
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+
  if(k1<=0.) {cout<<"k1 must be >0"<<endl; return -1;}
- if(npt<=0) npt = 100;
+ if(npt<=0) npt = 1000;
+ if(k2<=k1) {k1=1.e-4; k2=1.e+2;}
  cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
- if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%lf",&h100,&Om0,&Ob0,&tcmb);
  cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
+ cout<<"fcmbfast="<<fcmbfast<<endl;
 
  cout<<"TransfertEisenstein with baryon"<<endl;
@@ -43 +69 @@
  TransfertEisenstein tfnob(h100,Om0-Ob0,Ob0,tcmb,true);  // No baryons
 
- const int n = 5;
- char *vname[n] = {"k","t","tnosc1","tnosc2","tnob"};
+ TransfertTabulate tfcmbfast(h100,Om0-Ob0,Ob0);
+ bool cmbfastOK = false;
+ if(fcmbfast.size()>0) {
+   cout<<endl<<"TransfertTabulate for CMBfast"<<endl;
+   tfcmbfast.ReadCMBFast(fcmbfast);
+   tfcmbfast.SetInterpTyp(1);
+   if(tfcmbfast.NPoints()>0) cmbfastOK = true;
+ }
+
+ const int n = 6;
+ char *vname[n] = {"k","t","tnosc1","tnosc2","tnob","tcmbf"};
  NTuple nt(n,vname);
  double xnt[n];
+ for(int i=0;i<n;i++) xnt[i]=0.;
 
  double lnk1 = log10(k1), lnk2=log10(k2), dlnk=(lnk2-lnk1)/npt;
@@ -56 +92 @@
    xnt[3]  = tfnosc2(k);
    xnt[4]  = tfnob(k);
+   if(cmbfastOK) xnt[5]  = tfcmbfast(k);
    nt.Fill(xnt);
  }
@@ -73 +110 @@
 set k log10(k)
 
-n/plot nt.t%$k ! ! "nsta crossmarker3"
-n/plot nt.tnob%$k ! ! "nsta red crossmarker3 same"
-n/plot nt.tnosc1%$k ! ! "nsta blue circlemarker3 same"
-n/plot nt.tnosc2%$k ! ! "nsta green trianglemarker3 same"
+# Eisenstein
+n/plot nt.t%$k ! ! "nsta connectpoints"
+n/plot nt.tnob%$k ! ! "nsta red same connectpoints"
+n/plot nt.tnosc2%$k ! ! "nsta blue same connectpoints"
+n/plot nt.tnosc1%$k ! ! "nsta green same connectpoints"
 
-n/plot nt.t/tnosc1%$k ! ! "nsta crossmarker3"
-n/plot nt.t/tnosc2%$k ! ! "nsta crossmarker3 red same"
+n/plot nt.t/tnob%$k   tnob>0   ! "nsta red connectpoints"
+n/plot nt.t/tnosc2%$k tnosc2>0 ! "nsta blue same connectpoints"
+n/plot nt.t/tnosc1%$k tnosc1>0 ! "nsta green same connectpoints"
+
+# CMBFast
+n/plot nt.t%$k ! ! "nsta connectpoints"
+n/plot nt.tcmbf%$k ! ! "nsta red same connectpoints"
+
+n/plot nt.(tcmbf-t)%$k tcmbf>0. ! "nsta connectpoints"
+n/plot nt.(tcmbf-t)/tcmbf%$k tcmbf>0. ! "nsta connectpoints"
 */

trunk/Cosmo/SimLSS/cmvtstpk.cc

@@ -28 +28 @@
 int main(int narg,char *arg[])
 {
- double Ob0 = 0.0444356;
  // -- WMAP
- double h100=0.71, Om0=0.267804, Ol0=0.73,w0=-1.;
- // -- ouvert matter only
- //double h100=0.71, Om0=0.3, Ol0=0.,w0=-1.;
- // -- plat matter only
- //double h100=0.71, Om0=1., Ol0=0.,w0=-1.;
-
+ double h100=0.71, Om0=0.267804, Ob0 = 0.0444356, Ol0=0.73,w0=-1.;
 
  double ns = 1., as = 1.;
-
  int npt = 10000;
  double lkmin = -3., lkmax=2.;

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -65 +65 @@
  zero_();
  Init_();
+}
+
+void TransfertEisenstein::zero_(void)
+{
+ th2p7_=zeq_=keq_=zd_=Req_=Rd_=s_=ksilk_=alphac_=betac_=bnode_
+       =alphab_=betab_=alphag_=sfit_=kpeak_=1.e99;
 }
 
@@ -275 +281 @@
 
 ///////////////////////////////////////////////////////////
+//******************* TransfertTabulate *****************//
+///////////////////////////////////////////////////////////
+
+TransfertTabulate::TransfertTabulate(double h100,double OmegaCDM0,double OmegaBaryon0)
+: Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , kmin_(1.) , kmax_(-1.)
+, interptyp_(0)
+{
+}
+
+TransfertTabulate::TransfertTabulate(TransfertTabulate& tf)
+: Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , kmin_(tf.kmin_) , kmax_(tf.kmax_)
+, interptyp_(tf.interptyp_) , k_(tf.k_) , tf_(tf.tf_)
+{
+}
+
+TransfertTabulate::~TransfertTabulate(void)
+{
+}
+
+void TransfertTabulate::SetInterpTyp(int typ)
+// see comment in InterpTab
+{
+ if(typ<0) typ=0; else if(typ>2) typ=2;
+ interptyp_ = typ;
+}
+
+double TransfertTabulate::operator() (double k)
+{
+ return InterpTab(k,k_,tf_,interptyp_);
+}
+
+int TransfertTabulate::ReadCMBFast(string filename)
+{
+ FILE *file = fopen(filename.c_str(),"r");
+ if(file==NULL) return -1;
+
+ const int lenline = 512;
+ char *line = new char[lenline];
+
+ int nread = 0;
+ double tmax = -1.;
+ while ( fgets(line,lenline,file) != NULL ) {
+   double k,tc,tb,tf;
+   sscanf(line,"%lf %lf %lf",&k,&tc,&tb);
+   k *= h_;     // convert h^-1 Mpc  ->  Mpc
+   tf = (Oc_*tc+Ob_*tb)/(Oc_+Ob_);
+   if(tf>tmax) tmax = tf;
+   k_.push_back(k);
+   tf_.push_back(tf);
+   nread++;
+ }
+
+ cout<<"TransfertTabulate::ReadCMBFast: nread="<<nread<<" tf_max="<<tmax<<endl;
+ delete [] line;
+ if(nread==0) return nread;
+
+ for(unsigned int i=0;i<tf_.size();i++) tf_[i] /= tmax;
+
+ return nread;
+}
+
+///////////////////////////////////////////////////////////
 //********************* GrowthFactor ********************//
 ///////////////////////////////////////////////////////////

trunk/Cosmo/SimLSS/pkspectrum.h

@@ -43 +43 @@
   double T0tild(double k,double alphac,double betac);
   void Init_(void);
-  inline void zero_(void)
-    {th2p7_=zeq_=keq_=zd_=Req_=Rd_=s_=ksilk_=alphac_=betac_=bnode_=
-     alphab_=betab_=alphag_=sfit_=kpeak_=0.;}
+  void zero_(void);
+};
+
+//-----------------------------------------------------------------------------------
+class TransfertTabulate : public GenericFunc {
+public:
+  TransfertTabulate(double h100,double OmegaCDM0,double OmegaBaryon0);
+  TransfertTabulate(TransfertTabulate& tf);
+  virtual ~TransfertTabulate(void);
+  virtual double operator() (double k);
+  int NPoints(void) {return k_.size();}
+  void SetInterpTyp(int typ=0);
+  int ReadCMBFast(string filename);
+protected:
+  double Oc_,Ob_,h_;
+  double kmin_,kmax_;
+  int interptyp_;
+  vector<double> k_, tf_;
 };