Changeset 3348 in Sophya for trunk/Cosmo/SimLSS/pkspectrum.cc


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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.