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

Last change on this file since 3615 was 3524, checked in by cmv, 17 years ago

intro du remplissage du NTuple de debug cmv 22/09/2008
File size: 6.6 KB
Line 
1#ifndef GENEFLUCT3D_SEEN
2#define GENEFLUCT3D_SEEN
3
4#include "machdefs.h"
5#include <math.h>
6#include "genericfunc.h"
7#include "tarray.h"
8#include "histerr.h"
9#include "hist2err.h"
10#include "perandom.h"
11
12#include "FFTW/fftw3.h"
13#include "FitsIO/fitsio.h"
14
15#include <vector>
16
17#include "cosmocalc.h"
18#include "pkspectrum.h"
19
20#define WITH_FFTW_THREAD
21
22//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE #define GEN3D_FLOAT
23
24#if defined(GEN3D_FLOAT)
25#define GEN3D_TYPE r_4
26#define GEN3D_FFTW_PLAN fftwf_plan
27#define GEN3D_FFTW_COMPLEX fftwf_complex
28#else
29#define GEN3D_TYPE r_8
30#define GEN3D_FFTW_PLAN fftw_plan
31#define GEN3D_FFTW_COMPLEX fftw_complex
32#endif
33
34namespace SOPHYA {
35
36//-----------------------------------------------------------------------------------
37class GeneFluct3D {
38public:
39 GeneFluct3D(long nx,long ny,long nz,double dx,double dy,double dz,unsigned short nthread=0,int lp=0); // Mpc
40 GeneFluct3D(unsigned short nthread=0);
41 virtual ~GeneFluct3D(void);
42
43 // Distance los comobile a l'observateur
44 void SetObservator(double redshref=0.,double kredshref=0.);
45 inline double DXcom(long i) {return i*Dx_ - xobs_[0];}
46 inline double DYcom(long j) {return j*Dy_ - xobs_[1];}
47 inline double DZcom(long k) {return k*Dz_ - xobs_[2];}
48 inline double Dcom(long i,long j,long k) {
49 double dx=DXcom(i), dy=DYcom(j), dz=DZcom(k);
50 return sqrt(dx*dx+dy*dy+dz*dz);
51 }
52 void SetCosmology(CosmoCalc& cosmo);
53 void SetGrowthFactor(GrowthFactor& growth);
54 long LosComRedshift(double zinc=0.001,long npoints=-1);
55
56 TArray< complex<GEN3D_TYPE> >& GetComplexArray(void) {return T_;}
57 GEN3D_FFTW_COMPLEX * GetComplexPointer(void) {return fdata_;}
58 TArray<GEN3D_TYPE>& GetRealArray(void) {return R_;}
59 GEN3D_TYPE* GetRealPointer(void) {return data_;}
60
61 // Pour adressage data_[ip]
62 inline int_8 IndexR(long i,long j,long k) {return (int_8)(k+NTz_*(j+Ny_*i));}
63 // Pour adressage fdata_[ip][0-1]
64 inline int_8 IndexC(long i,long j,long k) {return (int_8)(k+NCz_*(j+Ny_*i));}
65 // On peut aussi adresser:
66 // TArray< complex<r_8> >& pk = gf3d.GetComplexArray();
67 // pk(k,j,i) avec k=[0,NCz_[ j=[0,Ny_[ i=[0,Nx_[
68 // pk[IndexC(i,j,k)]
69 // TArray<r_8>& rgen = gf3d.GetRealArray();
70 // rgen(k,j,i) avec k=[0,NTz_[ j=[0,Ny_[ i=[0,Nx_[
71 // mais seul k=[0,Nz_[ est utile
72 // rgen[IndexR(i,j,k)]
73 // ATTENTION: TArray adresse en memoire a l'envers du C !
74 // Tarray(n1,n2,n3) == Carray[n3][n2][n1]
75
76 vector<long> GetNpix(void) {return N_;}
77 int_8 NPix(void) {return NRtot_;}
78 long GetNx(void) {return Nx_;}
79 long GetNy(void) {return Ny_;}
80 long GetNz(void) {return Nz_;}
81
82 // Return |K_i| module relative to pixel indices
83 inline r_8 Kx(long i) {long ii=(i>Nx_/2)? Nx_-i :i; return ii*Dkx_;}
84 inline r_8 Ky(long j) {long jj=(j>Ny_/2)? Ny_-j :j; return jj*Dky_;}
85 inline r_8 Kz(long l) {return l*Dkz_;}
86
87 vector<r_8> GetDinc(void) {return D_;}
88 double GetDVol(void) {return dVol_;}
89 double GetVol(void) {return Vol_;}
90
91 vector<r_8> GetKinc(void) {return Dk_;}
92 vector<r_8> GetKnyq(void) {return Knyq_;}
93 double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
94 double GetKTmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_);}
95 double GetKincMin(void)
96 {vector<r_8>::const_iterator it = min_element(Dk_.begin(), Dk_.end()); return *it;}
97 double GetKincMax(void)
98 {vector<r_8>::const_iterator it = max_element(Dk_.begin(), Dk_.end()); return *it;}
99 double GetKTincMin(void) {return min(Dk_[0],Dk_[1]);}
100 double GetKTincMax(void) {return max(Dk_[0],Dk_[1]);}
101
102 void ComputeFourier0(GenericFunc& pk_at_z);
103 void ComputeFourier(GenericFunc& pk_at_z);
104 void FilterByPixel(void);
105
106 void ComputeReal(void);
107 void ApplyGrowthFactor(int type_evol=1);
108
109 void ReComputeFourier(void);
110
111 int ComputeSpectrum(HistoErr& herr);
112 int ComputeSpectrum2D(Histo2DErr& herr);
113 int ComputeSpectrum(HistoErr& herr,double sigma,bool pixcor);
114 int ComputeSpectrum2D(Histo2DErr& herr,double sigma,bool pixcor);
115
116 int_8 VarianceFrReal(double R,double& var);
117 int_8 MeanSigma2(double& rm,double& rs2,double vmin=1.,double vmax=-1.
118 ,bool useout=false,double vout=0.);
119 int_8 MinMax(double& xmin,double& xmax,double vmin=1.,double vmax=-1.);
120 int_8 NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
121 int_8 SetToVal(double vmin, double vmax,double val0=0.);
122 void ScaleOffset(double scalecube=1.,double offsetcube=0.);
123
124 void TurnFluct2Mass(void);
125 double TurnFluct2MeanNumber(double val_by_mpc3);
126 double ApplyPoisson(void);
127 double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
128 double TurnNGal2MassQuick(SchechterMassDist& schmdist);
129 double TurnMass2Flux(void);
130 //void AddAGN(double lfjy,double lsigma,double powlaw=0.);
131 void AddNoise2Real(double snoise,int type_evol=0);
132
133 void WriteFits(string cfname,int bitpix=FLOAT_IMG);
134 void ReadFits(string cfname);
135
136 void WritePPF(string cfname,bool write_real=true);
137 void ReadPPF(string cfname);
138 void WriteSlicePPF(string cfname);
139 void NTupleCheck(POutPersist &pos,string ntname,unsigned long nent);
140
141 void SetPrtLevel(int lp=0) {lp_ = lp;}
142 void Print(void);
143
144//-------------------------------------------------------------------
145
146protected:
147 void init_default(void);
148 void setsize(long nx,long ny,long nz,double dx,double dy,double dz);
149 void setalloc(void);
150 void setpointers(bool from_real);
151 void init_fftw(void);
152 void delete_fftw(void);
153 long manage_coefficients(void);
154 double compute_power_carte(void);
155 void check_array_alloc(void);
156 inline double pixelfilter(double x)
157 {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
158
159 // valeurs dans l'espace reel
160 long Nx_,Ny_,Nz_; vector<long> N_;
161 long NCz_,NTz_;
162 int_8 NRtot_;
163
164 double Dx_,Dy_,Dz_; vector<double> D_;
165
166 // valeurs dans l'espace des K
167 double Dkx_,Dky_,Dkz_; vector<double> Dk_;
168 double Knyqx_,Knyqy_,Knyqz_; vector<double> Knyq_;
169 double Dk3_;
170 double dVol_, Vol_;
171
172 // la gestion de la FFT
173 bool is_set_fft_plan;
174 GEN3D_FFTW_PLAN pf_,pb_;
175 unsigned short nthread_;
176 int lp_;
177
178 // le stockage du Cube de donnees et les pointeurs
179 bool array_allocated_; // true if array has been allocated
180 unsigned short array_type; // 0=empty, 1=real, 2=complex
181 TArray< complex<GEN3D_TYPE> > T_;
182 GEN3D_FFTW_COMPLEX *fdata_;
183 TArray<GEN3D_TYPE> R_;
184 GEN3D_TYPE *data_;
185
186 // l'observateur
187 CosmoCalc *cosmo_;
188 GrowthFactor *growth_;
189 double redsh_ref_,kredsh_ref_,dred_ref_;
190 double loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
191 double nu_ref_, dnu_ref_ ;
192 double xobs_[3];
193 double loscom_min_, loscom_max_;
194 vector<double> zred_, loscom_;
195 double loscom2zred_min_, loscom2zred_max_;
196 vector<double> loscom2zred_;
197
198};
199
200} // Fin du namespace SOPHYA
201
202#endif
Note: See TracBrowser for help on using the repository browser.