Changeset 3770 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

May 7, 2010, 7:32:40 PM (15 years ago)

Author:

cmv

Message:

relecture des redshift-distorsion + calcul des spectres, cmv 07/05/2010

Location:

trunk/Cosmo/SimLSS

Files:

• Makefile (modified) (3 diffs)
• cmvginit3d.cc (modified) (2 diffs)
• cmvrvloscor.cc (added)
• genefluct3d.cc (modified) (14 diffs)
• genefluct3d.h (modified) (4 diffs)

trunk/Cosmo/SimLSS/Makefile

@@ -18 +18 @@
 PROGS = \
        $(EXE)cmvobserv3d $(EXE)cmvobserv3df \
-       $(EXE)cmvginit3d $(EXE)cmvginit3df \
+       $(EXE)cmvginit3d $(EXE)cmvginit3df $(EXE)cmvrvloscorf \
        $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec $(EXE)cmvdefsurv \
@@ -26 +26 @@
 PROGSOBJ = \
           $(OBJ)cmvobserv3d.o $(OBJ)cmvobserv3df.o \
-          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o \
+          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o $(OBJ)cmvrvloscorf.o \
           $(OBJ)cmvtuniv.o $(OBJ)cmvtransf.o $(OBJ)cmvtgrowth.o $(OBJ)cmvtstpk.o \
           $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
@@ -206 +206 @@
 
 ##############################################################################
+cmvrvloscorf: $(EXE)cmvrvloscorf
+        echo $@ " done"
+$(EXE)cmvrvloscorf: $(OBJ)cmvrvloscorf.o $(LIBR) $(LIBG4)
+        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvrvloscorf.o $(MYLIB4) 
+$(OBJ)cmvrvloscorf.o: cmvrvloscor.cc
+        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -DGEN3D_FLOAT -o $@ cmvrvloscor.cc
+
+##############################################################################
 cmvtintfun: $(EXE)cmvtintfun
         echo $@ " done"

trunk/Cosmo/SimLSS/cmvginit3d.cc

@@ -509 +509 @@
    //-----------------------------------------------------------------
    cout<<"\n--- Modifying cube for Transverse velocity"<<endl;
-   fluct3d.ToVelTrans();
+   fluct3d.ToVelLoS();
    fluct3d.NTupleCheck(posobs,string("ntpvgenf"),ntnent);
    PrtTim(">>>> End Modifying cube for Transverse velocity");
@@ -601 +601 @@
  //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
  } catch (PException& exc) {
-   cerr<<"cmvobserv3d.cc catched PException"<<exc.Msg()<<endl;
+   cerr<<"cmvginit3d.cc catched PException"<<exc.Msg()<<endl;
    return 77;
  } catch (std::exception& sex) {
-   cerr << "cmvobserv3d.cc std::exception :" 
+   cerr << "cmvginit3d.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    return 78;
  } catch (...) {
-   cerr << "cmvobserv3d.cc catched unknown (...) exception  " << endl; 
+   cerr << "cmvginit3d.cc catched unknown (...) exception  " << endl; 
    return 79;
  }

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -92 +92 @@
  kredsh_ref_ = 0.;
  dred_ref_ = -999.;
+ h_ref_ = -999.;
  loscom_ref_ = -999.;
  dtrc_ref_ = dlum_ref_ = dang_ref_ = -999.;
@@ -277 +278 @@
  dlum_ref_ = cosmo_->Dlum(redsh_ref_);
  dang_ref_ = cosmo_->Dang(redsh_ref_);
- h_ref = cosmo_->H(redsh_ref_);
+ h_ref_ = cosmo_->H(redsh_ref_);
  growth_ref_ = (*growth_)(redsh_ref_);
  dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
@@ -287 +288 @@
        <<", nuref="<<nu_ref_ <<" GHz"
        <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
-       <<"   H="<<h_ref<<" km/s/Mpc"
+       <<"   H="<<h_ref_<<" km/s/Mpc"
        <<", D="<<growth_ref_
        <<", dD/dz="<<dsdz_growth_ref_<<endl
@@ -445 +446 @@
    fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
    fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
+   fwrt.WriteKey("HREF",h_ref_," ref H(z)");
    fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
    fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
@@ -503 +505 @@
    R_.Info()["ZREF"] = (r_8)redsh_ref_;
    R_.Info()["KZREF"] = (r_8)kredsh_ref_;
+   R_.Info()["HREF"] = (r_8)h_ref_;
    R_.Info()["Growth"] = (r_8)Dref();
    R_.Info()["GrowthPk"] = (r_8)DrefPk();
@@ -691 +694 @@
  long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);  kx *= kx;
-   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<kx<<"\t pk="<<pk_at_z(kx)<<endl;
+   double kx = AbsKx(i);  kx *= kx;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);  kz *= kz;
+       double kz = AbsKz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -734 +737 @@
  long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);  kx *= kx;
-   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<kx<<"\t pk="<<pk_at_z(kx)<<endl;
+   double kx = AbsKx(i);  kx *= kx;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);  kz *= kz;
+       double kz = AbsKz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -877 +880 @@
 
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i) *Dx_/2;
+   double kx = AbsKx(i) *Dx_/2;
    double pk_x = pixelfilter(kx);
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j) *Dy_/2;
+     double ky = AbsKy(j) *Dy_/2;
      double pk_y = pixelfilter(ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l) *Dz_/2;
+       double kz = AbsKz(l) *Dz_/2;
        double pk_z =  pixelfilter(kz);
        T_(l,j,i) *= pk_x*pk_y*pk_z;
@@ -893 +896 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ToVelTrans(void)
-{
-  cout<<"-------------------- TO BE VERIFIED ToVelTrans"<<endl;
-
+void GeneFluct3D::ToVelLoS(void)
+{
  double zpk = compute_pk_redsh_ref_;
  double dpsd = -cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
- if(lp_>0) cout<<"--- ToVelTrans --- at z="<<zpk<<", D'/D="<<dpsd
+ if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", D'/D="<<dpsd
                <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
  check_array_alloc();
@@ -963 +964 @@
 // Apply Conformal derivative of Growth to real space for transverse velocity cube
 {
-  cout<<"-------------------- TO BE VERIFIED ApplyDerGrowthFactor"<<endl;
-
  if(lp_>0) cout<<"--- ApplyDerGrowthFactor: evol="<<type_evol<<endl;
  check_array_alloc();
@@ -1044 +1043 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);  kx *= kx;
+   double kx = AbsKx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i));
@@ -1074 +1073 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);  kx *= kx;
+   double kx = AbsKx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      double kt = sqrt(kx+ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double pk = MODULE2(T_(l,j,i));
        herr.Add(kt,kz,pk);
@@ -1118 +1117 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);
+   double kx = AbsKx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);
+     double ky = AbsKy(j);
      double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
      ky *= ky; fy *= fy;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
@@ -1163 +1162 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);
+   double kx = AbsKx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);
+     double ky = AbsKy(j);
      double fy = (pixcor) ? pixelfilter(ky*Dy_/2) : 1.;
      ky *= ky; fy *= fy;
      double kt = sqrt(kx+ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
        double fz = (pixcor) ? vfz(l): 1.;

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -22 +22 @@
 #define WITH_FFTW_THREAD
 
-//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE ("g++ -DGEN3D_FLOAT ...")
+//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE ("g++ -c -DGEN3D_FLOAT ...")
 #if defined(GEN3D_FLOAT)
 #define GEN3D_TYPE r_4
@@ -85 +85 @@
 
   // Return |K_i| module relative to pixel indices
-  inline r_8 Kx(long i) {long ii=(i>Nx_/2)? Nx_-i :i; return ii*Dkx_;}
-  inline r_8 Ky(long j) {long jj=(j>Ny_/2)? Ny_-j :j; return jj*Dky_;}
+  inline r_8 AbsKx(long i) {long ii=(i>Nx_/2)? Nx_-i :i; return ii*Dkx_;}
+  inline r_8 AbsKy(long j) {long jj=(j>Ny_/2)? Ny_-j :j; return jj*Dky_;}
+  inline r_8 AbsKz(long l) {return l*Dkz_;}
+  inline r_8 Kx(long i) {long ii=(i>Nx_/2)? i-Nx_ :i; return ii*Dkx_;}
+  inline r_8 Ky(long j) {long jj=(j>Ny_/2)? j-Ny_ :j; return jj*Dky_;}
   inline r_8 Kz(long l) {return l*Dkz_;}
 
@@ -109 +112 @@
   double DrefPk(void)
     {return (growth_!=NULL && compute_pk_redsh_ref_>=0.) ? (*growth_)(compute_pk_redsh_ref_): -999.;}
+  double Href(void) {return h_ref_;}
 
   void ComputeFourier0(PkSpectrumZ& pk_at_z);
   void ComputeFourier(PkSpectrumZ& pk_at_z);
   void FilterByPixel(void);
-  void ToVelTrans(void);
+  void ToVelLoS(void);
 
   void ComputeReal(void);
@@ -202 +206 @@
   double redsh_ref_,kredsh_ref_,dred_ref_;
   double compute_pk_redsh_ref_;
-  double h_ref, loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
+  double h_ref_, loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
   double growth_ref_, dsdz_growth_ref_;
   double nu_ref_, dnu_ref_ ;