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

Timestamp:

Jan 19, 2007, 5:19:49 PM (19 years ago)

Author:

cmv

Message:

modifs d'options de cmvobserv3d.cc
creation de init_fftw() et deplacement a l'initialisation
creation de SetObservator pour positionner l'observateur + ecriture FITS et PPF

(mais le codage des distance/redshift reste a coder)

cmv 19/01/2007

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -32 +32 @@
 //-------------------------------------------------------
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
-  : T_(T) , array_allocated_(false) , lp_(0)
+  : T_(T) , Nx_(0) , Ny_(0) , Nz_(0) , array_allocated_(false) , lp_(0)
+  , redshref_(-999.) , kredshref_(-999.)
 {
  SetNThread();
@@ -52 +53 @@
   setalloc();
   setpointers(false);
+  init_fftw();
+}
+
+void GeneFluct3D::SetObservator(double redshref,double kredshref)
+// L'observateur est au redshift z=0
+//               est situe sur la "perpendiculaire" a la face x,y
+//                         issue du centre de cette face
+// Il faut positionner le cube sur l'axe des z cad des redshifts:
+//     redshref  = redshift de reference
+//                 Si redshref<0 alors redshref=0
+//     kredshref = indice (en double) correspondant a ce redshift
+//                 Si kredshref<0 alors kredshref=0
+{
+ if(redshref<0.) redshref = 0.;
+ if(kredshref<0.) kredshref = 0.;
+ redshref_ = redshref;
+ kredshref_ = kredshref;
 }
 
@@ -126 +144 @@
 }
 
+void GeneFluct3D::init_fftw(void)
+{
+ // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
+ if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: before fftw_plan ---");
+#ifdef FFTW_THREAD
+ if(nthread_>0) {
+   cout<<"Computing with "<<nthread_<<" threads"<<endl;
+   fftw_init_threads();
+   fftw_plan_with_nthreads(nthread_);
+ }
+#endif
+ pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
+ pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
+ if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: after fftw_plan ---");
+}
+
 
 //-------------------------------------------------------
@@ -142 +176 @@
    fwrt.WriteKey("DKY",Dky_," Mpc^-1");
    fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
+   fwrt.WriteKey("ZREF",redshref_," reference redshift");
+   fwrt.WriteKey("KZREF",kredshref_," reference redshift on z axe");
    fwrt.Write(R_);
  } catch (PThrowable & exc) {
@@ -165 +201 @@
    double dy = fimg.ReadKey("DY");
    double dz = fimg.ReadKey("DZ");
+   double zref = fimg.ReadKey("ZREF");
+   double kzref = fimg.ReadKey("KZREF");
    setsize(nx,ny,nz,dx,dy,dz);
    setpointers(true);
+   init_fftw();
+   SetObservator(zref,kzref);
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -188 +228 @@
    R_.Info()["DY"] = (r_8)Dy_;
    R_.Info()["DZ"] = (r_8)Dz_;
+   R_.Info()["ZREF"] = (r_8)redshref_;
+   R_.Info()["KZREF"] = (r_8)kredshref_;
    POutPersist pos(cfname.c_str());
    if(write_real) pos << PPFNameTag("rgen")  << R_;
@@ -217 +259 @@
      pis >> PPFNameTag("rgen")  >> R_;
    }
+   setpointers(from_real);  // a mettre ici pour relire les DVInfo
    int_8 nx = R_.Info()["NX"];
    int_8 ny = R_.Info()["NY"];
@@ -223 +266 @@
    r_8 dy = R_.Info()["DY"];
    r_8 dz = R_.Info()["DZ"];
+   r_8 zref = R_.Info()["ZREF"];
+   r_8 kzref = R_.Info()["KZREF"];
    setsize(nx,ny,nz,dx,dy,dz);
-   setpointers(from_real);
+   init_fftw();
+   SetObservator(zref,kzref);
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -250 +296 @@
      <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
  cout<<"          (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
+ cout<<"Redshift "<<redshref_<<" for z axe at k="<<kredshref_<<endl;
 }
 
@@ -257 +304 @@
 //    (attention on fait une fft pour realiser le spectre)
 {
- // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp_>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
-#ifdef FFTW_THREAD
- if(nthread_>0) {
-   cout<<"Computing with "<<nthread_<<" threads"<<endl;
-   fftw_init_threads();
-   fftw_plan_with_nthreads(nthread_);
- }
-#endif
- pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
- pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
- if(lp_>1) PrtTim("--- ComputeFourier0: after fftw_plan ---");
 
  // --- realisation d'un tableau de tirage gaussiens
@@ -326 +361 @@
 //                                          sum(fb(x_i)^2) = S*N^2
 {
- // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp_>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
-#ifdef FFTW_THREAD
- if(nthread_>0) {
-   cout<<"Computing with "<<nthread_<<" threads"<<endl;
-   fftw_init_threads();
-   fftw_plan_with_nthreads(nthread_);
- }
-#endif
- pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
- pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
- //fftw_print_plan(pb_); cout<<endl;
- if(lp_>1) PrtTim("--- ComputeFourier: after fftw_plan ---");
-
  // --- RaZ du tableau
  T_ = complex<r_8>(0.);