Index: /trunk/Cosmo/SimLSS/cmvdefsurv.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3805)
@@ -183,5 +183,5 @@
  univ.SetDynParam(h100,om0,or0,ol0,w0);
 
- GrowthFactor growth(om0,ol0);
+ GrowthEisenstein growth(om0,ol0);
 
  // --- 
Index: /trunk/Cosmo/SimLSS/cmvfitpk.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvfitpk.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvfitpk.cc	(revision 3805)
@@ -26,5 +26,5 @@
 class MyFCN : public FCNBase {
 public:
-  MyFCN(HistoErr& herr,PkSpectrumZ& pk);
+  MyFCN(HistoErr& herr,PkEisenstein& pk);
   void Init(const vector<double>& par) const;
   double Func(double x,const vector<double>& par) const;
@@ -33,5 +33,5 @@
 private:
   HistoErr& herr_;
-  PkSpectrumZ& pk_;
+  PkEisenstein& pk_;
 };
 } }  // namespace ROOT + Minuit2
@@ -84,11 +84,10 @@
  // --- Create spectrum
  cout<<endl<<"\n--- Create Spectrum"<<endl;
- InitialSpectrum pkini(ns,as);
+ InitialPowerLaw pkini(ns,as);
  bool nobaryon = false;
  double ocdm0 = om0-ob0;
  TransfertEisenstein tf(h100,ocdm0,ob0,T_CMB_Par,nobaryon);
- GrowthFactor growth(om0,ol0); cout<<"...Growth="<<growth(zref)<<endl;
- PkSpectrum0 pk0(pkini,tf);
- PkSpectrumZ pkz(pk0,growth,zref);
+ GrowthEisenstein growth(om0,ol0); cout<<"...Growth="<<growth(zref)<<endl;
+ PkEisenstein pkz(pkini,tf,growth,zref);
 
  // --- Compute variance and normalize spectrum
@@ -256,5 +255,5 @@
 
 //--------------------------------------------------
-MyFCN::MyFCN(HistoErr& herr,PkSpectrumZ& pk)
+MyFCN::MyFCN(HistoErr& herr,PkEisenstein& pk)
   : herr_(herr) , pk_(pk)
 {
@@ -264,7 +263,7 @@
 {
  double A=par[1], Ocdm0=par[2], Ob0=par[3], Ol0=par[4], h100=par[5], ns=par[6];
- pk_.GetPk0().GetPkIni().SetSlope(ns);
- pk_.GetPk0().GetPkIni().SetNorm(A);
- pk_.GetPk0().GetTransfert().SetParTo(h100,Ocdm0,Ob0);
+ pk_.GetPkIni().SetSlope(ns);
+ pk_.GetPkIni().SetNorm(A);
+ pk_.GetTransfert().SetParTo(h100,Ocdm0,Ob0);
  pk_.GetGrowthFactor().SetParTo(Ocdm0+Ob0,Ol0);
 }
Index: /trunk/Cosmo/SimLSS/cmvginit3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvginit3d.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvginit3d.cc	(revision 3805)
@@ -218,5 +218,5 @@
  cout<<endl<<"\n--- Create Spectrum"<<endl;
 
- InitialSpectrum pkini(ns,as);
+ InitialPowerLaw pkini(ns,as);
 
  TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
@@ -224,14 +224,11 @@
    tfnos.SetNoOscEnv(2);
 
- GrowthFactor growth(om0,ol0);
- // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
+ GrowthEisenstein growth(om0,ol0);
+ // GrowthEisenstein growth(1.,0.); // D(z) = 1/(1+z)
  double growth_at_z = growth(zref);
  cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
 
- PkSpectrum0 pk0(pkini,tf);
- PkSpectrum0 pk0nos(pkini,tfnos);
-
- PkSpectrumZ pkz(pk0,growth,zref);
- PkSpectrumZ pkznos(pk0nos,growth,zref);
+ PkSpecCalc pkz(pkini,tf,growth,zref);
+ PkSpecCalc pkznos(pkini,tfnos,growth,zref);
 
  //-----------------------------------------------------------------
Index: /trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3805)
@@ -299,17 +299,15 @@
  cout<<endl<<"\n--- Create Spectrum"<<endl;
 
- InitialSpectrum pkini(ns,as);
+ InitialPowerLaw pkini(ns,as);
 
  TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
  //tf.SetNoOscEnv(2);
 
- GrowthFactor growth(om0,ol0);
- // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
+ GrowthEisenstein growth(om0,ol0);
+ // GrowthEisenstein growth(1.,0.); // D(z) = 1/(1+z)
  double growth_at_z = growth(zref);
  cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
 
- PkSpectrum0 pk0(pkini,tf);
-
- PkSpectrumZ pkz(pk0,growth,zref);
+ PkSpecCalc pkz(pkini,tf,growth,zref);
 
  //-----------------------------------------------------------------
Index: /trunk/Cosmo/SimLSS/cmvtgrowth.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtgrowth.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvtgrowth.cc	(revision 3805)
@@ -30,5 +30,5 @@
  cout<<"Om0="<<om0<<"  Ol0="<<ol0<<endl;
 
- GrowthFactor growth(om0,ol0);
+ GrowthEisenstein growth(om0,ol0);
  cout<<"D1(z=0) = "<<growth(0.)<<endl;
 
Index: /trunk/Cosmo/SimLSS/cmvtintfun.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtintfun.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvtintfun.cc	(revision 3805)
@@ -44,6 +44,4 @@
   unsigned short glorder=4;
 
-  if(narg>1) sscanf(arg[1],"%lf",&n);
-  if(narg>2) sscanf(arg[2],"%lf,%lf,%d",&xmin,&xmax,&npt);
   cout<<"n="<<n<<" xmin="<<xmin<<" xmax="<<xmax<<" npt="<<npt<<endl;
 
Index: /trunk/Cosmo/SimLSS/cmvtransf.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtransf.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvtransf.cc	(revision 3805)
@@ -17,7 +17,7 @@
 cout
 <<"cmvtransf [options]"<<endl
-<<"  -k k1,k2,npt : k range in Mpc^-1"<<endl
+<<"  -k npt,k1,k2 : k range in Mpc^-1"<<endl
 <<"  -U h100,Om0,Ob0,tcmb : cosmology"<<endl
-<<"  -F filename : read also CMBFast file"<<endl
+<<"  -F filename : read also CAMB file"<<endl
 <<endl;
 }
@@ -27,5 +27,5 @@
  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 = "";
+ string fcmbfile = "";
 
  char c;
@@ -39,5 +39,5 @@
     break;
   case 'F' :
-    fcmbfast = optarg;
+    fcmbfile = optarg;
     break;
   case 'h' :
@@ -52,5 +52,5 @@
  cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
  cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
- cout<<"fcmbfast="<<fcmbfast<<endl;
+ cout<<"fcmbfile="<<fcmbfile<<endl;
 
  cout<<"TransfertEisenstein with baryon"<<endl;
@@ -69,11 +69,11 @@
  TransfertEisenstein tfnob(h100,Om0-Ob0,Ob0,tcmb,true);  // No baryons
 
- TransfertTabulate tfcmbfast;
- bool cmbfastOK = false;
- if(fcmbfast.size()>0) {
-   cout<<endl<<"TransfertTabulate for CMBfast"<<endl;
-   tfcmbfast.ReadCMBFast(fcmbfast,h100,Om0-Ob0,Ob0);
-   tfcmbfast.SetInterpTyp(1);
-   if(tfcmbfast.NPoints()>0) cmbfastOK = true;
+ TransfertTabulate tf_cmbfile;
+ bool cmbfileOK = false;
+ if(fcmbfile.size()>0) {
+   cout<<endl<<"TransfertTabulate for CAMB"<<endl;
+   tf_cmbfile.ReadCAMB(fcmbfile,h100);
+   tf_cmbfile.SetInterpTyp(1);
+   if(tf_cmbfile.NPoints()>0) cmbfileOK = true;
  }
 
@@ -92,5 +92,5 @@
    xnt[3]  = tfnosc2(k);
    xnt[4]  = tfnob(k);
-   if(cmbfastOK) xnt[5]  = tfcmbfast(k);
+   if(cmbfileOK) xnt[5]  = tf_cmbfile(k);
    nt.Fill(xnt);
  }
@@ -107,22 +107,22 @@
 openppf cmvtransf.ppf
 
-set k k
-set k log10(k)
+# Eisenstein
+n/plot nt.t%k t>0&&k>0 ! "nsta connectpoints logx logy"
+n/plot nt.tnob%k tnob>0&&k>0 ! "nsta red same connectpoints logx logy"
+n/plot nt.tnosc2%k tnosc2>0&&k>0 ! "nsta blue same connectpoints logx logy"
+n/plot nt.tnosc1%k tnosc1>0&&k>0 ! "nsta green same connectpoints logx logy"
 
-# 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/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"
 
-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"
+# CAMB
+n/plot nt.t%k t>0&&k>0 ! "nsta connectpoints logx logy"
+n/plot nt.tcmbf%k tcmbf>0&&k>0 ! "nsta red same connectpoints logx logy"
 
-# 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"
 
-n/plot nt.(tcmbf-t)%$k tcmbf>0. ! "nsta connectpoints"
-n/plot nt.(tcmbf-t)/tcmbf%$k tcmbf>0. ! "nsta connectpoints"
+n/plot nt.t/tnosc2%k tnosc2>0&&t>0&&k>0 ! "nsta connectpoints logx logy"
+n/plot nt.tcmbf/tnosc2%k tnosc2>0&&tcmbf>0&&k>0 ! "nsta connectpoints same red logx logy"
 */
Index: /trunk/Cosmo/SimLSS/cmvtstpk.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3805)
@@ -19,9 +19,11 @@
 void usage(void);
 void usage(void) {
- cout<<"cmvtstpk [options] z_redshift"<<endl
-     <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
-     <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
-     <<" -s scale : on normalise le spectre avec scale"<<endl
-     <<" -w : write spectra on ASCII file"<<endl;
+cout
+<<"cmvtstpk [options] z_redshift"<<endl
+<<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
+<<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
+<<" -s scale : on normalise le spectre avec scale"<<endl
+<<" -w : write spectra on ASCII file"<<endl
+<<"  -F filename : read also CAMB file"<<endl;
 }
 
@@ -36,7 +38,8 @@
  double scale = 2.54499e+07;  // normalisation du spectre a z=0 selon SDSS
  bool wrascii = false;
+ string fcmbfile = "";
 
  char c;
-  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:")) != -1) {
+  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:F:")) != -1) {
   switch (c) {
   case 'H' :
@@ -63,4 +66,7 @@
     wrascii = true;
     break;
+  case 'F' :
+    fcmbfile = optarg;
+    break;
   case 'h' :
   default :
@@ -78,7 +84,8 @@
  cout<<"scale="<<scale<<endl;
  cout<<"zval="<<zval<<endl;
-
- //--------------------------
- InitialSpectrum pkini(ns,as);
+ cout<<"fcmbfile="<<fcmbfile<<endl;
+
+ //--------------------------
+ InitialPowerLaw pkini(ns,as);
  pkini.SetNorm(scale);
 
@@ -89,16 +96,25 @@
  TransfertEisenstein tfnob(h100,Om0,0.,T_CMB_Par,true);
 
- GrowthFactor d1(Om0,Ol0);
+ GrowthEisenstein d1(Om0,Ol0);
  cout<<"GrowthFactor: "<<d1(zval)<<endl;
 
- PkSpectrum0 pk0(pkini,tf);
- PkSpectrum0 pk0nosc2(pkini,tfnosc2);
- PkSpectrum0 pk0nosc1(pkini,tfnosc1);
- PkSpectrum0 pk0nob(pkini,tfnob);
-
- PkSpectrumZ pkz(pk0,d1,zval);
- PkSpectrumZ pkznosc2(pk0nosc2,d1,zval);
- PkSpectrumZ pkznosc1(pk0nosc1,d1,zval);
- PkSpectrumZ pkznob(pk0nob,d1,zval);
+ PkSpecCalc pk0(pkini,tf,d1,0.);
+ PkSpecCalc pk0nosc2(pkini,tfnosc2,d1,0.);
+ PkSpecCalc pk0nosc1(pkini,tfnosc1,d1,0.);
+ PkSpecCalc pk0nob(pkini,tfnob,d1,0.);
+
+ PkSpecCalc pkz(pkini,tf,d1,zval);
+ PkSpecCalc pkznosc2(pkini,tfnosc2,d1,zval);
+ PkSpecCalc pkznosc1(pkini,tfnosc1,d1,zval);
+ PkSpecCalc pkznob(pkini,tfnob,d1,zval);
+
+ PkTabulate* pkcamb = NULL;
+ if(fcmbfile.size()>0) {
+   cout<<endl<<"PkTabulate for CAMB"<<endl;
+   pkcamb = new PkTabulate;
+   pkcamb->ReadCAMB(fcmbfile,h100,0.);
+   pkcamb->SetInterpTyp(1);
+   if(pkcamb->NPoints()==0) {delete pkcamb; pkcamb = NULL;}
+ }
 
  //--------------------------
@@ -118,5 +134,5 @@
 
  //--------------------------
- const int n = 14;
+ const int n = 15;
  const char *vname[n] = {
    "k","pkini",
@@ -124,5 +140,5 @@
    "tfnosc2","pk0nosc2","pknosc2",
    "tfnosc1","pk0nosc1","pknosc1",
-   "tfnob","pk0nob","pknob"
+   "tfnob","pk0nob","pknob","pktab"
  };
  NTuple nt(n,vname);
@@ -145,6 +161,9 @@
    xnt[12] = pk0nob(k);
    xnt[13] = pkznob(k,zval);
+   if(pkcamb != NULL) xnt[14] = (*pkcamb)(k); else xnt[14] = 0;
    nt.Fill(xnt);
  }
+
+   if(pkcamb != NULL) {delete pkcamb; pkcamb = NULL;}
 
  //--------------------------
@@ -190,49 +209,45 @@
 n/plot nt.log10(pkini)%log10(k) pkini>0.&&k>0. ! "nsta connectpoints"
 
-#### Gestion des abscisses
-set xk log10(k)
-set xk k
-
 #### fct de transfert
 zone
-n/plot nt.tf%$xk ! ! "nsta connectpoints"
-n/plot nt.tfnosc2%$xk ! ! "nsta connectpoints same red"
-n/plot nt.tfnosc1%$xk ! ! "nsta connectpoints same blue"
-n/plot nt.tfnob%$xk ! ! "nsta connectpoints same green"
-
-n/plot nt.tf/tfnosc2%$xk ! ! "nsta connectpoints red"
-n/plot nt.tf/tfnosc1%$xk ! ! "nsta connectpoints same blue"
-n/plot nt.tf/tfnob%$xk ! ! "nsta connectpoints same green"
+n/plot nt.tf%k k>0&&tf>0 ! "nsta connectpoints logx logy"
+n/plot nt.tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints same red logx logy"
+n/plot nt.tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx logy"
+n/plot nt.tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx logy"
+
+n/plot nt.tf/tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints red logx"
+n/plot nt.tf/tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx"
+n/plot nt.tf/tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx"
 addline -10 1 10 1
 
 #### Spectre a z=0
 zone
-n/plot nt.pk0%$xk ! ! "nsta connectpoints"
-n/plot nt.pk0nosc2%$xk ! ! "nsta connectpoints same red"
-n/plot nt.pk0nosc1%$xk ! ! "nsta connectpoints same blue"
-n/plot nt.pk0nob%$xk ! ! "nsta connectpoints same green"
+n/plot nt.pk0%k k>0&&pk0>0 ! "nsta connectpoints logx logy"
+n/plot nt.pk0nosc2%k k>0&&pk0nosc2>0 ! "nsta connectpoints same red logx logy"
+n/plot nt.pk0nosc1%k k>0&&pk0nosc1>0 ! "nsta connectpoints same blue logx logy"
+n/plot nt.pk0nob%k k>0&&pk0nob>0 ! "nsta connectpoints same green logx logy"
 
 # Check
 zone 2 2
-n/plot nt.pk0/pkini-tf*tf%$xk pkini>0 ! "nsta crossmarker3"
-n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%$xk pkini>0 ! "nsta crossmarker3"
-n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%$xk pkini>0 ! "nsta crossmarker3"
-n/plot nt.pk0nob/pkini-tfnob*tfnob%$xk pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0/pkini-tf*tf%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
+n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
+n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
+n/plot nt.pk0nob/pkini-tfnob*tfnob%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
 
 #### Spectre a z
 zone
-n/plot nt.pk%$xk ! ! "nsta connectpoints"
-n/plot nt.pknosc2%$xk ! ! "nsta connectpoints same red"
-n/plot nt.pknosc1%$xk ! ! "nsta connectpoints same blue"
-n/plot nt.pknob%$xk ! ! "nsta connectpoints same green"
-
-n/plot nt.pk/pknosc2%$xk ! ! "nsta connectpoints red"
-n/plot nt.pk/pknosc1%$xk ! ! "nsta connectpoints same blue"
-n/plot nt.pk/pknob%$xk ! ! "nsta connectpoints same green"
+n/plot nt.pk%k k>0&&pk>0 ! "nsta connectpoints logx logy"
+n/plot nt.pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx logy"
+n/plot nt.pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx logy"
+n/plot nt.pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx logy"
+
+n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints red logx"
+n/plot nt.pk/pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx"
+n/plot nt.pk/pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx"
 addline -10 1 10 1
 
 #### Le spectre version Delta^2
 set D2 k*k*k*pk/(2*M_PI*M_PI)
-n/plot nt.$D2%$xk  !  ! "nsta crossmarker3 connectpoints"
+n/plot nt.$D2%k  k>0  ! "nsta crossmarker3 connectpoints logx"
 
 #### Test des transferts dans Histo et TVector
@@ -248,3 +263,13 @@
 disp halea "same red"
 
+#### Le spectre tabule
+n/plot nt.pk%k pk>0&&k>0 ! "nsta connectpoints logx logy"
+n/plot nt.pktab%k pktab>0&&k>0 ! "nsta connectpoints same red logx logy"
+
+n/plot nt.(pktab-pk)%k k>0&&pk>0 ! "nsta connectpoints"
+n/plot nt.(pktab-pk)/pk%k k>0&&pk>0 ! "nsta connectpoints"
+
+n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints logx"
+n/plot nt.pktab/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx"
+
 */
Index: /trunk/Cosmo/SimLSS/cmvtvarspec.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtvarspec.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/cmvtvarspec.cc	(revision 3805)
@@ -75,14 +75,12 @@
 
  //-----------------------------------------------------------------
- InitialSpectrum pkini(ns,as);
+ InitialPowerLaw pkini(ns,as);
 
  TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
  //tf.SetNoOscEnv(2);
 
- GrowthFactor d1(om0,ol0);
+ GrowthEisenstein d1(om0,ol0);
 
- PkSpectrum0 pk0(pkini,tf);
-
- PkSpectrumZ pkz(pk0,d1,zref);
+ PkSpecCalc pkz(pkini,tf,d1,zref);
 
  //-----------------------------------------------------------------
Index: /trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3805)
@@ -678,5 +678,5 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier0(PkSpectrumZ& pk_at_z)
+void GeneFluct3D::ComputeFourier0(PkSpectrum& pk_at_z)
 // cf ComputeFourier() mais avec autre methode de realisation du spectre
 //    (attention on fait une fft pour realiser le spectre)
@@ -727,5 +727,5 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier(PkSpectrumZ& pk_at_z)
+void GeneFluct3D::ComputeFourier(PkSpectrum& pk_at_z)
 // Calcule une realisation du spectre "pk_at_z"
 // Attention: dans TArray le premier indice varie le + vite
Index: /trunk/Cosmo/SimLSS/genefluct3d.h
===================================================================
--- /trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3804)
+++ /trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3805)
@@ -115,6 +115,6 @@
   double Href(void) {return h_ref_;}
 
-  void ComputeFourier0(PkSpectrumZ& pk_at_z);
-  void ComputeFourier(PkSpectrumZ& pk_at_z);
+  void ComputeFourier0(PkSpectrum& pk_at_z);
+  void ComputeFourier(PkSpectrum& pk_at_z);
   void FilterByPixel(void);
   void ToRedshiftSpace(void);
Index: /trunk/Cosmo/SimLSS/pkspectrum.cc
===================================================================
--- /trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3804)
+++ /trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3805)
@@ -16,18 +16,18 @@
 
 ///////////////////////////////////////////////////////////
-//******************** 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,4 +298,6 @@
 : kmin_(1.) , kmax_(-1.) , interptyp_(0)
 {
+  k_.resize(0);
+  tf_.resize(0);
 }
 
@@ -331,5 +333,5 @@
  char *line = new char[lenline];
 
- int nread = 0;
+ k_.resize(0); tf_.resize(0);
  double tmax = -1.;
  while ( fgets(line,lenline,file) != NULL ) {
@@ -341,41 +343,80 @@
    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,5 +440,5 @@
 }
 
-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,6 +446,6 @@
 {
  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,5 +472,5 @@
 }
 
-void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
+void GrowthEisenstein::SetParTo(double OmegaMatter0,double OmegaLambda0)
 {
  if(OmegaMatter0>0.) O0_ = OmegaMatter0;
@@ -437,5 +478,5 @@
 }
 
-bool GrowthFactor::SetParTo(double OmegaMatter0)
+bool GrowthEisenstein::SetParTo(double OmegaMatter0)
 // idem precedent sans changer OmegaLambda0
 {
@@ -447,65 +488,156 @@
 
 ///////////////////////////////////////////////////////////
-//************** 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;
 }
-
 
 
Index: /trunk/Cosmo/SimLSS/pkspectrum.h
===================================================================
--- /trunk/Cosmo/SimLSS/pkspectrum.h	(revision 3804)
+++ /trunk/Cosmo/SimLSS/pkspectrum.h	(revision 3805)
@@ -10,7 +10,15 @@
 class InitialSpectrum : public GenericFunc {
 public:
-  InitialSpectrum(double n,double a=1.);
-  InitialSpectrum(InitialSpectrum& pkinf);
-  virtual ~InitialSpectrum(void);
+  InitialSpectrum(void) {};
+  InitialSpectrum(InitialSpectrum& pkinf) {};
+  virtual ~InitialSpectrum(void) {};
+};
+
+//-----------------------------------------------------------------------------------
+class InitialPowerLaw : public InitialSpectrum {
+public:
+  InitialPowerLaw(double n,double a=1.);
+  InitialPowerLaw(InitialPowerLaw& pkinf);
+  virtual ~InitialPowerLaw(void);
   virtual double operator() (double k) {return A_ * pow(k,n_);}
   void SetNorm(double a) {A_ = a;}
@@ -21,5 +29,12 @@
 
 //-----------------------------------------------------------------------------------
-class TransfertEisenstein : public GenericFunc {
+class TransfertFunction : public GenericFunc {
+public:
+  TransfertFunction(void) {};
+  virtual ~TransfertFunction(void) {};
+};
+
+//-----------------------------------------------------------------------------------
+class TransfertEisenstein : public TransfertFunction {
 public:
 
@@ -53,5 +68,5 @@
 
 //-----------------------------------------------------------------------------------
-class TransfertTabulate : public GenericFunc {
+class TransfertTabulate : public TransfertFunction {
 public:
   TransfertTabulate(void);
@@ -62,4 +77,5 @@
   void SetInterpTyp(int typ=0);
   int ReadCMBFast(string filename,double h100,double OmegaCDM0,double OmegaBaryon0);
+  int ReadCAMB(string filename, double h100=0.71);
 protected:
   double kmin_,kmax_;
@@ -72,7 +88,15 @@
 class GrowthFactor : public GenericFunc {
 public:
-  GrowthFactor(double OmegaMatter0,double OmegaLambda0);
-  GrowthFactor(GrowthFactor& d1);
-  virtual ~GrowthFactor(void);
+  GrowthFactor(void) {};
+  virtual ~GrowthFactor(void) {};
+  virtual double DsDz(double z, double);
+};
+
+//-----------------------------------------------------------------------------------
+class GrowthEisenstein : public GrowthFactor {
+public:
+  GrowthEisenstein(double OmegaMatter0,double OmegaLambda0);
+  GrowthEisenstein(GrowthEisenstein& d1);
+  virtual ~GrowthEisenstein(void);
   virtual double operator() (double z);
   virtual double DsDz(double z,double dzinc=0.01);
@@ -85,36 +109,20 @@
 
 //-----------------------------------------------------------------------------------
-class PkSpectrum0 : public GenericFunc {
-public:
-  PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf);
-  PkSpectrum0(PkSpectrum0& pk0);
-  virtual ~PkSpectrum0(void);
-  virtual double operator() (double z);
-  InitialSpectrum& GetPkIni(void) {return pkinf_;}
-  TransfertEisenstein& GetTransfert(void) {return tf_;}
-protected:
-  InitialSpectrum& pkinf_;
-  TransfertEisenstein& tf_;
-};
-
-//-----------------------------------------------------------------------------------
-class PkSpectrumZ : public GenericFunc {
-public:
+class PkSpectrum : public GenericFunc {
+public:
+  // typsec = PK : compute Pk(k)
+  //        = DELTA : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
   typedef enum {PK=0, DELTA=1} ReturnSpectrum;
-  PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref=0.);
-  PkSpectrumZ(PkSpectrumZ& pkz);
-  virtual ~PkSpectrumZ(void);
-  virtual double operator() (double k);
-  virtual double operator() (double k,double z);
-  void   SetZ(double z) {zref_ = z;}
-  double GetZ(void) {return zref_;}
-  void SetTypSpec(ReturnSpectrum typspec=PK);
-  void SetScale(double scale=1.) {scale_=scale;}
-  double GetScale(void) {return scale_;}
-  PkSpectrum0& GetPk0(void) {return pk0_;}
-  GrowthFactor& GetGrowthFactor(void) {return d1_;}
-protected:
-  PkSpectrum0& pk0_;
-  GrowthFactor& d1_;
+
+  PkSpectrum(void);
+  PkSpectrum(PkSpectrum& pk);
+  virtual ~PkSpectrum(void) {};
+  virtual void   SetZ(double z) {zref_ = z;}
+  virtual double GetZ(void) {return zref_;}
+  virtual void SetScale(double scale=1.) {scale_ = scale;}
+  virtual double GetScale(void) {return scale_;}
+  virtual void SetTypSpec(ReturnSpectrum typspec=PK) {typspec_ = typspec;}
+  virtual ReturnSpectrum GetTypSpec(void) {return typspec_ ;}
+protected:
   double zref_, scale_;
   ReturnSpectrum typspec_;
@@ -122,4 +130,55 @@
 
 //-----------------------------------------------------------------------------------
+class PkSpecCalc : public PkSpectrum {
+public:
+  PkSpecCalc(InitialSpectrum& pkinf,TransfertFunction& tf,GrowthFactor& d1,double zref=0.);
+  PkSpecCalc(PkSpecCalc& pkz);
+  virtual ~PkSpecCalc(void);
+  virtual double operator() (double k);
+  virtual double operator() (double k,double z);
+  InitialSpectrum& GetPkIni(void) {return pkinf_;}
+  TransfertFunction& GetTransfert(void) {return tf_;}
+  GrowthFactor& GetGrowthFactor(void) {return d1_;}
+protected:
+  InitialSpectrum& pkinf_;
+  TransfertFunction& tf_;
+  GrowthFactor& d1_;
+};
+
+//-----------------------------------------------------------------------------------
+class PkTabulate : public PkSpectrum {
+public:
+  PkTabulate(void);
+  PkTabulate(PkTabulate& pkz);
+  virtual ~PkTabulate(void);
+  virtual double operator() (double k);
+  virtual double operator() (double k,double z);
+  int NPoints(void) {return k_.size();}
+  void SetInterpTyp(int typ=0);
+  int ReadCAMB(string filename, double h100=0.71, double zreftab = 0.);
+protected:
+  double kmin_,kmax_;
+  int interptyp_;
+  vector<double> k_, pk_;
+};
+
+//-----------------------------------------------------------------------------------
+class PkEisenstein : public PkSpectrum {
+public:
+  PkEisenstein(InitialPowerLaw& pkinf,TransfertEisenstein& tf,GrowthEisenstein& d1,double zref=0.);
+  PkEisenstein(PkEisenstein& pkz);
+  virtual ~PkEisenstein(void);
+  virtual double operator() (double k);
+  virtual double operator() (double k,double z);
+  InitialPowerLaw& GetPkIni(void) {return pkinf_;}
+  TransfertEisenstein& GetTransfert(void) {return tf_;}
+  GrowthEisenstein& GetGrowthFactor(void) {return d1_;}
+protected:
+  InitialPowerLaw& pkinf_;
+  TransfertEisenstein& tf_;
+  GrowthEisenstein& d1_;
+};
+
+//-----------------------------------------------------------------------------------
 class VarianceSpectrum : public GenericFunc {
 public:
@@ -128,5 +187,5 @@
 
   VarianceSpectrum(GenericFunc& pk,double R,TypeFilter typfilter);
-  VarianceSpectrum(VarianceSpectrum& pkinf);
+  VarianceSpectrum(VarianceSpectrum& vpk);
   virtual ~VarianceSpectrum(void);
 
