Changeset 3348 in Sophya for trunk/Cosmo/SimLSS


Ignore:
Timestamp:
Oct 11, 2007, 4:37:03 PM (18 years ago)
Author:
cmv
Message:
  • definition des options par enum
  • mise en variable privee du rayon R du filtre de VarianceSpectrum Elle disparait des arguments des methodes: Variance FindMaximum FindLimits

cmv 11/10/2007

Location:
trunk/Cosmo/SimLSS
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvchkwhu.cc

    r3317 r3348  
    5959   double k = pow(10.,lnk);
    6060   //if(fabs(lnk-log10(0.0186209))>1e-5) continue;
    61    Tfull.SetReturnPart(0); double tf = Tfull(k);
     61   Tfull.SetReturnPart(TransfertEisenstein::ALL); double tf = Tfull(k);
    6262   double tfnw = Tnowig(k);
    6363   double tfnb = Tnob(k);
    64    Tfull.SetReturnPart(1); double tfc = Tfull(k);
    65    Tfull.SetReturnPart(2); double tfb = Tfull(k);
     64   Tfull.SetReturnPart(TransfertEisenstein::CDM); double tfc = Tfull(k);
     65   Tfull.SetReturnPart(TransfertEisenstein::BARYON); double tfb = Tfull(k);
    6666   double hutfc, hutfb;
    6767   double hutf = TFfit_onek(k,&hutfb,&hutfc);

  • trunk/Cosmo/SimLSS/cmvtvarspec.cc

    r3329 r3348  
    9191 //----
    9292 cout<<endl<<"Filtrage top_hat"<<endl;
    93  VarianceSpectrum varpk_th(pkz,0);
    94  double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
     93 VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
     94 double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
    9595 double pkmax_th = varpk_th(kfind_th);
    9696 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
    9797 k1=kmin, k2=kmax;
    98  rc = varpk_th.FindLimits(R,pkmax_th*fracmax,k1,k2,eps);
     98 rc = varpk_th.FindLimits(pkmax_th*fracmax,k1,k2,eps);
    9999 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
    100100
    101101 varpk_th.SetInteg(0.01,dlk,-1.,4);
    102  cout<<"varpk_th="<<varpk_th.Variance(R,k1,k2)<<endl;
     102 cout<<"varpk_th="<<varpk_th.Variance(k1,k2)<<endl;
    103103
    104104 //----
    105105 cout<<endl<<"Filtrage gaussien"<<endl;
    106  VarianceSpectrum varpk_ga(pkz,1);
    107  double kfind_ga = varpk_ga.FindMaximum(Rg,kmin,kmax,eps);
     106 VarianceSpectrum varpk_ga(pkz,Rg,VarianceSpectrum::GAUSSIAN);
     107 double kfind_ga = varpk_ga.FindMaximum(kmin,kmax,eps);
    108108 double pkmax_ga = varpk_ga(kfind_ga);
    109109 cout<<"kfind_ga = "<<kfind_ga<<" -> "<<log10(kfind_ga)<<", integrand="<<pkmax_ga<<endl;
    110110 k1=kmin, k2=kmax;
    111  rc = varpk_ga.FindLimits(Rg,pkmax_ga*fracmax,k1,k2,eps);
     111 rc = varpk_ga.FindLimits(pkmax_ga*fracmax,k1,k2,eps);
    112112 cout<<"limit_ga: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
    113113
    114114 varpk_ga.SetInteg(0.01,dlk,-1.,4);
    115  cout<<"varpk_ga="<<varpk_ga.Variance(Rg,k1,k2)<<endl;
     115 cout<<"varpk_ga="<<varpk_ga.Variance(k1,k2)<<endl;
    116116 
    117117 //----
    118118  cout<<endl<<"Filtrage 1 (integrale du spectre)"<<endl;
    119  VarianceSpectrum varpk_int(pkz,2);
    120  double kfind_int = varpk_int.FindMaximum(Rg,kmin,kmax,eps);
     119 VarianceSpectrum varpk_int(pkz,Rg,VarianceSpectrum::NOFILTER);
     120 double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
    121121 double pkmax_int = varpk_int(kfind_int);
    122122 cout<<"kfind_int = "<<kfind_int<<" -> "<<log10(kfind_int)<<", integrand="<<pkmax_int<<endl;
    123123 k1=kmin, k2=kmax;
    124  rc = varpk_int.FindLimits(Rg,pkmax_int*fracmax,k1,k2,eps);
     124 rc = varpk_int.FindLimits(pkmax_int*fracmax,k1,k2,eps);
    125125 cout<<"limit_int: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
    126126
    127127 varpk_int.SetInteg(0.01,dlk,-1.,4);
    128  cout<<"varpk_int="<<varpk_int.Variance(Rg,k1,k2)<<endl;
     128 cout<<"varpk_int="<<varpk_int.Variance(k1,k2)<<endl;
    129129 
    130130//-----------------------------------------------------------------

  • trunk/Cosmo/SimLSS/pkspectrum.cc

    r3329 r3348  
    5252  : lp_(lp)
    5353  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
    54   , nobaryon_(nobaryon) , nooscenv_(0), retpart_(0)
     54  , nobaryon_(nobaryon) , nooscenv_(0), retpart_(ALL)
    5555{
    5656 zero_();
    … …  
    192192}
    193193
    194 void TransfertEisenstein::SetReturnPart(unsigned short retpart)
     194void TransfertEisenstein::SetReturnPart(ReturnPart retpart)
    195195// To return only baryon or CDM part part of transfert function
    196 // retpart = 1 : return only CDM part of transfert function
    197 // retpart = 2 : return only Baryon part of transfert function
    198 // retpart = anything else: return only full transfert function
     196// retpart = ALL: return full transfert function
     197//         = CDM : return only CDM part of transfert function
     198//         = BARYON : return only Baryon part of transfert function
    199199// WARNING: only relevant for nobaryon_=false AND nooscenv!=2
    200200{
    201  if(retpart!=1 && retpart!=2) retpart = 0;
    202201 retpart_ = retpart;
    203202}
    … …  
    244243  double f = 1. / (1. + pow(k*s_/5.4,4.));
    245244  double Tc = f*T0tild(k,1.,betac_) + (1.-f)*T0tild(k,alphac_,betac_);
    246   if(retpart_ == 1) return Tc;
     245  if(retpart_ == CDM) return Tc;
    247246
    248247  // --- Baryons
    … …  
    264263  Tb += alphab_/(1.+pow(betab_/(k*s_),3.)) * exp(-pow(k/ksilk_,1.4));
    265264  Tb *= j0kst;
    266   if(retpart_ == 2) return Tb;
     265  if(retpart_ == BARYON) return Tb;
    267266
    268267  // --- Total
    … …  
    413412//------------------------------------
    414413PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
    415   : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(0)
     414  : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(PK)
    416415  , zold_(-1.) , d1old_(1.)
    417416{
    … …  
    419418
    420419PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
    421   : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(0)
     420  : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(PK)
    422421  , zold_(pkz.zold_) , d1old_(pkz.d1old_)
    423422{
    … …  
    428427}
    429428
    430 void PkSpectrumZ::SetTypSpec(unsigned short typspec)
    431   // typsec = 0 : compute Pk(k)
    432   //        = 1 : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
    433 {
    434   if(typspec>1) {
    435     cout<<"PkSpectrumZ::SetTypSpec: Error bad typspec value: "<<typspec<<endl;
    436     throw ParmError("PkSpectrumZ::SetTypSpec: Error bad typspec value");
    437   }
     429void PkSpectrumZ::SetTypSpec(ReturnSpectrum typspec)
     430// typsec = PK : compute Pk(k)
     431//        = DELTA : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
     432{
    438433  typspec_ = typspec;
    439434}
    … …  
    451446
    452447  double v = pk0_(k) * d1*d1;
    453   if(typspec_) v *= k*k*k/(2.*M_PI*M_PI);
     448  if(typspec_==DELTA) v *= k*k*k/(2.*M_PI*M_PI);
    454449
    455450  return scale_ * v;
    … …  
    462457///////////////////////////////////////////////////////////
    463458
    464 VarianceSpectrum::VarianceSpectrum(GenericFunc& pk,unsigned short typfilter=0)
    465   : pk_(pk) , R_(0.)
    466 {
     459VarianceSpectrum::VarianceSpectrum(GenericFunc& pk,double R,TypeFilter typfilter)
     460  : pk_(pk)
     461{
     462  SetRadius(R);
    467463  SetFilter(typfilter);
    468464}
    … …  
    478474}
    479475
     476void VarianceSpectrum::SetRadius(double R)
     477// R = taille du filter top-hat ou gaussien
     478{
     479 if(R<=0.) {
     480    cout<<"VarianceSpectrum::SetRadius: Error R<=0"<<endl;
     481    throw ParmError("VarianceSpectrum::SetRadius: Error R<=0");
     482  }
     483  R_ = R;
     484}
     485
    480486//------------------------------------
    481 void VarianceSpectrum::SetFilter(unsigned short typfilter)
    482 // typfilter = 0 : spherical 3D top-hat
    483 //           = 1 : spherical 3D gaussian
    484 //           = 2 : no filter juste integrate spectrum)
    485 {
    486   if(typfilter>2) {
    487     cout<<"VarianceSpectrum::SetFilter: Error bad value for type of filter"<<endl;
    488     throw ParmError("VarianceSpectrum::SetFilter: Error bad value for type of filter");
    489   }
     487void VarianceSpectrum::SetFilter(TypeFilter typfilter)
     488// typfilter = TOPHAT : spherical 3D top-hat
     489//           = GAUSSIAN : spherical 3D gaussian
     490//           = NOFILTER : no filter juste integrate spectrum)
     491// Remarque:
     492//   la meilleure approximation du filtre top-hat (R) est un filtre gaussien avec (Rg=R/sqrt(5))
     493{
    490494  typfilter_ = typfilter;
    491495}
    … …  
    507511{
    508512  // Just integrate the spectrum without filtering
    509   if(typfilter_ == 2) return 1.;
     513  if(typfilter_ == NOFILTER) return 1.;
    510514
    511515  double x2 = x*x;
    … …  
    515519  //             DL(x) = 1-x^2*(1-x^2/2)
    516520  //             pour x<0.01  |DL(x)-G(X)^2|<2.0e-13
    517   if(typfilter_ == 1)
     521  if(typfilter_ == GAUSSIAN)
    518522    if(x<0.01) return 1.-x2*(1.-x2/2.); else return exp(-x2);
    519523
    … …  
    534538}
    535539
    536 double VarianceSpectrum::Variance(double R,double kmin,double kmax)
     540double VarianceSpectrum::Variance(double kmin,double kmax)
    537541// Compute variance of spectrum pk_ by integration
    538542// Input:
    539 //     R = taille du filter top-hat ou gaussien
    540543//     kmin,kmax = bornes en k de l'integrale pour calculer la variance
    541544// Return:
    542545//     valeur de la variance (sigma^2)
    543546// Remarque:
    544 //   la meilleure approximation du filtre top-hat (R) est un filtre gaussien avec (Rg=R/sqrt(5))
    545547//   la variance renvoyee est la variance de la masse
    546548{
    547   if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
    548     cout<<"VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
    549     throw ParmError("VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax");
    550   }
    551 
    552   R_ = R;
     549  if(kmin<=0 || kmax<=0. || kmin>=kmax) {
     550    cout<<"VarianceSpectrum::Variance: Error kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
     551    throw ParmError("VarianceSpectrum::Variance: Error kmin<=0 or kmax<=0 or kmin>=kmax");
     552  }
     553
    553554  double lkmin = log10(kmin), lkmax = log10(kmax);
    554555
    … …  
    559560
    560561//------------------------------------
    561 double VarianceSpectrum::FindMaximum(double R,double kmin,double kmax,double eps)
     562double VarianceSpectrum::FindMaximum(double kmin,double kmax,double eps)
    562563// Retourne le maximum de la fonction a integrer
    563564// La recherche a lieu entre [kmin,kmax] par pas logarithmiques
    564565// Input:
    565 //     R : taille du filter top-hat ou gaussien
    566566//     kmin,kmax : intervalle de recherche
    567567//     eps : precision requise sur les valeurs
    … …  
    569569//     position (en k) du maximum
    570570{
    571   if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
    572     cout<<"VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
    573     throw ParmError("VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
    574   }
    575 
    576   R_ = R;
     571  if(kmin<=0 || kmax<=0. || kmin>=kmax) {
     572    cout<<"VarianceSpectrum::FindMaximum: Error kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
     573    throw ParmError("VarianceSpectrum::FindMaximum: Error kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
     574  }
    577575
    578576  int n = 10; // toujours >2
    … …  
    600598}
    601599
    602 int VarianceSpectrum::FindLimits(double R,double high,double &kmin,double &kmax,double eps)
     600int VarianceSpectrum::FindLimits(double high,double &kmin,double &kmax,double eps)
    603601// Retourne "[kmin,kmax]" tel que la fonction a integrer soit "f(k) <= high"
    604602//   La recherche a lieu entre [kmin,kmax] par pas logarithmiques
    605603// Input:
    606 //     R : taille du filter top-hat ou gaussien
    607604//     kmin,kmax : intervalle de recherche
    608605//     eps : precision requise sur les valeurs kmin et kmax
    … …  
    615612//     rc |= 4 "f(k) < high pour tout k"   (bit2 =1)
    616613{
    617   if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax  || eps<=0.) {
    618     cout<<"VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
    619     throw ParmError("VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
    620   }
    621 
    622   R_ = R;
     614  if(kmin<=0 || kmax<=0. || kmin>=kmax  || eps<=0.) {
     615    cout<<"VarianceSpectrum::FindLimits: Error kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
     616    throw ParmError("VarianceSpectrum::FindLimits: Error kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
     617  }
     618
    623619  int n = 10; // toujours >2
    624620

  • trunk/Cosmo/SimLSS/pkspectrum.h

    r3325 r3348  
    2323class TransfertEisenstein : public GenericFunc {
    2424public:
     25
     26  typedef enum{ALL=0, CDM=1, BARYON=2} ReturnPart;
     27
    2528  TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon=false,int lp=0);
    2629  TransfertEisenstein(TransfertEisenstein& tf);
    … …  
    2932  double KPeak(void);
    3033  void SetNoOscEnv(unsigned short nooscenv=0);
    31   void SetReturnPart(unsigned short retpart=0);
     34  void SetReturnPart(ReturnPart retpart=ALL);
    3235protected:
    3336  int lp_;
    … …  
    3942
    4043  bool nobaryon_;
    41   unsigned short nooscenv_, retpart_;
     44  unsigned short nooscenv_;
     45  ReturnPart retpart_;
    4246
    4347  double T0tild(double k,double alphac,double betac);
    … …  
    9498class PkSpectrumZ : public GenericFunc {
    9599public:
     100  typedef enum {PK=0, DELTA=1} ReturnSpectrum;
    96101  PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref=0.);
    97102  PkSpectrumZ(PkSpectrumZ& pkz);
    … …  
    101106  inline void   SetZ(double z) {zref_ = z;}
    102107  inline double GetZ(void) {return zref_;}
    103   void SetTypSpec(unsigned short typspec=0);
     108  void SetTypSpec(ReturnSpectrum typspec=PK);
    104109  inline void SetScale(double scale=1.) {scale_=scale; zold_=-1.;}
    105110  inline double GetScale(void) {return scale_;}
    … …  
    110115  GrowthFactor& d1_;
    111116  double zref_, scale_;
    112   unsigned short typspec_;
     117  ReturnSpectrum typspec_;
    113118  mutable double zold_, d1old_;
    114119};
    … …  
    117122class VarianceSpectrum : public GenericFunc {
    118123public:
    119   VarianceSpectrum(GenericFunc& pk,unsigned short typfilter);
     124
     125  typedef enum {TOPHAT=0, GAUSSIAN=1, NOFILTER=2} TypeFilter;
     126
     127  VarianceSpectrum(GenericFunc& pk,double R,TypeFilter typfilter);
    120128  VarianceSpectrum(VarianceSpectrum& pkinf);
    121129  virtual ~VarianceSpectrum(void);
    122130
    123   void SetFilter(unsigned short typfilter=0);
     131  void SetRadius(double R);
     132  void SetFilter(TypeFilter typfilter=TOPHAT);
    124133  void SetInteg(double dperc=0.1,double dlogkinc=-1.,double dlogkmax=-1.,unsigned short glorder=4);
    125134
    126   double Variance(double R,double kmin,double kmax);
     135  double Variance(double kmin,double kmax);
    127136
    128137  // ATTENTION: La fonction a integrer est : f(k)dk = k^3*Pk(k)/(2Pi^2) *filter2(k*R) *dk/k
    … …  
    131140
    132141  // Aide a l'integration
    133   double FindMaximum(double R,double kmin,double kmax,double eps=1.e-3);
    134   int FindLimits(double R,double high,double &kmin,double &kmax,double eps=1.e-3);
     142  double FindMaximum(double kmin,double kmax,double eps=1.e-3);
     143  int FindLimits(double high,double &kmin,double &kmax,double eps=1.e-3);
    135144
    136145protected:
    137146
    138147  GenericFunc& pk_;
    139   unsigned short typfilter_;
     148  TypeFilter typfilter_;
    140149  double R_;
    141150
Note: See TracChangeset for help on using the changeset viewer.