Changeset 3924 in Sophya for trunk/Cosmo/SimLSS/genefluct3d.cc

Timestamp:

Dec 10, 2010, 4:38:23 PM (15 years ago)

Author:

cmv

Message:

intro choix dabs calcul de beta pour les redshift-distorsions, cmv 10/12/2010

File:

• trunk/Cosmo/SimLSS/genefluct3d.cc (modified) (6 diffs)

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -904 +904 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ToRedshiftSpace(void)
+void GeneFluct3D::ToRedshiftSpace(int type_beta,double fixed_beta)
 // Apply redshift distortion corrections
-// ---- ATTENTION: Le spectre Pk doit etre (dRho/Rho)(k)
+// P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
+// type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
+//            >0 compute factor with fix beta=Omega^fixed_beta
+//            <0 compute factor with fix beta=fixed_beta
+// ---- ATTENTION:
+// 1- Le spectre Pk doit etre (dRho/Rho)(k)
+// 2- pas d'evolution en redhsuft possible,
+//    le spectre doit etre calcule au redshift final de la boite
 {
  double zpk = compute_pk_redsh_ref_;
- double beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
- if(lp_>0) cout<<"--- ToRedshiftSpace --- at z="<<zpk<<", Beta="<<beta
-               <<" (comp. for dz="<<good_dzinc_<<")"<<endl;
+ if(lp_>0) cout<<"--- ToRedshiftSpace --- at z="<<zpk<<", type_beta="<<type_beta<<endl;
+
+ double beta = 0.;
+ if(type_beta>0) {
+   double omegam = cosmo_->Omatter(zpk); beta = pow(omegam,fixed_beta);
+   if(lp_>0) cout<<"...beta = OmegaM^expo = "<<omegam<<"^"<<fixed_beta<<" = "<<beta<<endl;
+ } else if(type_beta<0) {
+   beta = fixed_beta;
+   if(lp_>0) cout<<"...beta = "<<beta<<endl;
+ } else {
+   beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+   if(lp_>0) cout<<"...beta = -(1+z)*D'/D = "<<beta<<endl;
+ }
+
  check_array_alloc();
 
@@ -929 +947 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ToVelLoS(void)
+void GeneFluct3D::ToVelLoS(int type_beta,double fixed_beta)
 /*
+// P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
+// type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
+//            >0 compute factor with fix beta=Omega^fixed_beta
+//            <0 compute factor with fix beta=fixed_beta
 ---- Vitesse particuliere
 Un atome au redshift "z" emet a la longueur d'onde "Le" si il est au repos
@@ -948 +970 @@
 {
  double zpk = compute_pk_redsh_ref_;
- double dpsd = -cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
- if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", -H*(1+z)*D'/D="<<dpsd
-               <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
+ if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", type_beta="<<type_beta<<endl;
+
+ double beta = 0;
+ if(type_beta>0) {
+   double omegam = cosmo_->Omatter(zpk);
+   beta = pow(omegam,fixed_beta);
+   if(lp_>0) cout<<"...beta = OmegaM^expo = "<<omegam<<"^"<<fixed_beta<<"="<<beta<<endl;
+ } else if(type_beta<0) {
+   beta = fixed_beta;
+   if(lp_>0) cout<<"...beta="<<beta<<endl;
+ } else {
+   beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+   if(lp_>0) cout<<"...beta = -(1+z)*D'/D = "<<beta<<endl;
+ }
+ double dpsd = cosmo_->H(zpk) * beta;
+ if(lp_>0) cout<<"......H*beta = "<<dpsd<<" (km/s)/Mpc"<<endl;
+
  check_array_alloc();
 
@@ -1019 +1055 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ApplyVelLosGrowthFactor(int type_evol)
+void GeneFluct3D::ApplyVelLosGrowthFactor(int type_evol,int type_beta,double fixed_beta)
 // Apply Growth(z) to transverse velocity real space cube
-{
- if(lp_>0) cout<<"--- ApplyVelLosGrowthFactor: evol="<<type_evol<<endl;
+// P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
+// type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
+//            >0 compute factor with fix beta=Omega^fixed_beta
+//            <0 compute factor with fix beta=fixed_beta
+{
+ if(lp_>0) cout<<"--- ApplyVelLosGrowthFactor: evol="<<type_evol<<", type_beta="<<type_beta<<endl;
+
  check_array_alloc();
 
@@ -1036 +1077 @@
  double zpk = compute_pk_redsh_ref_;
  double grw_orig = (*growth_)(zpk);
- double dpsd_orig = - cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / grw_orig;
- if(lp_>0) cout<<"    original growth="<<grw_orig<<" dpsd="<<dpsd_orig<<" computed at z="<<zpk<<endl;
+ double beta_orig = 0.;
+ if(type_beta>0) {double omegam = cosmo_->Omatter(zpk); beta_orig = pow(omegam,fixed_beta);}
+ else if(type_beta<0) beta_orig = fixed_beta;
+ else beta_orig = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+ double dpsd_orig = cosmo_->H(zpk) * beta_orig;
+ if(lp_>0) cout<<"...original at z="<<zpk<<": growth="<<grw_orig<<" beta="<<beta_orig
+               <<" H*beta="<<dpsd_orig<<" (km/s)/Mpc"<<endl;
+
+ if(type_beta<0) cout<<"***WARNING: type_beta<0 , no evolution used for beta !!!"<<endl;
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
@@ -1052 +1100 @@
        double z = interpinv(dz);   // interpolation par morceau
        double grw = (*growth_)(z);
-       double dpsd = interdpd(z);
+       double dpsd = 0;
+       if(type_beta>0) {
+         dpsd = cosmo_->H(z) * pow(cosmo_->Omatter(z),fixed_beta);
+       } else if(type_beta<0) {
+         dpsd = cosmo_->H(z) * fixed_beta;
+       } else {
+         dpsd = interdpd(z); // interpole -H(z)*(1+z)*D'(z)/D(z)
+       }
        int_8 ip = IndexR(i,j,l);
        // on remet le beta au bon z