Changeset 3805 in Sophya for trunk/Cosmo/SimLSS/cmvtstpk.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/cmvtstpk.cc (modified) (10 diffs)

trunk/Cosmo/SimLSS/cmvtstpk.cc

@@ -19 +19 @@
 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 +38 @@
  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 +66 @@
     wrascii = true;
     break;
+  case 'F' :
+    fcmbfile = optarg;
+    break;
   case 'h' :
   default :
@@ -78 +84 @@
  cout<<"scale="<<scale<<endl;
  cout<<"zval="<<zval<<endl;
-
- //--------------------------
- InitialSpectrum pkini(ns,as);
+ cout<<"fcmbfile="<<fcmbfile<<endl;
+
+ //--------------------------
+ InitialPowerLaw pkini(ns,as);
  pkini.SetNorm(scale);
 
@@ -89 +96 @@
  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 +134 @@
 
  //--------------------------
- const int n = 14;
+ const int n = 15;
  const char *vname[n] = {
    "k","pkini",
@@ -124 +140 @@
    "tfnosc2","pk0nosc2","pknosc2",
    "tfnosc1","pk0nosc1","pknosc1",
-   "tfnob","pk0nob","pknob"
+   "tfnob","pk0nob","pknob","pktab"
  };
  NTuple nt(n,vname);
@@ -145 +161 @@
    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 +209 @@
 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 +263 @@
 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"
+
 */