Index: trunk/Cosmo/SimLSS/genefluct3d.h
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3134)
+++ trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3141)
@@ -5,6 +5,10 @@
 #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,5 +16,4 @@
 #include "pkspectrum.h"
 
-#include "FFTW/fftw3.h"
 
 namespace SOPHYA {
@@ -19,5 +22,5 @@
 class GeneFluct3D {
 public:
-  GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz);
+  GeneFluct3D(TArray< complex<r_8 > >& T);
   virtual ~GeneFluct3D(void);
 
@@ -25,25 +28,42 @@
   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,36 +79,43 @@
   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_;
 };
 
