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

Last change on this file since 3283 was 3283, checked in by cmv, 18 years ago

scale du spectre final cmv 24/07/2007
File size: 5.0 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 "cosmocalc.h"
17#include "pkspectrum.h"
18
19
20namespace SOPHYA {
21
22//-----------------------------------------------------------------------------------
23class GeneFluct3D {
24public:
25 GeneFluct3D(TArray< complex<r_8 > >& T);
26 virtual ~GeneFluct3D(void);
27
28 void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
29 void SetSize(long nx,long ny,long nz,double dx,double dy,double dz); // Mpc
30 // Distance los comobile a l'observateur
31 void SetObservator(double redshref=0.,double kredshref=0.);
32 inline double DXcom(long i) {return i*Dx_ - xobs_[0];}
33 inline double DYcom(long j) {return j*Dy_ - xobs_[1];}
34 inline double DZcom(long k) {return k*Dz_ - xobs_[2];}
35 inline double Dcom(long i,long j,long k) {
36 double dx=DXcom(i), dy=DYcom(j), dz=DZcom(k);
37 return sqrt(dx*dx+dy*dy+dz*dz);
38 }
39 void SetCosmology(CosmoCalc& cosmo);
40 void SetGrowthFactor(GrowthFactor& growth);
41 long LosComRedshift(double zinc=0.001,long npoints=-1);
42
43 TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
44 fftw_complex* GetComplexPointer(void) {return fdata_;}
45 TArray<r_8>& GetRealArray(void) {return R_;}
46 r_8* GetRealPointer(void) {return data_;}
47
48 // Pour adressage data_[ip]
49 inline int_8 IndexR(long i,long j,long k) {return (int_8)(k+NTz_*(j+Ny_*i));}
50 // Pour adressage fdata_[ip][0-1]
51 inline int_8 IndexC(long i,long j,long k) {return (int_8)(k+NCz_*(j+Ny_*i));}
52 // - On peut aussi adresser:
53 // TArray< complex<r_8> >& pk = gf3d.GetComplexArray(); pk(i,j,k) = ...;
54 // TArray<r_8>& rgen = gf3d.GetRealArray(); rgen(i,j,k) = ...;
55
56 vector<long> GetNpix(void) {return N_;}
57 int_8 NPix(void) {return NRtot_;}
58
59 vector<r_8> GetDinc(void) {return D_;}
60 double GetDVol(void) {return dVol_;}
61 double GetVol(void) {return Vol_;}
62
63 vector<r_8> GetKinc(void) {return Dk_;}
64 vector<r_8> GetKnyq(void) {return Knyq_;}
65 double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
66 double GetKTmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_);}
67 double GetKincMin(void)
68 {vector<r_8>::const_iterator it = min_element(Dk_.begin(), Dk_.end()); return *it;}
69 double GetKTincMin(void) {return min(Dk_[0],Dk_[1]);}
70
71 void ComputeFourier0(GenericFunc& pk_at_z);
72 void ComputeFourier(GenericFunc& pk_at_z);
73 void FilterByPixel(void);
74
75 void ComputeReal(void);
76
77 void ReComputeFourier(void);
78
79 int ComputeSpectrum(HistoErr& herr);
80 int ComputeSpectrum2D(Histo2DErr& herr);
81 void ApplyGrowthFactor(void);
82
83 int_8 VarianceFrReal(double R,double& var);
84 int_8 MeanSigma2(double& rm,double& rs2,double vmin=1.,double vmax=-1.
85 ,bool useout=false,double vout=0.);
86 int_8 NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
87 int_8 SetToVal(double vmin, double vmax,double val0=0.);
88 void ScaleOffset(double scalecube=1.,double offsetcube=0.);
89
90 void TurnFluct2Mass(void);
91 double TurnMass2MeanNumber(double n_by_mpc3);
92 double ApplyPoisson(void);
93 double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
94 double TurnMass2Flux(void);
95 void AddAGN(double lfjy,double lsigma,double powlaw=0.);
96 void AddNoise2Real(double snoise,bool with_evol=false);
97
98 void WriteFits(string cfname,int bitpix=FLOAT_IMG);
99 void ReadFits(string cfname);
100
101 void WritePPF(string cfname,bool write_real=true);
102 void ReadPPF(string cfname);
103 void WriteSlicePPF(string cfname);
104
105 void SetPrtLevel(int lp=0) {lp_ = lp;}
106 void Print(void);
107
108//-------------------------------------------------------------------
109
110protected:
111 void setsize(long nx,long ny,long nz,double dx,double dy,double dz);
112 void setalloc(void);
113 void setpointers(bool from_real);
114 void init_fftw(void);
115 long manage_coefficients(void);
116 double compute_power_carte(void);
117 void check_array_alloc(void);
118 inline double pixelfilter(double x)
119 {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
120
121 // valeurs dans l'espace reel
122 long Nx_,Ny_,Nz_; vector<long> N_;
123 long NCz_,NTz_;
124 int_8 NRtot_;
125
126 double Dx_,Dy_,Dz_; vector<double> D_;
127
128 // valeurs dans l'espace des K
129 double Dkx_,Dky_,Dkz_; vector<double> Dk_;
130 double Knyqx_,Knyqy_,Knyqz_; vector<double> Knyq_;
131 double Dk3_;
132 double dVol_, Vol_;
133
134 // la gestion de la FFT
135 fftw_plan pf_,pb_;
136 unsigned short nthread_;
137 int lp_;
138
139 // le stockage du Cube de donnees et les pointeurs
140 bool array_allocated_; // true if array has been allocated
141 TArray< complex<r_8> >& T_;
142 fftw_complex *fdata_;
143 TArray<r_8> R_;
144 double *data_;
145
146 // l'observateur
147 CosmoCalc *cosmo_;
148 GrowthFactor *growth_;
149 double redsh_ref_,kredsh_ref_,dred_ref_;
150 double loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
151 double nu_ref_, dnu_ref_ ;
152 double xobs_[3];
153 double loscom_min_, loscom_max_;
154 vector<double> zred_, loscom_;
155 double loscom2zred_min_, loscom2zred_max_;
156 vector<double> loscom2zred_;
157
158};
159
160} // Fin du namespace
161
162#endif
Note: See TracBrowser for help on using the repository browser.