Changeset 3799 in Sophya

Timestamp:

Jun 30, 2010, 6:23:54 PM (15 years ago)

Author:

cmv

Message:

Vlos * by (1+z), Dlos(com) do not * by (1+z), cmv 30/06/2010

Location:

trunk/Cosmo/SimLSS

Files:

• cmvrvloscor.cc (modified) (7 diffs)
• genefluct3d.cc (modified) (3 diffs)

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -19 +19 @@
 #include "genefluct3d.h"
 // set simul = 6_0p0_780
-// cmvrvloscorf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
-// cmvrvloscorf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
-// cmvrvloscorf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -o rvlos_${simul}.ppf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -o rvlos_${simul}.ppf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
 
 void usage(void);
@@ -34 +34 @@
 <<"-S: compute cross-power spectrum of V*conj(R) (def: no)"<<endl
 <<"-N: do not create 3D cube and recompute 1D and 2D spectra (def: no do-it !)"<<endl
+<<"-p: do not distribute dRho/Rho in pixel, just take the nearest (def: no, distribute rho)"<<endl
 <<"-2: compute 2D projection fpr dRho and dRho(corrected) (def=no)"<<endl
+<<"-o: output ppf file name (def=\"cmvrvloscor.ppf\")"<<endl
 <<endl;
 }
@@ -41 +43 @@
 {
  int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
- bool docube=true, dopk = false, do2d = false;
+ bool docube=true, dopk = false, do2d = false, dodistrib=true;
+ string fnppf = "cmvrvloscor.ppf";
  
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"hn:K:SN2")) != -1) {
+ while((c = getopt(narg,arg,"hn:K:SNp2o:")) != -1) {
   switch (c) {
   case 'n' :
@@ -63 +66 @@
   case '2' :
     do2d = true;
+    break;
+  case 'p' :
+    dodistrib = false;
+    break;
+  case 'o' :
+    fnppf = optarg;
     break;
   case 'h' :
@@ -100 +109 @@
  Histo hmpc(-dmin*nmax,dmin*nmax,4*nmax);
 
- POutPersist pos("cmvrvloscor.ppf");
+ POutPersist pos(fnppf.c_str());
  DVList dvlcor;
 
@@ -173 +182 @@
    // Calcul du champ R redshift distordu
    for(int l=0;l<Nz;l++) {
-     double d = (1.+Zref) / Href * V(l);
+     double d = V(l) / Href;
      if(fhis) hmpc.Add(d);
      double lpd = (double)l + d/Dz; // valeur du deplacee
-     // on repartit proportionellement au recouvrement sur 2 pixels
-     long l1 = long(lpd); // pixel de gauche
-     long l2 = l1 + 1;  // pixel de droite
-     lpd -= (double)l1;  // recouvrement du pixel du dessus
-     if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
-     if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
+     if(dodistrib) { 
+       // on repartit proportionellement au recouvrement sur 2 pixels
+       long l1 = long(lpd); // pixel de gauche
+       long l2 = l1 + 1;  // pixel de droite
+       lpd -= (double)l1;  // recouvrement du pixel du dessus
+       if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
+       if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
+     } else { // take nearest pixel
+       long l1 = long(lpd + 0.5);
+       if(l1>=0 && l1<Nz) Rdis(l1) += R(l);
+     }
    }
    // On remplit eventuellement les matrices 2D
@@ -351 +365 @@
 n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
 
+defscript tom
+del m2
+c++exec TMatrix<r_8> m2(hpkrec2.NBinX(),hpkrec2.NBinY()); \
+for(int i=0;i<hpkrec2.NBinX();i++) \
+for(int j=0;j<hpkrec2.NBinY();j++) m2(i,j) = hpkrec2(i,j); \
+KeepObj(m2); cout<<"OK"<<endl;
+endscript
+
 # les matrices 2D
 set n 0

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -421 +421 @@
  for(int i=0;i<nptd;i++) {
    double z = zredmin_dpsd_ + i*dz;
-   double v = -cosmo_->H(z) * growth_->DsDz(z,good_dzinc_) / (*growth_)(z);
+   double om = cosmo_->Omatter(z);
+   double h = cosmo_->H(z);
+   double dsdz = growth_->DsDz(z,good_dzinc_);
+   double d = (*growth_)(z);
+   double v = -h * (1.+z) * dsdz / d;
    dpsdfrzred_.push_back(v);
-   if(lp_ && (i%nmod==0 || i==nptd-1)) cout<<"    z="<<z<<"  D'/D="<<v<<" km/s/Mpc"<<endl;
+   if(lp_ && (i%nmod==0 || i==nptd-1))
+     cout<<"    z="<<z<<" H="<<h<<" Beta="<<-(1.+z)*dsdz/d
+         <<" (Om^0.6="<<pow(om,0.6)
+         <<")  -H*(1+z)*D'/D="<<v<<" km/s/Mpc"<<endl;
  }
 
@@ -901 +908 @@
 {
  double zpk = compute_pk_redsh_ref_;
- double dpsd = -cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
- if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", D'/D="<<dpsd
+ 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;
  check_array_alloc();
@@ -980 +987 @@
 
  double zpk = compute_pk_redsh_ref_;
- double dpsd_orig = - cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+ double dpsd_orig = - cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);