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

Timestamp:

Nov 14, 2007, 7:25:34 PM (18 years ago)

Author:

cmv

Message:

1-/ pkspectrum: GrowthFactor + PkSpectrum on enleve les memorisations

des calculs precedents, tout est recalcule

2-/ cmvfitpk: mise de Ob+Ol+Oc+h100+ns dans le fit

cmv 14/11/2007

File:

• trunk/Cosmo/SimLSS/pkspectrum.cc (modified) (9 diffs)

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -31 +31 @@
 InitialSpectrum::~InitialSpectrum(void)
 {
-}
-
-void InitialSpectrum::SetNorm(double a)
-{
-  A_ = a;
-}
-
-void InitialSpectrum::SetSlope(double n)
-{
-  n_ = n;
 }
 
@@ -177 +167 @@
 }
 
-bool TransfertEisenstein::SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon)
+bool TransfertEisenstein::SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0)
 // Changement des valeurs des parametres (suivi de re-init eventuel)
+// Si h100,Omega...<=0. alors pas de changement, on garde l'ancienne valeur
 {
  bool haschanged = false;
@@ -185 +176 @@
  if(OmegaCDM0>0.) {Oc_ = OmegaCDM0; haschanged = true;}
  if(OmegaBaryon0>0.) {Ob_ = OmegaBaryon0; haschanged = true;}
- if(tcmb>0.) {tcmb_ = tcmb; haschanged = true;}
- if(nobaryon!=nobaryon_) {nobaryon_ = nobaryon; haschanged = true;}
 
  // et on recalcule les initialisations
@@ -371 +360 @@
 // Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
 GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
-  : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)
+  : O0_(OmegaMatter0) , Ol_(OmegaLambda0)
 {
   if(OmegaMatter0==0.) {
@@ -377 +366 @@
     throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
   }
-  norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
-  norm_ = (*this)(0.);
-  cout<<"GrowthFactor::GrowthFactor : norm="<<norm_<<endl;
 }
 
 GrowthFactor::GrowthFactor(GrowthFactor& d1)
-  : O0_(d1.O0_) , Ol_(d1.Ol_) , Ok_(d1.Ok_) , norm_(d1.norm_)
+  : O0_(d1.O0_) , Ol_(d1.Ol_)
 {
 }
@@ -396 +382 @@
  z += 1.;
  double z2 = z*z, z3 = z2*z;
- double den = Ol_ + Ok_*z2 + O0_*z3;
+
+ // Calcul de la normalisation (pour z=0 -> growth=1.)
+ double D1z0 = pow(O0_,4./7.) - Ol_ + (1.+O0_/2.)*(1.+Ol_/70.);
+ D1z0 = 2.5*O0_ / D1z0;
+
+ // Calcul du growthfactor pour z
+ double Ok = 1. - O0_ - Ol_;
+ double den = Ol_ + Ok*z2 + O0_*z3;
  double o0z = O0_ *z3 / den;
  double olz = Ol_ / den;
 
- // 4./7. = 0.571429
- double D1z = pow(o0z,0.571429) - olz + (1.+o0z/2.)*(1.+olz/70.);
+ double D1z = pow(o0z,4./7.) - olz + (1.+o0z/2.)*(1.+olz/70.);
  D1z = 2.5*o0z / z / D1z;
 
- return D1z / norm_;
-}
-
-bool GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
-{
- bool haschanged = false;
-
- if(OmegaMatter0>0.) {O0_ = OmegaMatter0; haschanged = true;}
- if(fabs(OmegaLambda0+12345.)>1e-6) {Ol_ = OmegaLambda0; haschanged = true;}
-
- // et on recalcule les initialisations
- if(haschanged) {
-   Ok_ = 1. - O0_ - Ol_;
-   norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
-   norm_ = (*this)(0.);
- }
-
- return haschanged;
+ return D1z / D1z0;
+}
+
+void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
+{
+ if(OmegaMatter0>0.) O0_ = OmegaMatter0;
+ Ol_ = OmegaLambda0;
+}
+
+bool GrowthFactor::SetParTo(double OmegaMatter0)
+// idem precedent sans changer OmegaLambda0
+{
+ if(OmegaMatter0<=0.) return false;
+ O0_ = OmegaMatter0;
+ return true;
 }
 
@@ -453 +442 @@
 PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
   : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(PK)
-  , zold_(-1.) , d1old_(1.)
 {
 }
@@ -459 +447 @@
 PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
   : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(PK)
-  , zold_(pkz.zold_) , d1old_(pkz.d1old_)
 {
 }
@@ -481 +468 @@
 double PkSpectrumZ::operator() (double k,double z)
 {
-  double d1;
-  if(z == zold_) d1 = d1old_;
-    else {d1 = d1old_ = d1_(z); zold_ = z;}
+  double d1 = d1_(z);
 
   double v = pk0_(k) * d1*d1;