//  Class examples to generate mass distribution 
// R.A. for A. Abate , Nov. 2008

#ifndef SPECPK_SEEN
#define SPECPK_SEEN

#include "machdefs.h"      
#include "sopnamsp.h"       
#include <math.h>
#include <iostream>
#include <vector>
#include <string> 

#include "genericfunc.h"  
#include "array.h"         
#include "histats.h"       
#include "fftwserver.h"    
#include "randinterf.h"      

#include "mdish.h"       

#define DeuxPI 2.*M_PI 

// -- SpectralShape class : test P(k) class
class SpectralShape : public GenericFunc {
public:
  SpectralShape(int typ);
// Return the value of power spectrum for wave number wk   
  virtual double operator() (double wk); 
  inline  double Value(double wk) { return((*this)(wk)); }
// Return a vector representing the power spectrum (for checking) 
  Histo GetPk(int n=256);
  double Sommek2Pk(double kmax=1000., int n=5000);
  inline void SetRenormFac(double f=1.) { renorm_fac=f; } 
  int typ_;
  double renorm_fac;
};


#define TF  r_4 

// -- Four3DPk class :  3D fourier amplitudes and power spectrum 
class Four3DPk {
public:
// Constructor
  Four3DPk(TArray< complex<TF> > & fourcoedd, RandomGeneratorInterface& rg);
  Four3DPk(RandomGeneratorInterface& rg, sa_size_t szx=128, sa_size_t szy=256, sa_size_t szz=128);
  virtual ~Four3DPk(); 

  inline void SetCellSize(double dkx=DeuxPI, double dky=DeuxPI, double dkz=DeuxPI) 
  { dkx_=dkx;  dky_=dky;  dkz_=dkz; } 
  inline int SetPrtLevel(int lev=0, int prtmod=10) 
  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
  void ComputeFourierAmp(SpectralShape& pk);
// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
// typically = ComovRadialDistance 
  void ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale=1., bool crmask=false);
  void ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales, Vector& noisp);

// Return the array size 
  inline sa_size_t NCells() { return fourAmp.Size(); }
  inline sa_size_t SizeX() { return fourAmp.SizeX(); }
  inline sa_size_t SizeY() { return fourAmp.SizeY(); }
  inline sa_size_t SizeZ() { return fourAmp.SizeZ(); }

// Set the cell size/step in Fourier Space
// Return the fourier amplitude matrix  
  TArray< complex<TF> > GetFourierAmp()
    { return fourAmp; }
// Return the mass density matrix    
  TArray<TF> ComputeMassDens(); 

// Return the reconstructed power spectrum as a profile histogram   
  HProf ComputePk(double s2cut=0., int nbin=256, double kmin=0., double kmax=-1., bool fgmodcnt=false);
  void  ComputePkCumul();

// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
// typically = ComovRadialDistance 
  HProf ComputeNoisePk(Four2DResponse& resp, double angscale=1., double s2cut=0., 
		       int nbin=256, double kmin=0., double kmax=-1.);

  // Fills a data table from the computed P(k) profile histogram and mode count 
  Histo FillPkDataTable(DataTable& dt);
  inline HProf& GetPk() { return *hp_pk_p_; }

protected:
  // member attribute
  RandomGeneratorInterface& rg_;
  TArray< complex<TF> > fourAmp;  // complex array of fourier coefficients
  double dkx_, dky_, dkz_;
  int prtlev_;
  int prtmodulo_;
  // Profile histograms for power spectrum and number of modes 
  HProf* hp_pk_p_;
  Histo* hmcnt_p_;
  Histo* hmcntok_p_;
  double s2cut_;
};

// --- PkNoiseCalculator : 
// - Classe de calcul du spectre de bruit PNoise(k) determine par une reponse
//   2D de l'instrument 
class PkNoiseCalculator
{
public:
  PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut=100., int ngen=1, const char* tit="PkNoise");
  
  inline void SetFreqRange(double freq0=835.,double dfreq=0.5)
  { freq0_=freq0;  dfreq_=dfreq; }
  inline void SetAngScaleConversion(double angscale=1.)
  { angscales_=angscale; }
  inline void SetAngScaleConversion(Vector& angscs)
  { angscales_=angscs; }
  inline void SetPNoiseFactor(double pnoisef=1.)
  { pnoisefac_=pnoisef; }
  inline void SetPNoiseFactor(Vector& pnoisefac)
  { pnoisefac_=pnoisefac; }
  inline void SetS2Cut(double s2cut=100.)
  {  S2CUT=s2cut; }
  inline double GetS2Cut() { return S2CUT; }
  HProf Compute(int nbin=256, double kmin=0., double kmax=-1.);
  inline int SetPrtLevel(int lev=0, int prtmod=10) 
  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }

protected:
  Four3DPk& pkn3d;
  Four2DResponse& frep;
  double freq0_,dfreq_;
  Vector angscales_; 
  Vector pnoisefac_; 
  double S2CUT;
  int NGEN;
  string title;
  int prtlev_;
  int prtmodulo_;
};



// -- MassDist2D class :  2D mass distribution 
class MassDist2D {
public:
// Constructor
  MassDist2D(GenericFunc& pk, int size=1024, double meandens=1.);
// Do the computation 
  void Compute();
// Return the array size 
  inline sa_size_t ArrSize() { return sizeA; }
// Return the fourier amplitude matrix  
  TMatrix< complex<r_8> > GetFourierAmp()
    { if (!fg_fourAmp) ComputeFourierAmp(); return fourAmp; }
// Return the mass density matrix    
  Matrix GetMassDens()
    { if (!fg_massDens) ComputeMassDens(); return massDens; }

// Return the reconstructed power spectrum as a profile histogram   
  HProf ReconstructPk(int nbin=0);
protected:
  void ComputeFourierAmp();
  void ComputeMassDens();
  
// member attribute
  GenericFunc& pkSpec;   // The spectralShape
  sa_size_t sizeA;       // 2D array size
  double meanRho;       // Mean Density
  bool fg_fourAmp;  // true -> fourAmp computed 
  TMatrix< complex<r_8> > fourAmp;  // complex array of fourier coefficients
  bool fg_massDens;  // true -> MassDens computed 
  TMatrix< r_8 > massDens;  	// real array of d rho/rho 
};


#endif