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

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

File:

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

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 ********************//
 ///////////////////////////////////////////////////////////