Changeset 3141 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

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

Author:

cmv

Message:

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

Location:

trunk/Cosmo/SimLSS

Files:

• Makefile (modified) (3 diffs)
• arrctcast.h (modified) (4 diffs)
• cmvconcherr.cc (added)
• cmvconchprof.cc (deleted)
• cmvobserv3d.cc (modified) (8 diffs)
• genefluct3d.cc (modified) (42 diffs)
• genefluct3d.h (modified) (5 diffs)
• genericfunc.h (modified) (2 diffs)

trunk/Cosmo/SimLSS/Makefile

@@ -19 +19 @@
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
        $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
-       $(EXE)cmvtpoisson $(EXE)cmvconchprof
+       $(EXE)cmvtpoisson $(EXE)cmvconcherr
 #$(EXE)cmvtluc
  
@@ -26 +26 @@
           $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
           $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o \
-          $(OBJ)cmvtpoisson.o $(OBJ)cmvconchprof.o $(OBJ)cmvtluc.o
+          $(OBJ)cmvtpoisson.o $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o
  
 LIBROBJ = \
@@ -148 +148 @@
 
 ##############################################################################
-cmvconchprof: $(EXE)cmvconchprof
+cmvconcherr: $(EXE)cmvconcherr
         echo $@ " done"
-$(EXE)cmvconchprof: $(OBJ)cmvconchprof.o $(LIB)libcmvsimbao.a
-        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvconchprof.o $(MYLIB) 
-$(OBJ)cmvconchprof.o: cmvconchprof.cc
-        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvconchprof.cc
+$(EXE)cmvconcherr: $(OBJ)cmvconcherr.o $(LIB)libcmvsimbao.a
+        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvconcherr.o $(MYLIB) 
+$(OBJ)cmvconcherr.o: cmvconcherr.cc
+        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvconcherr.cc
 
 

trunk/Cosmo/SimLSS/arrctcast.h

@@ -69 +69 @@
 TArray<T> ArrCastC2R(TArray< complex<T> > & a)  
 {
-  T x;
+  T x = 0;
   return ArrayCast(a, x);
 }
@@ -78 +78 @@
 TArray< complex<T> > ArrCastR2C(TArray< T > & a)  
 {
-  complex<T> x;
+  complex<T> x = 0;
   //  return ArrayCast< TArray< T > , complex<T> > (a, x);
   return ArrayCast(a, x);
@@ -89 +89 @@
 TArray<T> SDRealPart(TArray< complex<T> > & a)  
 {
-  T x;
+  T x = 0;
   return ArrayCast(a, x, 0, 2);
 }
@@ -98 +98 @@
 TArray<T> SDImagPart(TArray< complex<T> > & a)  
 {
-  T x;
+  T x = 0;
   return ArrayCast(a, x, 1, 2);
 }

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -10 +10 @@
 #include "ntuple.h"
 #include "matharr.h"
-#include "srandgen.h"
-#include "perandom.h"
-#include "fabtwriter.h"
-
-#include "arrctcast.h"
 
 #include "constcosmo.h"
@@ -208 +203 @@
  pkz.SetZ(zref);
  TArray< complex<r_8> > pkgen;
- GeneFluct3D fluct3d(pkgen,pkz);
+ GeneFluct3D fluct3d(pkgen);
  fluct3d.SetNThread(nthread);
  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ TArray<r_8>& rgen = fluct3d.GetRealArray();
  fluct3d.Print();
+
+ double dkmin = fluct3d.GetKincMin();
  double knyqmax = fluct3d.GetKmax();
- double dkmin = fluct3d.GetKinc()[0];
+ long nherr = long(knyqmax/dkmin+0.5);
+ cout<<"For HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
+
+ double dktmin = fluct3d.GetKTincMin();
+ double ktnyqmax = fluct3d.GetKTmax();
+ long nherrt = long(ktnyqmax/dktmin+0.5);
+ double dkzmin = fluct3d.GetKinc()[2];
+ double kznyqmax = fluct3d.GetKnyq()[2];
+ long nherrz = long(kznyqmax/dkzmin+0.5);
+ cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
+     <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
 
  cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
@@ -219 +227 @@
  // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
  // Vc/Vs = 0.675   ->  keff = kmax * (0.675)^(1/3) = kmax * 0.877
- knyqmax *= 0.877;
- ldlkint = (log10(knyqmax)-log10(k1int))/npt;
+ double knyqmax_mod = 0.877*knyqmax;
+ ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
  varpk_int.SetInteg(0.01,ldlkint,-1.,4);
- double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax);
- cout<<"varpk_int(<"<<knyqmax<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
+ double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax_mod);
+ cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
 
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
- if(computefourier0) fluct3d.ComputeFourier0();
-   else              fluct3d.ComputeFourier();
-
- cout<<"\n--- Checking realization spectra"<<endl;
- long nhprof = long(fluct3d.GetKmax()/dkmin+0.5);
- HProf hpkgen(0.,fluct3d.GetKmax(),nhprof);
- hpkgen.ReCenterBin();
- fluct3d.ComputeSpectrum(hpkgen);
- {
- tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
+ if(computefourier0) fluct3d.ComputeFourier0(pkz);
+   else              fluct3d.ComputeFourier(pkz);
+
+ if(1) {
+   cout<<"\n--- Checking realization spectra"<<endl;
+   HistoErr hpkgen(0.,knyqmax,nherr);
+   hpkgen.ReCenterBin(); hpkgen.Zero();
+   hpkgen.Show();
+   fluct3d.ComputeSpectrum(hpkgen);
+   {
+   tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
+   }
+ }
+
+ if(1) {
+   cout<<"\n--- Checking realization 2D spectra"<<endl;
+   Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkgen2.ReCenterBin(); hpkgen2.Zero();
+   hpkgen2.Show();
+   fluct3d.ComputeSpectrum2D(hpkgen2);
+   {
+   tagobs = "hpkgen2"; posobs.PutObject(hpkgen2,tagobs);
+   }
  }
 
@@ -241 +262 @@
    cout<<"\n--- Computing convolution by pixel shape"<<endl;
    fluct3d.FilterByPixel();
-
+ }
+
+ if(1) fluct3d.WriteFits("!cmvobserv3d_k.fits");
+ if(1) fluct3d.WritePPF("cmvobserv3d_k.ppf",false);
+
+ if(1) {
    cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
-   HProf hpkgenf(hpkgen); hpkgenf.Zero();
+   HistoErr hpkgenf(0.,knyqmax,nherr);
+   hpkgenf.ReCenterBin(); hpkgenf.Zero();
+   hpkgenf.Show();
    fluct3d.ComputeSpectrum(hpkgenf);
    {
@@ -250 +278 @@
  }
 
- if(0) {
-   cout<<"\n--- Writing cmvobserv3dk.ppf"<<endl;
-   string tag = "cmvobserv3dk.ppf";
-   POutPersist pos(tag);
-   tag = "pkgen"; pos.PutObject(pkgen,tag);
+ if(1) {
+   cout<<"\n--- Checking realization 2D spectra after pixel shape convol."<<endl;
+   Histo2DErr hpkgenf2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkgenf2.ReCenterBin(); hpkgenf2.Zero();
+   hpkgenf2.Show();
+   fluct3d.ComputeSpectrum2D(hpkgenf2);
+   {
+   tagobs = "hpkgenf2"; posobs.PutObject(hpkgenf2,tagobs);
+   }
  }
 
  cout<<"\n--- Computing a realization in real space"<<endl;
  fluct3d.ComputeReal();
- r_8 undouble=0.;
- TArray<r_8> rgen = ArrayCast(pkgen,undouble);
  double rmin,rmax; rgen.MinMax(rmin,rmax);
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
 
- cout<<"\n--- Check mean and variance in real space"<<endl;
- int_8 nlowone = fluct3d.NumberOfBad(-1.,1e+200);
- double rm,rs2; int_8 nm;
- nm = fluct3d.MeanSigma2(rm,rs2);
- cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<",  n(<-1)="<<nlowone<<endl;
+ if(1) fluct3d.WriteFits("!cmvobserv3d_r0.fits");
+ if(1) fluct3d.WritePPF("cmvobserv3d_r0.ppf",true);
+
+ int_8 nm;
+ double rm,rs2;
+ if(1) {
+   cout<<"\n--- Check mean and variance in real space"<<endl;
+   int_8 nlowone = fluct3d.NumberOfBad(-1.,1e+200);
+   nm = fluct3d.MeanSigma2(rm,rs2);
+   cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+       <<rs2<<" -> "<<sqrt(rs2)<<",  n(<-1)="<<nlowone<<endl;
+ }
 
  if(1) {
@@ -331 +367 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- if(1) {
-   cout<<"\n---Writing Cube to Fits file"<<endl;
-   FitsImg3DWriter fwrt("!cmvobserv3dr.fits",FLOAT_IMG,5);
-   fwrt.WriteKey("ZREF",zref," redshift");
-   vector<long> n = fluct3d.GetNpix();
-   fwrt.WriteKey("NX",n[0]," axe transverse 1");
-   fwrt.WriteKey("NY",n[1]," axe transverse 2");
-   fwrt.WriteKey("NZ",n[2]," axe longitudinal (redshift)");
-   vector<r_8> d;
-   d = fluct3d.GetDinc();
-   fwrt.WriteKey("DX",d[0]," Mpc");
-   fwrt.WriteKey("DY",d[1]," Mpc");
-   fwrt.WriteKey("DZ",d[2]," Mpc");
-   d = fluct3d.GetKinc();
-   fwrt.WriteKey("DKX",d[0]," Mpc^-1");
-   fwrt.WriteKey("DKY",d[1]," Mpc^-1");
-   fwrt.WriteKey("DKZ",d[2]," Mpc^-1");
-   fwrt.WriteKey("SNOISE",snoise," Msol");
-   fwrt.Write(rgen);
- }
+ if(1) fluct3d.WriteFits("!cmvobserv3d_r.fits");
+ if(1) fluct3d.WritePPF("cmvobserv3d_r.ppf",true);
 
  cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
@@ -358 +376 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- if(0) {
-   cout<<"\n--- Writing cmvobserv3dr.ppf"<<endl;
-   string tag = "cmvobserv3dr.ppf";
-   POutPersist pos(tag);
-   tag = "rgen"; pos.PutObject(rgen,tag);
- }
-
  //-----------------------------------------------------------------
  // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
+ 
  if(1) {
    cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
    fluct3d.ReComputeFourier();
-   HProf hpkrec(0.,fluct3d.GetKmax(),nhprof);
+ }
+
+ if(1) {
+   cout<<endl<<"\n--- Computing final spectrum"<<endl;
+   HistoErr hpkrec(0.,knyqmax,nherr);
    hpkrec.ReCenterBin();
+   hpkrec.Show();
    fluct3d.ComputeSpectrum(hpkrec);
    tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
  }
 
+ if(1) {
+   cout<<"\n--- Computing final 2D spectrum"<<endl;
+   Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkrec2.ReCenterBin(); hpkrec2.Zero();
+   hpkrec2.Show();
+   fluct3d.ComputeSpectrum2D(hpkrec2);
+   {
+   tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
+   }
+ }
+
  return 0;
 }
 
 /*
-openppf cmvobserv3dk.ppf
-openppf cmvobserv3dr.ppf
+######################################################
+readfits cmvobserv3d_k.fits
+readfits cmvobserv3d_r0.fits
+readfits cmvobserv3d_r.fits
+
+openppf cmvobserv3d_k.ppf
+openppf cmvobserv3d_r0.ppf
+openppf cmvobserv3d_r.ppf
+
+# pour le plot 2D d'une slice en Z du 3D: to2d nom_obj3D num_slice
+defscript to2d
+ objaoper $1 sliceyz $2
+ mv sliceyz_${2} ${1}_$2
+ disp ${1}_$2
+ echo display slice $2 of $1
+endscript 
+
+to2d $cobj 0
+
+######################################################
 openppf cmvobserv3d.ppf
 
-objaoper pkgen sliceyz 0
-
+zone
 n/plot hpkz.val%x ! ! "nsta connectpoints"
 n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
@@ -394 +439 @@
 n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
 
-defscript rgensl
- objaoper rgen sliceyz $1
- disp sliceyz_$1 
-endscript 
-
-rgensl 0
-
+zone 2 2
+imag hpkgen2
+imag hpkgenf2
+imag hpkrec2
+
+zone
 disp hmdndm
 disp tirhmdndm

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

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -5 +5 @@
 #include "genericfunc.h"
 #include "tarray.h"
-#include "hisprof.h"
+#include "histerr.h"
+#include "hist2err.h"
 #include "perandom.h"
+
+#include "FFTW/fftw3.h"
+#include "FitsIO/fitsio.h"
 
 #include <vector>
@@ -12 +16 @@
 #include "pkspectrum.h"
 
-#include "FFTW/fftw3.h"
 
 namespace SOPHYA {
@@ -19 +22 @@
 class GeneFluct3D {
 public:
-  GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz);
+  GeneFluct3D(TArray< complex<r_8 > >& T);
   virtual ~GeneFluct3D(void);
 
@@ -25 +28 @@
   void SetSize(long nx,long ny,long nz,double dx,double dy,double dz);  // Mpc
 
-  inline void SetZ(double z) {pkz_.SetZ(z);}
-  double GetZ(void) {return pkz_.GetZ();}
+  TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
+  fftw_complex* GetComplexPointer(void) {return fdata_;}
+  TArray<r_8>& GetRealArray(void) {return R_;}
+  r_8* GetRealPointer(void) {return data_;}
+
+  // Pour adressage data_[ip]
+  inline int_8 IndexR(long i,long j,long k) {return (int_8)(k+NTz_*(j+Ny_*i));}
+  // Pour adressage fdata_[ip][0-1]
+  inline int_8 IndexC(long i,long j,long k) {return (int_8)(k+NCz_*(j+Ny_*i));}
+  // - On peut aussi adresser:
+  // TArray< complex<r_8> >& pk = gf3d.GetComplexArray(); pk(i,j,k) = ...;
+  // TArray<r_8>& rgen = gf3d.GetRealArray(); rgen(i,j,k) = ...;
+
+  vector<long> GetNpix(void) {return N_;}
+  int_8 NPix(void) {return NRtot_;}
+
+  vector<r_8> GetDinc(void) {return D_;}
+  double GetDVol(void) {return dVol_;}
   double GetVol(void) {return Vol_;}
-  double GetDVol(void) {return dVol_;}
-  int_8 NPix(void) {return NRtot_;}
+
+  vector<r_8> GetKinc(void) {return Dk_;}
+  vector<r_8> GetKnyq(void) {return Knyq_;}
   double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
-  vector<r_8> GetKinc(void)
-    {vector<r_8> dk; dk.push_back(Dkx_); dk.push_back(Dky_); dk.push_back(Dkz_); return dk;}
-  vector<r_8> GetKnyq(void)
-    {vector<r_8> kn; kn.push_back(Knyqx_); kn.push_back(Knyqy_); kn.push_back(Knyqz_); return kn;}
-  vector<r_8> GetDinc(void)
-    {vector<r_8> d; d.push_back(Dx_); d.push_back(Dy_); d.push_back(Dz_); return d;}
-  vector<long> GetNpix(void)
-    {vector<long> d; d.push_back(Nx_); d.push_back(Ny_); d.push_back(Nz_); return d;}
+  double GetKTmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_);}
+  double GetKincMin(void)
+    {vector<r_8>::const_iterator it = min_element(Dk_.begin(), Dk_.end()); return *it;}
+  double GetKTincMin(void) {return min(Dk_[0],Dk_[1]);}
 
-  void ComputeFourier0(void);
-  void ComputeFourier(void);
+  void ComputeFourier0(GenericFunc& pk_at_z);
+  void ComputeFourier(GenericFunc& pk_at_z);
   void FilterByPixel(void);
-  int  ComputeSpectrum(HProf& hp);
+
   void ComputeReal(void);
+
   void ReComputeFourier(void);
+
+  int  ComputeSpectrum(HistoErr& herr);
+  int  ComputeSpectrum2D(Histo2DErr& herr);
 
   int_8 VarianceFrReal(double R,double& var);
@@ -59 +79 @@
   void AddNoise2Real(double snoise);
 
+  void WriteFits(string cfname,int bitpix=FLOAT_IMG);
+  void ReadFits(string cfname);
+
+  void WritePPF(string cfname,bool write_real=true);
+  void ReadPPF(string cfname);
+
   void Print(void);
 
 protected:
+  void setsize(long nx,long ny,long nz,double dx,double dy,double dz);
+  void setalloc(void);
+  void setpointers(bool from_real);
   long manage_coefficients(void);
   double compute_power_carte(void);
+  void check_array_alloc(void);
   inline double pixelfilter(double x)
     {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
 
 
-  long Nx_,Ny_,Nz_;
-  sa_size_t SzK_[3];
+  long Nx_,Ny_,Nz_;  vector<long> N_;
   long NCz_,NTz_;
   int_8 NRtot_;
-  double Dx_,Dy_,Dz_;
+
+  double Dx_,Dy_,Dz_;  vector<double> D_;
+
+  double Dkx_,Dky_,Dkz_;  vector<double> Dk_;
+  double Knyqx_,Knyqy_,Knyqz_;  vector<double> Knyq_;
+  double Dk3_;
   double dVol_, Vol_;
-  double Dkx_,Dky_,Dkz_;
-  double Knyqx_,Knyqy_,Knyqz_;
-  double Dk3_;
 
   fftw_plan pf_,pb_;
   unsigned short nthread_;
 
-  //   = 0 : T_ vide
-  //   = 1 : T_ contient une realisation du spectre (espace fourier)
-  //   = 2 : T_ contient une realisation dans l'espace reel
-  //   = 3 : T_ contient la TF d'une realisation dans l'espace reel (d_rho/rho)
-  //   = 4 : T_ contient (1 + d_rho/rho)
-  //   = ... etc ...
-  unsigned short Tcontent_;
+  bool array_allocated_;  // true if array has been allocated
 
-  TArray< complex<r_8 > >& T_;
-  PkSpectrumZ& pkz_;
+  TArray< complex<r_8> >& T_;
+  fftw_complex *fdata_;
+  TArray<r_8> R_;
+  double *data_;
 };
 

trunk/Cosmo/SimLSS/genericfunc.h

@@ -3 +3 @@
 
 #include "pexceptions.h"
+#include <vector>
 
 namespace SOPHYA {
@@ -32 +33 @@
   }
 
+  virtual double operator()(vector<double>& x) {
+    cout<<"GenericFunc::operator(vector<double>&) not implemented"<<endl;
+    throw NotAvailableOperation("GenericFunc::operator(vector<double>&) not implemented");
+  }
+
 };