1 | //#include "sopnamsp.h"
|
---|
2 | #include "machdefs.h"
|
---|
3 | #include <iostream>
|
---|
4 | #include <stdlib.h>
|
---|
5 | #include <stdio.h>
|
---|
6 | #include <string.h>
|
---|
7 | #include <math.h>
|
---|
8 | #include <unistd.h>
|
---|
9 |
|
---|
10 | #include "tarray.h"
|
---|
11 | #include "pexceptions.h"
|
---|
12 | #include "perandom.h"
|
---|
13 | #include "srandgen.h"
|
---|
14 |
|
---|
15 | #include "fabtcolread.h"
|
---|
16 | #include "fabtwriter.h"
|
---|
17 | #include "fioarr.h"
|
---|
18 | #include "ntuple.h"
|
---|
19 |
|
---|
20 | #include "arrctcast.h"
|
---|
21 |
|
---|
22 | #include "constcosmo.h"
|
---|
23 | #include "geneutils.h"
|
---|
24 | #include "schechter.h"
|
---|
25 |
|
---|
26 | #include "genefluct3d.h"
|
---|
27 |
|
---|
28 | #if defined(GEN3D_FLOAT)
|
---|
29 | #define GEN3D_FFTW_INIT_THREADS fftwf_init_threads
|
---|
30 | #define GEN3D_FFTW_CLEANUP_THREADS fftwf_cleanup_threads
|
---|
31 | #define GEN3D_FFTW_PLAN_WITH_NTHREADS fftwf_plan_with_nthreads
|
---|
32 | #define GEN3D_FFTW_PLAN_DFT_R2C_3D fftwf_plan_dft_r2c_3d
|
---|
33 | #define GEN3D_FFTW_PLAN_DFT_C2R_3D fftwf_plan_dft_c2r_3d
|
---|
34 | #define GEN3D_FFTW_DESTROY_PLAN fftwf_destroy_plan
|
---|
35 | #define GEN3D_FFTW_EXECUTE fftwf_execute
|
---|
36 | #else
|
---|
37 | #define GEN3D_FFTW_INIT_THREADS fftw_init_threads
|
---|
38 | #define GEN3D_FFTW_CLEANUP_THREADS fftw_cleanup_threads
|
---|
39 | #define GEN3D_FFTW_PLAN_WITH_NTHREADS fftw_plan_with_nthreads
|
---|
40 | #define GEN3D_FFTW_PLAN_DFT_R2C_3D fftw_plan_dft_r2c_3d
|
---|
41 | #define GEN3D_FFTW_PLAN_DFT_C2R_3D fftw_plan_dft_c2r_3d
|
---|
42 | #define GEN3D_FFTW_DESTROY_PLAN fftw_destroy_plan
|
---|
43 | #define GEN3D_FFTW_EXECUTE fftw_execute
|
---|
44 | #endif
|
---|
45 |
|
---|
46 | #define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
|
---|
47 |
|
---|
48 | namespace SOPHYA {
|
---|
49 |
|
---|
50 | //-------------------------------------------------------
|
---|
51 | GeneFluct3D::GeneFluct3D(long nx,long ny,long nz,double dx,double dy,double dz
|
---|
52 | ,unsigned short nthread,int lp)
|
---|
53 | {
|
---|
54 | init_default();
|
---|
55 |
|
---|
56 | lp_ = lp;
|
---|
57 | nthread_ = nthread;
|
---|
58 |
|
---|
59 | setsize(nx,ny,nz,dx,dy,dz);
|
---|
60 | setalloc();
|
---|
61 | setpointers(false);
|
---|
62 | init_fftw();
|
---|
63 | }
|
---|
64 |
|
---|
65 | GeneFluct3D::GeneFluct3D(unsigned short nthread)
|
---|
66 | {
|
---|
67 | init_default();
|
---|
68 | setsize(2,2,2,1.,1.,1.);
|
---|
69 | nthread_ = nthread;
|
---|
70 | setalloc();
|
---|
71 | setpointers(false);
|
---|
72 | init_fftw();
|
---|
73 | }
|
---|
74 |
|
---|
75 | GeneFluct3D::~GeneFluct3D(void)
|
---|
76 | {
|
---|
77 | delete_fftw();
|
---|
78 | }
|
---|
79 |
|
---|
80 | //-------------------------------------------------------
|
---|
81 | void GeneFluct3D::init_default(void)
|
---|
82 | {
|
---|
83 | Nx_ = Ny_ = Nz_ = 0;
|
---|
84 | is_set_fft_plan = false;
|
---|
85 | nthread_ = 0;
|
---|
86 | lp_ = 0;
|
---|
87 | array_allocated_ = false; array_type = 0;
|
---|
88 | cosmo_ = NULL;
|
---|
89 | growth_ = NULL;
|
---|
90 | good_dzinc_ = 0.01;
|
---|
91 | compute_pk_redsh_ref_ = -999.;
|
---|
92 | redsh_ref_ = -999.;
|
---|
93 | kredsh_ref_ = 0.;
|
---|
94 | dred_ref_ = -999.;
|
---|
95 | h_ref_ = -999.;
|
---|
96 | loscom_ref_ = -999.;
|
---|
97 | dtrc_ref_ = dlum_ref_ = dang_ref_ = -999.;
|
---|
98 | nu_ref_ = dnu_ref_ = -999.;
|
---|
99 | loscom_min_ = loscom_max_ = -999.;
|
---|
100 | loscom2zred_min_ = loscom2zred_max_ = 0.;
|
---|
101 | xobs_[0] = xobs_[1] = xobs_[2] = 0.;
|
---|
102 | zred_.resize(0);
|
---|
103 | loscom_.resize(0);
|
---|
104 | loscom2zred_.resize(0);
|
---|
105 | }
|
---|
106 |
|
---|
107 | void GeneFluct3D::setsize(long nx,long ny,long nz,double dx,double dy,double dz)
|
---|
108 | {
|
---|
109 | if(lp_>1) cout<<"--- GeneFluct3D::setsize: N="<<nx<<","<<ny<<","<<nz
|
---|
110 | <<" D="<<dx<<","<<dy<<","<<dz<<endl;
|
---|
111 | if(nx<=0 || dx<=0.) {
|
---|
112 | const char *bla = "GeneFluct3D::setsize_Error: bad value(s) for nn/dx";
|
---|
113 | cout<<bla<<endl; throw ParmError(bla);
|
---|
114 | }
|
---|
115 |
|
---|
116 | // Les tailles des tableaux
|
---|
117 | Nx_ = nx;
|
---|
118 | Ny_ = ny; if(Ny_ <= 0) Ny_ = Nx_;
|
---|
119 | Nz_ = nz; if(Nz_ <= 0) Nz_ = Nx_;
|
---|
120 | N_.resize(0); N_.push_back(Nx_); N_.push_back(Ny_); N_.push_back(Nz_);
|
---|
121 | NRtot_ = (int_8)Nx_*(int_8)Ny_*(int_8)Nz_; // nombre de pixels dans le survey
|
---|
122 | NCz_ = Nz_/2 +1;
|
---|
123 | NFcoef_ = (int_8)Nx_*(int_8)Ny_*(int_8)NCz_; // nombre de coeff de Fourier dans le survey
|
---|
124 | NTz_ = 2*NCz_;
|
---|
125 |
|
---|
126 | // le pas dans l'espace (Mpc)
|
---|
127 | Dx_ = dx;
|
---|
128 | Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
|
---|
129 | Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
|
---|
130 | D_.resize(0); D_.push_back(Dx_); D_.push_back(Dy_); D_.push_back(Dz_);
|
---|
131 | dVol_ = Dx_*Dy_*Dz_;
|
---|
132 | Vol_ = (Nx_*Dx_)*(Ny_*Dy_)*(Nz_*Dz_);
|
---|
133 |
|
---|
134 | // Le pas dans l'espace de Fourier (Mpc^-1)
|
---|
135 | Dkx_ = 2.*M_PI/(Nx_*Dx_);
|
---|
136 | Dky_ = 2.*M_PI/(Ny_*Dy_);
|
---|
137 | Dkz_ = 2.*M_PI/(Nz_*Dz_);
|
---|
138 | Dk_.resize(0); Dk_.push_back(Dkx_); Dk_.push_back(Dky_); Dk_.push_back(Dkz_);
|
---|
139 | Dk3_ = Dkx_*Dky_*Dkz_;
|
---|
140 |
|
---|
141 | // La frequence de Nyquist en k (Mpc^-1)
|
---|
142 | Knyqx_ = M_PI/Dx_;
|
---|
143 | Knyqy_ = M_PI/Dy_;
|
---|
144 | Knyqz_ = M_PI/Dz_;
|
---|
145 | Knyq_.resize(0); Knyq_.push_back(Knyqx_); Knyq_.push_back(Knyqy_); Knyq_.push_back(Knyqz_);
|
---|
146 | }
|
---|
147 |
|
---|
148 | void GeneFluct3D::setalloc(void)
|
---|
149 | {
|
---|
150 | #if defined(GEN3D_FLOAT)
|
---|
151 | if(lp_>1) cout<<"--- GeneFluct3D::setalloc FLOAT ---"<<endl;
|
---|
152 | #else
|
---|
153 | if(lp_>1) cout<<"--- GeneFluct3D::setalloc DOUBLE ---"<<endl;
|
---|
154 | #endif
|
---|
155 | // Dimensionnement du tableau complex<r_8>
|
---|
156 | // ATTENTION: TArray adresse en memoire a l'envers du C
|
---|
157 | // Tarray(n1,n2,n3) == Carray[n3][n2][n1]
|
---|
158 | sa_size_t SzK_[3] = {NCz_,Ny_,Nx_}; // a l'envers
|
---|
159 | try {
|
---|
160 | T_.ReSize(3,SzK_);
|
---|
161 | array_allocated_ = true; array_type=0;
|
---|
162 | if(lp_>1) cout<<" allocating: "<<T_.Size()*sizeof(complex<GEN3D_TYPE>)/1.e6<<" Mo"<<endl;
|
---|
163 | } catch (...) {
|
---|
164 | cout<<"GeneFluct3D::setalloc_Error: Problem allocating T_"<<endl;
|
---|
165 | }
|
---|
166 | }
|
---|
167 |
|
---|
168 | void GeneFluct3D::setpointers(bool from_real)
|
---|
169 | {
|
---|
170 | if(lp_>1) cout<<"--- GeneFluct3D::setpointers ---"<<endl;
|
---|
171 | if(from_real) T_ = ArrCastR2C(R_);
|
---|
172 | else R_ = ArrCastC2R(T_);
|
---|
173 | // On remplit les pointeurs
|
---|
174 | fdata_ = (GEN3D_FFTW_COMPLEX *) (&T_(0,0,0));
|
---|
175 | data_ = (GEN3D_TYPE *) (&R_(0,0,0));
|
---|
176 | }
|
---|
177 |
|
---|
178 | void GeneFluct3D::init_fftw(void)
|
---|
179 | {
|
---|
180 | if( is_set_fft_plan ) delete_fftw();
|
---|
181 |
|
---|
182 | // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
|
---|
183 | if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl;
|
---|
184 | #ifdef WITH_FFTW_THREAD
|
---|
185 | if(nthread_>0) {
|
---|
186 | cout<<"...Computing with "<<nthread_<<" threads"<<endl;
|
---|
187 | GEN3D_FFTW_INIT_THREADS();
|
---|
188 | GEN3D_FFTW_PLAN_WITH_NTHREADS(nthread_);
|
---|
189 | }
|
---|
190 | #endif
|
---|
191 | if(lp_>1) cout<<"...forward plan"<<endl;
|
---|
192 | pf_ = GEN3D_FFTW_PLAN_DFT_R2C_3D(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
|
---|
193 | if(lp_>1) cout<<"...backward plan"<<endl;
|
---|
194 | pb_ = GEN3D_FFTW_PLAN_DFT_C2R_3D(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
|
---|
195 | is_set_fft_plan = true;
|
---|
196 | }
|
---|
197 |
|
---|
198 | void GeneFluct3D::delete_fftw(void)
|
---|
199 | {
|
---|
200 | if( !is_set_fft_plan ) return;
|
---|
201 | GEN3D_FFTW_DESTROY_PLAN(pf_);
|
---|
202 | GEN3D_FFTW_DESTROY_PLAN(pb_);
|
---|
203 | #ifdef WITH_FFTW_THREAD
|
---|
204 | if(nthread_>0) GEN3D_FFTW_CLEANUP_THREADS();
|
---|
205 | #endif
|
---|
206 | is_set_fft_plan = false;
|
---|
207 | }
|
---|
208 |
|
---|
209 | void GeneFluct3D::check_array_alloc(void)
|
---|
210 | // Pour tester si le tableau T_ est alloue
|
---|
211 | {
|
---|
212 | if(array_allocated_) return;
|
---|
213 | char bla[90];
|
---|
214 | sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
|
---|
215 | cout<<bla<<endl; throw ParmError(bla);
|
---|
216 | }
|
---|
217 |
|
---|
218 | //-------------------------------------------------------
|
---|
219 | void GeneFluct3D::SetObservator(double redshref,double kredshref)
|
---|
220 | // L'observateur est au redshift z=0
|
---|
221 | // est situe sur la "perpendiculaire" a la face x,y
|
---|
222 | // issue du centre de cette face
|
---|
223 | // Il faut positionner le cube sur l'axe des z cad des redshifts:
|
---|
224 | // redshref = redshift de reference
|
---|
225 | // Si redshref<0 alors redshref=0
|
---|
226 | // kredshref = indice (en double) correspondant a ce redshift
|
---|
227 | // Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
|
---|
228 | // Exemple: redshref=1.5 kredshref=250.75
|
---|
229 | // -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
|
---|
230 | {
|
---|
231 | if(redshref<0.) redshref = 0.;
|
---|
232 | if(kredshref<0.) {
|
---|
233 | if(Nz_<=0) {
|
---|
234 | const char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 define cube geometry first";
|
---|
235 | cout<<bla<<endl; throw ParmError(bla);
|
---|
236 | }
|
---|
237 | kredshref = Nz_/2.;
|
---|
238 | }
|
---|
239 | redsh_ref_ = redshref;
|
---|
240 | kredsh_ref_ = kredshref;
|
---|
241 | if(lp_>0)
|
---|
242 | cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
|
---|
243 | }
|
---|
244 |
|
---|
245 | void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
|
---|
246 | {
|
---|
247 | cosmo_ = &cosmo;
|
---|
248 | if(lp_>1) cosmo_->Print();
|
---|
249 | }
|
---|
250 |
|
---|
251 | void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
|
---|
252 | {
|
---|
253 | growth_ = &growth;
|
---|
254 | }
|
---|
255 |
|
---|
256 | long GeneFluct3D::LosComRedshift(double zinc,long npoints)
|
---|
257 | // Given a position of the cube relative to the observer
|
---|
258 | // and a cosmology
|
---|
259 | // (SetObservator() and SetCosmology() should have been called !)
|
---|
260 | // This routine filled:
|
---|
261 | // the vector "zred_" of scanned redshift (by zinc increments)
|
---|
262 | // the vector "loscom_" of corresponding los comoving distance
|
---|
263 | // -- Input:
|
---|
264 | // zinc : redshift increment for computation
|
---|
265 | // npoints : number of points required for inverting loscom -> zred
|
---|
266 | //
|
---|
267 | {
|
---|
268 | if(zinc<=0.) zinc = good_dzinc_;
|
---|
269 | if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
|
---|
270 |
|
---|
271 | if(cosmo_ == NULL || growth_==NULL || redsh_ref_<0.) {
|
---|
272 | const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator, Cosmology and Growth first";
|
---|
273 | cout<<bla<<endl; throw ParmError(bla);
|
---|
274 | }
|
---|
275 |
|
---|
276 | // La distance angulaire/luminosite/Dnu au pixel de reference
|
---|
277 | dred_ref_ = Dz_/(cosmo_->Dhubble()/cosmo_->E(redsh_ref_));
|
---|
278 | loscom_ref_ = cosmo_->Dloscom(redsh_ref_);
|
---|
279 | dtrc_ref_ = cosmo_->Dtrcom(redsh_ref_);
|
---|
280 | dlum_ref_ = cosmo_->Dlum(redsh_ref_);
|
---|
281 | dang_ref_ = cosmo_->Dang(redsh_ref_);
|
---|
282 | h_ref_ = cosmo_->H(redsh_ref_);
|
---|
283 | growth_ref_ = (*growth_)(redsh_ref_);
|
---|
284 | dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
|
---|
285 | nu_ref_ = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
|
---|
286 | dnu_ref_ = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
|
---|
287 | if(lp_>0) {
|
---|
288 | cout<<"...reference pixel redshref="<<redsh_ref_
|
---|
289 | <<", dredref="<<dred_ref_
|
---|
290 | <<", nuref="<<nu_ref_ <<" GHz"
|
---|
291 | <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
|
---|
292 | <<" H="<<h_ref_<<" km/s/Mpc"
|
---|
293 | <<", D="<<growth_ref_
|
---|
294 | <<", dD/dz="<<dsdz_growth_ref_<<endl
|
---|
295 | <<" dlosc="<<loscom_ref_<<" Mpc com"
|
---|
296 | <<", dtrc="<<dtrc_ref_<<" Mpc com"
|
---|
297 | <<", dlum="<<dlum_ref_<<" Mpc"
|
---|
298 | <<", dang="<<dang_ref_<<" Mpc"<<endl;
|
---|
299 | }
|
---|
300 |
|
---|
301 | // On calcule les coordonnees de l'observateur dans le repere du cube
|
---|
302 | // cad dans le repere ou l'origine est au centre du pixel i=j=l=0.
|
---|
303 | // L'observateur est sur un axe centre sur le milieu de la face Oxy
|
---|
304 | xobs_[0] = Nx_/2.*Dx_;
|
---|
305 | xobs_[1] = Ny_/2.*Dy_;
|
---|
306 | xobs_[2] = kredsh_ref_*Dz_ - loscom_ref_;
|
---|
307 |
|
---|
308 | // L'observateur est-il dans le cube?
|
---|
309 | bool obs_in_cube = false;
|
---|
310 | if(xobs_[2]>=0. && xobs_[2]<=Nz_*Dz_) obs_in_cube = true;
|
---|
311 |
|
---|
312 | // Find MINIMUM los com distance to the observer:
|
---|
313 | // c'est le centre de la face a k=0
|
---|
314 | // (ou zero si l'observateur est dans le cube)
|
---|
315 | loscom_min_ = 0.;
|
---|
316 | if(!obs_in_cube) loscom_min_ = -xobs_[2];
|
---|
317 |
|
---|
318 | // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
|
---|
319 | if(loscom_min_<=1.e-50)
|
---|
320 | for(int i=0;i<10;i++)
|
---|
321 | cout<<"ATTENTION TOUTES LES PARTIES DU CODE NE MARCHENT PAS POUR UN OBSERVATEUR DANS LE CUBE"<<endl;
|
---|
322 | // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
|
---|
323 |
|
---|
324 |
|
---|
325 | // Find MAXIMUM los com distance to the observer:
|
---|
326 | // ou que soit positionne l'observateur, la distance
|
---|
327 | // maximal est sur un des coins du cube
|
---|
328 | loscom_max_ = 0.;
|
---|
329 | for(long i=0;i<=1;i++) {
|
---|
330 | double dx2 = DXcom(i*(Nx_-1)); dx2 *= dx2;
|
---|
331 | for(long j=0;j<=1;j++) {
|
---|
332 | double dy2 = DYcom(j*(Ny_-1)); dy2 *= dy2;
|
---|
333 | for(long k=0;k<=1;k++) {
|
---|
334 | double dz2 = DZcom(k*(Nz_-1)); dz2 *= dz2;
|
---|
335 | dz2 = sqrt(dx2+dy2+dz2);
|
---|
336 | if(dz2>loscom_max_) loscom_max_ = dz2;
|
---|
337 | }
|
---|
338 | }
|
---|
339 | }
|
---|
340 | if(lp_>0) {
|
---|
341 | cout<<"...zref="<<redsh_ref_<<" kzref="<<kredsh_ref_<<" losref="<<loscom_ref_<<" Mpc\n"
|
---|
342 | <<" xobs="<<xobs_[0]<<" , "<<xobs_[1]<<" , "<<xobs_[2]<<" Mpc "
|
---|
343 | <<" in_cube="<<obs_in_cube
|
---|
344 | <<" loscom_min="<<loscom_min_<<" loscom_max="<<loscom_max_<<" Mpc (com)"<<endl;
|
---|
345 | }
|
---|
346 |
|
---|
347 | // Fill the corresponding vectors for loscom and zred
|
---|
348 | // Be shure to have one dlc<loscom_min and one dlc>loscom_max
|
---|
349 | double zmin = 0., dlcmin=0.;
|
---|
350 | while(1) {
|
---|
351 | if(lp_>0)
|
---|
352 | cout<<"...Filling zred starting at zmin="<<zmin<<" with zinc="<<zinc
|
---|
353 | <<", loscom_min-max=["<<loscom_min_<<","<<loscom_max_<<"]"<<endl;
|
---|
354 | zred_.resize(0); loscom_.resize(0);
|
---|
355 | for(double z=zmin; ; z+=zinc) {
|
---|
356 | double dlc = cosmo_->Dloscom(z);
|
---|
357 | if(dlc<loscom_min_) {
|
---|
358 | dlcmin = dlc;
|
---|
359 | zmin = z;
|
---|
360 | zred_.resize(0); loscom_.resize(0);
|
---|
361 | }
|
---|
362 | zred_.push_back(z);
|
---|
363 | loscom_.push_back(dlc);
|
---|
364 | z += zinc;
|
---|
365 | if(dlc>loscom_max_) {
|
---|
366 | if(lp_>0)
|
---|
367 | cout<<" Min: z="<<zmin<<" dlc="<<dlcmin<<", Max: z="<<z<<" dlc="<<dlc<<endl;
|
---|
368 | break; // on sort apres avoir stoque un dlc>dlcmax
|
---|
369 | }
|
---|
370 | }
|
---|
371 | if(zred_.size()>=10) break;
|
---|
372 | zinc /= 10.;
|
---|
373 | cout<<" not enough points ("<<zred_.size()<<") for zref, retry with zinc="<<zinc<<endl;
|
---|
374 | }
|
---|
375 |
|
---|
376 | if(lp_>0) {
|
---|
377 | long n = zred_.size();
|
---|
378 | cout<<"...zred/loscom tables[zinc="<<zinc<<"]: n="<<n;
|
---|
379 | if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
|
---|
380 | if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
|
---|
381 | cout<<endl;
|
---|
382 | }
|
---|
383 |
|
---|
384 | // Compute the parameters and tables needed for inversion loscom->zred
|
---|
385 | if(npoints<3) npoints = zred_.size();
|
---|
386 | InverseFunc invfun(zred_,loscom_);
|
---|
387 | invfun.ComputeParab(npoints,loscom2zred_);
|
---|
388 | loscom2zred_min_ = invfun.YMin();
|
---|
389 | loscom2zred_max_ = invfun.YMax();
|
---|
390 |
|
---|
391 | if(lp_>0) {
|
---|
392 | long n = loscom2zred_.size();
|
---|
393 | cout<<"...loscom -> zred[npoints="<<npoints<<"]: n="<<n
|
---|
394 | <<" los_min="<<loscom2zred_min_
|
---|
395 | <<" los_max="<<loscom2zred_max_
|
---|
396 | <<" -> zred=[";
|
---|
397 | if(n>0) cout<<loscom2zred_[0];
|
---|
398 | cout<<",";
|
---|
399 | if(n>1) cout<<loscom2zred_[n-1];
|
---|
400 | cout<<"]"<<endl;
|
---|
401 | if(lp_>1 && n>0)
|
---|
402 | for(int i=0;i<n;i++)
|
---|
403 | if(i<2 || abs(i-n/2)<2 || i>=n-2)
|
---|
404 | cout<<" i="<<i
|
---|
405 | <<" d="<<loscom2zred_min_+i*(loscom2zred_max_-loscom2zred_min_)/(n-1.)
|
---|
406 | <<" Mpc z="<<loscom2zred_[i]<<endl;
|
---|
407 | }
|
---|
408 |
|
---|
409 |
|
---|
410 | // Compute the table for D'(z)/D(z) where D'(z)=dD/dEta (Eta is conformal time)
|
---|
411 | zredmin_dpsd_ = zred_[0];
|
---|
412 | zredmax_dpsd_ = zred_[zred_.size()-1];
|
---|
413 | long nptd = long(sqrt(Nx_*Nx_ + Ny_*Ny_ +Nz_*Nz_));
|
---|
414 | if(nptd<10) nptd = 10;
|
---|
415 | good_dzinc_ = (zredmax_dpsd_ - zredmin_dpsd_)/nptd;
|
---|
416 | if(lp_>0) cout<<"...good_dzinc changed to "<<good_dzinc_<<endl;
|
---|
417 | dpsdfrzred_.resize(0);
|
---|
418 | double dz = (zredmax_dpsd_ - zredmin_dpsd_)/(nptd-1);
|
---|
419 | int nmod = nptd/5; if(nmod==0) nmod = 1;
|
---|
420 | if(lp_>0) cout<<"...Compute table D'/D on "<<nptd<<" pts for z["
|
---|
421 | <<zredmin_dpsd_<<","<<zredmax_dpsd_<<"]"<<endl;
|
---|
422 | for(int i=0;i<nptd;i++) {
|
---|
423 | double z = zredmin_dpsd_ + i*dz;
|
---|
424 | double om = cosmo_->Omatter(z);
|
---|
425 | double h = cosmo_->H(z);
|
---|
426 | double dsdz = growth_->DsDz(z,good_dzinc_);
|
---|
427 | double d = (*growth_)(z);
|
---|
428 | double v = -h * (1.+z) * dsdz / d;
|
---|
429 | dpsdfrzred_.push_back(v);
|
---|
430 | if(lp_ && (i%nmod==0 || i==nptd-1))
|
---|
431 | cout<<" z="<<z<<" H="<<h<<" Beta="<<-(1.+z)*dsdz/d
|
---|
432 | <<" (Om^0.6="<<pow(om,0.6)
|
---|
433 | <<") -H*(1+z)*D'/D="<<v<<" km/s/Mpc"<<endl;
|
---|
434 | }
|
---|
435 |
|
---|
436 | return zred_.size();
|
---|
437 | }
|
---|
438 |
|
---|
439 | //-------------------------------------------------------
|
---|
440 | void GeneFluct3D::WriteFits(string cfname,int bitpix)
|
---|
441 | {
|
---|
442 | cout<<"--- GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
|
---|
443 | try {
|
---|
444 | FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
|
---|
445 | fwrt.WriteKey("NX",Nx_," axe transverse 1"); // NAXIS3
|
---|
446 | fwrt.WriteKey("NY",Ny_," axe transverse 2"); // NAXIS2
|
---|
447 | fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)"); // NAXIS1
|
---|
448 | fwrt.WriteKey("DX",Dx_," Mpc");
|
---|
449 | fwrt.WriteKey("DY",Dy_," Mpc");
|
---|
450 | fwrt.WriteKey("DZ",Dz_," Mpc");
|
---|
451 | fwrt.WriteKey("DKX",Dkx_," Mpc^-1");
|
---|
452 | fwrt.WriteKey("DKY",Dky_," Mpc^-1");
|
---|
453 | fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
|
---|
454 | fwrt.WriteKey("ZREFPK",compute_pk_redsh_ref_," Pk computed redshift");
|
---|
455 | fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
|
---|
456 | fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
|
---|
457 | fwrt.WriteKey("HREF",h_ref_," ref H(z)");
|
---|
458 | fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
|
---|
459 | fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
|
---|
460 | fwrt.Write(R_);
|
---|
461 | } catch (PThrowable & exc) {
|
---|
462 | cout<<"Exception : "<<(string)typeid(exc).name()
|
---|
463 | <<" - Msg= "<<exc.Msg()<<endl;
|
---|
464 | return;
|
---|
465 | } catch (...) {
|
---|
466 | cout<<" some other exception was caught !"<<endl;
|
---|
467 | return;
|
---|
468 | }
|
---|
469 | }
|
---|
470 |
|
---|
471 | void GeneFluct3D::ReadFits(string cfname)
|
---|
472 | {
|
---|
473 | cout<<"--- GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
|
---|
474 | try {
|
---|
475 | FitsImg3DRead fimg(cfname.c_str(),0,5);
|
---|
476 | fimg.Read(R_);
|
---|
477 | long nx = fimg.ReadKeyL("NX"); // NAXIS3
|
---|
478 | long ny = fimg.ReadKeyL("NY"); // NAXIS2
|
---|
479 | long nz = fimg.ReadKeyL("NZ"); // NAXIS1
|
---|
480 | double dx = fimg.ReadKey("DX");
|
---|
481 | double dy = fimg.ReadKey("DY");
|
---|
482 | double dz = fimg.ReadKey("DZ");
|
---|
483 | double pkzref = fimg.ReadKey("ZREFPK");
|
---|
484 | double zref = fimg.ReadKey("ZREF");
|
---|
485 | double kzref = fimg.ReadKey("KZREF");
|
---|
486 | setsize(nx,ny,nz,dx,dy,dz);
|
---|
487 | setpointers(true);
|
---|
488 | init_fftw();
|
---|
489 | SetObservator(zref,kzref);
|
---|
490 | array_allocated_ = true;
|
---|
491 | compute_pk_redsh_ref_ = pkzref;
|
---|
492 | } catch (PThrowable & exc) {
|
---|
493 | cout<<"Exception : "<<(string)typeid(exc).name()
|
---|
494 | <<" - Msg= "<<exc.Msg()<<endl;
|
---|
495 | return;
|
---|
496 | } catch (...) {
|
---|
497 | cout<<" some other exception was caught !"<<endl;
|
---|
498 | return;
|
---|
499 | }
|
---|
500 | }
|
---|
501 |
|
---|
502 | void GeneFluct3D::WritePPF(string cfname,bool write_real)
|
---|
503 | // On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> >
|
---|
504 | {
|
---|
505 | cout<<"--- GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
|
---|
506 | try {
|
---|
507 | R_.Info()["NX"] = (int_8)Nx_;
|
---|
508 | R_.Info()["NY"] = (int_8)Ny_;
|
---|
509 | R_.Info()["NZ"] = (int_8)Nz_;
|
---|
510 | R_.Info()["DX"] = (r_8)Dx_;
|
---|
511 | R_.Info()["DY"] = (r_8)Dy_;
|
---|
512 | R_.Info()["DZ"] = (r_8)Dz_;
|
---|
513 | R_.Info()["ZREFPK"] = (r_8)compute_pk_redsh_ref_;
|
---|
514 | R_.Info()["ZREF"] = (r_8)redsh_ref_;
|
---|
515 | R_.Info()["KZREF"] = (r_8)kredsh_ref_;
|
---|
516 | R_.Info()["HREF"] = (r_8)h_ref_;
|
---|
517 | R_.Info()["Growth"] = (r_8)Dref();
|
---|
518 | R_.Info()["GrowthPk"] = (r_8)DrefPk();
|
---|
519 | POutPersist pos(cfname.c_str());
|
---|
520 | if(write_real) pos << PPFNameTag("rgen") << R_;
|
---|
521 | else pos << PPFNameTag("pkgen") << T_;
|
---|
522 | } catch (PThrowable & exc) {
|
---|
523 | cout<<"Exception : "<<(string)typeid(exc).name()
|
---|
524 | <<" - Msg= "<<exc.Msg()<<endl;
|
---|
525 | return;
|
---|
526 | } catch (...) {
|
---|
527 | cout<<" some other exception was caught !"<<endl;
|
---|
528 | return;
|
---|
529 | }
|
---|
530 | }
|
---|
531 |
|
---|
532 | void GeneFluct3D::ReadPPF(string cfname)
|
---|
533 | {
|
---|
534 | cout<<"--- GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
|
---|
535 | try {
|
---|
536 | bool from_real = true;
|
---|
537 | PInPersist pis(cfname.c_str());
|
---|
538 | string name_tag_k = "pkgen";
|
---|
539 | bool found_tag_k = pis.GotoNameTag("pkgen");
|
---|
540 | if(found_tag_k) {
|
---|
541 | cout<<" ...reading spectrum into TArray<complex <r_8> >"<<endl;
|
---|
542 | pis >> PPFNameTag("pkgen") >> T_;
|
---|
543 | from_real = false;
|
---|
544 | } else {
|
---|
545 | cout<<" ...reading space into TArray<r_8>"<<endl;
|
---|
546 | pis >> PPFNameTag("rgen") >> R_;
|
---|
547 | }
|
---|
548 | setpointers(from_real); // a mettre ici pour relire les DVInfo
|
---|
549 | int_8 nx = R_.Info()["NX"];
|
---|
550 | int_8 ny = R_.Info()["NY"];
|
---|
551 | int_8 nz = R_.Info()["NZ"];
|
---|
552 | r_8 dx = R_.Info()["DX"];
|
---|
553 | r_8 dy = R_.Info()["DY"];
|
---|
554 | r_8 dz = R_.Info()["DZ"];
|
---|
555 | r_8 pkzref = R_.Info()["ZREFPK"];
|
---|
556 | r_8 zref = R_.Info()["ZREF"];
|
---|
557 | r_8 kzref = R_.Info()["KZREF"];
|
---|
558 | setsize(nx,ny,nz,dx,dy,dz);
|
---|
559 | init_fftw();
|
---|
560 | SetObservator(zref,kzref);
|
---|
561 | array_allocated_ = true;
|
---|
562 | compute_pk_redsh_ref_ = pkzref;
|
---|
563 | } catch (PThrowable & exc) {
|
---|
564 | cout<<"Exception : "<<(string)typeid(exc).name()
|
---|
565 | <<" - Msg= "<<exc.Msg()<<endl;
|
---|
566 | return;
|
---|
567 | } catch (...) {
|
---|
568 | cout<<" some other exception was caught !"<<endl;
|
---|
569 | return;
|
---|
570 | }
|
---|
571 | }
|
---|
572 |
|
---|
573 | void GeneFluct3D::WriteSlicePPF(string cfname)
|
---|
574 | // On ecrit 3 tranches du cube selon chaque axe
|
---|
575 | {
|
---|
576 | cout<<"--- GeneFluct3D::WriteSlicePPF: Writing Cube Slices "<<cfname<<endl;
|
---|
577 | try {
|
---|
578 |
|
---|
579 | POutPersist pos(cfname.c_str());
|
---|
580 | TMatrix<r_4> S;
|
---|
581 | char str[16];
|
---|
582 | long i,j,l;
|
---|
583 |
|
---|
584 | // Tranches en Z
|
---|
585 | for(int s=0;s<3;s++) {
|
---|
586 | S.ReSize(Nx_,Ny_);
|
---|
587 | if(s==0) l=0; else if(s==1) l=(Nz_+1)/2; else l=Nz_-1;
|
---|
588 | sprintf(str,"z%ld",l);
|
---|
589 | for(i=0;i<Nx_;i++) for(j=0;j<Ny_;j++) S(i,j)=data_[IndexR(i,j,l)];
|
---|
590 | pos<<PPFNameTag(str)<<S; S.RenewObjId();
|
---|
591 | }
|
---|
592 |
|
---|
593 | // Tranches en Y
|
---|
594 | for(int s=0;s<3;s++) {
|
---|
595 | S.ReSize(Nz_,Nx_);
|
---|
596 | if(s==0) j=0; else if(s==1) j=(Ny_+1)/2; else j=Ny_-1;
|
---|
597 | sprintf(str,"y%ld",j);
|
---|
598 | for(i=0;i<Nx_;i++) for(l=0;l<Nz_;l++) S(l,i)=data_[IndexR(i,j,l)];
|
---|
599 | pos<<PPFNameTag(str)<<S; S.RenewObjId();
|
---|
600 | }
|
---|
601 |
|
---|
602 | // Tranches en X
|
---|
603 | for(int s=0;s<3;s++) {
|
---|
604 | S.ReSize(Nz_,Ny_);
|
---|
605 | if(s==0) i=0; else if(s==1) i=(Nx_+1)/2; else i=Nx_-1;
|
---|
606 | sprintf(str,"x%ld",i);
|
---|
607 | for(j=0;j<Ny_;j++) for(l=0;l<Nz_;l++) S(l,j)=data_[IndexR(i,j,l)];
|
---|
608 | pos<<PPFNameTag(str)<<S; S.RenewObjId();
|
---|
609 | }
|
---|
610 |
|
---|
611 | } catch (PThrowable & exc) {
|
---|
612 | cout<<"Exception : "<<(string)typeid(exc).name()
|
---|
613 | <<" - Msg= "<<exc.Msg()<<endl;
|
---|
614 | return;
|
---|
615 | } catch (...) {
|
---|
616 | cout<<" some other exception was caught !"<<endl;
|
---|
617 | return;
|
---|
618 | }
|
---|
619 | }
|
---|
620 |
|
---|
621 | //-------------------------------------------------------
|
---|
622 | void GeneFluct3D::NTupleCheck(POutPersist &pos,string ntname,unsigned long nent)
|
---|
623 | // Remplit le NTuple "ntname" avec "nent" valeurs du cube (reel ou complex) et l'ecrit dans "pos"
|
---|
624 | {
|
---|
625 | if(ntname.size()<=0 || nent==0) return;
|
---|
626 | int nvar = 0;
|
---|
627 | if(array_type==1) nvar = 3;
|
---|
628 | else if(array_type==2) nvar = 4;
|
---|
629 | else return;
|
---|
630 | const char *vname[4] = {"t","z","re","im"};
|
---|
631 | float xnt[4];
|
---|
632 | NTuple nt(nvar,vname);
|
---|
633 |
|
---|
634 | if(array_type==1) {
|
---|
635 | unsigned long nmod = Nx_*Ny_*Nz_/nent; if(nmod==0) nmod=1;
|
---|
636 | unsigned long n=0;
|
---|
637 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
638 | if(n==nmod) {
|
---|
639 | int_8 ip = IndexR(i,j,l);
|
---|
640 | xnt[0]=sqrt(i*i+j*j); xnt[1]=l; xnt[2]=data_[ip];
|
---|
641 | nt.Fill(xnt);
|
---|
642 | n=0;
|
---|
643 | }
|
---|
644 | n++;
|
---|
645 | }
|
---|
646 | } else {
|
---|
647 | unsigned long nmod = Nx_*Ny_*NCz_/nent; if(nmod==0) nmod=1;
|
---|
648 | unsigned long n=0;
|
---|
649 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<NCz_;l++) {
|
---|
650 | if(n==nmod) {
|
---|
651 | xnt[0]=sqrt(i*i+j*j); xnt[1]=l; xnt[2]=T_(l,j,i).real(); xnt[3]=T_(l,j,i).imag();
|
---|
652 | nt.Fill(xnt);
|
---|
653 | n=0;
|
---|
654 | }
|
---|
655 | n++;
|
---|
656 | }
|
---|
657 | }
|
---|
658 |
|
---|
659 | pos.PutObject(nt,ntname);
|
---|
660 | }
|
---|
661 |
|
---|
662 | //-------------------------------------------------------
|
---|
663 | void GeneFluct3D::Print(void)
|
---|
664 | {
|
---|
665 | cout<<"GeneFluct3D(T_alloc="<<array_allocated_<<"):"<<endl;
|
---|
666 | cout<<"Space Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<Nz_<<" ("<<NTz_<<") size="
|
---|
667 | <<NRtot_<<endl;
|
---|
668 | cout<<" Resol: dx="<<Dx_<<" dy="<<Dy_<<" dz="<<Dz_<<" Mpc"
|
---|
669 | <<", dVol="<<dVol_<<", Vol="<<Vol_<<" Mpc^3"<<endl;
|
---|
670 | cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<" ncoeff="<<NFcoef_<<endl;
|
---|
671 | cout<<" Resol: dkx="<<Dkx_<<" dky="<<Dky_<<" dkz="<<Dkz_<<" Mpc^-1"
|
---|
672 | <<", Dk3="<<Dk3_<<" Mpc^-3"<<endl;
|
---|
673 | cout<<" (2Pi/k: "<<2.*M_PI/Dkx_<<" "<<2.*M_PI/Dky_<<" "<<2.*M_PI/Dkz_<<" Mpc)"<<endl;
|
---|
674 | cout<<" Nyquist: kx="<<Knyqx_<<" ky="<<Knyqy_<<" kz="<<Knyqz_<<" Mpc^-1"
|
---|
675 | <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
|
---|
676 | cout<<" (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
|
---|
677 | cout<<"Redshift "<<redsh_ref_<<" for z axe at k="<<kredsh_ref_<<endl;
|
---|
678 | }
|
---|
679 |
|
---|
680 | //-------------------------------------------------------
|
---|
681 | void GeneFluct3D::ComputeFourier0(PkSpectrum& pk_at_z)
|
---|
682 | // cf ComputeFourier() mais avec autre methode de realisation du spectre
|
---|
683 | // (attention on fait une fft pour realiser le spectre)
|
---|
684 | {
|
---|
685 | compute_pk_redsh_ref_ = pk_at_z.GetZ();
|
---|
686 | // --- realisation d'un tableau de tirage gaussiens
|
---|
687 | if(lp_>0) cout<<"--- ComputeFourier0 at z="<<compute_pk_redsh_ref_
|
---|
688 | <<": before gaussian filling ---"<<endl;
|
---|
689 | // On tient compte du pb de normalisation de FFTW3
|
---|
690 | double sntot = sqrt((double)NRtot_);
|
---|
691 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
692 | int_8 ip = IndexR(i,j,l);
|
---|
693 | data_[ip] = NorRand()/sntot;
|
---|
694 | }
|
---|
695 |
|
---|
696 | // --- realisation d'un tableau de tirage gaussiens
|
---|
697 | if(lp_>0) cout<<"...before fft real ---"<<endl;
|
---|
698 | GEN3D_FFTW_EXECUTE(pf_);
|
---|
699 |
|
---|
700 | // --- On remplit avec une realisation
|
---|
701 | if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
|
---|
702 | T_(0,0,0) = complex<GEN3D_TYPE>(0.); // on coupe le continue et on l'initialise
|
---|
703 | long lmod = Nx_/20; if(lmod<1) lmod=1;
|
---|
704 | for(long i=0;i<Nx_;i++) {
|
---|
705 | double kx = AbsKx(i); kx *= kx;
|
---|
706 | if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
|
---|
707 | for(long j=0;j<Ny_;j++) {
|
---|
708 | double ky = AbsKy(j); ky *= ky;
|
---|
709 | for(long l=0;l<NCz_;l++) {
|
---|
710 | double kz = AbsKz(l); kz *= kz;
|
---|
711 | if(i==0 && j==0 && l==0) continue; // Suppression du continu
|
---|
712 | double k = sqrt(kx+ky+kz);
|
---|
713 | // cf normalisation: Peacock, Cosmology, formule 16.38 p504
|
---|
714 | double pk = pk_at_z(k)/Vol_;
|
---|
715 | // ici pas de "/2" a cause de la remarque ci-dessus
|
---|
716 | T_(l,j,i) *= sqrt(pk);
|
---|
717 | }
|
---|
718 | }
|
---|
719 | }
|
---|
720 |
|
---|
721 | array_type = 2;
|
---|
722 |
|
---|
723 | if(lp_>0) cout<<"...computing power"<<endl;
|
---|
724 | double p = compute_power_carte();
|
---|
725 | if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;
|
---|
726 |
|
---|
727 | }
|
---|
728 |
|
---|
729 | //-------------------------------------------------------
|
---|
730 | void GeneFluct3D::ComputeFourier(PkSpectrum& pk_at_z)
|
---|
731 | // Calcule une realisation du spectre "pk_at_z"
|
---|
732 | // Attention: dans TArray le premier indice varie le + vite
|
---|
733 | // Explication normalisation: see Coles & Lucchin, Cosmology, p264-265
|
---|
734 | // FFTW3: on note N=Nx*Ny*Nz
|
---|
735 | // f --(FFT)--> F = TF(f) --(FFT^-1)--> fb = TF^-1(F) = TF^-1(TF(f))
|
---|
736 | // sum(f(x_i)^2) = S
|
---|
737 | // sum(F(nu_i)^2) = S*N
|
---|
738 | // sum(fb(x_i)^2) = S*N^2
|
---|
739 | {
|
---|
740 | // --- RaZ du tableau
|
---|
741 | T_ = complex<GEN3D_TYPE>(0.);
|
---|
742 | compute_pk_redsh_ref_ = pk_at_z.GetZ();
|
---|
743 |
|
---|
744 | // --- On remplit avec une realisation
|
---|
745 | if(lp_>0) cout<<"--- ComputeFourier at z="<<compute_pk_redsh_ref_<<" ---"<<endl;
|
---|
746 | long lmod = Nx_/20; if(lmod<1) lmod=1;
|
---|
747 | for(long i=0;i<Nx_;i++) {
|
---|
748 | double kx = AbsKx(i); kx *= kx;
|
---|
749 | if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
|
---|
750 | for(long j=0;j<Ny_;j++) {
|
---|
751 | double ky = AbsKy(j); ky *= ky;
|
---|
752 | for(long l=0;l<NCz_;l++) {
|
---|
753 | double kz = AbsKz(l); kz *= kz;
|
---|
754 | if(i==0 && j==0 && l==0) continue; // Suppression du continu
|
---|
755 | double k = sqrt(kx+ky+kz);
|
---|
756 | // cf normalisation: Peacock, Cosmology, formule 16.38 p504
|
---|
757 | double pk = pk_at_z(k)/Vol_;
|
---|
758 | // Explication de la division par 2: voir perandom.cc
|
---|
759 | // ou egalement Coles & Lucchin, Cosmology formula 13.7.2 p279
|
---|
760 | T_(l,j,i) = ComplexGaussianRand(sqrt(pk/2.));
|
---|
761 | }
|
---|
762 | }
|
---|
763 | }
|
---|
764 |
|
---|
765 | array_type = 2;
|
---|
766 | manage_coefficients(); // gros effet pour les spectres que l'on utilise !
|
---|
767 |
|
---|
768 | if(lp_>0) cout<<"...computing power"<<endl;
|
---|
769 | double p = compute_power_carte();
|
---|
770 | if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;
|
---|
771 |
|
---|
772 | }
|
---|
773 |
|
---|
774 | long GeneFluct3D::manage_coefficients(void)
|
---|
775 | // Take into account the real and complexe conjugate coefficients
|
---|
776 | // because we want a realization of a real data in real space
|
---|
777 | // On ecrit que: conj(P(k_x,k_y,k_z)) = P(-k_x,-k_y,-k_z)
|
---|
778 | // avec k_x = i, -k_x = N_x - i etc...
|
---|
779 | {
|
---|
780 | if(lp_>1) cout<<"...managing coefficients"<<endl;
|
---|
781 | check_array_alloc();
|
---|
782 |
|
---|
783 | // 1./ Le Continu et Nyquist sont reels
|
---|
784 | int_8 nreal = 0;
|
---|
785 | for(long kk=0;kk<2;kk++) {
|
---|
786 | long k=0; // continu
|
---|
787 | if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
|
---|
788 | for(long jj=0;jj<2;jj++) {
|
---|
789 | long j=0;
|
---|
790 | if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
|
---|
791 | for(long ii=0;ii<2;ii++) {
|
---|
792 | long i=0;
|
---|
793 | if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
|
---|
794 | int_8 ip = IndexC(i,j,k);
|
---|
795 | //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata_[ip][0]<<","<<fdata_[ip][1]<<")"<<endl;
|
---|
796 | fdata_[ip][1] = 0.; fdata_[ip][0] *= M_SQRT2;
|
---|
797 | nreal++;
|
---|
798 | }
|
---|
799 | }
|
---|
800 | }
|
---|
801 | if(lp_>1) cout<<"Number of forced real number ="<<nreal<<endl;
|
---|
802 |
|
---|
803 | // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist)
|
---|
804 |
|
---|
805 | // a./ les lignes et colonnes du continu et de nyquist
|
---|
806 | int_8 nconj1 = 0;
|
---|
807 | for(long kk=0;kk<2;kk++) {
|
---|
808 | long k=0; // continu
|
---|
809 | if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
|
---|
810 | for(long jj=0;jj<2;jj++) { // selon j
|
---|
811 | long j=0;
|
---|
812 | if(jj==1) {if(Ny_%2!=0) continue; else j = Ny_/2;}
|
---|
813 | for(long i=1;i<(Nx_+1)/2;i++) {
|
---|
814 | int_8 ip = IndexC(i,j,k);
|
---|
815 | int_8 ip1 = IndexC(Nx_-i,j,k);
|
---|
816 | fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
|
---|
817 | nconj1++;
|
---|
818 | }
|
---|
819 | }
|
---|
820 | for(long ii=0;ii<2;ii++) {
|
---|
821 | long i=0;
|
---|
822 | if(ii==1) {if(Nx_%2!=0) continue; else i = Nx_/2;}
|
---|
823 | for(long j=1;j<(Ny_+1)/2;j++) {
|
---|
824 | int_8 ip = IndexC(i,j,k);
|
---|
825 | int_8 ip1 = IndexC(i,Ny_-j,k);
|
---|
826 | fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
|
---|
827 | nconj1++;
|
---|
828 | }
|
---|
829 | }
|
---|
830 | }
|
---|
831 | if(lp_>1) cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;
|
---|
832 |
|
---|
833 | // b./ les lignes et colonnes hors continu et de nyquist
|
---|
834 | int_8 nconj2 = 0;
|
---|
835 | for(long kk=0;kk<2;kk++) {
|
---|
836 | long k=0; // continu
|
---|
837 | if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
|
---|
838 | for(long j=1;j<(Ny_+1)/2;j++) {
|
---|
839 | if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
|
---|
840 | for(long i=1;i<Nx_;i++) {
|
---|
841 | if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
|
---|
842 | int_8 ip = IndexC(i,j,k);
|
---|
843 | int_8 ip1 = IndexC(Nx_-i,Ny_-j,k);
|
---|
844 | fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
|
---|
845 | nconj2++;
|
---|
846 | }
|
---|
847 | }
|
---|
848 | }
|
---|
849 | if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
|
---|
850 |
|
---|
851 | if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(NFcoef_-nconj1-nconj2)-nreal<<endl;
|
---|
852 |
|
---|
853 | return nreal+nconj1+nconj2;
|
---|
854 | }
|
---|
855 |
|
---|
856 | double GeneFluct3D::compute_power_carte(void)
|
---|
857 | // Calcul la puissance de la realisation du spectre Pk
|
---|
858 | {
|
---|
859 | check_array_alloc();
|
---|
860 |
|
---|
861 | double s2 = 0.;
|
---|
862 | for(long l=0;l<NCz_;l++)
|
---|
863 | for(long j=0;j<Ny_;j++)
|
---|
864 | for(long i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));
|
---|
865 |
|
---|
866 | double s20 = 0.;
|
---|
867 | for(long j=0;j<Ny_;j++)
|
---|
868 | for(long i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));
|
---|
869 |
|
---|
870 | double s2n = 0.;
|
---|
871 | if(Nz_%2==0)
|
---|
872 | for(long j=0;j<Ny_;j++)
|
---|
873 | for(long i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));
|
---|
874 |
|
---|
875 | return 2.*s2 -s20 -s2n;
|
---|
876 | }
|
---|
877 |
|
---|
878 | //-------------------------------------------------------------------
|
---|
879 | void GeneFluct3D::FilterByPixel(void)
|
---|
880 | // Filtrage par la fonction fenetre du pixel (parallelepipede)
|
---|
881 | // TF = 1/(dx*dy*dz)*Int[{-dx/2,dx/2},{-dy/2,dy/2},{-dz/2,dz/2}]
|
---|
882 | // e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
|
---|
883 | // = 2/(k_x*dx) * sin(k_x*dx/2) * (idem y) * (idem z)
|
---|
884 | // Gestion divergence en 0: sin(y)/y = 1 - y^2/6*(1-y^2/20)
|
---|
885 | // avec y = k_x*dx/2
|
---|
886 | //--- Ex a 1D: TF(k) = sinc(k D/2) = sin(k D/2) / (k D/2)
|
---|
887 | // Le premier zero du sinc, cad k1 tq TF(k1)=0 est a: k1 D/2 = Pi -> k1 = 2Pi/D
|
---|
888 | // Avec un echantillonage de pas D, on a: k_nyq = Pi/D = k1/2
|
---|
889 | // donc pour la frequence de Nyquist: TF(k_nyq) = sinc(Pi/2) = 2/Pi ~= 0.6366
|
---|
890 | //---
|
---|
891 | {
|
---|
892 | if(lp_>0) cout<<"--- FilterByPixel ---"<<endl;
|
---|
893 | check_array_alloc();
|
---|
894 |
|
---|
895 | for(long i=0;i<Nx_;i++) {
|
---|
896 | double kx = AbsKx(i) *Dx_/2;
|
---|
897 | double pk_x = pixelfilter(kx);
|
---|
898 | for(long j=0;j<Ny_;j++) {
|
---|
899 | double ky = AbsKy(j) *Dy_/2;
|
---|
900 | double pk_y = pixelfilter(ky);
|
---|
901 | for(long l=0;l<NCz_;l++) {
|
---|
902 | double kz = AbsKz(l) *Dz_/2;
|
---|
903 | double pk_z = pixelfilter(kz);
|
---|
904 | T_(l,j,i) *= pk_x*pk_y*pk_z;
|
---|
905 | }
|
---|
906 | }
|
---|
907 | }
|
---|
908 |
|
---|
909 | }
|
---|
910 |
|
---|
911 | //-------------------------------------------------------------------
|
---|
912 | void GeneFluct3D::ToRedshiftSpace(int type_beta,double fixed_beta)
|
---|
913 | // Apply redshift distortion corrections
|
---|
914 | // P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
|
---|
915 | // type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
|
---|
916 | // >0 compute factor with fix beta=Omega^fixed_beta
|
---|
917 | // <0 compute factor with fix beta=fixed_beta
|
---|
918 | // ---- ATTENTION:
|
---|
919 | // 1- Le spectre Pk doit etre (dRho/Rho)(k)
|
---|
920 | // 2- pas d'evolution en redhsuft possible,
|
---|
921 | // le spectre doit etre calcule au redshift final de la boite
|
---|
922 | {
|
---|
923 | double zpk = compute_pk_redsh_ref_;
|
---|
924 | if(lp_>0) cout<<"--- ToRedshiftSpace --- at z="<<zpk<<", type_beta="<<type_beta<<endl;
|
---|
925 |
|
---|
926 | double beta = 0.;
|
---|
927 | if(type_beta>0) {
|
---|
928 | double omegam = cosmo_->Omatter(zpk); beta = pow(omegam,fixed_beta);
|
---|
929 | if(lp_>0) cout<<"...beta = OmegaM^expo = "<<omegam<<"^"<<fixed_beta<<" = "<<beta<<endl;
|
---|
930 | } else if(type_beta<0) {
|
---|
931 | beta = fixed_beta;
|
---|
932 | if(lp_>0) cout<<"...beta = "<<beta<<endl;
|
---|
933 | } else {
|
---|
934 | beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
|
---|
935 | if(lp_>0) cout<<"...beta = -(1+z)*D'/D = "<<beta<<endl;
|
---|
936 | }
|
---|
937 |
|
---|
938 | check_array_alloc();
|
---|
939 |
|
---|
940 | for(long i=0;i<Nx_;i++) {
|
---|
941 | double kx = Kx(i);
|
---|
942 | for(long j=0;j<Ny_;j++) {
|
---|
943 | double ky = Ky(j);
|
---|
944 | double kt2 = kx*kx + ky*ky;
|
---|
945 | for(long l=0;l<NCz_;l++) {
|
---|
946 | double kz = Kz(l);
|
---|
947 | if(l==0) continue;
|
---|
948 | T_(l,j,i) *= (1. + beta*kz*kz/(kt2 + kz*kz));
|
---|
949 | }
|
---|
950 | }
|
---|
951 | }
|
---|
952 | }
|
---|
953 |
|
---|
954 | //-------------------------------------------------------------------
|
---|
955 | void GeneFluct3D::ToVelLoS(int type_beta,double fixed_beta)
|
---|
956 | /*
|
---|
957 | // P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
|
---|
958 | // type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
|
---|
959 | // >0 compute factor with fix beta=Omega^fixed_beta
|
---|
960 | // <0 compute factor with fix beta=fixed_beta
|
---|
961 | ---- Vitesse particuliere
|
---|
962 | Un atome au redshift "z" emet a la longueur d'onde "Le" si il est au repos
|
---|
963 | Il a une vitesse (particuliere) "V" dans le referentiel comobile au redshift "z"
|
---|
964 | Dans le referentiel comobile au redshift "z", il emet a la longueur d'onde "Le*(1+V/C)"
|
---|
965 | L'observateur voit donc une longueur d'onde "Le*(1+V/C)*(1+z) = Le*[1+z+V/C*(1+z)]"
|
---|
966 | et croit que l'atome est au redshift "z' = z + V/C*(1+z)"
|
---|
967 | L'observateur observe donc une vitesse particuliere (1+z)*V/c
|
---|
968 | --> vitesse particuliere dans le repere comobile en z: V/C
|
---|
969 | vitesse particuliere dans le repere comobile de l'observateur: (1+z)*V/C
|
---|
970 | ---- ATTENTION: Le spectre Pk doit etre (dRho/Rho)(k)
|
---|
971 | ---- On genere la vitesse particuliere (dans le referentiel comobile au redshift "z")
|
---|
972 | Vlos(k) = i*Beta(z)*k(los)/k^2*(dRho/Rho)(k) (unite Mpc)
|
---|
973 | .avec Beta(z) = dln(D)/da = -(1+z)/D*dD/dz
|
---|
974 | .puis on convertit en vitesse en faisant Vlos(k) * H(z) (unite lm/s)
|
---|
975 | */
|
---|
976 | {
|
---|
977 | double zpk = compute_pk_redsh_ref_;
|
---|
978 | if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", type_beta="<<type_beta<<endl;
|
---|
979 |
|
---|
980 | double beta = 0;
|
---|
981 | if(type_beta>0) {
|
---|
982 | double omegam = cosmo_->Omatter(zpk);
|
---|
983 | beta = pow(omegam,fixed_beta);
|
---|
984 | if(lp_>0) cout<<"...beta = OmegaM^expo = "<<omegam<<"^"<<fixed_beta<<"="<<beta<<endl;
|
---|
985 | } else if(type_beta<0) {
|
---|
986 | beta = fixed_beta;
|
---|
987 | if(lp_>0) cout<<"...beta="<<beta<<endl;
|
---|
988 | } else {
|
---|
989 | beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
|
---|
990 | if(lp_>0) cout<<"...beta = -(1+z)*D'/D = "<<beta<<endl;
|
---|
991 | }
|
---|
992 | double dpsd = cosmo_->H(zpk) * beta;
|
---|
993 | if(lp_>0) cout<<"......H*beta = "<<dpsd<<" (km/s)/Mpc"<<endl;
|
---|
994 |
|
---|
995 | check_array_alloc();
|
---|
996 |
|
---|
997 | for(long i=0;i<Nx_;i++) {
|
---|
998 | double kx = Kx(i);
|
---|
999 | for(long j=0;j<Ny_;j++) {
|
---|
1000 | double ky = Ky(j);
|
---|
1001 | double kt2 = kx*kx + ky*ky;
|
---|
1002 | for(long l=0;l<NCz_;l++) {
|
---|
1003 | double kz = Kz(l);
|
---|
1004 | double k2 = kt2 + kz*kz;
|
---|
1005 | if(l==0) k2 = 0.; else k2 = dpsd*kz/k2;
|
---|
1006 | T_(l,j,i) *= complex<double>(0.,k2);
|
---|
1007 | }
|
---|
1008 | }
|
---|
1009 | }
|
---|
1010 |
|
---|
1011 | }
|
---|
1012 |
|
---|
1013 | //-------------------------------------------------------------------
|
---|
1014 | void GeneFluct3D::ApplyGrowthFactor(int type_evol)
|
---|
1015 | // Apply Growth to real space d(Rho)/Rho cube
|
---|
1016 | // Using the correspondance between redshift and los comoving distance
|
---|
1017 | // describe in vector "zred_" "loscom_"
|
---|
1018 | // type_evol = 1 : evolution avec la distance a l'observateur
|
---|
1019 | // 2 : evolution avec la distance du plan Z
|
---|
1020 | // (tous les pixels d'un plan Z sont mis au meme redshift z que celui du milieu)
|
---|
1021 | {
|
---|
1022 | if(lp_>0) cout<<"--- ApplyGrowthFactor: evol="<<type_evol<<" (zpk="<<compute_pk_redsh_ref_<<")"<<endl;
|
---|
1023 | check_array_alloc();
|
---|
1024 |
|
---|
1025 | if(growth_ == NULL) {
|
---|
1026 | const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
|
---|
1027 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1028 | }
|
---|
1029 | if(type_evol<1 || type_evol>2) {
|
---|
1030 | const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
|
---|
1031 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1032 | }
|
---|
1033 |
|
---|
1034 | InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
|
---|
1035 |
|
---|
1036 | for(long i=0;i<Nx_;i++) {
|
---|
1037 | double dx2 = DXcom(i); dx2 *= dx2;
|
---|
1038 | for(long j=0;j<Ny_;j++) {
|
---|
1039 | double dy2 = DYcom(j); dy2 *= dy2;
|
---|
1040 | for(long l=0;l<Nz_;l++) {
|
---|
1041 | double dz = DZcom(l);
|
---|
1042 | if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
|
---|
1043 | else dz = fabs(dz); // tous les plans Z au meme redshift
|
---|
1044 | double z = interpinv(dz); // interpolation par morceau
|
---|
1045 | double dzgr = (*growth_)(z);
|
---|
1046 | int_8 ip = IndexR(i,j,l);
|
---|
1047 | data_[ip] *= dzgr;
|
---|
1048 | }
|
---|
1049 | }
|
---|
1050 | }
|
---|
1051 |
|
---|
1052 | }
|
---|
1053 |
|
---|
1054 | //-------------------------------------------------------------------
|
---|
1055 | void GeneFluct3D::ApplyRedshiftSpaceGrowthFactor(void)
|
---|
1056 | // Apply Growth(z) to redshift-distorted real space cube
|
---|
1057 | {
|
---|
1058 | const char *bla = "GeneFluct3D::ApplyRedshiftSpaceGrowthFactor_Error: redshift evolution impossible because factor depends on k";
|
---|
1059 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1060 | }
|
---|
1061 |
|
---|
1062 | //-------------------------------------------------------------------
|
---|
1063 | void GeneFluct3D::ApplyVelLosGrowthFactor(int type_evol,int type_beta,double fixed_beta)
|
---|
1064 | // Apply Growth(z) to transverse velocity real space cube
|
---|
1065 | // P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
|
---|
1066 | // type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
|
---|
1067 | // >0 compute factor with fix beta=Omega^fixed_beta
|
---|
1068 | // <0 compute factor with fix beta=fixed_beta
|
---|
1069 | {
|
---|
1070 | if(lp_>0) cout<<"--- ApplyVelLosGrowthFactor: evol="<<type_evol<<", type_beta="<<type_beta<<endl;
|
---|
1071 |
|
---|
1072 | check_array_alloc();
|
---|
1073 |
|
---|
1074 | if(growth_ == NULL) {
|
---|
1075 | const char *bla = "GeneFluct3D::ApplyVelLosGrowthFactor_Error: set GrowthFactor first";
|
---|
1076 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1077 | }
|
---|
1078 | if(type_evol<1 || type_evol>2) {
|
---|
1079 | const char *bla = "GeneFluct3D::ApplyVelLosGrowthFactor_Error: bad type_evol value";
|
---|
1080 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1081 | }
|
---|
1082 |
|
---|
1083 | double zpk = compute_pk_redsh_ref_;
|
---|
1084 | double grw_orig = (*growth_)(zpk);
|
---|
1085 | double beta_orig = 0.;
|
---|
1086 | if(type_beta>0) {double omegam = cosmo_->Omatter(zpk); beta_orig = pow(omegam,fixed_beta);}
|
---|
1087 | else if(type_beta<0) beta_orig = fixed_beta;
|
---|
1088 | else beta_orig = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
|
---|
1089 | double dpsd_orig = cosmo_->H(zpk) * beta_orig;
|
---|
1090 | if(lp_>0) cout<<"...original at z="<<zpk<<": growth="<<grw_orig<<" beta="<<beta_orig
|
---|
1091 | <<" H*beta="<<dpsd_orig<<" (km/s)/Mpc"<<endl;
|
---|
1092 |
|
---|
1093 | if(type_beta<0) cout<<"***WARNING: type_beta<0 , no evolution used for beta !!!"<<endl;
|
---|
1094 |
|
---|
1095 | InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
|
---|
1096 | InterpFunc interdpd(zredmin_dpsd_,zredmax_dpsd_,dpsdfrzred_);
|
---|
1097 |
|
---|
1098 | for(long i=0;i<Nx_;i++) {
|
---|
1099 | double dx2 = DXcom(i); dx2 *= dx2;
|
---|
1100 | for(long j=0;j<Ny_;j++) {
|
---|
1101 | double dy2 = DYcom(j); dy2 *= dy2;
|
---|
1102 | for(long l=0;l<Nz_;l++) {
|
---|
1103 | double dz = DZcom(l);
|
---|
1104 | if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
|
---|
1105 | else dz = fabs(dz); // tous les plans Z au meme redshift
|
---|
1106 | double z = interpinv(dz); // interpolation par morceau
|
---|
1107 | double grw = (*growth_)(z);
|
---|
1108 | double dpsd = 0;
|
---|
1109 | if(type_beta>0) {
|
---|
1110 | dpsd = cosmo_->H(z) * pow(cosmo_->Omatter(z),fixed_beta);
|
---|
1111 | } else if(type_beta<0) {
|
---|
1112 | dpsd = cosmo_->H(z) * fixed_beta;
|
---|
1113 | } else {
|
---|
1114 | dpsd = interdpd(z); // interpole -H(z)*(1+z)*D'(z)/D(z)
|
---|
1115 | }
|
---|
1116 | int_8 ip = IndexR(i,j,l);
|
---|
1117 | // on remet le beta au bon z
|
---|
1118 | // on corrige du growth factor car data_ a ete calcule avec pk(zpk)
|
---|
1119 | data_[ip] *= (dpsd / dpsd_orig) * (grw / grw_orig);
|
---|
1120 | }
|
---|
1121 | }
|
---|
1122 | }
|
---|
1123 |
|
---|
1124 | }
|
---|
1125 |
|
---|
1126 | //-------------------------------------------------------------------
|
---|
1127 | void GeneFluct3D::ComputeReal(void)
|
---|
1128 | // Calcule une realisation dans l'espace reel
|
---|
1129 | {
|
---|
1130 | if(lp_>0) cout<<"--- ComputeReal --- from spectrum at z="<<compute_pk_redsh_ref_<<endl;
|
---|
1131 | check_array_alloc();
|
---|
1132 |
|
---|
1133 | // On fait la FFT
|
---|
1134 | GEN3D_FFTW_EXECUTE(pb_);
|
---|
1135 | array_type = 1;
|
---|
1136 | }
|
---|
1137 |
|
---|
1138 | //-------------------------------------------------------------------
|
---|
1139 | void GeneFluct3D::ReComputeFourier(void)
|
---|
1140 | {
|
---|
1141 | if(lp_>0) cout<<"--- ReComputeFourier ---"<<endl;
|
---|
1142 | check_array_alloc();
|
---|
1143 |
|
---|
1144 | // On fait la FFT
|
---|
1145 | GEN3D_FFTW_EXECUTE(pf_);
|
---|
1146 | array_type = 2;
|
---|
1147 |
|
---|
1148 | // On corrige du pb de la normalisation de FFTW3
|
---|
1149 | complex<r_8> v((r_8)NRtot_,0.);
|
---|
1150 | for(long i=0;i<Nx_;i++)
|
---|
1151 | for(long j=0;j<Ny_;j++)
|
---|
1152 | for(long l=0;l<NCz_;l++) T_(l,j,i) /= v;
|
---|
1153 | }
|
---|
1154 |
|
---|
1155 | //-------------------------------------------------------------------
|
---|
1156 | int GeneFluct3D::ComputeSpectrum(HistoErr& herr)
|
---|
1157 | // Compute spectrum from "T" and fill HistoErr "herr"
|
---|
1158 | // T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
|
---|
1159 | // cad T(kz,ky,kx) avec 0<kz<kz_nyq -ky_nyq<ky<ky_nyq -kx_nyq<kx<kx_nyq
|
---|
1160 | {
|
---|
1161 | if(lp_>0) cout<<"--- ComputeSpectrum ---"<<endl;
|
---|
1162 | check_array_alloc();
|
---|
1163 |
|
---|
1164 | if(herr.NBins()<0) return -1;
|
---|
1165 | herr.Zero();
|
---|
1166 |
|
---|
1167 | // Attention a l'ordre
|
---|
1168 | for(long i=0;i<Nx_;i++) {
|
---|
1169 | double kx = AbsKx(i); kx *= kx;
|
---|
1170 | for(long j=0;j<Ny_;j++) {
|
---|
1171 | double ky = AbsKy(j); ky *= ky;
|
---|
1172 | for(long l=0;l<NCz_;l++) {
|
---|
1173 | double kz = AbsKz(l);
|
---|
1174 | double k = sqrt(kx+ky+kz*kz);
|
---|
1175 | double pk = MODULE2(T_(l,j,i));
|
---|
1176 | herr.Add(k,pk);
|
---|
1177 | }
|
---|
1178 | }
|
---|
1179 | }
|
---|
1180 | herr.ToVariance();
|
---|
1181 |
|
---|
1182 | // renormalize to directly compare to original spectrum
|
---|
1183 | double norm = Vol_;
|
---|
1184 | herr *= norm;
|
---|
1185 |
|
---|
1186 | return 0;
|
---|
1187 | }
|
---|
1188 |
|
---|
1189 | int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
|
---|
1190 | {
|
---|
1191 | if(lp_>0) cout<<"--- ComputeSpectrum2D ---"<<endl;
|
---|
1192 | check_array_alloc();
|
---|
1193 |
|
---|
1194 | if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
|
---|
1195 | herr.Zero();
|
---|
1196 |
|
---|
1197 | // Attention a l'ordre
|
---|
1198 | for(long i=0;i<Nx_;i++) {
|
---|
1199 | double kx = AbsKx(i); kx *= kx;
|
---|
1200 | for(long j=0;j<Ny_;j++) {
|
---|
1201 | double ky = AbsKy(j); ky *= ky;
|
---|
1202 | double kt = sqrt(kx+ky);
|
---|
1203 | for(long l=0;l<NCz_;l++) {
|
---|
1204 | double kz = AbsKz(l);
|
---|
1205 | double pk = MODULE2(T_(l,j,i));
|
---|
1206 | herr.Add(kt,kz,pk);
|
---|
1207 | }
|
---|
1208 | }
|
---|
1209 | }
|
---|
1210 | herr.ToVariance();
|
---|
1211 |
|
---|
1212 | // renormalize to directly compare to original spectrum
|
---|
1213 | double norm = Vol_;
|
---|
1214 | herr *= norm;
|
---|
1215 |
|
---|
1216 | return 0;
|
---|
1217 | }
|
---|
1218 |
|
---|
1219 | //-------------------------------------------------------------------
|
---|
1220 | int GeneFluct3D::ComputeSpectrum(HistoErr& herr,double sigma,bool pixcor)
|
---|
1221 | // Compute spectrum from "T" and fill HistoErr "herr"
|
---|
1222 | // AVEC la soustraction du niveau de bruit et la correction par filterpixel()
|
---|
1223 | // Si on ne fait pas ca, alors on obtient un spectre non-isotrope!
|
---|
1224 | //
|
---|
1225 | // T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
|
---|
1226 | // cad T(kz,ky,kx) avec 0<kz<kz_nyq -ky_nyq<ky<ky_nyq -kx_nyq<kx<kx_nyq
|
---|
1227 | {
|
---|
1228 | if(lp_>0) cout<<"--- ComputeSpectrum: sigma="<<sigma<<endl;
|
---|
1229 | check_array_alloc();
|
---|
1230 |
|
---|
1231 | if(sigma<=0.) sigma = 0.;
|
---|
1232 | double sigma2 = sigma*sigma / (double)NRtot_;
|
---|
1233 |
|
---|
1234 | if(herr.NBins()<0) return -1;
|
---|
1235 | herr.Zero();
|
---|
1236 |
|
---|
1237 | TVector<r_8> vfz(NCz_);
|
---|
1238 | if(pixcor) // kz = l*Dkz_
|
---|
1239 | for(long l=0;l<NCz_;l++) {vfz(l)=pixelfilter(l*Dkz_ *Dz_/2); vfz(l)*=vfz(l);}
|
---|
1240 |
|
---|
1241 | // Attention a l'ordre
|
---|
1242 | for(long i=0;i<Nx_;i++) {
|
---|
1243 | double kx = AbsKx(i);
|
---|
1244 | double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
|
---|
1245 | kx *= kx; fx *= fx;
|
---|
1246 | for(long j=0;j<Ny_;j++) {
|
---|
1247 | double ky = AbsKy(j);
|
---|
1248 | double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
|
---|
1249 | ky *= ky; fy *= fy;
|
---|
1250 | for(long l=0;l<NCz_;l++) {
|
---|
1251 | double kz = AbsKz(l);
|
---|
1252 | double k = sqrt(kx+ky+kz*kz);
|
---|
1253 | double pk = MODULE2(T_(l,j,i)) - sigma2;
|
---|
1254 | double fz = (pixcor) ? vfz(l): 1.;
|
---|
1255 | double f = fx*fy*fz;
|
---|
1256 | if(f>0.) herr.Add(k,pk/f);
|
---|
1257 | }
|
---|
1258 | }
|
---|
1259 | }
|
---|
1260 | herr.ToVariance();
|
---|
1261 | for(int i=0;i<herr.NBins();i++) herr(i) += sigma2;
|
---|
1262 |
|
---|
1263 | // renormalize to directly compare to original spectrum
|
---|
1264 | double norm = Vol_;
|
---|
1265 | herr *= norm;
|
---|
1266 |
|
---|
1267 | return 0;
|
---|
1268 | }
|
---|
1269 |
|
---|
1270 | int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr,double sigma,bool pixcor)
|
---|
1271 | // AVEC la soustraction du niveau de bruit et la correction par filterpixel()
|
---|
1272 | {
|
---|
1273 | if(lp_>0) cout<<"--- ComputeSpectrum2D: sigma="<<sigma<<endl;
|
---|
1274 | check_array_alloc();
|
---|
1275 |
|
---|
1276 | if(sigma<=0.) sigma = 0.;
|
---|
1277 | double sigma2 = sigma*sigma / (double)NRtot_;
|
---|
1278 |
|
---|
1279 | if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
|
---|
1280 | herr.Zero();
|
---|
1281 |
|
---|
1282 | TVector<r_8> vfz(NCz_);
|
---|
1283 | if(pixcor) // kz = l*Dkz_
|
---|
1284 | for(long l=0;l<NCz_;l++) {vfz(l)=pixelfilter(l*Dkz_ *Dz_/2); vfz(l)*=vfz(l);}
|
---|
1285 |
|
---|
1286 | // Attention a l'ordre
|
---|
1287 | for(long i=0;i<Nx_;i++) {
|
---|
1288 | double kx = AbsKx(i);
|
---|
1289 | double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
|
---|
1290 | kx *= kx; fx *= fx;
|
---|
1291 | for(long j=0;j<Ny_;j++) {
|
---|
1292 | double ky = AbsKy(j);
|
---|
1293 | double fy = (pixcor) ? pixelfilter(ky*Dy_/2) : 1.;
|
---|
1294 | ky *= ky; fy *= fy;
|
---|
1295 | double kt = sqrt(kx+ky);
|
---|
1296 | for(long l=0;l<NCz_;l++) {
|
---|
1297 | double kz = AbsKz(l);
|
---|
1298 | double pk = MODULE2(T_(l,j,i)) - sigma2;
|
---|
1299 | double fz = (pixcor) ? vfz(l): 1.;
|
---|
1300 | double f = fx*fy*fz;
|
---|
1301 | if(f>0.) herr.Add(kt,kz,pk/f);
|
---|
1302 | }
|
---|
1303 | }
|
---|
1304 | }
|
---|
1305 | herr.ToVariance();
|
---|
1306 | for(int i=0;i<herr.NBinX();i++)
|
---|
1307 | for(int j=0;j<herr.NBinY();j++) herr(i,j) += sigma2;
|
---|
1308 |
|
---|
1309 | // renormalize to directly compare to original spectrum
|
---|
1310 | double norm = Vol_;
|
---|
1311 | herr *= norm;
|
---|
1312 |
|
---|
1313 | return 0;
|
---|
1314 | }
|
---|
1315 |
|
---|
1316 | //-------------------------------------------------------
|
---|
1317 | int_8 GeneFluct3D::VarianceFrReal(double R,double& var)
|
---|
1318 | // Recompute MASS variance in spherical top-hat (rayon=R)
|
---|
1319 | // Par definition: SigmaR^2 = <(M-<M>)^2>/<M>^2
|
---|
1320 | // ou M = masse dans sphere de rayon R
|
---|
1321 | // --- ATTENTION: la variance calculee a une tres grande dispersion
|
---|
1322 | // (surtout si le volume du cube est petit). Pour verifier
|
---|
1323 | // que le sigmaR calcule par cette methode est en accord avec
|
---|
1324 | // le sigmaR en input, il faut faire plusieurs simulations (~100)
|
---|
1325 | // et regarder la moyenne des sigmaR reconstruits
|
---|
1326 | {
|
---|
1327 | if(lp_>0) cout<<"--- VarianceFrReal R="<<R<<endl;
|
---|
1328 | check_array_alloc();
|
---|
1329 |
|
---|
1330 | long dnx = long(R/Dx_)+1; if(dnx<=0) dnx = 1;
|
---|
1331 | long dny = long(R/Dy_)+1; if(dny<=0) dny = 1;
|
---|
1332 | long dnz = long(R/Dz_)+1; if(dnz<=0) dnz = 1;
|
---|
1333 | if(lp_>0) cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
|
---|
1334 |
|
---|
1335 | double sum=0., sum2=0., sn=0., r2 = R*R;
|
---|
1336 | int_8 nsum=0;
|
---|
1337 |
|
---|
1338 | for(long i=dnx;i<Nx_-dnx;i+=2*dnx) {
|
---|
1339 | for(long j=dny;j<Ny_-dny;j+=2*dny) {
|
---|
1340 | for(long l=dnz;l<Nz_-dnz;l+=2*dnz) {
|
---|
1341 | double m=0.; int_8 n=0;
|
---|
1342 | for(long ii=i-dnx;ii<=i+dnx;ii++) {
|
---|
1343 | double x = (ii-i)*Dx_; x *= x;
|
---|
1344 | for(long jj=j-dny;jj<=j+dny;jj++) {
|
---|
1345 | double y = (jj-j)*Dy_; y *= y;
|
---|
1346 | for(long ll=l-dnz;ll<=l+dnz;ll++) {
|
---|
1347 | double z = (ll-l)*Dz_; z *= z;
|
---|
1348 | if(x+y+z>r2) continue;
|
---|
1349 | int_8 ip = IndexR(ii,jj,ll);
|
---|
1350 | m += 1.+data_[ip]; // 1+drho/rho
|
---|
1351 | n++;
|
---|
1352 | }
|
---|
1353 | }
|
---|
1354 | }
|
---|
1355 | if(n>0) {sum += m; sum2 += m*m; nsum++; sn += n;}
|
---|
1356 | //cout<<i<<","<<j<<","<<l<<" n="<<n<<" m="<<m<<" sum="<<sum<<" sum2="<<sum2<<endl;
|
---|
1357 | }
|
---|
1358 | }
|
---|
1359 | }
|
---|
1360 |
|
---|
1361 | if(nsum<=1) {var=0.; return nsum;}
|
---|
1362 | sum /= nsum;
|
---|
1363 | sum2 = sum2/nsum - sum*sum;
|
---|
1364 | sn /= nsum;
|
---|
1365 | if(lp_>0) cout<<"...<n>="<<sn<<", nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
|
---|
1366 | var = sum2/(sum*sum); // <dM^2>/<M>^2
|
---|
1367 | if(lp_>0) cout<<"...sigmaR^2 = <(M-<M>)^2>/<M>^2 = "<<var
|
---|
1368 | <<" -> sigmaR = "<<sqrt(var)<<endl;
|
---|
1369 |
|
---|
1370 | return nsum;
|
---|
1371 | }
|
---|
1372 |
|
---|
1373 | //-------------------------------------------------------
|
---|
1374 | int_8 GeneFluct3D::NumberOfBad(double vmin,double vmax)
|
---|
1375 | // number of pixels outside of ]vmin,vmax[ extremites exclues
|
---|
1376 | // -> vmin and vmax are considered as bad
|
---|
1377 | {
|
---|
1378 | check_array_alloc();
|
---|
1379 |
|
---|
1380 | int_8 nbad = 0;
|
---|
1381 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1382 | int_8 ip = IndexR(i,j,l);
|
---|
1383 | double v = data_[ip];
|
---|
1384 | if(v<=vmin || v>=vmax) nbad++;
|
---|
1385 | }
|
---|
1386 |
|
---|
1387 | if(lp_>0) cout<<"--- NumberOfBad "<<nbad<<" px out of ]"<<vmin<<","<<vmax
|
---|
1388 | <<"[ i.e. frac="<<nbad/(double)NRtot_<<endl;
|
---|
1389 | return nbad;
|
---|
1390 | }
|
---|
1391 |
|
---|
1392 | int_8 GeneFluct3D::MinMax(double& xmin,double& xmax,double vmin,double vmax)
|
---|
1393 | // Calcul des valeurs xmin et xmax dans le cube reel avec valeurs ]vmin,vmax[ extremites exclues
|
---|
1394 | {
|
---|
1395 | bool tstval = (vmax>vmin)? true: false;
|
---|
1396 | if(lp_>0) {
|
---|
1397 | cout<<"--- MinMax";
|
---|
1398 | if(tstval) cout<<" range=]"<<vmin<<","<<vmax<<"[";
|
---|
1399 | cout<<endl;
|
---|
1400 | }
|
---|
1401 | check_array_alloc();
|
---|
1402 |
|
---|
1403 | int_8 n = 0;
|
---|
1404 | xmin = xmax = data_[0];
|
---|
1405 |
|
---|
1406 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1407 | int_8 ip = IndexR(i,j,l);
|
---|
1408 | double x = data_[ip];
|
---|
1409 | if(tstval && (x<=vmin || x>=vmax)) continue;
|
---|
1410 | if(x<xmin) xmin = x;
|
---|
1411 | if(x>xmax) xmax = x;
|
---|
1412 | n++;
|
---|
1413 | }
|
---|
1414 |
|
---|
1415 | if(lp_>0) cout<<" n="<<n<<" min="<<xmin<<" max="<<xmax<<endl;
|
---|
1416 |
|
---|
1417 | return n;
|
---|
1418 | }
|
---|
1419 |
|
---|
1420 | int_8 GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax
|
---|
1421 | ,bool useout,double vout)
|
---|
1422 | // Calcul de mean,sigma2 dans le cube reel avec valeurs ]vmin,vmax[ extremites exclues
|
---|
1423 | // useout = false: ne pas utiliser les pixels hors limites pour calculer mean,sigma2
|
---|
1424 | // true : utiliser les pixels hors limites pour calculer mean,sigma2
|
---|
1425 | // en remplacant leurs valeurs par "vout"
|
---|
1426 | {
|
---|
1427 | bool tstval = (vmax>vmin)? true: false;
|
---|
1428 | if(lp_>0) {
|
---|
1429 | cout<<"--- MeanSigma2";
|
---|
1430 | if(tstval) cout<<" range=]"<<vmin<<","<<vmax<<"[";
|
---|
1431 | if(useout) cout<<", useout="<<useout<<" vout="<<vout;
|
---|
1432 | cout<<endl;
|
---|
1433 | }
|
---|
1434 | check_array_alloc();
|
---|
1435 |
|
---|
1436 | int_8 n = 0;
|
---|
1437 | rm = rs2 = 0.;
|
---|
1438 |
|
---|
1439 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1440 | int_8 ip = IndexR(i,j,l);
|
---|
1441 | double v = data_[ip];
|
---|
1442 | if(tstval) {
|
---|
1443 | if(v<=vmin || v>=vmax) {if(useout) v=vout; else continue;}
|
---|
1444 | }
|
---|
1445 | rm += v;
|
---|
1446 | rs2 += v*v;
|
---|
1447 | n++;
|
---|
1448 | }
|
---|
1449 |
|
---|
1450 | if(n>1) {
|
---|
1451 | rm /= (double)n;
|
---|
1452 | rs2 = rs2/(double)n - rm*rm;
|
---|
1453 | }
|
---|
1454 |
|
---|
1455 | if(lp_>0) cout<<" n="<<n<<" m="<<rm<<" s2="<<rs2<<" s="<<sqrt(fabs(rs2))<<endl;
|
---|
1456 |
|
---|
1457 | return n;
|
---|
1458 | }
|
---|
1459 |
|
---|
1460 | int_8 GeneFluct3D::SetToVal(double vmin, double vmax,double val0)
|
---|
1461 | // set to "val0" if out of range ]vmin,vmax[ extremites exclues
|
---|
1462 | // cad set to "val0" if in [vmin,vmax] -> vmin and vmax are set to val0
|
---|
1463 | {
|
---|
1464 | check_array_alloc();
|
---|
1465 |
|
---|
1466 | int_8 nbad = 0;
|
---|
1467 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1468 | int_8 ip = IndexR(i,j,l);
|
---|
1469 | double v = data_[ip];
|
---|
1470 | if(v<=vmin || v>=vmax) {data_[ip] = val0; nbad++;}
|
---|
1471 | }
|
---|
1472 |
|
---|
1473 | if(lp_>0) cout<<"--- SetToVal "<<nbad<<" px set to="<<val0
|
---|
1474 | <<" because out of range=]"<<vmin<<","<<vmax<<"["<<endl;
|
---|
1475 | return nbad;
|
---|
1476 | }
|
---|
1477 |
|
---|
1478 | void GeneFluct3D::ScaleOffset(double scalecube,double offsetcube)
|
---|
1479 | // Replace "V" by "scalecube * ( V + offsetcube )"
|
---|
1480 | {
|
---|
1481 | if(lp_>0) cout<<"--- ScaleCube scale="<<scalecube<<" offset="<<offsetcube<<endl;
|
---|
1482 |
|
---|
1483 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1484 | int_8 ip = IndexR(i,j,l);
|
---|
1485 | data_[ip] = scalecube * ( data_[ip] + offsetcube );
|
---|
1486 | }
|
---|
1487 |
|
---|
1488 | return;
|
---|
1489 | }
|
---|
1490 |
|
---|
1491 | //-------------------------------------------------------
|
---|
1492 | void GeneFluct3D::TurnFluct2Mass(void)
|
---|
1493 | // d_rho/rho -> Mass (add one!)
|
---|
1494 | {
|
---|
1495 | if(lp_>0) cout<<"--- TurnFluct2Mass ---"<<endl;
|
---|
1496 | check_array_alloc();
|
---|
1497 |
|
---|
1498 |
|
---|
1499 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1500 | int_8 ip = IndexR(i,j,l);
|
---|
1501 | data_[ip] += 1.;
|
---|
1502 | }
|
---|
1503 | }
|
---|
1504 |
|
---|
1505 | double GeneFluct3D::TurnFluct2MeanNumber(double val_by_mpc3)
|
---|
1506 | // ATTENTION: la gestion des pixels<0 proposee ici induit une perte de variance
|
---|
1507 | // sur la carte, le spectre Pk reconstruit sera plus faible!
|
---|
1508 | // L'effet sera d'autant plus grand que le nombre de pixels<0 sera grand.
|
---|
1509 | {
|
---|
1510 | if(lp_>0) cout<<"--- TurnFluct2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;
|
---|
1511 |
|
---|
1512 | // First convert dRho/Rho into 1+dRho/Rho
|
---|
1513 | int_8 nball = 0; double sumall = 0., sumall2 = 0.;
|
---|
1514 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1515 | int_8 ip = IndexR(i,j,l);
|
---|
1516 | data_[ip] += 1.;
|
---|
1517 | nball++; sumall += data_[ip]; sumall2 += data_[ip]*data_[ip];
|
---|
1518 | }
|
---|
1519 | if(nball>2) {
|
---|
1520 | sumall /= (double)nball;
|
---|
1521 | sumall2 = sumall2/(double)nball - sumall*sumall;
|
---|
1522 | if(lp_>0) cout<<"1+dRho/Rho: mean="<<sumall<<" variance="<<sumall2
|
---|
1523 | <<" -> "<<sqrt(fabs(sumall2))<<endl;
|
---|
1524 | }
|
---|
1525 |
|
---|
1526 | // Find contribution for positive pixels
|
---|
1527 | int_8 nbpos = 0; double sumpos = 0. , sumpos2 = 0.;
|
---|
1528 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1529 | int_8 ip = IndexR(i,j,l);
|
---|
1530 | double v = data_[ip];
|
---|
1531 | if(data_[ip]>0.) {nbpos++; sumpos += v; sumpos2 += v*v;}
|
---|
1532 | }
|
---|
1533 | if(nbpos<1) {
|
---|
1534 | cout<<"TurnFluct2MeanNumber_Error: nbpos<1"<<endl;
|
---|
1535 | throw RangeCheckError("TurnFluct2MeanNumber_Error: nbpos<1");
|
---|
1536 | }
|
---|
1537 | sumpos2 = sumpos2/nball - sumpos*sumpos/(nball*nball);
|
---|
1538 | if(lp_>0)
|
---|
1539 | cout<<"1+dRho/Rho with v<0 set to zero: mean="<<sumpos/nball
|
---|
1540 | <<" variance="<<sumpos2<<" -> "<<sqrt(fabs(sumpos2))<<endl;
|
---|
1541 | cout<<"Sum of positive values: sumpos="<<sumpos
|
---|
1542 | <<" (n(v>0) = "<<nbpos<<" frac(v>0)="<<nbpos/(double)NRtot_<<")"<<endl;
|
---|
1543 |
|
---|
1544 | // - Mettre exactement val_by_mpc3*Vol galaxies (ou Msol) dans notre survey
|
---|
1545 | // - Uniquement dans les pixels de masse >0.
|
---|
1546 | // - Mise a zero des pixels <0
|
---|
1547 | double dn = val_by_mpc3 * Vol_ / sumpos;
|
---|
1548 | if(lp_>0) cout<<"...density move from "
|
---|
1549 | <<val_by_mpc3*dVol_<<" to "<<dn<<" / pixel"<<endl;
|
---|
1550 |
|
---|
1551 | double sum = 0.;
|
---|
1552 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1553 | int_8 ip = IndexR(i,j,l);
|
---|
1554 | if(data_[ip]<=0.) data_[ip] = 0.;
|
---|
1555 | else {
|
---|
1556 | data_[ip] *= dn;
|
---|
1557 | sum += data_[ip];
|
---|
1558 | }
|
---|
1559 | }
|
---|
1560 |
|
---|
1561 | if(lp_>0) cout<<"...quantity put into survey "<<sum<<" / "<<val_by_mpc3*Vol_<<endl;
|
---|
1562 |
|
---|
1563 | return sum;
|
---|
1564 | }
|
---|
1565 |
|
---|
1566 | double GeneFluct3D::ApplyPoisson(void)
|
---|
1567 | // do NOT treate negative or nul mass -> let it as it is
|
---|
1568 | {
|
---|
1569 | if(lp_>0) cout<<"--- ApplyPoisson ---"<<endl;
|
---|
1570 | check_array_alloc();
|
---|
1571 |
|
---|
1572 | double sum = 0.;
|
---|
1573 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1574 | int_8 ip = IndexR(i,j,l);
|
---|
1575 | double v = data_[ip];
|
---|
1576 | if(v>0.) {
|
---|
1577 | double dn = PoissonRand(v,10.);
|
---|
1578 | data_[ip] = dn;
|
---|
1579 | sum += dn;
|
---|
1580 | }
|
---|
1581 | }
|
---|
1582 | if(lp_>0) cout<<sum<<" galaxies put into survey"<<endl;
|
---|
1583 |
|
---|
1584 | return sum;
|
---|
1585 | }
|
---|
1586 |
|
---|
1587 | double GeneFluct3D::TurnNGal2Mass(FunRan& massdist,bool axeslog)
|
---|
1588 | // do NOT treate negative or nul mass -> let it as it is
|
---|
1589 | // INPUT:
|
---|
1590 | // massdist : distribution de masse (m*dn/dm)
|
---|
1591 | // axeslog = false : retourne la masse
|
---|
1592 | // = true : retourne le log10(mass)
|
---|
1593 | // RETURN la masse totale
|
---|
1594 | {
|
---|
1595 | if(lp_>0) cout<<"--- TurnNGal2Mass ---"<<endl;
|
---|
1596 | check_array_alloc();
|
---|
1597 |
|
---|
1598 | double sum = 0.;
|
---|
1599 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1600 | int_8 ip = IndexR(i,j,l);
|
---|
1601 | double v = data_[ip];
|
---|
1602 | if(v>0.) {
|
---|
1603 | long ngal = long(v+0.1);
|
---|
1604 | data_[ip] = 0.;
|
---|
1605 | for(long i=0;i<ngal;i++) {
|
---|
1606 | double m = massdist.RandomInterp(); // massdist.Random();
|
---|
1607 | if(axeslog) m = pow(10.,m);
|
---|
1608 | data_[ip] += m;
|
---|
1609 | }
|
---|
1610 | sum += data_[ip];
|
---|
1611 | }
|
---|
1612 | }
|
---|
1613 | if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;
|
---|
1614 |
|
---|
1615 | return sum;
|
---|
1616 | }
|
---|
1617 |
|
---|
1618 | double GeneFluct3D::TurnNGal2MassQuick(SchechterMassDist& schmdist)
|
---|
1619 | // idem TurnNGal2Mass mais beaucoup plus rapide
|
---|
1620 | {
|
---|
1621 | if(lp_>0) cout<<"--- TurnNGal2MassQuick ---"<<endl;
|
---|
1622 | check_array_alloc();
|
---|
1623 |
|
---|
1624 | double sum = 0.;
|
---|
1625 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
|
---|
1626 | int_8 ip = IndexR(i,j,l);
|
---|
1627 | double v = data_[ip];
|
---|
1628 | if(v>0.) {
|
---|
1629 | long ngal = long(v+0.1);
|
---|
1630 | data_[ip] = schmdist.TirMass(ngal);
|
---|
1631 | sum += data_[ip];
|
---|
1632 | }
|
---|
1633 | }
|
---|
1634 | if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;
|
---|
1635 |
|
---|
1636 | return sum;
|
---|
1637 | }
|
---|
1638 |
|
---|
1639 | void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
|
---|
1640 | // add noise to every pixels (meme les <=0 !)
|
---|
1641 | // type_evol = 0 : pas d'evolution de la puissance du bruit
|
---|
1642 | // 1 : evolution de la puissance du bruit avec la distance a l'observateur
|
---|
1643 | // 2 : evolution de la puissance du bruit avec la distance du plan Z
|
---|
1644 | // (tous les plans Z sont mis au meme redshift z de leur milieu)
|
---|
1645 | {
|
---|
1646 | if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
|
---|
1647 | check_array_alloc();
|
---|
1648 |
|
---|
1649 | if(type_evol<0) type_evol = 0;
|
---|
1650 | if(type_evol>2) {
|
---|
1651 | const char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
|
---|
1652 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1653 | }
|
---|
1654 |
|
---|
1655 | vector<double> correction;
|
---|
1656 | InterpFunc *intercor = NULL;
|
---|
1657 |
|
---|
1658 | if(type_evol>0) {
|
---|
1659 | // Sigma_Noise(en mass) :
|
---|
1660 | // Slim ~ 1/sqrt(DNu) * sqrt(nlobe) en W/m^2Hz
|
---|
1661 | // Flim ~ sqrt(DNu) * sqrt(nlobe) en W/m^2
|
---|
1662 | // Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe) en Msol
|
---|
1663 | // nlobe ~ 1/Dtrcom^2
|
---|
1664 | // Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
|
---|
1665 | if(cosmo_ == NULL || redsh_ref_<0. || loscom2zred_.size()<1) {
|
---|
1666 | const char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
|
---|
1667 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1668 | }
|
---|
1669 | InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
|
---|
1670 | long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
|
---|
1671 | for(long i=0;i<nsz;i++) {
|
---|
1672 | double d = interpinv.X(i);
|
---|
1673 | double zred = interpinv(d);
|
---|
1674 | double dtrc = cosmo_->Dtrcom(zred); // pour variation angle solide
|
---|
1675 | double dlum = cosmo_->Dlum(zred); // pour variation conversion mass HI
|
---|
1676 | double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
|
---|
1677 | double dnu = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
|
---|
1678 | double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
|
---|
1679 | if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
|
---|
1680 | cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
|
---|
1681 | <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
|
---|
1682 | correction.push_back(corr);
|
---|
1683 | }
|
---|
1684 | intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
|
---|
1685 | }
|
---|
1686 |
|
---|
1687 | double corrlim[2] = {1.,1.};
|
---|
1688 | for(long i=0;i<Nx_;i++) {
|
---|
1689 | double dx2 = DXcom(i); dx2 *= dx2;
|
---|
1690 | for(long j=0;j<Ny_;j++) {
|
---|
1691 | double dy2 = DYcom(j); dy2 *= dy2;
|
---|
1692 | for(long l=0;l<Nz_;l++) {
|
---|
1693 | double corr = 1.;
|
---|
1694 | if(type_evol>0) {
|
---|
1695 | double dz = DZcom(l);
|
---|
1696 | if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
|
---|
1697 | else dz = fabs(dz); // tous les plans Z au meme redshift
|
---|
1698 | corr = (*intercor)(dz);
|
---|
1699 | if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
|
---|
1700 | }
|
---|
1701 | int_8 ip = IndexR(i,j,l);
|
---|
1702 | data_[ip] += snoise*corr*NorRand();
|
---|
1703 | }
|
---|
1704 | }
|
---|
1705 | }
|
---|
1706 | if(type_evol>0)
|
---|
1707 | cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;
|
---|
1708 |
|
---|
1709 | if(intercor!=NULL) delete intercor;
|
---|
1710 | }
|
---|
1711 |
|
---|
1712 | } // Fin namespace SOPHYA
|
---|
1713 |
|
---|
1714 |
|
---|
1715 |
|
---|
1716 |
|
---|
1717 | /*********************************************************************
|
---|
1718 | void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
|
---|
1719 | // Add AGN flux into simulation:
|
---|
1720 | // --- Procedure:
|
---|
1721 | // 1. lancer "cmvdefsurv" avec les parametres du survey
|
---|
1722 | // (au redshift de reference du survey)
|
---|
1723 | // et recuperer l'angle solide "angsol sr" du pixel elementaire
|
---|
1724 | // au centre du cube.
|
---|
1725 | // 2. lancer "cmvtstagn" pour cet angle solide -> cmvtstagn.ppf
|
---|
1726 | // 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
|
---|
1727 | // cad une moyenne <log10(S)> et un sigma "sig"
|
---|
1728 | // Attention: la distribution n'est pas gaussienne les "mean,sigma"
|
---|
1729 | // de l'histo ne sont pas vraiment ce que l'on veut
|
---|
1730 | // --- Limitations actuelle du code:
|
---|
1731 | // . les AGN sont supposes evoluer avec la meme loi de puissance pour tout theta,phi
|
---|
1732 | // . le flux des AGN est mis dans une colonne Oz (indice k) et pas sur la ligne de visee
|
---|
1733 | // . la distribution est approximee a une gaussienne
|
---|
1734 | // ... C'est une approximation pour un observateur loin du centre du cube
|
---|
1735 | // et pour un cube peu epais / distance observateur
|
---|
1736 | // --- Parametres de la routine:
|
---|
1737 | // llfy : c'est le <log10(S)> du flux depose par les AGN
|
---|
1738 | // dans l'angle solide du pixel elementaire de reference du cube
|
---|
1739 | // lsigma : c'est le sigma de la distribution des log10(S)
|
---|
1740 | // powlaw : c'est la pente de la distribution cad que le flux "lmsol"
|
---|
1741 | // et considere comme le flux a 1.4GHz et qu'on suppose une loi
|
---|
1742 | // F(nu) = (1.4GHz/nu)^powlaw * F(1.4GHz)
|
---|
1743 | // - Comme on est en echelle log10():
|
---|
1744 | // on tire log10(Msol) + X
|
---|
1745 | // ou X est une realisation sur une gaussienne de variance "sig^2"
|
---|
1746 | // La masse realisee est donc: Msol*10^X
|
---|
1747 | // - Pas de probleme de pixel negatif car on a une multiplication!
|
---|
1748 | {
|
---|
1749 | if(lp_>0) cout<<"--- AddAGN: <log10(S Jy)> = "<<lfjy<<" , sigma = "<<lsigma<<endl;
|
---|
1750 | check_array_alloc();
|
---|
1751 |
|
---|
1752 | if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
|
---|
1753 | char *bla = "GeneFluct3D::AddAGN_Error: set Observator and Cosmology first";
|
---|
1754 | cout<<bla<<endl; throw ParmError(bla);
|
---|
1755 | }
|
---|
1756 |
|
---|
1757 | // Le flux des AGN en Jy et en mass solaire
|
---|
1758 | double fagnref = pow(10.,lfjy)*(dnu_ref_*1.e9); // Jy.Hz = W/m^2
|
---|
1759 | double magnref = FluxHI2Msol(fagnref*Jansky2Watt_cst,dlum_ref_); // Msol
|
---|
1760 | if(lp_>0)
|
---|
1761 | cout<<"Au pixel de ref: fagnref="<<fagnref
|
---|
1762 | <<" Jy.Hz (a 1.4GHz), magnref="<<magnref<<" Msol"<<endl;
|
---|
1763 |
|
---|
1764 | if(powlaw!=0.) {
|
---|
1765 | // F(nu) = F(1.4GHz)*(nu GHz/1.4 Ghz)^p = F(1.4GHz)*(1/(1+z))^p , car nu = 1.4 GHz/(1+z)
|
---|
1766 | magnref *= pow(1/(1.+redsh_ref_),powlaw);
|
---|
1767 | if(lp_>0) cout<<" powlaw="<<powlaw<<" -> change magnref to "<<magnref<<" Msol"<<endl;
|
---|
1768 | }
|
---|
1769 |
|
---|
1770 | // Les infos en fonction de l'indice "l" selon Oz
|
---|
1771 | vector<double> correction;
|
---|
1772 | InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
|
---|
1773 | long nzmod = ((Nz_>10)?Nz_/10:1);
|
---|
1774 | for(long l=0;l<Nz_;l++) {
|
---|
1775 | double z = fabs(DZcom(l));
|
---|
1776 | double zred = interpinv(z);
|
---|
1777 | double dtrc = cosmo_->Dtrcom(zred); // pour variation angle solide
|
---|
1778 | double dlum = cosmo_->Dlum(zred); // pour variation conversion mass HI
|
---|
1779 | double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
|
---|
1780 | double dnu = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
|
---|
1781 | // on a: Mass ~ DNu * Dlum^2 / Dtrcom^2
|
---|
1782 | double corr = dnu/dnu_ref_*pow(dtrc_ref_/dtrc*dlum/dlum_ref_,2.);
|
---|
1783 | // F(nu) = F(1.4GHz)*(nu GHz/1.4 Ghz)^p = F(1.4GHz)*(1/(1+z))^p , car nu = 1.4 GHz/(1+z)
|
---|
1784 | if(powlaw!=0.) corr *= pow((1.+redsh_ref_)/(1.+zred),powlaw);
|
---|
1785 | correction.push_back(corr);
|
---|
1786 | if(lp_>0 && (l==0 || l==Nz_-1 || l%nzmod==0)) {
|
---|
1787 | cout<<"l="<<l<<" z="<<z<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
|
---|
1788 | <<" dtrc="<<dtrc<<" dlum="<<dlum
|
---|
1789 | <<" -> cor="<<corr<<endl;
|
---|
1790 | }
|
---|
1791 | }
|
---|
1792 |
|
---|
1793 | double sum=0., sum2=0., nsum=0.;
|
---|
1794 | for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
|
---|
1795 | double a = lsigma*NorRand();
|
---|
1796 | a = magnref*pow(10.,a);
|
---|
1797 | // On met le meme tirage le long de Oz (indice k)
|
---|
1798 | for(long l=0;l<Nz_;l++) {
|
---|
1799 | int_8 ip = IndexR(i,j,l);
|
---|
1800 | data_[ip] += a*correction[l];
|
---|
1801 | }
|
---|
1802 | sum += a; sum2 += a*a; nsum += 1.;
|
---|
1803 | }
|
---|
1804 |
|
---|
1805 | if(lp_>0 && nsum>1.) {
|
---|
1806 | sum /= nsum;
|
---|
1807 | sum2 = sum2/nsum - sum*sum;
|
---|
1808 | cout<<"...Mean mass="<<sum<<" Msol , s^2="<<sum2<<" s="<<sqrt(fabs(sum2))<<endl;
|
---|
1809 | }
|
---|
1810 |
|
---|
1811 | }
|
---|
1812 | *********************************************************************/
|
---|