Changeset 3806 in Sophya

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

Location:

trunk/Cosmo/SimLSS

Files:

• cmvginit3d.cc (modified) (18 diffs)
• genefluct3d.cc (modified) (4 diffs)
• pkspectrum.cc (modified) (10 diffs)
• pkspectrum.h (modified) (1 diff)

trunk/Cosmo/SimLSS/cmvginit3d.cc

@@ -20 +20 @@
 #include "geneutils.h"
 #include "genefluct3d.h"
+
+string decodecambarg(string arg,double& h, double& z);
 
 void usage(void);
@@ -36 +38 @@
      <<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
      <<" -Z zref : redshift for the center of the simulation cube"<<endl
+     <<" -f cambspec.dat,h100tab,ztab : use CAMB file (def: Eisenstein spectrum)"<<endl
      <<" -1 : compute 1D spectrum (default: no)"<<endl
      <<" -2 : compute 2D spectrum (default: no)"<<endl
@@ -69 +72 @@
  unsigned short flat = 0;
  double ob0 = 0.0444356;
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73, w0=-1.;
  double zref = 0.5;
  double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
@@ -78 +81 @@
  double sigmaR = 1.;
 
- double kmin=1e-5,kmax=1000.;
- int npt = 10000;
- double lkmin=log10(kmin), lkmax=log10(kmax);
+ double kmin=1e-4,kmax=100.;
+ int npt = 1000;
  double eps=1.e-3;
+
+ // *** Spectrum read from CAMB file
+ string cambfread = "";
 
  // *** type de generation
@@ -115 +120 @@
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"haG:F:x:y:z:Z:128:v:n:CO:So:VT:")) != -1) {
+ while((c = getopt(narg,arg,"haG:F:x:y:z:Z:128:v:n:CO:So:VT:f:")) != -1) {
   int nth = 0;
   switch (c) {
@@ -175 +180 @@
     nthread = (nth<1)? 0: nth;
     break;
+  case 'f' :
+    cambfread = optarg;
+    break;
   case 'h' :
   default :
@@ -189 +197 @@
  cout<<"zref="<<zref<<endl;
  cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
- cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"<<", npt="<<npt<<endl;
+ cout<<"kmin="<<kmin<<" ("<<log10(kmin)<<"), kmax="<<kmax<<" ("<<log10(kmax)<<") Mpc^-1"<<", npt="<<npt<<endl;
  cout<<"Filter by pixel = "<<filter_by_pixel<<endl;
  if(comptransveloc) cout<<"Tansverse velocity generation requested"<<endl;
+ if(cambfread.size()>0) cout<<"use CAMB file: cambfread="<<cambfread<<endl;
  cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
  cout<<"Use_growth_factor = "<<use_growth_factor<<endl;
@@ -216 +225 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Create Spectrum"<<endl;
-
- InitialPowerLaw pkini(ns,as);
-
- TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
- TransfertEisenstein tfnos(h100,om0-ob0,ob0,T_CMB_Par,false);
-   tfnos.SetNoOscEnv(2);
+ cout<<endl<<"\n--- Create GrowthFactor"<<endl;
 
  GrowthEisenstein growth(om0,ol0);
@@ -229 +232 @@
  cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
 
- PkSpecCalc pkz(pkini,tf,growth,zref);
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Spectrum"<<endl;
+
+ InitialPowerLaw pkini(ns,as);
+
+ TransfertEisenstein tfnos(h100,om0-ob0,ob0,T_CMB_Par,false);
+   tfnos.SetNoOscEnv(2);
+ TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+
  PkSpecCalc pkznos(pkini,tfnos,growth,zref);
-
- //-----------------------------------------------------------------
- pkz.SetZ(0.);
+ PkSpectrum* pkz = NULL;
+ if(cambfread.size()>0) {  // read pk in CAMB file
+   double h100tab, ztab;
+   string fn = decodecambarg(cambfread,h100tab,ztab);
+   if(h100tab<=0.) h100tab = h100;
+   PkTabulate* pkzt = new PkTabulate;
+   pkzt->ReadCAMB(fn,h100tab,ztab);
+   pkzt->SetGrowthFactor(&growth);
+   pkzt->SetZ(zref);
+   pkz = pkzt;
+   // change k limits
+   kmin = pkzt->KMin();
+   kmax = pkzt->KMax();
+   cout<<"Change k limits to kmin="<<kmin<<" kmax="<<kmax<<endl;
+ } else {   // use Eisenstein pk
+   PkSpecCalc* pkze = new PkSpecCalc(pkini,tf,growth,zref);
+   pkz = pkze;
+ }
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Normalize Spectrum"<<endl;
+
+ pkz->SetZ(0.);
  pkznos.SetZ(0.);
- cout<<endl<<"\n--- Compute variance for top-hat R="<<R<<" at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
- double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
- double pkmax_th = varpk_th(kfind_th);
- cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
- double k1=kmin, k2=kmax;
- int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
- cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
-     <<k2<<" ("<<log10(k2)<<")"<<endl;
-
- double ldlk = (log10(k2)-log10(k1))/npt;
- varpk_th.SetInteg(0.01,ldlk,-1.,4);
- double sr2 = varpk_th.Variance(k1,k2);
- cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
-
- double normpkz = sigmaR*sigmaR/sr2;
- pkz.SetScale(normpkz);
- pkznos.SetScale(normpkz);
- cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
+
+ double normpkz[2] = {0.,0.};
+ for(int i=0;i<2;i++) {
+   PkSpectrum* pkvar = NULL; string nam = "";
+   if(i==0) pkvar = pkz; else {pkvar = &pkznos; nam = "NoOsc";}
+   cout<<endl<<"\n--- Compute variance for Pk "<<nam<<" with top-hat R="<<R<<" at z="<<pkvar->GetZ()<<endl;
+   VarianceSpectrum varpk_th(*pkvar,R,VarianceSpectrum::TOPHAT);
+   double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
+   double pkmax_th = varpk_th(kfind_th);
+   cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
+   double k1=kmin, k2=kmax;
+   int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
+   cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
+       <<k2<<" ("<<log10(k2)<<")"<<endl;
+   double ldlk = (log10(k2)-log10(k1))/npt;
+   varpk_th.SetInteg(0.01,ldlk,-1.,4);
+   double sr2 = varpk_th.Variance(k1,k2);
+   normpkz[i] = sigmaR*sigmaR/sr2;
+   cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<" normpkz="<<normpkz[i]<<endl;
+ }
+ cout<<endl;
+
+ if(cambfread.size()>0) {
+   pkz->SetScale(1.);
+   pkznos.SetScale(normpkz[1]);
+   cout<<"Warning: CAMB spectrum, no normalisation applied, scale="<<pkz->GetScale()<<endl;
+ } else {
+   pkz->SetScale(normpkz[0]);
+   pkznos.SetScale(normpkz[0]);
+   cout<<"Spectrum normalisation (osc+noosc), scale = "<<pkz->GetScale()<<endl;
+ }
+
  {
- Histo hpkz0(lkmin,lkmax,npt); hpkz0.ReCenterBin();
- FuncToHisto(pkz,hpkz0,true);
+ Histo hpkz0(log10(kmin),log10(kmax),npt); hpkz0.ReCenterBin();
+ FuncToHisto(*pkz,hpkz0,true);
  posobs.PutObject(hpkz0,"hpkz0");
  Histo hpkz0nos(hpkz0);
@@ -263 +307 @@
  }
 
- pkz.SetZ(zref);
+ pkz->SetZ(zref);
  pkznos.SetZ(zref);
+
  {
- Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
- FuncToHisto(pkz,hpkz,true);
+ Histo hpkz(log10(kmin),log10(kmax),npt); hpkz.ReCenterBin();
+ FuncToHisto(*pkz,hpkz,true);
  posobs.PutObject(hpkz,"hpkz");
  Histo hpkznos(hpkz);
@@ -275 +320 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
+ cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz->GetZ()<<endl;
+ VarianceSpectrum varpk_int(*pkz,R,VarianceSpectrum::NOFILTER);
  double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
  double pkmax_int = varpk_int(kfind_int);
@@ -325 +370 @@
 
  //-----------------------------------------------------------------
- cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
+ cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz->GetZ()<<endl;
  // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
  // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
@@ -339 +384 @@
  //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
- if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
- cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
- fluct3d.ComputeFourier(pkz);
+ if(use_growth_factor>0) pkz->SetZ(0.); else pkz->SetZ(zref);
+ cout<<"Power spectrum set at redshift: "<<pkz->GetZ()<<endl;
+ fluct3d.ComputeFourier(*pkz);
  fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
  PrtTim(">>>> End Computing a realization in Fourier space");
@@ -411 +456 @@
  //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in real space"<<endl;
+ double vdum1,vdum2;
  fluct3d.ComputeReal();
- double rmin,rmax; fluct3d.MinMax(rmin,rmax);
- cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+ fluct3d.MinMax(vdum1,vdum2);
+ fluct3d.MeanSigma2(vdum1,vdum2);
  fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
  PrtTim(">>>> End Computing a realization in real space");
@@ -421 +467 @@
    cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
    fluct3d.ApplyGrowthFactor(use_growth_factor);
-   fluct3d.MinMax(rmin,rmax);
-   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.MinMax(vdum1,vdum2);
+   fluct3d.MeanSigma2(vdum1,vdum2);
    fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
    PrtTim(">>>> End Applying growth factor");
  }
-
- int_8 nm;
- double rmref,rs2ref;
- cout<<"\n--- Computing reference variance in real space"<<endl;
- nm = fluct3d.MeanSigma2(rmref,rs2ref);
- cout<<" rs2ref= "<<rs2ref<<" , rmref="<<rmref<<" ("<<nm<<")"<<endl;
- PrtTim(">>>> End Computing reference variance in real space");
 
  if(whattowrt[1]&1) {
@@ -561 +600 @@
    cout<<"\n--- Computing a realization in real space for Transverse velocity"<<endl;
    fluct3d.ComputeReal();
-   double rmin,rmax; fluct3d.MinMax(rmin,rmax);
-   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.MinMax(vdum1,vdum2);
+   fluct3d.MeanSigma2(vdum1,vdum2);
    fluct3d.NTupleCheck(posobs,string("nvreal"),ntnent);
    PrtTim(">>>> End Computing a realization in real space");
 
     if(use_growth_factor>0) {
-      cout<<"\n--- Apply Growth factor"<<endl;
+      cout<<"\n--- Apply Growth factor for transverse velocity"<<endl;
       cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
       fluct3d.ApplyVelLosGrowthFactor(use_growth_factor);
-      fluct3d.MinMax(rmin,rmax);
-      cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+      fluct3d.MinMax(vdum1,vdum2);
+      fluct3d.MeanSigma2(vdum1,vdum2);
       fluct3d.NTupleCheck(posobs,string("nvgrow"),ntnent);
       PrtTim(">>>> End Applying growth factor");
     }
-
-    int_8 nmv;
-    double vmref,vs2ref;
-    cout<<"\n--- Computing reference variance in real space"<<endl;
-    nmv = fluct3d.MeanSigma2(vmref,vs2ref);
-    cout<<" vs2ref= "<<vs2ref<<" , vmref="<<vmref<<" ("<<nmv<<")"<<endl;
-    PrtTim(">>>> End Computing reference variance in real space");
 
     if(whattowrt[1]&1) {
@@ -610 +642 @@
  dvl("s8") = sigmaR;
  dvl("Dtype") = (int_4)use_growth_factor; dvl("Pxfilt") = (int_4)filter_by_pixel; 
+ if(cambfread.size()>0) dvl("fCAMB") = cambfread;
  posobs.PutObject(dvl,"siminfo");
  }
 
  delete RandGen;
+ if(pkz != NULL) delete pkz;
  PrtTim(">>>> End Of Job");
 
@@ -634 +668 @@
 }
 
+#include "strutilxx.h"
+string decodecambarg(string arg,double& h, double& z)
+// - decode argument for CAMB generation:
+// input: arg = "cambfilename,h100tab,ztab"
+// output: h = h100 use to generate CAMB file
+//         z = redshift use to generate CAMB file
+// return: CAMB file name
+{
+  h = -1.;
+  z = 0.;
+  vector<string> vs;
+  FillVStringFrString(arg,vs,',');
+  if(vs.size()>1) h = atof(vs[1].c_str());
+  if(vs.size()>2) z = atof(vs[2].c_str());
+  cout<<"decodecambarg: "<<arg<<" fn=\""<<vs[0]<<"\" h="<<h<<" z="<<z<<endl;
+  return vs[0];
+}

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -978 +978 @@
 //             (tous les pixels d'un plan Z sont mis au meme redshift z que celui du milieu)
 {
- if(lp_>0) cout<<"--- ApplyGrowthFactor: evol="<<type_evol<<endl;
+  if(lp_>0) cout<<"--- ApplyGrowthFactor: evol="<<type_evol<<" (zpk="<<compute_pk_redsh_ref_<<")"<<endl;
  check_array_alloc();
 
@@ -1035 +1035 @@
 
  double zpk = compute_pk_redsh_ref_;
- double dpsd_orig = - cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+ double grw_orig = (*growth_)(zpk);
+ double dpsd_orig = - cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / grw_orig;
+ if(lp_>0) cout<<"    original growth="<<grw_orig<<" dpsd="<<dpsd_orig<<" computed at z="<<zpk<<endl;
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
@@ -1049 +1051 @@
          else dz = fabs(dz); // tous les plans Z au meme redshift
        double z = interpinv(dz);   // interpolation par morceau
+       double grw = (*growth_)(z);
        double dpsd = interdpd(z);
        int_8 ip = IndexR(i,j,l);
-       data_[ip] *= dpsd / dpsd_orig;
+       // on remet le beta au bon z
+       // on corrige du growth factor car data_ a ete calcule avec pk(zpk)
+       data_[ip] *= (dpsd / dpsd_orig) * (grw / grw_orig);
      }
    }
@@ -1062 +1067 @@
 // Calcule une realisation dans l'espace reel
 {
- if(lp_>0) cout<<"--- ComputeReal ---"<<endl;
- check_array_alloc();
-
- // On fait la FFT
- GEN3D_FFTW_EXECUTE(pb_);
- array_type = 1;
+  if(lp_>0) cout<<"--- ComputeReal --- from spectrum at z="<<compute_pk_redsh_ref_<<endl;
+  check_array_alloc();
+
+  // On fait la FFT
+  GEN3D_FFTW_EXECUTE(pb_);
+  array_type = 1;
 }
 

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");
   }
 

trunk/Cosmo/SimLSS/pkspectrum.h

@@ -154 +154 @@
   virtual double operator() (double k);
   virtual double operator() (double k,double z);
+  virtual void SetZ(double z);
   int NPoints(void) {return k_.size();}
   void SetInterpTyp(int typ=0);
-  int ReadCAMB(string filename, double h100=0.71, double zreftab = 0.);
+  int ReadCAMB(string filename, double h100tab=0.71, double zreftab=0.);
+  void SetGrowthFactor(GrowthFactor* d1) {d1_ = d1;}
+  GrowthFactor* GetGrowthFactor(void) {return d1_;}
+  double KMin(void) {if(k_.size()!=0) return k_[0]; else return 0.;}
+  double KMax(void) {if(k_.size()!=0) return k_[k_.size()-1]; else return 0.;}
 protected:
   double kmin_,kmax_;
   int interptyp_;
   vector<double> k_, pk_;
+  GrowthFactor* d1_;
 };