Changeset 3805 in Sophya for trunk/Cosmo


Ignore:
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

Location:
trunk/Cosmo/SimLSS
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvdefsurv.cc

    r3749 r3805  
    183183 univ.SetDynParam(h100,om0,or0,ol0,w0);
    184184
    185  GrowthFactor growth(om0,ol0);
     185 GrowthEisenstein growth(om0,ol0);
    186186
    187187 // ---

  • trunk/Cosmo/SimLSS/cmvfitpk.cc

    r3572 r3805  
    2626class MyFCN : public FCNBase {
    2727public:
    28   MyFCN(HistoErr& herr,PkSpectrumZ& pk);
     28  MyFCN(HistoErr& herr,PkEisenstein& pk);
    2929  void Init(const vector<double>& par) const;
    3030  double Func(double x,const vector<double>& par) const;
    … …  
    3333private:
    3434  HistoErr& herr_;
    35   PkSpectrumZ& pk_;
     35  PkEisenstein& pk_;
    3636};
    3737} }  // namespace ROOT + Minuit2
    … …  
    8484 // --- Create spectrum
    8585 cout<<endl<<"\n--- Create Spectrum"<<endl;
    86  InitialSpectrum pkini(ns,as);
     86 InitialPowerLaw pkini(ns,as);
    8787 bool nobaryon = false;
    8888 double ocdm0 = om0-ob0;
    8989 TransfertEisenstein tf(h100,ocdm0,ob0,T_CMB_Par,nobaryon);
    90  GrowthFactor growth(om0,ol0); cout<<"...Growth="<<growth(zref)<<endl;
    91  PkSpectrum0 pk0(pkini,tf);
    92  PkSpectrumZ pkz(pk0,growth,zref);
     90 GrowthEisenstein growth(om0,ol0); cout<<"...Growth="<<growth(zref)<<endl;
     91 PkEisenstein pkz(pkini,tf,growth,zref);
    9392
    9493 // --- Compute variance and normalize spectrum
    … …  
    256255
    257256//--------------------------------------------------
    258 MyFCN::MyFCN(HistoErr& herr,PkSpectrumZ& pk)
     257MyFCN::MyFCN(HistoErr& herr,PkEisenstein& pk)
    259258  : herr_(herr) , pk_(pk)
    260259{
    … …  
    264263{
    265264 double A=par[1], Ocdm0=par[2], Ob0=par[3], Ol0=par[4], h100=par[5], ns=par[6];
    266  pk_.GetPk0().GetPkIni().SetSlope(ns);
    267  pk_.GetPk0().GetPkIni().SetNorm(A);
    268  pk_.GetPk0().GetTransfert().SetParTo(h100,Ocdm0,Ob0);
     265 pk_.GetPkIni().SetSlope(ns);
     266 pk_.GetPkIni().SetNorm(A);
     267 pk_.GetTransfert().SetParTo(h100,Ocdm0,Ob0);
    269268 pk_.GetGrowthFactor().SetParTo(Ocdm0+Ob0,Ol0);
    270269}

  • trunk/Cosmo/SimLSS/cmvginit3d.cc

    r3800 r3805  
    218218 cout<<endl<<"\n--- Create Spectrum"<<endl;
    219219
    220  InitialSpectrum pkini(ns,as);
     220 InitialPowerLaw pkini(ns,as);
    221221
    222222 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
    … …  
    224224   tfnos.SetNoOscEnv(2);
    225225
    226  GrowthFactor growth(om0,ol0);
    227  // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
     226 GrowthEisenstein growth(om0,ol0);
     227 // GrowthEisenstein growth(1.,0.); // D(z) = 1/(1+z)
    228228 double growth_at_z = growth(zref);
    229229 cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
    230230
    231  PkSpectrum0 pk0(pkini,tf);
    232  PkSpectrum0 pk0nos(pkini,tfnos);
    233 
    234  PkSpectrumZ pkz(pk0,growth,zref);
    235  PkSpectrumZ pkznos(pk0nos,growth,zref);
     231 PkSpecCalc pkz(pkini,tf,growth,zref);
     232 PkSpecCalc pkznos(pkini,tfnos,growth,zref);
    236233
    237234 //-----------------------------------------------------------------

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3768 r3805  
    299299 cout<<endl<<"\n--- Create Spectrum"<<endl;
    300300
    301  InitialSpectrum pkini(ns,as);
     301 InitialPowerLaw pkini(ns,as);
    302302
    303303 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
    304304 //tf.SetNoOscEnv(2);
    305305
    306  GrowthFactor growth(om0,ol0);
    307  // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
     306 GrowthEisenstein growth(om0,ol0);
     307 // GrowthEisenstein growth(1.,0.); // D(z) = 1/(1+z)
    308308 double growth_at_z = growth(zref);
    309309 cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
    310310
    311  PkSpectrum0 pk0(pkini,tf);
    312 
    313  PkSpectrumZ pkz(pk0,growth,zref);
     311 PkSpecCalc pkz(pkini,tf,growth,zref);
    314312
    315313 //-----------------------------------------------------------------

  • trunk/Cosmo/SimLSS/cmvtgrowth.cc

    r3768 r3805  
    3030 cout<<"Om0="<<om0<<"  Ol0="<<ol0<<endl;
    3131
    32  GrowthFactor growth(om0,ol0);
     32 GrowthEisenstein growth(om0,ol0);
    3333 cout<<"D1(z=0) = "<<growth(0.)<<endl;
    3434

  • trunk/Cosmo/SimLSS/cmvtintfun.cc

    r3196 r3805  
    4444  unsigned short glorder=4;
    4545
    46   if(narg>1) sscanf(arg[1],"%lf",&n);
    47   if(narg>2) sscanf(arg[2],"%lf,%lf,%d",&xmin,&xmax,&npt);
    4846  cout<<"n="<<n<<" xmin="<<xmin<<" xmax="<<xmax<<" npt="<<npt<<endl;
    4947

  • trunk/Cosmo/SimLSS/cmvtransf.cc

    r3802 r3805  
    1717cout
    1818<<"cmvtransf [options]"<<endl
    19 <<"  -k k1,k2,npt : k range in Mpc^-1"<<endl
     19<<"  -k npt,k1,k2 : k range in Mpc^-1"<<endl
    2020<<"  -U h100,Om0,Ob0,tcmb : cosmology"<<endl
    21 <<"  -F filename : read also CMBFast file"<<endl
     21<<"  -F filename : read also CAMB file"<<endl
    2222<<endl;
    2323}
    … …  
    2727 double k1 = 1e-6, k2 = 10.; int_4 npt = 100;
    2828 double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;
    29  string fcmbfast = "";
     29 string fcmbfile = "";
    3030
    3131 char c;
    … …  
    3939    break;
    4040  case 'F' :
    41     fcmbfast = optarg;
     41    fcmbfile = optarg;
    4242    break;
    4343  case 'h' :
    … …  
    5252 cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
    5353 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
    54  cout<<"fcmbfast="<<fcmbfast<<endl;
     54 cout<<"fcmbfile="<<fcmbfile<<endl;
    5555
    5656 cout<<"TransfertEisenstein with baryon"<<endl;
    … …  
    6969 TransfertEisenstein tfnob(h100,Om0-Ob0,Ob0,tcmb,true);  // No baryons
    7070
    71  TransfertTabulate tfcmbfast;
    72  bool cmbfastOK = false;
    73  if(fcmbfast.size()>0) {
    74    cout<<endl<<"TransfertTabulate for CMBfast"<<endl;
    75    tfcmbfast.ReadCMBFast(fcmbfast,h100,Om0-Ob0,Ob0);
    76    tfcmbfast.SetInterpTyp(1);
    77    if(tfcmbfast.NPoints()>0) cmbfastOK = true;
     71 TransfertTabulate tf_cmbfile;
     72 bool cmbfileOK = false;
     73 if(fcmbfile.size()>0) {
     74   cout<<endl<<"TransfertTabulate for CAMB"<<endl;
     75   tf_cmbfile.ReadCAMB(fcmbfile,h100);
     76   tf_cmbfile.SetInterpTyp(1);
     77   if(tf_cmbfile.NPoints()>0) cmbfileOK = true;
    7878 }
    7979
    … …  
    9292   xnt[3]  = tfnosc2(k);
    9393   xnt[4]  = tfnob(k);
    94    if(cmbfastOK) xnt[5]  = tfcmbfast(k);
     94   if(cmbfileOK) xnt[5]  = tf_cmbfile(k);
    9595   nt.Fill(xnt);
    9696 }
    … …  
    107107openppf cmvtransf.ppf
    108108
    109 set k k
    110 set k log10(k)
     109# Eisenstein
     110n/plot nt.t%k t>0&&k>0 ! "nsta connectpoints logx logy"
     111n/plot nt.tnob%k tnob>0&&k>0 ! "nsta red same connectpoints logx logy"
     112n/plot nt.tnosc2%k tnosc2>0&&k>0 ! "nsta blue same connectpoints logx logy"
     113n/plot nt.tnosc1%k tnosc1>0&&k>0 ! "nsta green same connectpoints logx logy"
    111114
    112 # Eisenstein
    113 n/plot nt.t%$k ! ! "nsta connectpoints"
    114 n/plot nt.tnob%$k ! ! "nsta red same connectpoints"
    115 n/plot nt.tnosc2%$k ! ! "nsta blue same connectpoints"
    116 n/plot nt.tnosc1%$k ! ! "nsta green same connectpoints"
     115n/plot nt.t/tnob%k   tnob>0   ! "nsta red connectpoints"
     116n/plot nt.t/tnosc2%k tnosc2>0 ! "nsta blue same connectpoints"
     117n/plot nt.t/tnosc1%k tnosc1>0 ! "nsta green same connectpoints"
    117118
    118 n/plot nt.t/tnob%$k   tnob>0   ! "nsta red connectpoints"
    119 n/plot nt.t/tnosc2%$k tnosc2>0 ! "nsta blue same connectpoints"
    120 n/plot nt.t/tnosc1%$k tnosc1>0 ! "nsta green same connectpoints"
     119# CAMB
     120n/plot nt.t%k t>0&&k>0 ! "nsta connectpoints logx logy"
     121n/plot nt.tcmbf%k tcmbf>0&&k>0 ! "nsta red same connectpoints logx logy"
    121122
    122 # CMBFast
    123 n/plot nt.t%$k ! ! "nsta connectpoints"
    124 n/plot nt.tcmbf%$k ! ! "nsta red same connectpoints"
     123n/plot nt.(tcmbf-t)%k tcmbf>0. ! "nsta connectpoints"
     124n/plot nt.(tcmbf-t)/tcmbf%k tcmbf>0. ! "nsta connectpoints"
    125125
    126 n/plot nt.(tcmbf-t)%$k tcmbf>0. ! "nsta connectpoints"
    127 n/plot nt.(tcmbf-t)/tcmbf%$k tcmbf>0. ! "nsta connectpoints"
     126n/plot nt.t/tnosc2%k tnosc2>0&&t>0&&k>0 ! "nsta connectpoints logx logy"
     127n/plot nt.tcmbf/tnosc2%k tnosc2>0&&tcmbf>0&&k>0 ! "nsta connectpoints same red logx logy"
    128128*/

  • trunk/Cosmo/SimLSS/cmvtstpk.cc

    r3679 r3805  
    1919void usage(void);
    2020void usage(void) {
    21  cout<<"cmvtstpk [options] z_redshift"<<endl
    22      <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
    23      <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
    24      <<" -s scale : on normalise le spectre avec scale"<<endl
    25      <<" -w : write spectra on ASCII file"<<endl;
     21cout
     22<<"cmvtstpk [options] z_redshift"<<endl
     23<<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
     24<<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
     25<<" -s scale : on normalise le spectre avec scale"<<endl
     26<<" -w : write spectra on ASCII file"<<endl
     27<<"  -F filename : read also CAMB file"<<endl;
    2628}
    2729
    … …  
    3638 double scale = 2.54499e+07;  // normalisation du spectre a z=0 selon SDSS
    3739 bool wrascii = false;
     40 string fcmbfile = "";
    3841
    3942 char c;
    40   while((c = getopt(narg,arg,"hwk:s:H:M:B:L:")) != -1) {
     43  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:F:")) != -1) {
    4144  switch (c) {
    4245  case 'H' :
    … …  
    6366    wrascii = true;
    6467    break;
     68  case 'F' :
     69    fcmbfile = optarg;
     70    break;
    6571  case 'h' :
    6672  default :
    … …  
    7884 cout<<"scale="<<scale<<endl;
    7985 cout<<"zval="<<zval<<endl;
    80 
    81  //--------------------------
    82  InitialSpectrum pkini(ns,as);
     86 cout<<"fcmbfile="<<fcmbfile<<endl;
     87
     88 //--------------------------
     89 InitialPowerLaw pkini(ns,as);
    8390 pkini.SetNorm(scale);
    8491
    … …  
    8996 TransfertEisenstein tfnob(h100,Om0,0.,T_CMB_Par,true);
    9097
    91  GrowthFactor d1(Om0,Ol0);
     98 GrowthEisenstein d1(Om0,Ol0);
    9299 cout<<"GrowthFactor: "<<d1(zval)<<endl;
    93100
    94  PkSpectrum0 pk0(pkini,tf);
    95  PkSpectrum0 pk0nosc2(pkini,tfnosc2);
    96  PkSpectrum0 pk0nosc1(pkini,tfnosc1);
    97  PkSpectrum0 pk0nob(pkini,tfnob);
    98 
    99  PkSpectrumZ pkz(pk0,d1,zval);
    100  PkSpectrumZ pkznosc2(pk0nosc2,d1,zval);
    101  PkSpectrumZ pkznosc1(pk0nosc1,d1,zval);
    102  PkSpectrumZ pkznob(pk0nob,d1,zval);
     101 PkSpecCalc pk0(pkini,tf,d1,0.);
     102 PkSpecCalc pk0nosc2(pkini,tfnosc2,d1,0.);
     103 PkSpecCalc pk0nosc1(pkini,tfnosc1,d1,0.);
     104 PkSpecCalc pk0nob(pkini,tfnob,d1,0.);
     105
     106 PkSpecCalc pkz(pkini,tf,d1,zval);
     107 PkSpecCalc pkznosc2(pkini,tfnosc2,d1,zval);
     108 PkSpecCalc pkznosc1(pkini,tfnosc1,d1,zval);
     109 PkSpecCalc pkznob(pkini,tfnob,d1,zval);
     110
     111 PkTabulate* pkcamb = NULL;
     112 if(fcmbfile.size()>0) {
     113   cout<<endl<<"PkTabulate for CAMB"<<endl;
     114   pkcamb = new PkTabulate;
     115   pkcamb->ReadCAMB(fcmbfile,h100,0.);
     116   pkcamb->SetInterpTyp(1);
     117   if(pkcamb->NPoints()==0) {delete pkcamb; pkcamb = NULL;}
     118 }
    103119
    104120 //--------------------------
    … …  
    118134
    119135 //--------------------------
    120  const int n = 14;
     136 const int n = 15;
    121137 const char *vname[n] = {
    122138   "k","pkini",
    … …  
    124140   "tfnosc2","pk0nosc2","pknosc2",
    125141   "tfnosc1","pk0nosc1","pknosc1",
    126    "tfnob","pk0nob","pknob"
     142   "tfnob","pk0nob","pknob","pktab"
    127143 };
    128144 NTuple nt(n,vname);
    … …  
    145161   xnt[12] = pk0nob(k);
    146162   xnt[13] = pkznob(k,zval);
     163   if(pkcamb != NULL) xnt[14] = (*pkcamb)(k); else xnt[14] = 0;
    147164   nt.Fill(xnt);
    148165 }
     166
     167   if(pkcamb != NULL) {delete pkcamb; pkcamb = NULL;}
    149168
    150169 //--------------------------
    … …  
    190209n/plot nt.log10(pkini)%log10(k) pkini>0.&&k>0. ! "nsta connectpoints"
    191210
    192 #### Gestion des abscisses
    193 set xk log10(k)
    194 set xk k
    195 
    196211#### fct de transfert
    197212zone
    198 n/plot nt.tf%$xk ! ! "nsta connectpoints"
    199 n/plot nt.tfnosc2%$xk ! ! "nsta connectpoints same red"
    200 n/plot nt.tfnosc1%$xk ! ! "nsta connectpoints same blue"
    201 n/plot nt.tfnob%$xk ! ! "nsta connectpoints same green"
    202 
    203 n/plot nt.tf/tfnosc2%$xk ! ! "nsta connectpoints red"
    204 n/plot nt.tf/tfnosc1%$xk ! ! "nsta connectpoints same blue"
    205 n/plot nt.tf/tfnob%$xk ! ! "nsta connectpoints same green"
     213n/plot nt.tf%k k>0&&tf>0 ! "nsta connectpoints logx logy"
     214n/plot nt.tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints same red logx logy"
     215n/plot nt.tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx logy"
     216n/plot nt.tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx logy"
     217
     218n/plot nt.tf/tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints red logx"
     219n/plot nt.tf/tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx"
     220n/plot nt.tf/tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx"
    206221addline -10 1 10 1
    207222
    208223#### Spectre a z=0
    209224zone
    210 n/plot nt.pk0%$xk ! ! "nsta connectpoints"
    211 n/plot nt.pk0nosc2%$xk ! ! "nsta connectpoints same red"
    212 n/plot nt.pk0nosc1%$xk ! ! "nsta connectpoints same blue"
    213 n/plot nt.pk0nob%$xk ! ! "nsta connectpoints same green"
     225n/plot nt.pk0%k k>0&&pk0>0 ! "nsta connectpoints logx logy"
     226n/plot nt.pk0nosc2%k k>0&&pk0nosc2>0 ! "nsta connectpoints same red logx logy"
     227n/plot nt.pk0nosc1%k k>0&&pk0nosc1>0 ! "nsta connectpoints same blue logx logy"
     228n/plot nt.pk0nob%k k>0&&pk0nob>0 ! "nsta connectpoints same green logx logy"
    214229
    215230# Check
    216231zone 2 2
    217 n/plot nt.pk0/pkini-tf*tf%$xk pkini>0 ! "nsta crossmarker3"
    218 n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%$xk pkini>0 ! "nsta crossmarker3"
    219 n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%$xk pkini>0 ! "nsta crossmarker3"
    220 n/plot nt.pk0nob/pkini-tfnob*tfnob%$xk pkini>0 ! "nsta crossmarker3"
     232n/plot nt.pk0/pkini-tf*tf%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
     233n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
     234n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
     235n/plot nt.pk0nob/pkini-tfnob*tfnob%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
    221236
    222237#### Spectre a z
    223238zone
    224 n/plot nt.pk%$xk ! ! "nsta connectpoints"
    225 n/plot nt.pknosc2%$xk ! ! "nsta connectpoints same red"
    226 n/plot nt.pknosc1%$xk ! ! "nsta connectpoints same blue"
    227 n/plot nt.pknob%$xk ! ! "nsta connectpoints same green"
    228 
    229 n/plot nt.pk/pknosc2%$xk ! ! "nsta connectpoints red"
    230 n/plot nt.pk/pknosc1%$xk ! ! "nsta connectpoints same blue"
    231 n/plot nt.pk/pknob%$xk ! ! "nsta connectpoints same green"
     239n/plot nt.pk%k k>0&&pk>0 ! "nsta connectpoints logx logy"
     240n/plot nt.pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx logy"
     241n/plot nt.pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx logy"
     242n/plot nt.pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx logy"
     243
     244n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints red logx"
     245n/plot nt.pk/pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx"
     246n/plot nt.pk/pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx"
    232247addline -10 1 10 1
    233248
    234249#### Le spectre version Delta^2
    235250set D2 k*k*k*pk/(2*M_PI*M_PI)
    236 n/plot nt.$D2%$xk  !  ! "nsta crossmarker3 connectpoints"
     251n/plot nt.$D2%k  k>0  ! "nsta crossmarker3 connectpoints logx"
    237252
    238253#### Test des transferts dans Histo et TVector
    … …  
    248263disp halea "same red"
    249264
     265#### Le spectre tabule
     266n/plot nt.pk%k pk>0&&k>0 ! "nsta connectpoints logx logy"
     267n/plot nt.pktab%k pktab>0&&k>0 ! "nsta connectpoints same red logx logy"
     268
     269n/plot nt.(pktab-pk)%k k>0&&pk>0 ! "nsta connectpoints"
     270n/plot nt.(pktab-pk)/pk%k k>0&&pk>0 ! "nsta connectpoints"
     271
     272n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints logx"
     273n/plot nt.pktab/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx"
     274
    250275*/

  • trunk/Cosmo/SimLSS/cmvtvarspec.cc

    r3572 r3805  
    7575
    7676 //-----------------------------------------------------------------
    77  InitialSpectrum pkini(ns,as);
     77 InitialPowerLaw pkini(ns,as);
    7878
    7979 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
    8080 //tf.SetNoOscEnv(2);
    8181
    82  GrowthFactor d1(om0,ol0);
     82 GrowthEisenstein d1(om0,ol0);
    8383
    84  PkSpectrum0 pk0(pkini,tf);
    85 
    86  PkSpectrumZ pkz(pk0,d1,zref);
     84 PkSpecCalc pkz(pkini,tf,d1,zref);
    8785
    8886 //-----------------------------------------------------------------

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3800 r3805  
    678678
    679679//-------------------------------------------------------
    680 void GeneFluct3D::ComputeFourier0(PkSpectrumZ& pk_at_z)
     680void GeneFluct3D::ComputeFourier0(PkSpectrum& pk_at_z)
    681681// cf ComputeFourier() mais avec autre methode de realisation du spectre
    682682//    (attention on fait une fft pour realiser le spectre)
    … …  
    727727
    728728//-------------------------------------------------------
    729 void GeneFluct3D::ComputeFourier(PkSpectrumZ& pk_at_z)
     729void GeneFluct3D::ComputeFourier(PkSpectrum& pk_at_z)
    730730// Calcule une realisation du spectre "pk_at_z"
    731731// Attention: dans TArray le premier indice varie le + vite

  • trunk/Cosmo/SimLSS/genefluct3d.h

    r3800 r3805  
    115115  double Href(void) {return h_ref_;}
    116116
    117   void ComputeFourier0(PkSpectrumZ& pk_at_z);
    118   void ComputeFourier(PkSpectrumZ& pk_at_z);
     117  void ComputeFourier0(PkSpectrum& pk_at_z);
     118  void ComputeFourier(PkSpectrum& pk_at_z);
    119119  void FilterByPixel(void);
    120120  void ToRedshiftSpace(void);

  • trunk/Cosmo/SimLSS/pkspectrum.cc

    r3802 r3805  
    1616
    1717///////////////////////////////////////////////////////////
    18 //******************** InitialSpectrum ******************//
    19 ///////////////////////////////////////////////////////////
    20 
    21 InitialSpectrum::InitialSpectrum(double n,double a)
     18//******************** InitialPowerLaw ******************//
     19///////////////////////////////////////////////////////////
     20
     21InitialPowerLaw::InitialPowerLaw(double n,double a)
    2222  : n_(n), A_(a)
    2323{
    2424}
    2525
    26 InitialSpectrum::InitialSpectrum(InitialSpectrum& pkinf)
     26InitialPowerLaw::InitialPowerLaw(InitialPowerLaw& pkinf)
    2727  : n_(pkinf.n_), A_(pkinf.A_)
    2828{
    2929}
    3030
    31 InitialSpectrum::~InitialSpectrum(void)
     31InitialPowerLaw::~InitialPowerLaw(void)
    3232{
    3333}
    … …  
    298298: kmin_(1.) , kmax_(-1.) , interptyp_(0)
    299299{
     300  k_.resize(0);
     301  tf_.resize(0);
    300302}
    301303
    … …  
    331333 char *line = new char[lenline];
    332334
    333  int nread = 0;
     335 k_.resize(0); tf_.resize(0);
    334336 double tmax = -1.;
    335337 while ( fgets(line,lenline,file) != NULL ) {
    … …  
    341343   k_.push_back(k);
    342344   tf_.push_back(tf);
    343    nread++;
    344345 }
    345346
    346  cout<<"TransfertTabulate::ReadCMBFast: nread="<<nread<<" tf_max="<<tmax<<endl;
     347 cout<<"TransfertTabulate::ReadCMBFast: nread="<<tf_.size()<<" tf_max="<<tmax<<endl;
    347348 delete [] line;
    348  if(nread==0) return nread;
     349 if(tf_.size()==0) return (int)tf_.size();
    349350
    350351 for(unsigned int i=0;i<tf_.size();i++) tf_[i] /= tmax;
    351352
    352  return nread;
    353 }
     353 return (int)tf_.size();
     354}
     355
     356int TransfertTabulate::ReadCAMB(string filename, double h100)
     357{
     358 FILE *file = fopen(filename.c_str(),"r");
     359 if(file==NULL) return -1;
     360 cout<<"TransfertTabulate::ReadCAMB: fn="<<filename<<" h100="<<h100<<endl;
     361
     362 const int lenline = 512;
     363 char *line = new char[lenline];
     364
     365 k_.resize(0); tf_.resize(0);
     366 double tmax = -1.;
     367 while ( fgets(line,lenline,file) != NULL ) {
     368   double k,tcdm,tbar,tph,trel,tnu,ttot, tf;
     369   sscanf(line,"%lf %lf %lf %lf %lf %lf %lf",&k,&tcdm,&tbar,&tph,&trel,&tnu,&ttot);
     370   k *= h100;     // convert h Mpc^-1  ->  Mpc^-1
     371   tf = ttot;
     372   if(tf>tmax) tmax = tf;
     373   k_.push_back(k);
     374   tf_.push_back(tf);
     375 }
     376
     377 cout<<"TransfertTabulate::ReadCAMB nread="<<tf_.size()<<" tf_max="<<tmax<<endl;
     378 delete [] line;
     379 if(tf_.size()==0) return (int)tf_.size();
     380
     381 for(unsigned int i=0;i<tf_.size();i++) tf_[i] /= tmax;
     382
     383 return (int)tf_.size();
     384}
     385
    354386
    355387///////////////////////////////////////////////////////////
    356388//********************* GrowthFactor ********************//
     389///////////////////////////////////////////////////////////
     390double GrowthFactor::DsDz(double z, double)
     391{
     392  cout<<"Error: GrowthFactor::DsDz not implemented"<<endl;
     393  throw AllocationError("Error: GrowthFactor::DsDz not implemented");
     394}
     395
     396///////////////////////////////////////////////////////////
     397//********************* GrowthEisenstein ********************//
    357398///////////////////////////////////////////////////////////
    358399
    359400// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
    360401// Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
    361 GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
     402GrowthEisenstein::GrowthEisenstein(double OmegaMatter0,double OmegaLambda0)
    362403  : O0_(OmegaMatter0) , Ol_(OmegaLambda0)
    363404{
    364405  if(OmegaMatter0==0.) {
    365     cout<<"GrowthFactor::GrowthFactor: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
    366     throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
    367   }
    368 }
    369 
    370 GrowthFactor::GrowthFactor(GrowthFactor& d1)
     406    cout<<"GrowthEisenstein::GrowthEisenstein: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
     407    throw ParmError("GrowthEisenstein::GrowthEisenstein: Error badOmegaMatter0  value");
     408  }
     409}
     410
     411GrowthEisenstein::GrowthEisenstein(GrowthEisenstein& d1)
    371412  : O0_(d1.O0_) , Ol_(d1.Ol_)
    372413{
    373414}
    374415
    375 GrowthFactor::~GrowthFactor(void)
    376 {
    377 }
    378 
    379 double GrowthFactor::operator() (double z)
     416GrowthEisenstein::~GrowthEisenstein(void)
     417{
     418}
     419
     420double GrowthEisenstein::operator() (double z)
    380421// see Formulae A4 + A5 + A6 page 614
    381422{
    … …  
    399440}
    400441
    401 double GrowthFactor::DsDz(double z,double dzinc)
     442double GrowthEisenstein::DsDz(double z,double dzinc)
    402443// y-y0 = a*(x-x0)^2 + b*(x-x0)
    403444// dy = a*dx^2 + b*dx   avec  dx = x-x0
    … …  
    405446{
    406447 if(z<0. || dzinc<=0.) {
    407    cout<<"GrowthFactor::DsDz: z<0 or dzinc<=0. !"<<endl;
    408    throw ParmError("GrowthFactor::DsDz: z<0 or dzinc<=0. !");
     448   cout<<"GrowthEisenstein::DsDz: z<0 or dzinc<=0. !"<<endl;
     449   throw ParmError("GrowthEisenstein::DsDz: z<0 or dzinc<=0. !");
    409450 }
    410451
    … …  
    431472}
    432473
    433 void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
     474void GrowthEisenstein::SetParTo(double OmegaMatter0,double OmegaLambda0)
    434475{
    435476 if(OmegaMatter0>0.) O0_ = OmegaMatter0;
    … …  
    437478}
    438479
    439 bool GrowthFactor::SetParTo(double OmegaMatter0)
     480bool GrowthEisenstein::SetParTo(double OmegaMatter0)
    440481// idem precedent sans changer OmegaLambda0
    441482{
    … …  
    447488
    448489///////////////////////////////////////////////////////////
    449 //************** PkSpectrum0 et PkSpectrumZ *************//
    450 ///////////////////////////////////////////////////////////
    451 
    452 PkSpectrum0::PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf)
    453   : pkinf_(pkinf) , tf_(tf)
    454 {
    455 }
    456 
    457 PkSpectrum0::PkSpectrum0(PkSpectrum0& pk0)
    458   : pkinf_(pk0.pkinf_) , tf_(pk0.tf_)
    459 {
    460 }
    461 
    462 PkSpectrum0::~PkSpectrum0(void)
    463 {
    464 }
    465 
    466 double PkSpectrum0::operator() (double k)
     490//********************** PkSpectrum *********************//
     491///////////////////////////////////////////////////////////
     492
     493PkSpectrum::PkSpectrum(void)
     494  : zref_(0.) , scale_(1.) , typspec_(PK)
     495{
     496}
     497
     498PkSpectrum::PkSpectrum(PkSpectrum& pk)
     499  : zref_(pk.zref_) , scale_(pk.scale_) , typspec_(pk.typspec_)
     500{
     501}
     502
     503
     504///////////////////////////////////////////////////////////
     505//********************** PkSpecCalc *********************//
     506///////////////////////////////////////////////////////////
     507
     508PkSpecCalc::PkSpecCalc(InitialSpectrum& pkinf,TransfertFunction& tf,GrowthFactor& d1,double zref)
     509  : pkinf_(pkinf) , tf_(tf) , d1_(d1)
     510{
     511  zref_ = zref;
     512}
     513
     514PkSpecCalc::PkSpecCalc(PkSpecCalc& pkz)
     515  : pkinf_(pkz.pkinf_) , tf_(pkz.tf_) , d1_(pkz.d1_)
     516{
     517}
     518
     519PkSpecCalc::~PkSpecCalc(void)
     520{
     521}
     522
     523double PkSpecCalc::operator() (double k)
     524{
     525  return (*this)(k,zref_);
     526}
     527
     528double PkSpecCalc::operator() (double k,double z)
    467529{
    468530  double tf = tf_(k);
    469531  double pkinf = pkinf_(k);
    470   return pkinf *tf*tf;
    471 }
    472 
    473 //------------------------------------
    474 PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
    475   : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(PK)
    476 {
    477 }
    478 
    479 PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
    480   : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(PK)
    481 {
    482 }
    483 
    484 PkSpectrumZ::~PkSpectrumZ(void)
    485 {
    486 }
    487 
    488 void PkSpectrumZ::SetTypSpec(ReturnSpectrum typspec)
    489 // typsec = PK : compute Pk(k)
    490 //        = DELTA : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
    491 {
    492   typspec_ = typspec;
    493 }
    494 
    495 double PkSpectrumZ::operator() (double k)
     532  double d1 = d1_(z);
     533
     534  double v = pkinf * (tf*tf) * (d1*d1);
     535  if(typspec_ == DELTA) v *= k*k*k/(2.*M_PI*M_PI);
     536
     537  return scale_ * v;
     538}
     539
     540///////////////////////////////////////////////////////////
     541//********************** PkTabulate *********************//
     542///////////////////////////////////////////////////////////
     543
     544PkTabulate::PkTabulate(void)
     545: kmin_(1.) , kmax_(-1.) , interptyp_(0)
     546{
     547  k_.resize(0);
     548  pk_.resize(0);
     549}
     550
     551PkTabulate::PkTabulate(PkTabulate& pkz)
     552 : kmin_(pkz.kmin_) , kmax_(pkz.kmax_) , interptyp_(pkz.interptyp_) , k_(pkz.k_) , pk_(pkz.pk_)
     553{
     554}
     555
     556
     557PkTabulate::~PkTabulate(void)
     558{
     559}
     560
     561void PkTabulate ::SetInterpTyp(int typ)
     562// see comment in InterpTab
     563{
     564 if(typ<0) typ=0; else if(typ>2) typ=2;
     565 interptyp_ = typ;
     566}
     567
     568double PkTabulate::operator() (double k)
     569{
     570  double v = InterpTab(k,k_,pk_,interptyp_);
     571  if(typspec_ == DELTA) v *= k*k*k/(2.*M_PI*M_PI);
     572  return scale_ * v;
     573}
     574
     575double PkTabulate::operator() (double k,double z)
     576{
     577  cout<<"Error: PkTabulate::operator(double k,double z) not implemented"<<endl;
     578  throw AllocationError("Error: PkTabulate::operator(double k,double z) not implemented");
     579}
     580
     581int PkTabulate::ReadCAMB(string filename, double h100, double zreftab)
     582{
     583 FILE *file = fopen(filename.c_str(),"r");
     584 if(file==NULL) return -1;
     585 cout<<"PkTabulate::ReadCAMB: fn="<<filename<<" h100="<<h100<<" zreftab = "<<zreftab<<endl;
     586
     587 const int lenline = 512;
     588 char *line = new char[lenline];
     589
     590 k_.resize(0); pk_.resize(0);
     591 while ( fgets(line,lenline,file) != NULL ) {
     592   double k, pk;
     593   sscanf(line,"%lf %lf",&k,&pk);
     594   k *= h100;             // convert h    Mpc^-1  ->  Mpc^-1
     595   pk /= h100*h100*h100;  // convert h^-3 Mpc^3   ->  Mpc^3
     596   k_.push_back(k);
     597   pk_.push_back(pk);
     598 }
     599
     600 SetZ(zreftab);
     601 cout<<"PkTabulate::ReadCAMB nread="<<pk_.size()<<"zref="<<GetZ()<<endl;
     602 delete [] line;
     603
     604 return (int)pk_.size();
     605}
     606
     607///////////////////////////////////////////////////////////
     608//********************* PkEisenstein ********************//
     609///////////////////////////////////////////////////////////
     610
     611PkEisenstein::PkEisenstein(InitialPowerLaw& pkinf,TransfertEisenstein& tf,GrowthEisenstein& d1,double zref)
     612  : pkinf_(pkinf) , tf_(tf) , d1_(d1)
     613{
     614  zref_ = zref;
     615}
     616
     617PkEisenstein::PkEisenstein(PkEisenstein& pkz)
     618  : pkinf_(pkz.pkinf_) , tf_(pkz.tf_) , d1_(pkz.d1_)
     619{
     620}
     621
     622PkEisenstein::~PkEisenstein(void)
     623{
     624}
     625
     626double PkEisenstein::operator() (double k)
    496627{
    497628  return (*this)(k,zref_);
    498629}
    499630
    500 double PkSpectrumZ::operator() (double k,double z)
    501 {
     631double PkEisenstein::operator() (double k,double z)
     632{
     633  double tf = tf_(k);
     634  double pkinf = pkinf_(k);
    502635  double d1 = d1_(z);
    503636
    504   double v = pk0_(k) * d1*d1;
    505   if(typspec_==DELTA) v *= k*k*k/(2.*M_PI*M_PI);
     637  double v = pkinf * (tf*tf) * (d1*d1);
     638  if(typspec_ == DELTA) v *= k*k*k/(2.*M_PI*M_PI);
    506639
    507640  return scale_ * v;
    508641}
    509 
    510642
    511643

  • trunk/Cosmo/SimLSS/pkspectrum.h

    r3802 r3805  
    1010class InitialSpectrum : public GenericFunc {
    1111public:
    12   InitialSpectrum(double n,double a=1.);
    13   InitialSpectrum(InitialSpectrum& pkinf);
    14   virtual ~InitialSpectrum(void);
     12  InitialSpectrum(void) {};
     13  InitialSpectrum(InitialSpectrum& pkinf) {};
     14  virtual ~InitialSpectrum(void) {};
     15};
     16
     17//-----------------------------------------------------------------------------------
     18class InitialPowerLaw : public InitialSpectrum {
     19public:
     20  InitialPowerLaw(double n,double a=1.);
     21  InitialPowerLaw(InitialPowerLaw& pkinf);
     22  virtual ~InitialPowerLaw(void);
    1523  virtual double operator() (double k) {return A_ * pow(k,n_);}
    1624  void SetNorm(double a) {A_ = a;}
    … …  
    2129
    2230//-----------------------------------------------------------------------------------
    23 class TransfertEisenstein : public GenericFunc {
     31class TransfertFunction : public GenericFunc {
     32public:
     33  TransfertFunction(void) {};
     34  virtual ~TransfertFunction(void) {};
     35};
     36
     37//-----------------------------------------------------------------------------------
     38class TransfertEisenstein : public TransfertFunction {
    2439public:
    2540
    … …  
    5368
    5469//-----------------------------------------------------------------------------------
    55 class TransfertTabulate : public GenericFunc {
     70class TransfertTabulate : public TransfertFunction {
    5671public:
    5772  TransfertTabulate(void);
    … …  
    6277  void SetInterpTyp(int typ=0);
    6378  int ReadCMBFast(string filename,double h100,double OmegaCDM0,double OmegaBaryon0);
     79  int ReadCAMB(string filename, double h100=0.71);
    6480protected:
    6581  double kmin_,kmax_;
    … …  
    7288class GrowthFactor : public GenericFunc {
    7389public:
    74   GrowthFactor(double OmegaMatter0,double OmegaLambda0);
    75   GrowthFactor(GrowthFactor& d1);
    76   virtual ~GrowthFactor(void);
     90  GrowthFactor(void) {};
     91  virtual ~GrowthFactor(void) {};
     92  virtual double DsDz(double z, double);
     93};
     94
     95//-----------------------------------------------------------------------------------
     96class GrowthEisenstein : public GrowthFactor {
     97public:
     98  GrowthEisenstein(double OmegaMatter0,double OmegaLambda0);
     99  GrowthEisenstein(GrowthEisenstein& d1);
     100  virtual ~GrowthEisenstein(void);
    77101  virtual double operator() (double z);
    78102  virtual double DsDz(double z,double dzinc=0.01);
    … …  
    85109
    86110//-----------------------------------------------------------------------------------
    87 class PkSpectrum0 : public GenericFunc {
    88 public:
    89   PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf);
    90   PkSpectrum0(PkSpectrum0& pk0);
    91   virtual ~PkSpectrum0(void);
    92   virtual double operator() (double z);
    93   InitialSpectrum& GetPkIni(void) {return pkinf_;}
    94   TransfertEisenstein& GetTransfert(void) {return tf_;}
    95 protected:
    96   InitialSpectrum& pkinf_;
    97   TransfertEisenstein& tf_;
    98 };
    99 
    100 //-----------------------------------------------------------------------------------
    101 class PkSpectrumZ : public GenericFunc {
    102 public:
     111class PkSpectrum : public GenericFunc {
     112public:
     113  // typsec = PK : compute Pk(k)
     114  //        = DELTA : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
    103115  typedef enum {PK=0, DELTA=1} ReturnSpectrum;
    104   PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref=0.);
    105   PkSpectrumZ(PkSpectrumZ& pkz);
    106   virtual ~PkSpectrumZ(void);
    107   virtual double operator() (double k);
    108   virtual double operator() (double k,double z);
    109   void   SetZ(double z) {zref_ = z;}
    110   double GetZ(void) {return zref_;}
    111   void SetTypSpec(ReturnSpectrum typspec=PK);
    112   void SetScale(double scale=1.) {scale_=scale;}
    113   double GetScale(void) {return scale_;}
    114   PkSpectrum0& GetPk0(void) {return pk0_;}
    115   GrowthFactor& GetGrowthFactor(void) {return d1_;}
    116 protected:
    117   PkSpectrum0& pk0_;
    118   GrowthFactor& d1_;
     116
     117  PkSpectrum(void);
     118  PkSpectrum(PkSpectrum& pk);
     119  virtual ~PkSpectrum(void) {};
     120  virtual void   SetZ(double z) {zref_ = z;}
     121  virtual double GetZ(void) {return zref_;}
     122  virtual void SetScale(double scale=1.) {scale_ = scale;}
     123  virtual double GetScale(void) {return scale_;}
     124  virtual void SetTypSpec(ReturnSpectrum typspec=PK) {typspec_ = typspec;}
     125  virtual ReturnSpectrum GetTypSpec(void) {return typspec_ ;}
     126protected:
    119127  double zref_, scale_;
    120128  ReturnSpectrum typspec_;
    … …  
    122130
    123131//-----------------------------------------------------------------------------------
     132class PkSpecCalc : public PkSpectrum {
     133public:
     134  PkSpecCalc(InitialSpectrum& pkinf,TransfertFunction& tf,GrowthFactor& d1,double zref=0.);
     135  PkSpecCalc(PkSpecCalc& pkz);
     136  virtual ~PkSpecCalc(void);
     137  virtual double operator() (double k);
     138  virtual double operator() (double k,double z);
     139  InitialSpectrum& GetPkIni(void) {return pkinf_;}
     140  TransfertFunction& GetTransfert(void) {return tf_;}
     141  GrowthFactor& GetGrowthFactor(void) {return d1_;}
     142protected:
     143  InitialSpectrum& pkinf_;
     144  TransfertFunction& tf_;
     145  GrowthFactor& d1_;
     146};
     147
     148//-----------------------------------------------------------------------------------
     149class PkTabulate : public PkSpectrum {
     150public:
     151  PkTabulate(void);
     152  PkTabulate(PkTabulate& pkz);
     153  virtual ~PkTabulate(void);
     154  virtual double operator() (double k);
     155  virtual double operator() (double k,double z);
     156  int NPoints(void) {return k_.size();}
     157  void SetInterpTyp(int typ=0);
     158  int ReadCAMB(string filename, double h100=0.71, double zreftab = 0.);
     159protected:
     160  double kmin_,kmax_;
     161  int interptyp_;
     162  vector<double> k_, pk_;
     163};
     164
     165//-----------------------------------------------------------------------------------
     166class PkEisenstein : public PkSpectrum {
     167public:
     168  PkEisenstein(InitialPowerLaw& pkinf,TransfertEisenstein& tf,GrowthEisenstein& d1,double zref=0.);
     169  PkEisenstein(PkEisenstein& pkz);
     170  virtual ~PkEisenstein(void);
     171  virtual double operator() (double k);
     172  virtual double operator() (double k,double z);
     173  InitialPowerLaw& GetPkIni(void) {return pkinf_;}
     174  TransfertEisenstein& GetTransfert(void) {return tf_;}
     175  GrowthEisenstein& GetGrowthFactor(void) {return d1_;}
     176protected:
     177  InitialPowerLaw& pkinf_;
     178  TransfertEisenstein& tf_;
     179  GrowthEisenstein& d1_;
     180};
     181
     182//-----------------------------------------------------------------------------------
    124183class VarianceSpectrum : public GenericFunc {
    125184public:
    … …  
    128187
    129188  VarianceSpectrum(GenericFunc& pk,double R,TypeFilter typfilter);
    130   VarianceSpectrum(VarianceSpectrum& pkinf);
     189  VarianceSpectrum(VarianceSpectrum& vpk);
    131190  virtual ~VarianceSpectrum(void);
    132191
Note: See TracChangeset for help on using the changeset viewer.