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

Timestamp:

Jan 17, 2007, 6:44:04 PM (19 years ago)

Author:

cmv

Message:

chgt HProf->HistoErr + spectre 2D cmv 17/01/2007

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -14 +14 @@
 #include "srandgen.h"
 
+#include "fabtcolread.h"
+#include "fabtwriter.h"
+#include "fioarr.h"
+
+#include "arrctcast.h"
+
 #include "constcosmo.h"
 #include "integfunc.h"
@@ -25 +31 @@
 
 //-------------------------------------------------------
-GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz)
-  : T_(T) , pkz_(pkz)
+GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
+  : T_(T) , array_allocated_(false)
 {
  SetNThread();
- SetSize(1,-1,1,1.,-1.,1.);
 }
 
@@ -44 +49 @@
 void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
 {
- if(nx<=0 || dx<=0. ) {
+  setsize(nx,ny,nz,dx,dy,dz);
+  setalloc();
+  setpointers(false);
+}
+
+void GeneFluct3D::setsize(long nx,long ny,long nz,double dx,double dy,double dz)
+{
+ if(nx<=0 || dx<=0.) {
    cout<<"GeneFluct3D::SetSize_Error: bad value(s)"<<endl;
    throw ParmError("GeneFluct3D::SetSize_Error: bad value(s)");
  }
 
- Tcontent_ = 0;
-
- // Les taille des tableaux
+ // Les tailles des tableaux
  Nx_ = nx;
  Ny_ = ny;  if(Ny_ <= 0) Ny_ = Nx_;
  Nz_ = nz;  if(Nz_ <= 0) Nz_ = Nx_;
+ N_.resize(0); N_.push_back(Nx_); N_.push_back(Ny_); N_.push_back(Nz_);
  NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
  NCz_ =  Nz_/2 +1;
  NTz_ = 2*NCz_;
- SzK_[2] = Nx_; SzK_[1] = Ny_; SzK_[0] = NCz_; // a l'envers
 
  // le pas dans l'espace (Mpc)
@@ -64 +74 @@
  Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
  Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
+ D_.resize(0); D_.push_back(Dx_); D_.push_back(Dy_); D_.push_back(Dz_);
  dVol_ = Dx_*Dy_*Dz_;
  Vol_ = (Nx_*Dx_)*(Ny_*Dy_)*(Nz_*Dz_);
@@ -71 +82 @@
  Dky_ = 2.*M_PI/(Ny_*Dy_);
  Dkz_ = 2.*M_PI/(Nz_*Dz_);
+ Dk_.resize(0); Dk_.push_back(Dkx_); Dk_.push_back(Dky_); Dk_.push_back(Dkz_);
  Dk3_ = Dkx_*Dky_*Dkz_;
  
@@ -77 +89 @@
  Knyqy_ = M_PI/Dy_;
  Knyqz_ = M_PI/Dz_;
-
+ Knyq_.resize(0); Knyq_.push_back(Knyqx_); Knyq_.push_back(Knyqy_); Knyq_.push_back(Knyqz_);
+}
+
+void GeneFluct3D::setalloc(void)
+{
+ // Dimensionnement du tableau complex<r_8>
+ // ATTENTION: TArray adresse en memoire a l'envers du C
+ //            Tarray(n1,n2,n3) == Carray[n3][n2][n1]
+ sa_size_t SzK_[3] = {NCz_,Ny_,Nx_}; // a l'envers
+ try {
+   T_.ReSize(3,SzK_);
+   array_allocated_ = true;
+ } catch (...) {
+   cout<<"GeneFluct3D::SetSize_Error: Problem allocating T_"<<endl;
+ }
+ T_.SetMemoryMapping(BaseArray::CMemoryMapping);
+}
+
+void GeneFluct3D::setpointers(bool from_real)
+{
+ if(from_real) T_ = ArrCastR2C(R_);
+   else        R_ = ArrCastC2R(T_);
+ // On remplit les pointeurs
+ fdata_ = (fftw_complex *) (&T_(0,0,0));
+ data_ = (double *) (&R_(0,0,0));
+}
+
+void GeneFluct3D::check_array_alloc(void)
+// Pour tester si le tableau T_ est alloue
+{
+ if(array_allocated_) return;
+ char bla[90];
+ sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
+ cout<<bla<<endl;
+ throw ParmError(bla);
+}
+
+
+//-------------------------------------------------------
+void GeneFluct3D::WriteFits(string cfname,int bitpix)
+{
+ cout<<"GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
+ try {
+   FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
+   fwrt.WriteKey("NX",Nx_," axe transverse 1");
+   fwrt.WriteKey("NY",Ny_," axe transverse 2");
+   fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");
+   fwrt.WriteKey("DX",Dx_," Mpc");
+   fwrt.WriteKey("DY",Dy_," Mpc");
+   fwrt.WriteKey("DZ",Dz_," Mpc");
+   fwrt.WriteKey("DKX",Dkx_," Mpc^-1");
+   fwrt.WriteKey("DKY",Dky_," Mpc^-1");
+   fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
+   fwrt.Write(R_);
+ } catch (PThrowable & exc) {
+   cout<<"Exception : "<<(string)typeid(exc).name()
+       <<" - Msg= "<<exc.Msg()<<endl;
+   return;
+ } catch (...) {
+   cout<<" some other exception was caught !"<<endl;
+   return;
+ }
+}
+
+void GeneFluct3D::ReadFits(string cfname)
+{
+ cout<<"GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
+ try {
+   FitsImg3DRead fimg(cfname.c_str(),0,5);
+   fimg.Read(R_);
+   long nx = fimg.ReadKeyL("NX");
+   long ny = fimg.ReadKeyL("NY");
+   long nz = fimg.ReadKeyL("NZ");
+   double dx = fimg.ReadKey("DX");
+   double dy = fimg.ReadKey("DY");
+   double dz = fimg.ReadKey("DZ");
+   setsize(nx,ny,nz,dx,dy,dz);
+   setpointers(true);
+ } catch (PThrowable & exc) {
+   cout<<"Exception : "<<(string)typeid(exc).name()
+       <<" - Msg= "<<exc.Msg()<<endl;
+   return;
+ } catch (...) {
+   cout<<" some other exception was caught !"<<endl;
+   return;
+ }
+}
+
+void GeneFluct3D::WritePPF(string cfname,bool write_real)
+// On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> >
+{
+ cout<<"GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
+ try {
+   R_.Info()["NX"] = (int_8)Nx_;
+   R_.Info()["NY"] = (int_8)Ny_;
+   R_.Info()["NZ"] = (int_8)Nz_;
+   R_.Info()["DX"] = (r_8)Dx_;
+   R_.Info()["DY"] = (r_8)Dy_;
+   R_.Info()["DZ"] = (r_8)Dz_;
+   POutPersist pos(cfname.c_str());
+   if(write_real) pos << PPFNameTag("rgen")  << R_;
+     else         pos << PPFNameTag("pkgen") << T_;
+ } catch (PThrowable & exc) {
+   cout<<"Exception : "<<(string)typeid(exc).name()
+       <<" - Msg= "<<exc.Msg()<<endl;
+   return;
+ } catch (...) {
+   cout<<" some other exception was caught !"<<endl;
+   return;
+ }
+}
+
+void GeneFluct3D::ReadPPF(string cfname)
+{
+ cout<<"GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
+ try {
+   bool from_real = true;
+   PInPersist pis(cfname.c_str());
+   string name_tag_k = "pkgen";
+   bool found_tag_k = pis.GotoNameTag("pkgen");
+   if(found_tag_k) {
+     cout<<"           ...reading spectrun into TArray<complex <r_8> >"<<endl;
+     pis >> PPFNameTag("pkgen")  >> T_;
+     from_real = false;
+   } else {
+     cout<<"           ...reading space into TArray<r_8>"<<endl;
+     pis >> PPFNameTag("rgen")  >> R_;
+   }
+   int_8 nx = R_.Info()["NX"];
+   int_8 ny = R_.Info()["NY"];
+   int_8 nz = R_.Info()["NZ"];
+   r_8 dx = R_.Info()["DX"];
+   r_8 dy = R_.Info()["DY"];
+   r_8 dz = R_.Info()["DZ"];
+   setsize(nx,ny,nz,dx,dy,dz);
+   setpointers(from_real);
+ } catch (PThrowable & exc) {
+   cout<<"Exception : "<<(string)typeid(exc).name()
+       <<" - Msg= "<<exc.Msg()<<endl;
+   return;
+ } catch (...) {
+   cout<<" some other exception was caught !"<<endl;
+   return;
+ }
 }
 
@@ -83 +238 @@
 void GeneFluct3D::Print(void)
 {
- cout<<"GeneFluct3D(T_typ="<<Tcontent_<<"): z="<<pkz_.GetZ()<<endl;
+ cout<<"GeneFluct3D(T_alloc="<<array_allocated_<<"):"<<endl;
  cout<<"Space Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<Nz_<<" ("<<NTz_<<")  size="
      <<NRtot_<<endl;
@@ -98 +253 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier0(void)
+void GeneFluct3D::ComputeFourier0(GenericFunc& pk_at_z)
 // cf ComputeFourier() mais avec autre methode de realisation du spectre
 //    (attention on fait une fft pour realiser le spectre)
@@ -104 +259 @@
  int lp=2;
 
- T_.ReSize(3,SzK_);
- T_.SetMemoryMapping(BaseArray::CMemoryMapping);
-
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
  if(lp>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
- fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
- double *data = (double *) (&T_(0,0,0));
 #ifdef FFTW_THREAD
  if(nthread_>0) {
@@ -118 +268 @@
  }
 #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);
+ 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 ---");
 
@@ -127 +277 @@
  double sntot = sqrt((double)NRtot_);
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   data[ip] = NorRand()/sntot;
+   int_8 ip = IndexR(i,j,l);
+   data_[ip] = NorRand()/sntot;
  }
  if(lp>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
@@ -153 +303 @@
        double k = sqrt(kx+ky+kz);
        // cf normalisation: Peacock, Cosmology, formule 16.38 p504
-       double pk = pkz_(k)/Vol_;
+       double pk = pk_at_z(k)/Vol_;
        // ici pas de "/2" a cause de la remarque ci-dessus
        T_(l,j,i) *= sqrt(pk);
@@ -165 +315 @@
  if(lp>1) PrtTim("--- ComputeFourier0: after Computing power ---");
 
- Tcontent_ = 1;
-
 }
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier(void)
-// Calcule une realisation du spectre "pkz_"
+void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
+// Calcule une realisation du spectre "pk_at_z"
 // Attention: dans TArray le premier indice varie le + vite
 // Explication normalisation: see Coles & Lucchin, Cosmology, p264-265
@@ -181 +329 @@
 {
  int lp=2;
- 
- // --- Dimensionnement du tableau
- // ATTENTION: TArray adresse en memoire a l'envers du C
- //            Tarray(n1,n2,n3) == Carray[n3][n2][n1]
- T_.ReSize(3,SzK_);
- T_.SetMemoryMapping(BaseArray::CMemoryMapping);
 
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
  if(lp>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
- fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
- double *data = (double *) (&T_(0,0,0));
 #ifdef FFTW_THREAD
  if(nthread_>0) {
@@ -199 +339 @@
  }
 #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);
+ 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 ---");
@@ -222 +362 @@
        double k = sqrt(kx+ky+kz);
        // cf normalisation: Peacock, Cosmology, formule 16.38 p504
-       double pk = pkz_(k)/Vol_;
+       double pk = pk_at_z(k)/Vol_;
        // Explication de la division par 2: voir perandom.cc
        // ou egalement Coles & Lucchin, Cosmology formula 13.7.2 p279
@@ -238 +378 @@
  if(lp>1) PrtTim("--- ComputeFourier: after before Computing power ---");
 
- Tcontent_ = 1;
-
 }
 
@@ -246 +384 @@
 // because we want a realization of a real data in real space
 {
- fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
+ check_array_alloc();
 
  // 1./ Le Continu et Nyquist sont reels
@@ -259 +397 @@
        long i=0;
        if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
-       int_8 ip = k+NCz_*(j+Ny_*i);
-       //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata[ip][0]<<","<<fdata[ip][1]<<")"<<endl;
-       fdata[ip][1] = 0.; fdata[ip][0] *= M_SQRT2;
+       int_8 ip = IndexC(i,j,k);
+       //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata_[ip][0]<<","<<fdata_[ip][1]<<")"<<endl;
+       fdata_[ip][1] = 0.; fdata_[ip][0] *= M_SQRT2;
        nreal++;
      }
@@ -279 +417 @@
      if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
      for(long i=1;i<(Nx_+1)/2;i++) {
-       int_8 ip = k+NCz_*(j+Ny_*i);
-       int_8 ip1 = k+NCz_*(j+Ny_*(Nx_-i));
-       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
+       int_8 ip = IndexC(i,j,k);
+       int_8 ip1 = IndexC(Nx_-i,j,k);
+       fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
        nconj1++;
      }
@@ -289 +427 @@
      if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
      for(long j=1;j<(Ny_+1)/2;j++) {
-       int_8 ip = k+NCz_*(j+Ny_*i);
-       int_8 ip1 = k+NCz_*((Ny_-j)+Ny_*i);
-       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
+       int_8 ip = IndexC(i,j,k);
+       int_8 ip1 = IndexC(i,Ny_-j,k);
+       fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
        nconj1++;
      }
@@ -307 +445 @@
      for(long i=1;i<Nx_;i++) {
        if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
-       int_8 ip = k+NCz_*(j+Ny_*i);
-       int_8 ip1 = k+NCz_*((Ny_-j)+Ny_*(Nx_-i));
-       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
+       int_8 ip = IndexC(i,j,k);
+       int_8 ip1 = IndexC(Nx_-i,Ny_-j,k);
+       fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
        nconj2++;
      }
@@ -324 +462 @@
 // Calcul la puissance de la realisation du spectre Pk
 {
+ check_array_alloc();
+
  double s2 = 0.;
  for(long l=0;l<NCz_;l++)
@@ -351 +491 @@
 {
  int lp=2;
+ check_array_alloc();
+
  if(lp>1) PrtTim("--- FilterByPixel: before ---");
 
@@ -356 +498 @@
    long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_ *Dx_/2;
-   double pkx = pixelfilter(kx); 
+   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 pky =  pixelfilter(ky); 
+     double pk_y =  pixelfilter(ky); 
      for(long l=0;l<NCz_;l++) {
        double kz = l*Dkz_ *Dz_/2;
-       double pkz =  pixelfilter(kz);
-       T_(l,j,i) *= pkx*pky*pkz;
+       double pk_z =  pixelfilter(kz);
+       T_(l,j,i) *= pk_x*pk_y*pk_z;
      }
    }
@@ -377 +519 @@
 {
  int lp=2;
-
- if( Tcontent_==0 ) {
-   cout<<"GeneFluct3D::ComputeReal_Error: empty array"<<endl;
-   throw ParmError("GeneFluct3D::ComputeReal_Error: empty array");
- }
+ check_array_alloc();
 
  // On fait la FFT
@@ -387 +525 @@
  fftw_execute(pb_);
  if(lp>1) PrtTim("--- ComputeReal: after fftw backward ---");
-
- Tcontent_ = 2;
 }
 
@@ -395 +531 @@
 {
  int lp=2;
+ check_array_alloc();
 
  // On fait la FFT
@@ -405 +542 @@
      for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
 
- Tcontent_ = 3;
  if(lp>1) PrtTim("--- ComputeReal: after fftw forward ---");
 }
 
 //-------------------------------------------------------------------
-int GeneFluct3D::ComputeSpectrum(HProf& hp)
-// Compute spectrum from "T" and fill HProf "hp"
+int GeneFluct3D::ComputeSpectrum(HistoErr& herr)
+// Compute spectrum from "T" and fill HistoErr "herr"
 // T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
 // cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
 {
-
- if(hp.NBins()<0) return -1;
- hp.Zero();
- hp.SetErrOpt(false);
+ check_array_alloc();
+
+ if(herr.NBins()<0) return -1;
+ herr.Zero();
 
  // Attention a l'ordre
@@ -431 +567 @@
        double k = sqrt(kx+ky+kz);
        double pk = MODULE2(T_(l,j,i));
-       hp.Add(k,pk);
-     }
-   }
- }
- hp.UpdateHisto();
+       herr.Add(k,pk);
+     }
+   }
+ }
+ herr.ToCorrel();
 
  // renormalize to directly compare to original spectrum
  double norm = Vol_;
- hp *= norm;
+ herr *= norm;
+
+ return 0;
+}
+
+int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
+{
+ check_array_alloc();
+
+ if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
+ herr.Zero();
+
+ // 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;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;  ky *= ky;
+     double kt = sqrt(kx+ky);
+     for(long l=0;l<NCz_;l++) {
+       double kz = l*Dkz_;
+       double pk = MODULE2(T_(l,j,i));
+       herr.Add(kt,kz,pk);
+     }
+   }
+ }
+ herr.ToCorrel();
+
+ // renormalize to directly compare to original spectrum
+ double norm = Vol_;
+ herr *= norm;
 
  return 0;
@@ -449 +616 @@
 {
  int lp=2;
+ check_array_alloc();
+
  if(lp>1) PrtTim("--- VarianceFrReal: before computation ---");
 
- double *data = (double *) (&T_(0,0,0));
  long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
  long dny = long(R/Dy_+0.5); if(dny<=0) dny = 1;
@@ -470 +638 @@
              double z = (ll-l)*Dz_; z *= z;
              if(x+y+z>r2) continue;
-             int_8 ip = ll+NTz_*(jj+Ny_*ii);
-             s += 1.+data[ip];
+             int_8 ip = IndexR(ii,jj,ll);
+             s += 1.+data_[ip];
              n++;
            }
@@ -500 +668 @@
 //     ->  vmin and vmax are considered as bad
 {
- double *data = (double *) (&T_(0,0,0));
+ check_array_alloc();
 
  int_8 nbad = 0;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   double v = data[ip];
+   int_8 ip = IndexR(i,j,l);
+   double v = data_[ip];
    if(v<=vmin || v>=vmax) nbad++;
  }
@@ -516 +684 @@
 //   -> mean and sigma^2 are NOT computed with pixels values vmin and vmax
 {
- double *data = (double *) (&T_(0,0,0));
+ check_array_alloc();
 
  int_8 n = 0;
@@ -522 +690 @@
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   double v = data[ip];
+   int_8 ip = IndexR(i,j,l);
+   double v = data_[ip];
    if(v<=vmin || v>=vmax) continue;
    rm += v;
@@ -542 +710 @@
 //     ->  vmin and vmax are set to val0
 {
- double *data = (double *) (&T_(0,0,0));
+ check_array_alloc();
 
  int_8 nbad = 0;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   double v = data[ip];
-   if(v<=vmin || v>=vmax) {data[ip] = val0; nbad++;}
+   int_8 ip = IndexR(i,j,l);
+   double v = data_[ip];
+   if(v<=vmin || v>=vmax) {data_[ip] = val0; nbad++;}
  }
 
@@ -559 +727 @@
 {
  int lp=2;
+ check_array_alloc();
+
  if(lp>1) PrtTim("--- TurnFluct2Mass: before computation ---");
- double *data = (double *) (&T_(0,0,0));
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   data[ip] += 1.;
- }
-
- Tcontent_ = 4;
+   int_8 ip = IndexR(i,j,l);
+   data_[ip] += 1.;
+ }
+
  if(lp>1) PrtTim("--- TurnFluct2Mass: after computation ---");
 }
@@ -576 +744 @@
  int lp=2;
  if(lp>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
-
- double *data = (double *) (&T_(0,0,0));
 
  double m,s2;
@@ -596 +762 @@
  double sum = 0.;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   data[ip] *= dn; // par coherence on multiplie aussi les <=0
-   if(data[ip]>0.) sum += data[ip];  // mais on ne les compte pas
+   int_8 ip = IndexR(i,j,l);
+   data_[ip] *= dn; // par coherence on multiplie aussi les <=0
+   if(data_[ip]>0.) sum += data_[ip];  // mais on ne les compte pas
  }
  if(lp) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
 
- Tcontent_ = 6;
  if(lp>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
  return sum;
@@ -611 +776 @@
 {
  int lp=2;
+ check_array_alloc();
+
  if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
-
- double *data = (double *) (&T_(0,0,0));
 
  double sum = 0.;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   double v = data[ip];
+   int_8 ip = IndexR(i,j,l);
+   double v = data_[ip];
    if(v>0.) {
      unsigned long dn = PoissRandLimit(v,10.);
-     data[ip] = (double)dn;
+     data_[ip] = (double)dn;
      sum += (double)dn;
    }
  }
  if(lp) cout<<sum<<" galaxies put into survey"<<endl;
- Tcontent_ = 8;
 
  if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
@@ -641 +805 @@
 {
  int lp=2;
+ check_array_alloc();
+
  if(lp>1) PrtTim("--- TurnNGal2Mass: before computation ---");
-
-  double *data = (double *) (&T_(0,0,0));
 
  double sum = 0.;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   double v = data[ip];
+   int_8 ip = IndexR(i,j,l);
+   double v = data_[ip];
    if(v>0.) {
      long ngal = long(v+0.1);
-     data[ip] = 0.;
+     data_[ip] = 0.;
      for(long i=0;i<ngal;i++) {
        double m = massdist.RandomInterp();  // massdist.Random();
        if(axeslog) m = pow(10.,m);
-       data[ip] += m;
-     }
-     sum += data[ip];
+       data_[ip] += m;
+     }
+     sum += data_[ip];
    }
  }
  if(lp) cout<<sum<<" MSol HI mass put into survey"<<endl;
- Tcontent_ = 10;
 
  if(lp>1) PrtTim("--- TurnNGal2Mass: after computation ---");
@@ -671 +834 @@
 {
  int lp=2;
+ check_array_alloc();
+
  if(lp>1) PrtTim("--- AddNoise2Real: before computation ---");
 
- double *data = (double *) (&T_(0,0,0));
-
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = l+NTz_*(j+Ny_*i);
-   data[ip] += snoise*NorRand();
- }
- Tcontent_ = 12;
+   int_8 ip = IndexR(i,j,l);
+   data_[ip] += snoise*NorRand();
+ }
 
  if(lp>1) PrtTim("--- AddNoise2Real: after computation ---");