source: Sophya/trunk/Cosmo/SimLSS/genefluct3d.h@ 3115

#ifndef GENEFLUCT3D_SEEN
#define GENEFLUCT3D_SEEN

#include "machdefs.h"
#include "genericfunc.h"
#include "tarray.h"
#include "hisprof.h"
#include "perandom.h"

#include <vector>

#include "pkspectrum.h"

#include "FFTW/fftw3.h"

namespace SOPHYA {

//-----------------------------------------------------------------------------------
class GeneFluct3D {
public:
  GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz);
  virtual ~GeneFluct3D(void);

  void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
  void SetSize(int nx,int ny,int nz,double dx,double dy,double dz);  // Mpc

  inline void SetZ(double z) {pkz_.SetZ(z);}
  double GetZ(void) {return pkz_.GetZ();}
  double GetVol(void) {return Vol_;}
  double GetDVol(void) {return dVol_;}
  size_t NPix(void) {return NRtot_;}
  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;}

  void ComputeFourier0(void);
  void ComputeFourier(void);
  void FilterByPixel(void);
  int  ComputeSpectrum(HProf& hp);
  void ComputeReal(void);
  void ReComputeFourier(void);

  size_t VarianceFrReal(double R,double& var);
  size_t MeanSigma2(double& rm,double& rs2,double vmin=-1.e+150,double vmax=1.e+150);
  size_t NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
  size_t SetToVal(double vmin, double vmax,double val0=0.);

  void TurnFluct2Mass(void);
  double TurnMass2MeanNumber(double n_by_mpc3);
  double ApplyPoisson(void);
  double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
  double TurnMass2Flux(void);
  void AddNoise2Real(double snoise);

  void Print(void);

protected:
  int manage_coefficients(void);
  double compute_power_carte(void);
  inline double pixelfilter(double kx)
    {return (kx<0.05) ? 1.-kx*kx/12.*(1.-kx*kx/30.): 2.*(1.-cos(kx))/(kx*kx);}


  int Nx_,Ny_,Nz_,SzK_[3];
  int NCz_,NTz_;
  size_t NRtot_;
  double Dx_,Dy_,Dz_;
  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_;

  TArray< complex<r_8 > >& T_;
  PkSpectrumZ& pkz_;
};

} // Fin du namespace

#endif