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

Timestamp:

Jul 24, 2010, 12:01:34 AM (15 years ago)

Author:

cmv

Message:

Refonte totale de la structure des classes InitialSpectrum, TransfertFunction, GrowthFactor, PkSpectrum et classes tabulate pour lecture des fichiers CAMB, cmv 24/07/2010

File:

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

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -16 +16 @@
 
 ///////////////////////////////////////////////////////////
-//******************** InitialSpectrum ******************//
-///////////////////////////////////////////////////////////
-
-InitialSpectrum::InitialSpectrum(double n,double a)
+//******************** InitialPowerLaw ******************//
+///////////////////////////////////////////////////////////
+
+InitialPowerLaw::InitialPowerLaw(double n,double a)
   : n_(n), A_(a)
 {
 }
 
-InitialSpectrum::InitialSpectrum(InitialSpectrum& pkinf)
+InitialPowerLaw::InitialPowerLaw(InitialPowerLaw& pkinf)
   : n_(pkinf.n_), A_(pkinf.A_)
 {
 }
 
-InitialSpectrum::~InitialSpectrum(void)
+InitialPowerLaw::~InitialPowerLaw(void)
 {
 }
@@ -298 +298 @@
 : kmin_(1.) , kmax_(-1.) , interptyp_(0)
 {
+  k_.resize(0);
+  tf_.resize(0);
 }
 
@@ -331 +333 @@
  char *line = new char[lenline];
 
- int nread = 0;
+ k_.resize(0); tf_.resize(0);
  double tmax = -1.;
  while ( fgets(line,lenline,file) != NULL ) {
@@ -341 +343 @@
    k_.push_back(k);
    tf_.push_back(tf);
-   nread++;
  }
 
- cout<<"TransfertTabulate::ReadCMBFast: nread="<<nread<<" tf_max="<<tmax<<endl;
+ cout<<"TransfertTabulate::ReadCMBFast: nread="<<tf_.size()<<" tf_max="<<tmax<<endl;
  delete [] line;
- if(nread==0) return nread;
+ if(tf_.size()==0) return (int)tf_.size();
 
  for(unsigned int i=0;i<tf_.size();i++) tf_[i] /= tmax;
 
- return nread;
-}
+ return (int)tf_.size();
+}
+
+int TransfertTabulate::ReadCAMB(string filename, double h100)
+{
+ FILE *file = fopen(filename.c_str(),"r");
+ if(file==NULL) return -1;
+ cout<<"TransfertTabulate::ReadCAMB: fn="<<filename<<" h100="<<h100<<endl;
+
+ const int lenline = 512;
+ char *line = new char[lenline];
+
+ k_.resize(0); tf_.resize(0);
+ double tmax = -1.;
+ while ( fgets(line,lenline,file) != NULL ) {
+   double k,tcdm,tbar,tph,trel,tnu,ttot, tf;
+   sscanf(line,"%lf %lf %lf %lf %lf %lf %lf",&k,&tcdm,&tbar,&tph,&trel,&tnu,&ttot);
+   k *= h100;     // convert h Mpc^-1  ->  Mpc^-1
+   tf = ttot;
+   if(tf>tmax) tmax = tf;
+   k_.push_back(k);
+   tf_.push_back(tf);
+ }
+
+ cout<<"TransfertTabulate::ReadCAMB nread="<<tf_.size()<<" tf_max="<<tmax<<endl;
+ delete [] line;
+ if(tf_.size()==0) return (int)tf_.size();
+
+ for(unsigned int i=0;i<tf_.size();i++) tf_[i] /= tmax;
+
+ return (int)tf_.size();
+}
+
 
 ///////////////////////////////////////////////////////////
 //********************* GrowthFactor ********************//
+///////////////////////////////////////////////////////////
+double GrowthFactor::DsDz(double z, double)
+{
+  cout<<"Error: GrowthFactor::DsDz not implemented"<<endl;
+  throw AllocationError("Error: GrowthFactor::DsDz not implemented");
+}
+
+///////////////////////////////////////////////////////////
+//********************* GrowthEisenstein ********************//
 ///////////////////////////////////////////////////////////
 
 // From Eisenstein & Hu ApJ 496:605-614 1998 April 1
 // Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
-GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
+GrowthEisenstein::GrowthEisenstein(double OmegaMatter0,double OmegaLambda0)
   : O0_(OmegaMatter0) , Ol_(OmegaLambda0)
 {
   if(OmegaMatter0==0.) {
-    cout<<"GrowthFactor::GrowthFactor: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
-    throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
-  }
-}
-
-GrowthFactor::GrowthFactor(GrowthFactor& d1)
+    cout<<"GrowthEisenstein::GrowthEisenstein: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
+    throw ParmError("GrowthEisenstein::GrowthEisenstein: Error badOmegaMatter0  value");
+  }
+}
+
+GrowthEisenstein::GrowthEisenstein(GrowthEisenstein& d1)
   : O0_(d1.O0_) , Ol_(d1.Ol_)
 {
 }
 
-GrowthFactor::~GrowthFactor(void)
-{
-}
-
-double GrowthFactor::operator() (double z)
+GrowthEisenstein::~GrowthEisenstein(void)
+{
+}
+
+double GrowthEisenstein::operator() (double z)
 // see Formulae A4 + A5 + A6 page 614
 {
@@ -399 +440 @@
 }
 
-double GrowthFactor::DsDz(double z,double dzinc)
+double GrowthEisenstein::DsDz(double z,double dzinc)
 // y-y0 = a*(x-x0)^2 + b*(x-x0)
 // dy = a*dx^2 + b*dx   avec  dx = x-x0
@@ -405 +446 @@
 {
  if(z<0. || dzinc<=0.) {
-   cout<<"GrowthFactor::DsDz: z<0 or dzinc<=0. !"<<endl;
-   throw ParmError("GrowthFactor::DsDz: z<0 or dzinc<=0. !");
+   cout<<"GrowthEisenstein::DsDz: z<0 or dzinc<=0. !"<<endl;
+   throw ParmError("GrowthEisenstein::DsDz: z<0 or dzinc<=0. !");
  }
 
@@ -431 +472 @@
 }
 
-void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
+void GrowthEisenstein::SetParTo(double OmegaMatter0,double OmegaLambda0)
 {
  if(OmegaMatter0>0.) O0_ = OmegaMatter0;
@@ -437 +478 @@
 }
 
-bool GrowthFactor::SetParTo(double OmegaMatter0)
+bool GrowthEisenstein::SetParTo(double OmegaMatter0)
 // idem precedent sans changer OmegaLambda0
 {
@@ -447 +488 @@
 
 ///////////////////////////////////////////////////////////
-//************** PkSpectrum0 et PkSpectrumZ *************//
-///////////////////////////////////////////////////////////
-
-PkSpectrum0::PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf)
-  : pkinf_(pkinf) , tf_(tf)
-{
-}
-
-PkSpectrum0::PkSpectrum0(PkSpectrum0& pk0)
-  : pkinf_(pk0.pkinf_) , tf_(pk0.tf_)
-{
-}
-
-PkSpectrum0::~PkSpectrum0(void)
-{
-}
-
-double PkSpectrum0::operator() (double k)
+//********************** PkSpectrum *********************//
+///////////////////////////////////////////////////////////
+
+PkSpectrum::PkSpectrum(void)
+  : zref_(0.) , scale_(1.) , typspec_(PK)
+{
+}
+
+PkSpectrum::PkSpectrum(PkSpectrum& pk)
+  : zref_(pk.zref_) , scale_(pk.scale_) , typspec_(pk.typspec_)
+{
+}
+
+
+///////////////////////////////////////////////////////////
+//********************** PkSpecCalc *********************//
+///////////////////////////////////////////////////////////
+
+PkSpecCalc::PkSpecCalc(InitialSpectrum& pkinf,TransfertFunction& tf,GrowthFactor& d1,double zref)
+  : pkinf_(pkinf) , tf_(tf) , d1_(d1)
+{
+  zref_ = zref;
+}
+
+PkSpecCalc::PkSpecCalc(PkSpecCalc& pkz)
+  : pkinf_(pkz.pkinf_) , tf_(pkz.tf_) , d1_(pkz.d1_)
+{
+}
+
+PkSpecCalc::~PkSpecCalc(void) 
+{
+}
+
+double PkSpecCalc::operator() (double k)
+{
+  return (*this)(k,zref_);
+}
+
+double PkSpecCalc::operator() (double k,double z)
 {
   double tf = tf_(k);
   double pkinf = pkinf_(k);
-  return pkinf *tf*tf;
-}
-
-//------------------------------------
-PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
-  : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(PK)
-{
-}
-
-PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
-  : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(PK)
-{
-}
-
-PkSpectrumZ::~PkSpectrumZ(void) 
-{
-}
-
-void PkSpectrumZ::SetTypSpec(ReturnSpectrum typspec)
-// typsec = PK : compute Pk(k)
-//        = DELTA : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
-{
-  typspec_ = typspec;
-}
-
-double PkSpectrumZ::operator() (double k)
+  double d1 = d1_(z);
+
+  double v = pkinf * (tf*tf) * (d1*d1);
+  if(typspec_ == DELTA) v *= k*k*k/(2.*M_PI*M_PI);
+
+  return scale_ * v;
+}
+
+///////////////////////////////////////////////////////////
+//********************** PkTabulate *********************//
+///////////////////////////////////////////////////////////
+
+PkTabulate::PkTabulate(void)
+: kmin_(1.) , kmax_(-1.) , interptyp_(0)
+{
+  k_.resize(0);
+  pk_.resize(0);
+}
+
+PkTabulate::PkTabulate(PkTabulate& pkz)
+ : kmin_(pkz.kmin_) , kmax_(pkz.kmax_) , interptyp_(pkz.interptyp_) , k_(pkz.k_) , pk_(pkz.pk_)
+{
+}
+
+
+PkTabulate::~PkTabulate(void) 
+{
+}
+
+void PkTabulate ::SetInterpTyp(int typ)
+// see comment in InterpTab
+{
+ if(typ<0) typ=0; else if(typ>2) typ=2;
+ interptyp_ = typ;
+}
+
+double PkTabulate::operator() (double k)
+{
+  double v = InterpTab(k,k_,pk_,interptyp_);
+  if(typspec_ == DELTA) v *= k*k*k/(2.*M_PI*M_PI);
+  return scale_ * v;
+}
+
+double PkTabulate::operator() (double k,double z)
+{
+  cout<<"Error: PkTabulate::operator(double k,double z) not implemented"<<endl;
+  throw AllocationError("Error: PkTabulate::operator(double k,double z) not implemented");
+}
+
+int PkTabulate::ReadCAMB(string filename, double h100, double zreftab)
+{
+ FILE *file = fopen(filename.c_str(),"r");
+ if(file==NULL) return -1;
+ cout<<"PkTabulate::ReadCAMB: fn="<<filename<<" h100="<<h100<<" zreftab = "<<zreftab<<endl;
+
+ const int lenline = 512;
+ char *line = new char[lenline];
+
+ k_.resize(0); pk_.resize(0);
+ while ( fgets(line,lenline,file) != NULL ) {
+   double k, pk;
+   sscanf(line,"%lf %lf",&k,&pk);
+   k *= h100;             // convert h    Mpc^-1  ->  Mpc^-1
+   pk /= h100*h100*h100;  // convert h^-3 Mpc^3   ->  Mpc^3
+   k_.push_back(k);
+   pk_.push_back(pk);
+ }
+
+ SetZ(zreftab);
+ cout<<"PkTabulate::ReadCAMB nread="<<pk_.size()<<"zref="<<GetZ()<<endl;
+ delete [] line;
+
+ return (int)pk_.size();
+}
+
+///////////////////////////////////////////////////////////
+//********************* PkEisenstein ********************//
+///////////////////////////////////////////////////////////
+
+PkEisenstein::PkEisenstein(InitialPowerLaw& pkinf,TransfertEisenstein& tf,GrowthEisenstein& d1,double zref)
+  : pkinf_(pkinf) , tf_(tf) , d1_(d1)
+{
+  zref_ = zref;
+}
+
+PkEisenstein::PkEisenstein(PkEisenstein& pkz)
+  : pkinf_(pkz.pkinf_) , tf_(pkz.tf_) , d1_(pkz.d1_)
+{
+}
+
+PkEisenstein::~PkEisenstein(void) 
+{
+}
+
+double PkEisenstein::operator() (double k)
 {
   return (*this)(k,zref_);
 }
 
-double PkSpectrumZ::operator() (double k,double z)
-{
+double PkEisenstein::operator() (double k,double z)
+{
+  double tf = tf_(k);
+  double pkinf = pkinf_(k);
   double d1 = d1_(z);
 
-  double v = pk0_(k) * d1*d1;
-  if(typspec_==DELTA) v *= k*k*k/(2.*M_PI*M_PI);
+  double v = pkinf * (tf*tf) * (d1*d1);
+  if(typspec_ == DELTA) v *= k*k*k/(2.*M_PI*M_PI);
 
   return scale_ * v;
 }
-