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

Last change on this file since 3141 was 3141, checked in by cmv, 19 years ago

chgt HProf->HistoErr + spectre 2D cmv 17/01/2007
File size: 3.5 KB
Line 
1#ifndef GENEFLUCT3D_SEEN
2#define GENEFLUCT3D_SEEN
3
4#include "machdefs.h"
5#include "genericfunc.h"
6#include "tarray.h"
7#include "histerr.h"
8#include "hist2err.h"
9#include "perandom.h"
10
11#include "FFTW/fftw3.h"
12#include "FitsIO/fitsio.h"
13
14#include <vector>
15
16#include "pkspectrum.h"
17
18
19namespace SOPHYA {
20
21//-----------------------------------------------------------------------------------
22class GeneFluct3D {
23public:
24 GeneFluct3D(TArray< complex<r_8 > >& T);
25 virtual ~GeneFluct3D(void);
26
27 void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
28 void SetSize(long nx,long ny,long nz,double dx,double dy,double dz); // Mpc
29
30 TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
31 fftw_complex* GetComplexPointer(void) {return fdata_;}
32 TArray<r_8>& GetRealArray(void) {return R_;}
33 r_8* GetRealPointer(void) {return data_;}
34
35 // Pour adressage data_[ip]
36 inline int_8 IndexR(long i,long j,long k) {return (int_8)(k+NTz_*(j+Ny_*i));}
37 // Pour adressage fdata_[ip][0-1]
38 inline int_8 IndexC(long i,long j,long k) {return (int_8)(k+NCz_*(j+Ny_*i));}
39 // - On peut aussi adresser:
40 // TArray< complex<r_8> >& pk = gf3d.GetComplexArray(); pk(i,j,k) = ...;
41 // TArray<r_8>& rgen = gf3d.GetRealArray(); rgen(i,j,k) = ...;
42
43 vector<long> GetNpix(void) {return N_;}
44 int_8 NPix(void) {return NRtot_;}
45
46 vector<r_8> GetDinc(void) {return D_;}
47 double GetDVol(void) {return dVol_;}
48 double GetVol(void) {return Vol_;}
49
50 vector<r_8> GetKinc(void) {return Dk_;}
51 vector<r_8> GetKnyq(void) {return Knyq_;}
52 double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
53 double GetKTmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_);}
54 double GetKincMin(void)
55 {vector<r_8>::const_iterator it = min_element(Dk_.begin(), Dk_.end()); return *it;}
56 double GetKTincMin(void) {return min(Dk_[0],Dk_[1]);}
57
58 void ComputeFourier0(GenericFunc& pk_at_z);
59 void ComputeFourier(GenericFunc& pk_at_z);
60 void FilterByPixel(void);
61
62 void ComputeReal(void);
63
64 void ReComputeFourier(void);
65
66 int ComputeSpectrum(HistoErr& herr);
67 int ComputeSpectrum2D(Histo2DErr& herr);
68
69 int_8 VarianceFrReal(double R,double& var);
70 int_8 MeanSigma2(double& rm,double& rs2,double vmin=-1.e+150,double vmax=1.e+150);
71 int_8 NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
72 int_8 SetToVal(double vmin, double vmax,double val0=0.);
73
74 void TurnFluct2Mass(void);
75 double TurnMass2MeanNumber(double n_by_mpc3);
76 double ApplyPoisson(void);
77 double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
78 double TurnMass2Flux(void);
79 void AddNoise2Real(double snoise);
80
81 void WriteFits(string cfname,int bitpix=FLOAT_IMG);
82 void ReadFits(string cfname);
83
84 void WritePPF(string cfname,bool write_real=true);
85 void ReadPPF(string cfname);
86
87 void Print(void);
88
89protected:
90 void setsize(long nx,long ny,long nz,double dx,double dy,double dz);
91 void setalloc(void);
92 void setpointers(bool from_real);
93 long manage_coefficients(void);
94 double compute_power_carte(void);
95 void check_array_alloc(void);
96 inline double pixelfilter(double x)
97 {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
98
99
100 long Nx_,Ny_,Nz_; vector<long> N_;
101 long NCz_,NTz_;
102 int_8 NRtot_;
103
104 double Dx_,Dy_,Dz_; vector<double> D_;
105
106 double Dkx_,Dky_,Dkz_; vector<double> Dk_;
107 double Knyqx_,Knyqy_,Knyqz_; vector<double> Knyq_;
108 double Dk3_;
109 double dVol_, Vol_;
110
111 fftw_plan pf_,pb_;
112 unsigned short nthread_;
113
114 bool array_allocated_; // true if array has been allocated
115
116 TArray< complex<r_8> >& T_;
117 fftw_complex *fdata_;
118 TArray<r_8> R_;
119 double *data_;
120};
121
122} // Fin du namespace
123
124#endif
Note: See TracBrowser for help on using the repository browser.