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

Timestamp:

Jul 24, 2010, 6:15:07 PM (15 years ago)

Author:

cmv

Message:

suite de la mise au point pour lecture fichiers CAMB, cmv 24/07/2010

File:

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

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -390 +390 @@
 double GrowthFactor::DsDz(double z, double)
 {
-  cout<<"Error: GrowthFactor::DsDz not implemented"<<endl;
-  throw AllocationError("Error: GrowthFactor::DsDz not implemented");
+  cout<<"GrowthFactor::DsDz_Error not implemented"<<endl;
+  throw AllocationError(" GrowthFactor::DsDz_Error not implemented");
 }
 
@@ -404 +404 @@
 {
   if(OmegaMatter0==0.) {
-    cout<<"GrowthEisenstein::GrowthEisenstein: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
-    throw ParmError("GrowthEisenstein::GrowthEisenstein: Error badOmegaMatter0  value");
+    cout<<"GrowthEisenstein::GrowthEisenstein_Error: bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
+    throw ParmError("GrowthEisenstein::GrowthEisenstein_Error:  bad OmegaMatter0  value");
   }
 }
@@ -446 +446 @@
 {
  if(z<0. || dzinc<=0.) {
-   cout<<"GrowthEisenstein::DsDz: z<0 or dzinc<=0. !"<<endl;
-   throw ParmError("GrowthEisenstein::DsDz: z<0 or dzinc<=0. !");
+   cout<<"GrowthEisenstein::DsDz_Error: z<0 or dzinc<=0. !"<<endl;
+   throw ParmError("GrowthEisenstein::DsDz_Error: z<0 or dzinc<=0. !");
  }
 
@@ -543 +543 @@
 
 PkTabulate::PkTabulate(void)
-: kmin_(1.) , kmax_(-1.) , interptyp_(0)
+  : kmin_(1.) , kmax_(-1.) , interptyp_(0), d1_(NULL)
 {
   k_.resize(0);
@@ -550 +550 @@
 
 PkTabulate::PkTabulate(PkTabulate& pkz)
- : kmin_(pkz.kmin_) , kmax_(pkz.kmax_) , interptyp_(pkz.interptyp_) , k_(pkz.k_) , pk_(pkz.pk_)
-{
-}
-
+ : kmin_(pkz.kmin_) , kmax_(pkz.kmax_) , interptyp_(pkz.interptyp_)
+ , k_(pkz.k_) , pk_(pkz.pk_)
+ , d1_(pkz.d1_)
+{
+}
 
 PkTabulate::~PkTabulate(void) 
@@ -575 +576 @@
 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)
+  cout<<"PkTabulate::operator(double k,double z)_Error: not implemented"<<endl;
+  throw AllocationError("PkTabulate::operator(double k,double z)_Error: not implemented");
+}
+
+void PkTabulate::SetZ(double z)
+{
+  if(d1_ == NULL) {
+    cout<<"PkTabulate::SetZ_Error: d1==NULL, no possible redshift change for tabulated Pk"<<endl;
+    throw ParmError("PkTabulate::SetZ_Error: d1==NULL, no possible redshift change for tabulated Pk");
+  }
+  if(pk_.size() == 0) {
+    cout<<"PkTabulate::SetZ_Error: pk_.size()==0, no possible redshift change for tabulated Pk"<<endl;
+    throw ParmError("PkTabulate::SetZ_Error: pk_.size()==0, no possible redshift change for tabulated Pk");
+  }
+
+  double zold = zref_;
+  if(fabs(z-zold)<1.e-4) return;
+
+  zref_ = z;
+  double d0 = (*d1_)(zold);
+  double d1 = (*d1_)(zref_);
+  double conv = d1*d1 / (d0*d0);
+  cout<<"PkTabulate::SetZ: change redshift from "<<zold<<" (d="<<d0
+      <<") to "<<zref_<<" (d="<<d1<<") conv="<<conv<<endl;
+  for(unsigned int i=0;i<pk_.size();i++) pk_[i] *= conv;
+}
+
+int PkTabulate::ReadCAMB(string filename, double h100tab, double zreftab)
 {
  FILE *file = fopen(filename.c_str(),"r");
  if(file==NULL) return -1;
- cout<<"PkTabulate::ReadCAMB: fn="<<filename<<" h100="<<h100<<" zreftab = "<<zreftab<<endl;
+ cout<<"PkTabulate::ReadCAMB: fn="<<filename<<" h100="<<h100tab<<" zreftab = "<<zreftab<<endl;
 
  const int lenline = 512;
  char *line = new char[lenline];
+ double h = h100tab, h3 = pow(h100tab,3.);
 
  k_.resize(0); pk_.resize(0);
+ double kmax = 0., pkmax = 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 *= h;      // convert h    Mpc^-1  ->  Mpc^-1
+   pk /= h3;    // convert h^-3 Mpc^3   ->  Mpc^3
+   if(pk>pkmax) {pkmax = pk; kmax = k;}
    k_.push_back(k);
    pk_.push_back(pk);
  }
 
- SetZ(zreftab);
- cout<<"PkTabulate::ReadCAMB nread="<<pk_.size()<<"zref="<<GetZ()<<endl;
+ zref_ = zreftab;
+ cout<<"  nread="<<pk_.size()<<" zref="<<GetZ()<<" , k,pk: max="<<kmax<<","<<pkmax;
+ if(pk_.size()>0) cout<<" [0]="<<k_[0]<<","<<pk_[0]
+                      <<" [n]="<<k_[pk_.size()-1]<<","<<pk_[pk_.size()-1];
+ cout<<endl;
+
  delete [] line;
 
@@ -667 +698 @@
 {
  if(R<=0.) {
-    cout<<"VarianceSpectrum::SetRadius: Error R<=0"<<endl;
-    throw ParmError("VarianceSpectrum::SetRadius: Error R<=0");
+    cout<<"VarianceSpectrum::SetRadius_Error: R<=0"<<endl;
+    throw ParmError("VarianceSpectrum::SetRadius_Error: R<=0");
   }
   R_ = R;
@@ -737 +768 @@
 {
   if(kmin<=0 || kmax<=0. || kmin>=kmax) {
-    cout<<"VarianceSpectrum::Variance: Error kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
-    throw ParmError("VarianceSpectrum::Variance: Error kmin<=0 or kmax<=0 or kmin>=kmax");
+    cout<<"VarianceSpectrum::Variance_Error: kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
+    throw ParmError("VarianceSpectrum::Variance_Error: kmin<=0 or kmax<=0 or kmin>=kmax");
   }
 
@@ -759 +790 @@
 {
   if(kmin<=0 || kmax<=0. || kmin>=kmax) {
-    cout<<"VarianceSpectrum::FindMaximum: Error kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
-    throw ParmError("VarianceSpectrum::FindMaximum: Error kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
+    cout<<"VarianceSpectrum::FindMaximum_Error: kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
+    throw ParmError("VarianceSpectrum::FindMaximum_Error: kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
   }
 
@@ -802 +833 @@
 {
   if(kmin<=0 || kmax<=0. || kmin>=kmax  || eps<=0.) {
-    cout<<"VarianceSpectrum::FindLimits: Error kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
-    throw ParmError("VarianceSpectrum::FindLimits: Error kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
+    cout<<"VarianceSpectrum::FindLimits_Error: kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
+    throw ParmError("VarianceSpectrum::FindLimits_Error: kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
   }