Changeset 3154 in Sophya for trunk/Cosmo/SimLSS

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

Location:

trunk/Cosmo/SimLSS

Files:

• cmvconcherr.cc (modified) (3 diffs)
• cmvobserv3d.cc (modified) (15 diffs)
• genefluct3d.cc (modified) (11 diffs)
• genefluct3d.h (modified) (5 diffs)

trunk/Cosmo/SimLSS/cmvconcherr.cc

@@ -119 +119 @@
    nread++;
    for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
-   sum+=herr(itest); sum2+=herr(itest)*herr(itest), nsum++;
+   sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;
  }
  herrconc->ToVariance();
@@ -182 +182 @@
    nread++;
    for(int i=0;i<nherrx;i++)
-     for(int j=0;j<nherrx;j++) herrconc->AddBin(i,j,herr(i,j));
-   sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty), nsum++;
+     for(int j=0;j<nherry;j++) herrconc->AddBin(i,j,herr(i,j));
+   sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty); nsum++;
  }
  herrconc->ToVariance();
+ if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
  cout<<"Test bin "<<itestx<<","<<itesty
      <<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
@@ -214 +215 @@
 openppf cmvconcherr.ppf
 
+zone 2 2
 disp herrconc "hbincont err"
 disp herrconc "hbinerr"
 disp herrconc "hbinent"
 
+zone
 n/plot herrconc.val%log10(x) x>0 ! "connectpoints"
-n/plot herrconc.sqrt(err2)%log10(x) x>0&&err2>0. ! "connectpoints same red"
-
-n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0.&&val>0 ! "connectpoints"
+n/plot herrconc.sqrt(err2)%log10(x) x>0&&err2>0 ! "connectpoints same red"
+
+n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0&&val>0 ! "connectpoints"
 
 #### Histo2DErr 2D
 openppf cmvconcherr2.ppf
 
-disp herrconc2 "hbincont"
-disp herrconc2 "hbinerr"
-disp herrconc2 "hbinent"
+zone 2 2
+imag herrconc2 "hbincont"
+imag herrconc2 "hbinerr"
+imag herrconc2 "hbinent"
 */

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -28 +28 @@
      <<" -Z zref : redshift median"<<endl
      <<" -s snoise : sigma du bruit en Msol"<<endl
+     <<" -2 : compute 2D spectrum"<<endl
+     <<" -F : write cube in FITS format"<<endl
+     <<" -P : write cube in PPF format"<<endl
+     <<" -V : compute variance from real space"<<endl
      <<endl;
 }
@@ -73 +77 @@
  unsigned short nthread=4;
 
+ // *** What to do
+ bool comp2dspec = false;
+ bool wfits = false;
+ bool wppf = false;
+ bool compvarreal = false;
+
  // --- Decodage arguments
 
  char c;
- while((c = getopt(narg,arg,"ha0x:y:z:s:Z:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2x:y:z:s:Z:")) != -1) {
   switch (c) {
   case 'a' :
@@ -98 +108 @@
   case 'Z' :
     sscanf(optarg,"%lf",&zref);
+    break;
+  case '2' :
+    comp2dspec = true;
+    break;
+  case 'V' :
+    compvarreal = true;
+    break;
+  case 'W' :
+    wfits = true;
+    break;
+  case 'P' :
+    wppf = true;
     break;
   case 'h' :
@@ -185 +207 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Compute galaxies density number"<<endl;
+ cout<<endl<<"\n--- Compute galaxy density number"<<endl;
 
  sch.SetOutValue(0);
@@ -249 +271 @@
  }
 
- if(1) {
+ if(comp2dspec) {
    cout<<"\n--- Checking realization 2D spectra"<<endl;
    Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
@@ -265 +287 @@
  }
 
- if(1) fluct3d.WriteFits("!cmvobserv3d_k.fits");
- if(1) fluct3d.WritePPF("cmvobserv3d_k.ppf",false);
+ if(wfits) fluct3d.WriteFits("!cmvobserv3d_k0.fits");
+ if(wppf) fluct3d.WritePPF("cmvobserv3d_k0.ppf",false);
 
  if(1) {
@@ -279 +301 @@
  }
 
- if(1) {
+ if(comp2dspec) {
    cout<<"\n--- Checking realization 2D spectra after pixel shape convol."<<endl;
    Histo2DErr hpkgenf2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
@@ -295 +317 @@
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
 
- if(1) fluct3d.WriteFits("!cmvobserv3d_r0.fits");
- if(1) fluct3d.WritePPF("cmvobserv3d_r0.ppf",true);
+ if(wfits) fluct3d.WriteFits("!cmvobserv3d_r0.fits");
+ if(wppf) fluct3d.WritePPF("cmvobserv3d_r0.ppf",true);
 
  int_8 nm;
@@ -308 +330 @@
  }
 
- if(1) {
+ if(compvarreal) {
    cout<<"\n--- Check variance sigmaR in real space"<<endl;
    double varr;
@@ -336 +358 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- cout<<"\n--- Apply poisson on galaxie number"<<endl;
+ cout<<"\n--- Apply poisson on galaxy number"<<endl;
  nm = fluct3d.ApplyPoisson();
  cout<<nm<<" galaxies into survey after poisson"<<endl;
@@ -343 +365 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- cout<<"\n--- Convert Galaxies number to HI mass"<<endl;
+ cout<<"\n--- Convert Galaxy number to HI mass"<<endl;
  long nhmdndm = long( (lschmax-lschmin+1.)*100. + 0.5);
  Histo hmdndm(lschmin,lschmax,nhmdndm);
@@ -368 +390 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- if(1) fluct3d.WriteFits("!cmvobserv3d_r.fits");
- if(1) fluct3d.WritePPF("cmvobserv3d_r.ppf",true);
+ if(wfits) fluct3d.WriteFits("!cmvobserv3d_r.fits");
+ if(wppf) fluct3d.WritePPF("cmvobserv3d_r.ppf",true);
 
  cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
@@ -384 +406 @@
    fluct3d.ReComputeFourier();
  }
+
+ if(wfits) fluct3d.WriteFits("!cmvobserv3d_k.fits");
+ if(wppf) fluct3d.WritePPF("cmvobserv3d_k.ppf",false);
 
  if(1) {
@@ -394 +419 @@
  }
 
- if(1) {
+ if(comp2dspec) {
    cout<<"\n--- Computing final 2D spectrum"<<endl;
    Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
@@ -410 +435 @@
 /*
 ######################################################
+readfits cmvobserv3d_k0.fits
 readfits cmvobserv3d_k.fits
 readfits cmvobserv3d_r0.fits
 readfits cmvobserv3d_r.fits
 
+openppf cmvobserv3d_k0.ppf
 openppf cmvobserv3d_k.ppf
 openppf cmvobserv3d_r0.ppf

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.);

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -27 +27 @@
   void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
   void SetSize(long nx,long ny,long nz,double dx,double dy,double dz);  // Mpc
+  void SetObservator(double redshref=0.,double kredshref=-1.);
 
   TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
@@ -92 +93 @@
   void setalloc(void);
   void setpointers(bool from_real);
+  void init_fftw(void);
   long manage_coefficients(void);
   double compute_power_carte(void);
@@ -98 +100 @@
     {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
 
+  // valeurs dans l'espace reel
   long Nx_,Ny_,Nz_;  vector<long> N_;
   long NCz_,NTz_;
@@ -104 +107 @@
   double Dx_,Dy_,Dz_;  vector<double> D_;
 
+  // valeurs dans l'espace des K
   double Dkx_,Dky_,Dkz_;  vector<double> Dk_;
   double Knyqx_,Knyqy_,Knyqz_;  vector<double> Knyq_;
@@ -109 +113 @@
   double dVol_, Vol_;
 
+  // la gestion de la FFT
   fftw_plan pf_,pb_;
   unsigned short nthread_;
   int lp_;
 
+  // le stockage du Cube de donnees et les pointeurs
   bool array_allocated_;  // true if array has been allocated
-
   TArray< complex<r_8> >& T_;
   fftw_complex *fdata_;
   TArray<r_8> R_;
   double *data_;
+
+  // l'observateur
+  double redshref_,kredshref_;
 };