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

Timestamp:

May 3, 2010, 4:08:34 PM (15 years ago)

Author:

cmv

Message:

• refonte du code pour creer uniquement des conditions initiales
• introduction du tirage des vitesse LOS pour les redshift-distortion

cmv 03/05/2010

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -1 @@
-#include "sopnamsp.h"
 #include "machdefs.h"
 #include <iostream>
@@ -88 +87 @@
  cosmo_ = NULL;
  growth_ = NULL;
+ good_dzinc_ = 0.01;
+ compute_pk_redsh_ref_ = -999.;
  redsh_ref_ = -999.;
  kredsh_ref_ = 0.;
@@ -262 +263 @@
 // 
 {
+ if(zinc<=0.) zinc = good_dzinc_;
  if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
 
- if(cosmo_ == NULL || redsh_ref_<0.) {
-   const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
+ if(cosmo_ == NULL || growth_==NULL || redsh_ref_<0.) {
+   const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator, Cosmology and Growth first";
    cout<<bla<<endl; throw ParmError(bla);
  }
@@ -275 +277 @@
  dlum_ref_ = cosmo_->Dlum(redsh_ref_);
  dang_ref_ = cosmo_->Dang(redsh_ref_);
+ h_ref = cosmo_->H(redsh_ref_);
+ growth_ref_ = (*growth_)(redsh_ref_);
+ dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
  nu_ref_   = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
  dnu_ref_  = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
@@ -282 +287 @@
        <<", nuref="<<nu_ref_ <<" GHz"
        <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
+       <<"   H="<<h_ref<<" km/s/Mpc"
+       <<", D="<<growth_ref_
+       <<", dD/dz="<<dsdz_growth_ref_<<endl
        <<"   dlosc="<<loscom_ref_<<" Mpc com"
        <<", dtrc="<<dtrc_ref_<<" Mpc com"
@@ -335 +343 @@
 
  // Fill the corresponding vectors for loscom and zred
- // Be shure to have one dlc <loscom_min and one >loscom_max
- if(zinc<=0.) zinc = 0.01;
+ // Be shure to have one dlc<loscom_min and one dlc>loscom_max
  double zmin = 0., dlcmin=0.;
  while(1) {
@@ -397 +404 @@
  }
 
+
+ // Compute the table for D'(z)/D(z) where D'(z)=dD/dEta (Eta is conformal time)
+ zredmin_dpsd_ = zred_[0];
+ zredmax_dpsd_ = zred_[zred_.size()-1];
+ long nptd = long(sqrt(Nx_*Nx_ + Ny_*Ny_ +Nz_*Nz_));
+ if(nptd<10) nptd = 10;
+ good_dzinc_ = (zredmax_dpsd_ - zredmin_dpsd_)/nptd;
+ if(lp_>0) cout<<"...good_dzinc changed to "<<good_dzinc_<<endl;
+ dpsdfrzred_.resize(0);
+ double dz = (zredmax_dpsd_ - zredmin_dpsd_)/(nptd-1);
+ int nmod = nptd/5; if(nmod==0) nmod = 1;
+ if(lp_>0) cout<<"...Compute table D'/D on "<<nptd<<" pts for z["
+               <<zredmin_dpsd_<<","<<zredmax_dpsd_<<"]"<<endl;
+ 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);
+   dpsdfrzred_.push_back(v);
+   if(lp_ && (i%nmod==0 || i==nptd-1)) cout<<"    z="<<z<<"  D'/D="<<v<<" km/s/Mpc"<<endl;
+ }
+
  return zred_.size();
 }
@@ -415 +442 @@
    fwrt.WriteKey("DKY",Dky_," Mpc^-1");
    fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
+   fwrt.WriteKey("ZREFPK",compute_pk_redsh_ref_," Pk computed redshift");
    fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
    fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
+   fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
+   fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
    fwrt.Write(R_);
  } catch (PThrowable & exc) {
@@ -440 +470 @@
    double dy = fimg.ReadKey("DY");
    double dz = fimg.ReadKey("DZ");
+   double pkzref = fimg.ReadKey("ZREFPK");
    double zref = fimg.ReadKey("ZREF");
    double kzref = fimg.ReadKey("KZREF");
@@ -447 +478 @@
    SetObservator(zref,kzref);
    array_allocated_ = true;
+   compute_pk_redsh_ref_ = pkzref;
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -468 +500 @@
    R_.Info()["DY"] = (r_8)Dy_;
    R_.Info()["DZ"] = (r_8)Dz_;
+   R_.Info()["ZREFPK"] = (r_8)compute_pk_redsh_ref_;
    R_.Info()["ZREF"] = (r_8)redsh_ref_;
    R_.Info()["KZREF"] = (r_8)kredsh_ref_;
+   R_.Info()["Growth"] = (r_8)Dref();
+   R_.Info()["GrowthPk"] = (r_8)DrefPk();
    POutPersist pos(cfname.c_str());
    if(write_real) pos << PPFNameTag("rgen")  << R_;
@@ -506 +541 @@
    r_8 dy = R_.Info()["DY"];
    r_8 dz = R_.Info()["DZ"];
+   r_8 pkzref = R_.Info()["ZREFPK"];
    r_8 zref = R_.Info()["ZREF"];
    r_8 kzref = R_.Info()["KZREF"];
@@ -512 +548 @@
    SetObservator(zref,kzref);
    array_allocated_ = true;
+   compute_pk_redsh_ref_ = pkzref;
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -630 +667 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier0(GenericFunc& pk_at_z)
+void GeneFluct3D::ComputeFourier0(PkSpectrumZ& pk_at_z)
 // cf ComputeFourier() mais avec autre methode de realisation du spectre
 //    (attention on fait une fft pour realiser le spectre)
 {
-
+  compute_pk_redsh_ref_ = pk_at_z.GetZ();
  // --- realisation d'un tableau de tirage gaussiens
- if(lp_>0) cout<<"--- ComputeFourier0: before gaussian filling ---"<<endl;
+  if(lp_>0) cout<<"--- ComputeFourier0 at z="<<compute_pk_redsh_ref_
+                <<": before gaussian filling ---"<<endl;
  // On tient compte du pb de normalisation de FFTW3
  double sntot = sqrt((double)NRtot_);
@@ -651 +689 @@
  if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
  T_(0,0,0) = complex<GEN3D_TYPE>(0.);  // on coupe le continue et on l'initialise
- long lmod = Nx_/10; if(lmod<1) lmod=1;
+ long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
-   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   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;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;  kz *= kz;
+       double kz = Kz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -680 +716 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
+void GeneFluct3D::ComputeFourier(PkSpectrumZ& pk_at_z)
 // Calcule une realisation du spectre "pk_at_z"
 // Attention: dans TArray le premier indice varie le + vite
@@ -692 +728 @@
  // --- RaZ du tableau
  T_ = complex<GEN3D_TYPE>(0.);
+ compute_pk_redsh_ref_ = pk_at_z.GetZ();
 
  // --- On remplit avec une realisation
- if(lp_>0) cout<<"--- ComputeFourier ---"<<endl;
- long lmod = Nx_/10; if(lmod<1) lmod=1;
+ if(lp_>0) cout<<"--- ComputeFourier at z="<<compute_pk_redsh_ref_<<" ---"<<endl;
+ long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
-   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   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;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;  kz *= kz;
+       double kz = Kz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -728 +763 @@
 // Take into account the real and complexe conjugate coefficients
 // because we want a realization of a real data in real space
+// On ecrit que: P(k_x,k_y,k_z) = P(-k_x,-k_y,-k_z)
+//               avec k_x = i, -k_x = N_x - i  etc...
 {
  if(lp_>1) cout<<"...managing coefficients"<<endl;
@@ -800 +837 @@
  if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
 
- if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
+ if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-nreal<<endl;
 
  return nreal+nconj1+nconj2;
@@ -840 +877 @@
 
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_ *Dx_/2;
+   double kx = Kx(i) *Dx_/2;
    double pk_x = pixelfilter(kx);
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_ *Dy_/2;
+     double ky = Ky(j) *Dy_/2;
      double pk_y = pixelfilter(ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_ *Dz_/2;
+       double kz = Kz(l) *Dz_/2;
        double pk_z =  pixelfilter(kz);
        T_(l,j,i) *= pk_x*pk_y*pk_z;
@@ -855 +890 @@
  }
 
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::ToVelTrans(void)
+{
+  cout<<"-------------------- TO BE VERIFIED ToVelTrans"<<endl;
+
+ 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
+               <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
+ check_array_alloc();
+
+ for(long i=0;i<Nx_;i++) {
+   double kx = Kx(i);
+   for(long j=0;j<Ny_;j++) {
+     double ky = Ky(j);
+     double kt2 = kx*kx + ky*ky;
+     for(long l=0;l<NCz_;l++) {
+       double kz = Kz(l);
+       double k2 = kt2 + kz*kz;
+       if(k2<=0.) continue;
+       T_(l,j,i) *= complex<double>(0.,dpsd*kz/k2);
+     }
+   }
+ }
 }
 
@@ -879 +940 @@
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
- //unsigned short ok;
-
- //CHECK: Histo hgr(0.9*zred_[0],1.1*zred_[n-1],1000);
+
  for(long i=0;i<Nx_;i++) {
    double dx2 = DXcom(i); dx2 *= dx2;
@@ -890 +949 @@
        if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
          else dz = fabs(dz); // tous les plans Z au meme redshift
-       double z = interpinv(dz);
-       //CHECK: hgr.Add(z);
-       double dzgr = (*growth_)(z);   // interpolation par morceau
-       //double dzgr = growth_->Linear(z,ok);  // interpolation lineaire
-       //double dzgr = growth_->Parab(z,ok);  // interpolation parabolique
+       double z = interpinv(dz);   // interpolation par morceau
+       double dzgr = (*growth_)(z);
        int_8 ip = IndexR(i,j,l);
        data_[ip] *= dzgr;
@@ -901 +957 @@
  }
 
- //CHECK: {POutPersist pos("applygrowth.ppf"); string tag="hgr"; pos.PutObject(hgr,tag);}
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::ApplyDerGrowthFactor(int type_evol)
+// 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();
+
+ if(growth_ == NULL) {
+   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+ if(type_evol<1 || type_evol>2) {
+   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ double zpk = compute_pk_redsh_ref_;
+ double dpsd_orig = - cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+
+ InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
+ InterpFunc interdpd(zredmin_dpsd_,zredmax_dpsd_,dpsdfrzred_);
+
+ for(long i=0;i<Nx_;i++) {
+   double dx2 = DXcom(i); dx2 *= dx2;
+   for(long j=0;j<Ny_;j++) {
+     double dy2 = DYcom(j); dy2 *= dy2;
+     for(long l=0;l<Nz_;l++) {
+       double dz = DZcom(l);
+       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
+         else dz = fabs(dz); // tous les plans Z au meme redshift
+       double z = interpinv(dz);   // interpolation par morceau
+       double dpsd = interdpd(z);
+       int_8 ip = IndexR(i,j,l);
+       data_[ip] *= dpsd / dpsd_orig;
+     }
+   }
+ }
 
 }
@@ -948 +1044 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
+   double kx = Kx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i));
@@ -980 +1074 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
+   double kx = Kx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      double kt = sqrt(kx+ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double pk = MODULE2(T_(l,j,i));
        herr.Add(kt,kz,pk);
@@ -1026 +1118 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;
+   double kx = Kx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;
+     double ky = Ky(j);
      double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
      ky *= ky; fy *= fy;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
@@ -1073 +1163 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;
+   double kx = Kx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;
+     double ky = Ky(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 = l*Dkz_;
+       double kz = Kz(l);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
        double fz = (pixcor) ? vfz(l): 1.;