source: Sophya/trunk/Cosmo/SimLSS/genefluct3d.cc@ 3271

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

intro evolution bruit avec redshift cmv 20/06/2007
File size: 38.9 KB
Line 
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
19#include "arrctcast.h"
20
21#include "constcosmo.h"
22#include "geneutils.h"
23#include "schechter.h"
24
25#include "genefluct3d.h"
26
27//#define FFTW_THREAD
28
29#define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
30
31//-------------------------------------------------------
32GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
33 : Nx_(0) , Ny_(0) , Nz_(0)
34 , lp_(0)
35 , array_allocated_(false) , T_(T)
36 , cosmo_(NULL) , growth_(NULL)
37 , redsh_ref_(-999.), kredsh_ref_(0.), dred_ref_(-999.)
38 , loscom_ref_(-999.), dtrc_ref_(-999.), dlum_ref_(-999.), dang_ref_(-999.)
39 , nu_ref_(-999.), dnu_ref_ (-999.)
40 , loscom_min_(-999.), loscom_max_(-999.)
41 , loscom2zred_min_(0.) , loscom2zred_max_(0.)
42{
43 xobs_[0] = xobs_[1] = xobs_[2] = 0.;
44 zred_.resize(0);
45 loscom_.resize(0);
46 loscom2zred_.resize(0);
47 SetNThread();
48}
49
50GeneFluct3D::~GeneFluct3D(void)
51{
52 fftw_destroy_plan(pf_);
53 fftw_destroy_plan(pb_);
54#ifdef FFTW_THREAD
55 if(nthread_>0) fftw_cleanup_threads();
56#endif
57}
58
59//-------------------------------------------------------
60void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
61{
62 setsize(nx,ny,nz,dx,dy,dz);
63 setalloc();
64 setpointers(false);
65 init_fftw();
66}
67
68void GeneFluct3D::SetObservator(double redshref,double kredshref)
69// L'observateur est au redshift z=0
70// est situe sur la "perpendiculaire" a la face x,y
71// issue du centre de cette face
72// Il faut positionner le cube sur l'axe des z cad des redshifts:
73// redshref = redshift de reference
74// Si redshref<0 alors redshref=0
75// kredshref = indice (en double) correspondant a ce redshift
76// Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
77// Exemple: redshref=1.5 kredshref=250.75
78// -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
79{
80 if(redshref<0.) redshref = 0.;
81 if(kredshref<0.) {
82 if(Nz_<=0) {
83 char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 SetSize should be called first";
84 cout<<bla<<endl; throw ParmError(bla);
85 }
86 kredshref = Nz_/2.;
87 }
88 redsh_ref_ = redshref;
89 kredsh_ref_ = kredshref;
90 if(lp_>0)
91 cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
92}
93
94void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
95{
96 cosmo_ = &cosmo;
97 if(lp_>1) cosmo_->Print();
98}
99
100void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
101{
102 growth_ = &growth;
103}
104
105void GeneFluct3D::setsize(long nx,long ny,long nz,double dx,double dy,double dz)
106{
107 if(lp_>1) cout<<"--- GeneFluct3D::setsize: N="<<nx<<","<<ny<<","<<nz
108 <<" D="<<dx<<","<<dy<<","<<dz<<endl;
109 if(nx<=0 || dx<=0.) {
110 char *bla = "GeneFluct3D::setsize_Error: bad value(s) for nn/dx";
111 cout<<bla<<endl; throw ParmError(bla);
112 }
113
114 // Les tailles des tableaux
115 Nx_ = nx;
116 Ny_ = ny; if(Ny_ <= 0) Ny_ = Nx_;
117 Nz_ = nz; if(Nz_ <= 0) Nz_ = Nx_;
118 N_.resize(0); N_.push_back(Nx_); N_.push_back(Ny_); N_.push_back(Nz_);
119 NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
120 NCz_ = Nz_/2 +1;
121 NTz_ = 2*NCz_;
122
123 // le pas dans l'espace (Mpc)
124 Dx_ = dx;
125 Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
126 Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
127 D_.resize(0); D_.push_back(Dx_); D_.push_back(Dy_); D_.push_back(Dz_);
128 dVol_ = Dx_*Dy_*Dz_;
129 Vol_ = (Nx_*Dx_)*(Ny_*Dy_)*(Nz_*Dz_);
130
131 // Le pas dans l'espace de Fourier (Mpc^-1)
132 Dkx_ = 2.*M_PI/(Nx_*Dx_);
133 Dky_ = 2.*M_PI/(Ny_*Dy_);
134 Dkz_ = 2.*M_PI/(Nz_*Dz_);
135 Dk_.resize(0); Dk_.push_back(Dkx_); Dk_.push_back(Dky_); Dk_.push_back(Dkz_);
136 Dk3_ = Dkx_*Dky_*Dkz_;
137
138 // La frequence de Nyquist en k (Mpc^-1)
139 Knyqx_ = M_PI/Dx_;
140 Knyqy_ = M_PI/Dy_;
141 Knyqz_ = M_PI/Dz_;
142 Knyq_.resize(0); Knyq_.push_back(Knyqx_); Knyq_.push_back(Knyqy_); Knyq_.push_back(Knyqz_);
143}
144
145void GeneFluct3D::setalloc(void)
146{
147 if(lp_>1) cout<<"--- GeneFluct3D::setalloc ---"<<endl;
148 // Dimensionnement du tableau complex<r_8>
149 // ATTENTION: TArray adresse en memoire a l'envers du C
150 // Tarray(n1,n2,n3) == Carray[n3][n2][n1]
151 sa_size_t SzK_[3] = {NCz_,Ny_,Nx_}; // a l'envers
152 try {
153 T_.ReSize(3,SzK_);
154 array_allocated_ = true;
155 if(lp_>1) cout<<" allocating: "<<T_.Size()*sizeof(complex<r_8>)/1.e6<<" Mo"<<endl;
156 } catch (...) {
157 cout<<"GeneFluct3D::setalloc_Error: Problem allocating T_"<<endl;
158 }
159 T_.SetMemoryMapping(BaseArray::CMemoryMapping);
160}
161
162void GeneFluct3D::setpointers(bool from_real)
163{
164 if(lp_>1) cout<<"--- GeneFluct3D::setpointers ---"<<endl;
165 if(from_real) T_ = ArrCastR2C(R_);
166 else R_ = ArrCastC2R(T_);
167 // On remplit les pointeurs
168 fdata_ = (fftw_complex *) (&T_(0,0,0));
169 data_ = (double *) (&R_(0,0,0));
170}
171
172void GeneFluct3D::check_array_alloc(void)
173// Pour tester si le tableau T_ est alloue
174{
175 if(array_allocated_) return;
176 char bla[90];
177 sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
178 cout<<bla<<endl; throw ParmError(bla);
179}
180
181void GeneFluct3D::init_fftw(void)
182{
183 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
184 if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl;
185#ifdef FFTW_THREAD
186 if(nthread_>0) {
187 cout<<"...Computing with "<<nthread_<<" threads"<<endl;
188 fftw_init_threads();
189 fftw_plan_with_nthreads(nthread_);
190 }
191#endif
192 if(lp_>1) cout<<"...forward plan"<<endl;
193 pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
194 if(lp_>1) cout<<"...backward plan"<<endl;
195 pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
196}
197
198//-------------------------------------------------------
199long GeneFluct3D::LosComRedshift(double zinc,long npoints)
200// Given a position of the cube relative to the observer
201// and a cosmology
202// (SetObservator() and SetCosmology() should have been called !)
203// This routine filled:
204// the vector "zred_" of scanned redshift (by zinc increments)
205// the vector "loscom_" of corresponding los comoving distance
206// -- Input:
207// zinc : redshift increment for computation
208// npoints : number of points required for inverting loscom -> zred
209//
210{
211 if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
212
213 if(cosmo_ == NULL || redsh_ref_<0.) {
214 char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
215 cout<<bla<<endl; throw ParmError(bla);
216 }
217
218 // La distance angulaire/luminosite/Dnu au pixel de reference
219 dred_ref_ = Dz_/(cosmo_->Dhubble()/cosmo_->E(redsh_ref_));
220 loscom_ref_ = cosmo_->Dloscom(redsh_ref_);
221 dtrc_ref_ = cosmo_->Dtrcom(redsh_ref_);
222 dlum_ref_ = cosmo_->Dlum(redsh_ref_);
223 dang_ref_ = cosmo_->Dang(redsh_ref_);
224 nu_ref_ = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
225 dnu_ref_ = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
226 if(lp_>0) {
227 cout<<"...reference pixel redshref="<<redsh_ref_
228 <<", dredref="<<dred_ref_
229 <<", nuref="<<nu_ref_ <<" GHz"
230 <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
231 <<" dlosc="<<loscom_ref_<<" Mpc com"
232 <<", dtrc="<<dtrc_ref_<<" Mpc com"
233 <<", dlum="<<dlum_ref_<<" Mpc"
234 <<", dang="<<dang_ref_<<" Mpc"<<endl;
235 }
236
237 // On calcule les coordonnees de l'observateur dans le repere du cube
238 // cad dans le repere ou l'origine est au centre du pixel i=j=l=0.
239 // L'observateur est sur un axe centre sur le milieu de la face Oxy
240 xobs_[0] = Nx_/2.*Dx_;
241 xobs_[1] = Ny_/2.*Dy_;
242 xobs_[2] = kredsh_ref_*Dz_ - loscom_ref_;
243
244 // L'observateur est-il dans le cube?
245 bool obs_in_cube = false;
246 if(xobs_[2]>=0. && xobs_[2]<=Nz_*Dz_) obs_in_cube = true;
247
248 // Find MINIMUM los com distance to the observer:
249 // c'est le centre de la face a k=0
250 // (ou zero si l'observateur est dans le cube)
251 loscom_min_ = 0.;
252 if(!obs_in_cube) loscom_min_ = -xobs_[2];
253
254 // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
255 if(loscom_min_<=1.e-50)
256 for(int i=0;i<50;i++)
257 cout<<"ATTENTION TOUTES LES PARTIES DU CODE NE MARCHENT PAS POUR UN OBSERVATEUR DANS LE CUBE"<<endl;
258 // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
259
260
261 // Find MAXIMUM los com distance to the observer:
262 // ou que soit positionne l'observateur, la distance
263 // maximal est sur un des coins du cube
264 loscom_max_ = 0.;
265 for(long i=0;i<=1;i++) {
266 double dx2 = DXcom(i*(Nx_-1)); dx2 *= dx2;
267 for(long j=0;j<=1;j++) {
268 double dy2 = DYcom(j*(Ny_-1)); dy2 *= dy2;
269 for(long k=0;k<=1;k++) {
270 double dz2 = DZcom(k*(Nz_-1)); dz2 *= dz2;
271 dz2 = sqrt(dx2+dy2+dz2);
272 if(dz2>loscom_max_) loscom_max_ = dz2;
273 }
274 }
275 }
276 if(lp_>0) {
277 cout<<"...zref="<<redsh_ref_<<" kzref="<<kredsh_ref_<<" losref="<<loscom_ref_<<" Mpc\n"
278 <<" xobs="<<xobs_[0]<<" , "<<xobs_[1]<<" , "<<xobs_[2]<<" Mpc "
279 <<" in_cube="<<obs_in_cube
280 <<" loscom_min="<<loscom_min_<<" loscom_max="<<loscom_max_<<" Mpc (com)"<<endl;
281 }
282
283 // Fill the corresponding vectors for loscom and zred
284 // Be shure to have one dlc <loscom_min and one >loscom_max
285 if(zinc<=0.) zinc = 0.01;
286 for(double z=0.; ; z+=zinc) {
287 double dlc = cosmo_->Dloscom(z);
288 if(dlc<loscom_min_) {zred_.resize(0); loscom_.resize(0);}
289 zred_.push_back(z);
290 loscom_.push_back(dlc);
291 z += zinc;
292 if(dlc>loscom_max_) break; // on sort apres avoir stoque un dlc>dlcmax
293 }
294
295 if(lp_>0) {
296 long n = zred_.size();
297 cout<<"...zred/loscom tables[zinc="<<zinc<<"]: n="<<n;
298 if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
299 if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
300 cout<<endl;
301 }
302
303 // Compute the parameters and tables needed for inversion loscom->zred
304 if(npoints<3) npoints = zred_.size();
305 InverseFunc invfun(zred_,loscom_);
306 invfun.ComputeParab(npoints,loscom2zred_);
307 loscom2zred_min_ = invfun.YMin();
308 loscom2zred_max_ = invfun.YMax();
309
310 if(lp_>0) {
311 long n = loscom2zred_.size();
312 cout<<"...loscom -> zred[npoints="<<npoints<<"]: n="<<n
313 <<" los_min="<<loscom2zred_min_
314 <<" los_max="<<loscom2zred_max_
315 <<" -> zred=[";
316 if(n>0) cout<<loscom2zred_[0];
317 cout<<",";
318 if(n>1) cout<<loscom2zred_[n-1];
319 cout<<"]"<<endl;
320 if(lp_>1 && n>0)
321 for(int i=0;i<n;i++)
322 if(i==0 || abs(i-n/2)<2 || i==n-1)
323 cout<<" "<<i<<" "<<loscom2zred_[i]<<endl;
324 }
325
326 return zred_.size();
327}
328
329//-------------------------------------------------------
330void GeneFluct3D::WriteFits(string cfname,int bitpix)
331{
332 cout<<"--- GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
333 try {
334 FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
335 fwrt.WriteKey("NX",Nx_," axe transverse 1");
336 fwrt.WriteKey("NY",Ny_," axe transverse 2");
337 fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");
338 fwrt.WriteKey("DX",Dx_," Mpc");
339 fwrt.WriteKey("DY",Dy_," Mpc");
340 fwrt.WriteKey("DZ",Dz_," Mpc");
341 fwrt.WriteKey("DKX",Dkx_," Mpc^-1");
342 fwrt.WriteKey("DKY",Dky_," Mpc^-1");
343 fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
344 fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
345 fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
346 fwrt.Write(R_);
347 } catch (PThrowable & exc) {
348 cout<<"Exception : "<<(string)typeid(exc).name()
349 <<" - Msg= "<<exc.Msg()<<endl;
350 return;
351 } catch (...) {
352 cout<<" some other exception was caught !"<<endl;
353 return;
354 }
355}
356
357void GeneFluct3D::ReadFits(string cfname)
358{
359 cout<<"--- GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
360 try {
361 FitsImg3DRead fimg(cfname.c_str(),0,5);
362 fimg.Read(R_);
363 long nx = fimg.ReadKeyL("NX");
364 long ny = fimg.ReadKeyL("NY");
365 long nz = fimg.ReadKeyL("NZ");
366 double dx = fimg.ReadKey("DX");
367 double dy = fimg.ReadKey("DY");
368 double dz = fimg.ReadKey("DZ");
369 double zref = fimg.ReadKey("ZREF");
370 double kzref = fimg.ReadKey("KZREF");
371 setsize(nx,ny,nz,dx,dy,dz);
372 setpointers(true);
373 init_fftw();
374 SetObservator(zref,kzref);
375 } catch (PThrowable & exc) {
376 cout<<"Exception : "<<(string)typeid(exc).name()
377 <<" - Msg= "<<exc.Msg()<<endl;
378 return;
379 } catch (...) {
380 cout<<" some other exception was caught !"<<endl;
381 return;
382 }
383}
384
385void GeneFluct3D::WritePPF(string cfname,bool write_real)
386// On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> >
387{
388 cout<<"--- GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
389 try {
390 R_.Info()["NX"] = (int_8)Nx_;
391 R_.Info()["NY"] = (int_8)Ny_;
392 R_.Info()["NZ"] = (int_8)Nz_;
393 R_.Info()["DX"] = (r_8)Dx_;
394 R_.Info()["DY"] = (r_8)Dy_;
395 R_.Info()["DZ"] = (r_8)Dz_;
396 R_.Info()["ZREF"] = (r_8)redsh_ref_;
397 R_.Info()["KZREF"] = (r_8)kredsh_ref_;
398 POutPersist pos(cfname.c_str());
399 if(write_real) pos << PPFNameTag("rgen") << R_;
400 else pos << PPFNameTag("pkgen") << T_;
401 } catch (PThrowable & exc) {
402 cout<<"Exception : "<<(string)typeid(exc).name()
403 <<" - Msg= "<<exc.Msg()<<endl;
404 return;
405 } catch (...) {
406 cout<<" some other exception was caught !"<<endl;
407 return;
408 }
409}
410
411void GeneFluct3D::ReadPPF(string cfname)
412{
413 cout<<"--- GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
414 try {
415 bool from_real = true;
416 PInPersist pis(cfname.c_str());
417 string name_tag_k = "pkgen";
418 bool found_tag_k = pis.GotoNameTag("pkgen");
419 if(found_tag_k) {
420 cout<<" ...reading spectrum into TArray<complex <r_8> >"<<endl;
421 pis >> PPFNameTag("pkgen") >> T_;
422 from_real = false;
423 } else {
424 cout<<" ...reading space into TArray<r_8>"<<endl;
425 pis >> PPFNameTag("rgen") >> R_;
426 }
427 setpointers(from_real); // a mettre ici pour relire les DVInfo
428 int_8 nx = R_.Info()["NX"];
429 int_8 ny = R_.Info()["NY"];
430 int_8 nz = R_.Info()["NZ"];
431 r_8 dx = R_.Info()["DX"];
432 r_8 dy = R_.Info()["DY"];
433 r_8 dz = R_.Info()["DZ"];
434 r_8 zref = R_.Info()["ZREF"];
435 r_8 kzref = R_.Info()["KZREF"];
436 setsize(nx,ny,nz,dx,dy,dz);
437 init_fftw();
438 SetObservator(zref,kzref);
439 } catch (PThrowable & exc) {
440 cout<<"Exception : "<<(string)typeid(exc).name()
441 <<" - Msg= "<<exc.Msg()<<endl;
442 return;
443 } catch (...) {
444 cout<<" some other exception was caught !"<<endl;
445 return;
446 }
447}
448
449//-------------------------------------------------------
450void GeneFluct3D::Print(void)
451{
452 cout<<"GeneFluct3D(T_alloc="<<array_allocated_<<"):"<<endl;
453 cout<<"Space Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<Nz_<<" ("<<NTz_<<") size="
454 <<NRtot_<<endl;
455 cout<<" Resol: dx="<<Dx_<<" dy="<<Dy_<<" dz="<<Dz_<<" Mpc"
456 <<", dVol="<<dVol_<<", Vol="<<Vol_<<" Mpc^3"<<endl;
457 cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<endl;
458 cout<<" Resol: dkx="<<Dkx_<<" dky="<<Dky_<<" dkz="<<Dkz_<<" Mpc^-1"
459 <<", Dk3="<<Dk3_<<" Mpc^-3"<<endl;
460 cout<<" (2Pi/k: "<<2.*M_PI/Dkx_<<" "<<2.*M_PI/Dky_<<" "<<2.*M_PI/Dkz_<<" Mpc)"<<endl;
461 cout<<" Nyquist: kx="<<Knyqx_<<" ky="<<Knyqy_<<" kz="<<Knyqz_<<" Mpc^-1"
462 <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
463 cout<<" (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
464 cout<<"Redshift "<<redsh_ref_<<" for z axe at k="<<kredsh_ref_<<endl;
465}
466
467//-------------------------------------------------------
468void GeneFluct3D::ComputeFourier0(GenericFunc& pk_at_z)
469// cf ComputeFourier() mais avec autre methode de realisation du spectre
470// (attention on fait une fft pour realiser le spectre)
471{
472
473 // --- realisation d'un tableau de tirage gaussiens
474 if(lp_>0) cout<<"--- ComputeFourier0: before gaussian filling ---"<<endl;
475 // On tient compte du pb de normalisation de FFTW3
476 double sntot = sqrt((double)NRtot_);
477 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
478 int_8 ip = IndexR(i,j,l);
479 data_[ip] = NorRand()/sntot;
480 }
481
482 // --- realisation d'un tableau de tirage gaussiens
483 if(lp_>0) cout<<"...before fft real ---"<<endl;
484 fftw_execute(pf_);
485
486 // --- On remplit avec une realisation
487 if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
488 T_(0,0,0) = complex<r_8>(0.); // on coupe le continue et on l'initialise
489 long lmod = Nx_/10; if(lmod<1) lmod=1;
490 for(long i=0;i<Nx_;i++) {
491 long ii = (i>Nx_/2) ? Nx_-i : i;
492 double kx = ii*Dkx_; kx *= kx;
493 if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
494 for(long j=0;j<Ny_;j++) {
495 long jj = (j>Ny_/2) ? Ny_-j : j;
496 double ky = jj*Dky_; ky *= ky;
497 for(long l=0;l<NCz_;l++) {
498 double kz = l*Dkz_; kz *= kz;
499 if(i==0 && j==0 && l==0) continue; // Suppression du continu
500 double k = sqrt(kx+ky+kz);
501 // cf normalisation: Peacock, Cosmology, formule 16.38 p504
502 double pk = pk_at_z(k)/Vol_;
503 // ici pas de "/2" a cause de la remarque ci-dessus
504 T_(l,j,i) *= sqrt(pk);
505 }
506 }
507 }
508
509 if(lp_>0) cout<<"...computing power"<<endl;
510 double p = compute_power_carte();
511 if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;
512
513}
514
515//-------------------------------------------------------
516void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
517// Calcule une realisation du spectre "pk_at_z"
518// Attention: dans TArray le premier indice varie le + vite
519// Explication normalisation: see Coles & Lucchin, Cosmology, p264-265
520// FFTW3: on note N=Nx*Ny*Nz
521// f --(FFT)--> F = TF(f) --(FFT^-1)--> fb = TF^-1(F) = TF^-1(TF(f))
522// sum(f(x_i)^2) = S
523// sum(F(nu_i)^2) = S*N
524// sum(fb(x_i)^2) = S*N^2
525{
526 // --- RaZ du tableau
527 T_ = complex<r_8>(0.);
528
529 // --- On remplit avec une realisation
530 if(lp_>0) cout<<"--- ComputeFourier ---"<<endl;
531 long lmod = Nx_/10; if(lmod<1) lmod=1;
532 for(long i=0;i<Nx_;i++) {
533 long ii = (i>Nx_/2) ? Nx_-i : i;
534 double kx = ii*Dkx_; kx *= kx;
535 if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
536 for(long j=0;j<Ny_;j++) {
537 long jj = (j>Ny_/2) ? Ny_-j : j;
538 double ky = jj*Dky_; ky *= ky;
539 for(long l=0;l<NCz_;l++) {
540 double kz = l*Dkz_; kz *= kz;
541 if(i==0 && j==0 && l==0) continue; // Suppression du continu
542 double k = sqrt(kx+ky+kz);
543 // cf normalisation: Peacock, Cosmology, formule 16.38 p504
544 double pk = pk_at_z(k)/Vol_;
545 // Explication de la division par 2: voir perandom.cc
546 // ou egalement Coles & Lucchin, Cosmology formula 13.7.2 p279
547 T_(l,j,i) = ComplexGaussRan(sqrt(pk/2.));
548 }
549 }
550 }
551
552 manage_coefficients(); // gros effet pour les spectres que l'on utilise !
553
554 if(lp_>0) cout<<"...computing power"<<endl;
555 double p = compute_power_carte();
556 if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl;
557
558}
559
560long GeneFluct3D::manage_coefficients(void)
561// Take into account the real and complexe conjugate coefficients
562// because we want a realization of a real data in real space
563{
564 if(lp_>1) cout<<"...managing coefficients"<<endl;
565 check_array_alloc();
566
567 // 1./ Le Continu et Nyquist sont reels
568 long nreal = 0;
569 for(long kk=0;kk<2;kk++) {
570 long k=0; // continu
571 if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
572 for(long jj=0;jj<2;jj++) {
573 long j=0;
574 if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
575 for(long ii=0;ii<2;ii++) {
576 long i=0;
577 if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
578 int_8 ip = IndexC(i,j,k);
579 //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata_[ip][0]<<","<<fdata_[ip][1]<<")"<<endl;
580 fdata_[ip][1] = 0.; fdata_[ip][0] *= M_SQRT2;
581 nreal++;
582 }
583 }
584 }
585 if(lp_>1) cout<<"Number of forced real number ="<<nreal<<endl;
586
587 // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist)
588
589 // a./ les lignes et colonnes du continu et de nyquist
590 long nconj1 = 0;
591 for(long kk=0;kk<2;kk++) {
592 long k=0; // continu
593 if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
594 for(long jj=0;jj<2;jj++) { // selon j
595 long j=0;
596 if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
597 for(long i=1;i<(Nx_+1)/2;i++) {
598 int_8 ip = IndexC(i,j,k);
599 int_8 ip1 = IndexC(Nx_-i,j,k);
600 fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
601 nconj1++;
602 }
603 }
604 for(long ii=0;ii<2;ii++) {
605 long i=0;
606 if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
607 for(long j=1;j<(Ny_+1)/2;j++) {
608 int_8 ip = IndexC(i,j,k);
609 int_8 ip1 = IndexC(i,Ny_-j,k);
610 fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
611 nconj1++;
612 }
613 }
614 }
615 if(lp_>1) cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;
616
617 // b./ les lignes et colonnes hors continu et de nyquist
618 long nconj2 = 0;
619 for(long kk=0;kk<2;kk++) {
620 long k=0; // continu
621 if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
622 for(long j=1;j<(Ny_+1)/2;j++) {
623 if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
624 for(long i=1;i<Nx_;i++) {
625 if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
626 int_8 ip = IndexC(i,j,k);
627 int_8 ip1 = IndexC(Nx_-i,Ny_-j,k);
628 fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
629 nconj2++;
630 }
631 }
632 }
633 if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
634
635 if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
636
637 return nreal+nconj1+nconj2;
638}
639
640double GeneFluct3D::compute_power_carte(void)
641// Calcul la puissance de la realisation du spectre Pk
642{
643 check_array_alloc();
644
645 double s2 = 0.;
646 for(long l=0;l<NCz_;l++)
647 for(long j=0;j<Ny_;j++)
648 for(long i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));
649
650 double s20 = 0.;
651 for(long j=0;j<Ny_;j++)
652 for(long i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));
653
654 double s2n = 0.;
655 if(Nz_%2==0)
656 for(long j=0;j<Ny_;j++)
657 for(long i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));
658
659 return 2.*s2 -s20 -s2n;
660}
661
662//-------------------------------------------------------------------
663void GeneFluct3D::FilterByPixel(void)
664// Filtrage par la fonction fenetre du pixel (parallelepipede)
665// TF = 1/(dx*dy*dz)*Int[{-dx/2,dx/2},{-dy/2,dy/2},{-dz/2,dz/2}]
666// e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
667// = 2/(k_x*dx) * sin(k_x*dx/2) * (idem y) * (idem z)
668// Gestion divergence en 0: sin(y)/y = 1 - y^2/6*(1-y^2/20)
669// avec y = k_x*dx/2
670{
671 if(lp_>0) cout<<"--- FilterByPixel ---"<<endl;
672 check_array_alloc();
673
674 for(long i=0;i<Nx_;i++) {
675 long ii = (i>Nx_/2) ? Nx_-i : i;
676 double kx = ii*Dkx_ *Dx_/2;
677 double pk_x = pixelfilter(kx);
678 for(long j=0;j<Ny_;j++) {
679 long jj = (j>Ny_/2) ? Ny_-j : j;
680 double ky = jj*Dky_ *Dy_/2;
681 double pk_y = pixelfilter(ky);
682 for(long l=0;l<NCz_;l++) {
683 double kz = l*Dkz_ *Dz_/2;
684 double pk_z = pixelfilter(kz);
685 T_(l,j,i) *= pk_x*pk_y*pk_z;
686 }
687 }
688 }
689
690}
691
692//-------------------------------------------------------------------
693void GeneFluct3D::ApplyGrowthFactor(void)
694// Apply Growth to real space
695// Using the correspondance between redshift and los comoving distance
696// describe in vector "zred_" "loscom_"
697{
698 if(lp_>0) cout<<"--- ApplyGrowthFactor ---"<<endl;
699 check_array_alloc();
700
701 if(growth_ == NULL) {
702 char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
703 cout<<bla<<endl; throw ParmError(bla);
704 }
705
706 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
707 unsigned short ok;
708
709 //CHECK: Histo hgr(0.9*zred_[0],1.1*zred_[n-1],1000);
710 for(long i=0;i<Nx_;i++) {
711 double dx2 = xobs_[0] - i*Dx_; dx2 *= dx2;
712 for(long j=0;j<Ny_;j++) {
713 double dy2 = xobs_[1] - j*Dy_; dy2 *= dy2;
714 for(long l=0;l<Nz_;l++) {
715 double dz2 = xobs_[2] - l*Dz_; dz2 *= dz2;
716 dz2 = sqrt(dx2+dy2+dz2);
717 double z = interpinv(dz2);
718 //CHECK: hgr.Add(z);
719 double dzgr = (*growth_)(z); // interpolation par morceau
720 //double dzgr = growth_->Linear(z,ok); // interpolation lineaire
721 //double dzgr = growth_->Parab(z,ok); // interpolation parabolique
722 int_8 ip = IndexR(i,j,l);
723 data_[ip] *= dzgr;
724 }
725 }
726 }
727
728 //CHECK: {POutPersist pos("applygrowth.ppf"); string tag="hgr"; pos.PutObject(hgr,tag);}
729
730}
731
732//-------------------------------------------------------------------
733void GeneFluct3D::ComputeReal(void)
734// Calcule une realisation dans l'espace reel
735{
736 if(lp_>0) cout<<"--- ComputeReal ---"<<endl;
737 check_array_alloc();
738
739 // On fait la FFT
740 fftw_execute(pb_);
741}
742
743//-------------------------------------------------------------------
744void GeneFluct3D::ReComputeFourier(void)
745{
746 if(lp_>0) cout<<"--- ReComputeFourier ---"<<endl;
747 check_array_alloc();
748
749 // On fait la FFT
750 fftw_execute(pf_);
751 // On corrige du pb de la normalisation de FFTW3
752 double v = (double)NRtot_;
753 for(long i=0;i<Nx_;i++)
754 for(long j=0;j<Ny_;j++)
755 for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
756
757}
758
759//-------------------------------------------------------------------
760int GeneFluct3D::ComputeSpectrum(HistoErr& herr)
761// Compute spectrum from "T" and fill HistoErr "herr"
762// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
763// cad T(kz,ky,kx) avec 0<kz<kz_nyq -ky_nyq<ky<ky_nyq -kx_nyq<kx<kx_nyq
764{
765 if(lp_>0) cout<<"--- ComputeSpectrum ---"<<endl;
766 check_array_alloc();
767
768 if(herr.NBins()<0) return -1;
769 herr.Zero();
770
771 // Attention a l'ordre
772 for(long i=0;i<Nx_;i++) {
773 long ii = (i>Nx_/2) ? Nx_-i : i;
774 double kx = ii*Dkx_; kx *= kx;
775 for(long j=0;j<Ny_;j++) {
776 long jj = (j>Ny_/2) ? Ny_-j : j;
777 double ky = jj*Dky_; ky *= ky;
778 for(long l=0;l<NCz_;l++) {
779 double kz = l*Dkz_; kz *= kz;
780 double k = sqrt(kx+ky+kz);
781 double pk = MODULE2(T_(l,j,i));
782 herr.Add(k,pk);
783 }
784 }
785 }
786 herr.ToVariance();
787
788 // renormalize to directly compare to original spectrum
789 double norm = Vol_;
790 herr *= norm;
791
792 return 0;
793}
794
795int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
796{
797 if(lp_>0) cout<<"--- ComputeSpectrum2D ---"<<endl;
798 check_array_alloc();
799
800 if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
801 herr.Zero();
802
803 // Attention a l'ordre
804 for(long i=0;i<Nx_;i++) {
805 long ii = (i>Nx_/2) ? Nx_-i : i;
806 double kx = ii*Dkx_; kx *= kx;
807 for(long j=0;j<Ny_;j++) {
808 long jj = (j>Ny_/2) ? Ny_-j : j;
809 double ky = jj*Dky_; ky *= ky;
810 double kt = sqrt(kx+ky);
811 for(long l=0;l<NCz_;l++) {
812 double kz = l*Dkz_;
813 double pk = MODULE2(T_(l,j,i));
814 herr.Add(kt,kz,pk);
815 }
816 }
817 }
818 herr.ToVariance();
819
820 // renormalize to directly compare to original spectrum
821 double norm = Vol_;
822 herr *= norm;
823
824 return 0;
825}
826
827//-------------------------------------------------------
828int_8 GeneFluct3D::VarianceFrReal(double R,double& var)
829// Recompute MASS variance in spherical top-hat (rayon=R)
830{
831 if(lp_>0) cout<<"--- VarianceFrReal R="<<R<<endl;
832 check_array_alloc();
833
834 long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
835 long dny = long(R/Dy_+0.5); if(dny<=0) dny = 1;
836 long dnz = long(R/Dz_+0.5); if(dnz<=0) dnz = 1;
837 if(lp_>0) cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
838
839 double sum=0., sum2=0., r2 = R*R; int_8 nsum=0;
840
841 for(long i=dnx;i<Nx_-dnx;i+=dnx) {
842 for(long j=dny;j<Ny_-dny;j+=dny) {
843 for(long l=dnz;l<Nz_-dnz;l+=dnz) {
844 double s=0.; int_8 n=0;
845 for(long ii=i-dnx;ii<=i+dnx;ii++) {
846 double x = (ii-i)*Dx_; x *= x;
847 for(long jj=j-dny;jj<=j+dny;jj++) {
848 double y = (jj-j)*Dy_; y *= y;
849 for(long ll=l-dnz;ll<=l+dnz;ll++) {
850 double z = (ll-l)*Dz_; z *= z;
851 if(x+y+z>r2) continue;
852 int_8 ip = IndexR(ii,jj,ll);
853 s += 1.+data_[ip];
854 n++;
855 }
856 }
857 }
858 if(n>0) {sum += s; sum2 += s*s; nsum++;}
859 //cout<<i<<","<<j<<","<<l<<" n="<<n<<" s="<<s<<" sum="<<sum<<" sum2="<<sum2<<endl;
860 }
861 }
862 }
863
864 if(nsum<=1) {var=0.; return nsum;}
865
866 sum /= nsum;
867 sum2 = sum2/nsum - sum*sum;
868 if(lp_>0) cout<<"...nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
869
870 var = sum2/(sum*sum); // <dM>^2/<M>^2
871 if(lp_>0) cout<<"...sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
872
873 return nsum;
874}
875
876//-------------------------------------------------------
877int_8 GeneFluct3D::NumberOfBad(double vmin,double vmax)
878// number of pixels outside of ]vmin,vmax[ extremites exclues
879// -> vmin and vmax are considered as bad
880{
881 check_array_alloc();
882
883 int_8 nbad = 0;
884 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
885 int_8 ip = IndexR(i,j,l);
886 double v = data_[ip];
887 if(v<=vmin || v>=vmax) nbad++;
888 }
889
890 if(lp_>0) cout<<"--- NumberOfBad "<<nbad<<" px out of ]"<<vmin<<","<<vmax<<"["<<endl;
891 return nbad;
892}
893
894int_8 GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax
895 ,bool useout,double vout)
896// Calcul de mean,sigma2 dans le cube reel avec valeurs ]vmin,vmax[ extremites exclues
897// useout = false: ne pas utiliser les pixels hors limites pour calculer mean,sigma2
898// true : utiliser les pixels hors limites pour calculer mean,sigma2
899// en remplacant leurs valeurs par "vout"
900{
901 bool tstval = (vmax>vmin)? true: false;
902 if(lp_>0) {
903 cout<<"--- MeanSigma2";
904 if(tstval) cout<<" range=]"<<vmin<<","<<vmax<<"[";
905 if(useout) cout<<", useout="<<useout<<" vout="<<vout;
906 cout<<endl;
907 }
908 check_array_alloc();
909
910 int_8 n = 0;
911 rm = rs2 = 0.;
912
913 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
914 int_8 ip = IndexR(i,j,l);
915 double v = data_[ip];
916 if(tstval) {
917 if(v<=vmin || v>=vmax) {if(useout) v=vout; else continue;}
918 }
919 rm += v;
920 rs2 += v*v;
921 n++;
922 }
923
924 if(n>1) {
925 rm /= (double)n;
926 rs2 = rs2/(double)n - rm*rm;
927 }
928
929 if(lp_>0) cout<<" n="<<n<<" m="<<rm<<" s2="<<rs2<<" s="<<sqrt(fabs(rs2))<<endl;
930
931 return n;
932}
933
934int_8 GeneFluct3D::SetToVal(double vmin, double vmax,double val0)
935// set to "val0" if out of range ]vmin,vmax[ extremites exclues
936// cad set to "val0" if in [vmin,vmax] -> vmin and vmax are set to val0
937{
938 check_array_alloc();
939
940 int_8 nbad = 0;
941 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
942 int_8 ip = IndexR(i,j,l);
943 double v = data_[ip];
944 if(v<=vmin || v>=vmax) {data_[ip] = val0; nbad++;}
945 }
946
947 if(lp_>0) cout<<"--- SetToVal "<<nbad<<" px set to="<<val0
948 <<" because out of range=]"<<vmin<<","<<vmax<<"["<<endl;
949 return nbad;
950}
951
952//-------------------------------------------------------
953void GeneFluct3D::TurnFluct2Mass(void)
954// d_rho/rho -> Mass (add one!)
955{
956 if(lp_>0) cout<<"--- TurnFluct2Mass ---"<<endl;
957 check_array_alloc();
958
959
960 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
961 int_8 ip = IndexR(i,j,l);
962 data_[ip] += 1.;
963 }
964}
965
966double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
967// do NOT treate negative or nul values
968{
969 if(lp_>0) cout<<"--- TurnMass2MeanNumber ---"<<endl;
970
971 double mall=0., mgood=0.;
972 int_8 nall=0, ngood=0;
973 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
974 int_8 ip = IndexR(i,j,l);
975 mall += data_[ip]; nall++;
976 if(data_[ip]>0.) {mgood += data_[ip]; ngood++;}
977 }
978 if(ngood>0) mgood /= (double)ngood;
979 if(nall>0) mall /= (double)nall;
980 if(lp_>0) cout<<"...ngood="<<ngood<<" mgood="<<mgood
981 <<", nall="<<nall<<" mall="<<mall<<endl;
982 if(ngood<=0 || mall<=0.) {
983 cout<<"TurnMass2MeanNumber_Error: ngood="<<ngood<<" <=0 || mall="<<mall<<" <=0"<<endl;
984 throw RangeCheckError("TurnMass2MeanNumber_Error: ngood<=0 || mall<=0");
985 }
986
987 // On doit mettre n*Vol galaxies dans notre survey
988 // On en met uniquement dans les pixels de masse >0.
989 // On NE met PAS a zero les pixels <0
990 // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies
991 // comme on ne prend que les pixels >0, on doit normaliser
992 // a la moyenne de <1+d_rho/rho> sur ces pixels
993 // (rappel sur tout les pixels <1+d_rho/rho>=1)
994 // nb de gal a mettre ds 1 px:
995 double dn = n_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
996 if(lp_>0) cout<<"...galaxy density move from "
997 <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
998
999 double sum = 0.;
1000 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
1001 int_8 ip = IndexR(i,j,l);
1002 data_[ip] *= dn; // par coherence on multiplie aussi les <=0
1003 if(data_[ip]>0.) sum += data_[ip]; // mais on ne les compte pas
1004 }
1005
1006 if(lp_>0) cout<<sum<<"...galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
1007
1008 return sum;
1009}
1010
1011double GeneFluct3D::ApplyPoisson(void)
1012// do NOT treate negative or nul mass -> let it as it is
1013{
1014 if(lp_>0) cout<<"--- ApplyPoisson ---"<<endl;
1015 check_array_alloc();
1016
1017 double sum = 0.;
1018 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
1019 int_8 ip = IndexR(i,j,l);
1020 double v = data_[ip];
1021 if(v>0.) {
1022 unsigned long dn = PoissRandLimit(v,10.);
1023 data_[ip] = (double)dn;
1024 sum += (double)dn;
1025 }
1026 }
1027 if(lp_>0) cout<<sum<<" galaxies put into survey"<<endl;
1028
1029 return sum;
1030}
1031
1032double GeneFluct3D::TurnNGal2Mass(FunRan& massdist,bool axeslog)
1033// do NOT treate negative or nul mass -> let it as it is
1034// INPUT:
1035// massdist : distribution de masse (m*dn/dm)
1036// axeslog = false : retourne la masse
1037// = true : retourne le log10(mass)
1038// RETURN la masse totale
1039{
1040 if(lp_>0) cout<<"--- TurnNGal2Mass ---"<<endl;
1041 check_array_alloc();
1042
1043 double sum = 0.;
1044 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
1045 int_8 ip = IndexR(i,j,l);
1046 double v = data_[ip];
1047 if(v>0.) {
1048 long ngal = long(v+0.1);
1049 data_[ip] = 0.;
1050 for(long i=0;i<ngal;i++) {
1051 double m = massdist.RandomInterp(); // massdist.Random();
1052 if(axeslog) m = pow(10.,m);
1053 data_[ip] += m;
1054 }
1055 sum += data_[ip];
1056 }
1057 }
1058 if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;
1059
1060 return sum;
1061}
1062
1063void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
1064// Add AGN flux into simulation:
1065// --- Procedure:
1066// 1. lancer "cmvdefsurv" avec les parametres du survey
1067// (au redshift de reference du survey)
1068// et recuperer l'angle solide "angsol sr" du pixel elementaire
1069// au centre du cube.
1070// 2. lancer "cmvtstagn" pour cet angle solide -> cmvtstagn.ppf
1071// 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
1072// cad une moyenne <log10(S)> et un sigma "sig"
1073// Attention: la distribution n'est pas gaussienne les "mean,sigma"
1074// de l'histo ne sont pas vraiment ce que l'on veut
1075// --- Limitations actuelle du code:
1076// . les AGN sont supposes evoluer avec la meme loi de puissance pour tout theta,phi
1077// . le flux des AGN est mis dans une colonne Oz (indice k) et pas sur la ligne de visee
1078// . la distribution est approximee a une gaussienne
1079// ... C'est une approximation pour un observateur loin du centre du cube
1080// et pour un cube peu epais / distance observateur
1081// --- Parametres de la routine:
1082// llfy : c'est le <log10(S)> du flux depose par les AGN
1083// dans l'angle solide du pixel elementaire de reference du cube
1084// lsigma : c'est le sigma de la distribution des log10(S)
1085// powlaw : c'est la pente de la distribution cad que le flux "lmsol"
1086// et considere comme le flux a 1.4GHz et qu'on suppose une loi
1087// F(nu) = (1.4GHz/nu)^powlaw * F(1.4GHz)
1088// - Comme on est en echelle log10():
1089// on tire log10(Msol) + X
1090// ou X est une realisation sur une gaussienne de variance "sig^2"
1091// La masse realisee est donc: Msol*10^X
1092// - Pas de probleme de pixel negatif car on a une multiplication!
1093{
1094 if(lp_>0) cout<<"--- AddAGN: <log10(S Jy)> = "<<lfjy<<" , sigma = "<<lsigma<<endl;
1095 check_array_alloc();
1096
1097 if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
1098 char *bla = "GeneFluct3D::AddAGN_Error: set Observator and Cosmology first";
1099 cout<<bla<<endl; throw ParmError(bla);
1100 }
1101
1102 // Le flux des AGN en Jy et en mass solaire
1103 double fagnref = pow(10.,lfjy)*(dnu_ref_*1.e9); // Jy.Hz = W/m^2
1104 double magnref = FluxHI2Msol(fagnref*Jansky2Watt_cst,dlum_ref_); // Msol
1105 if(lp_>0)
1106 cout<<"Au pixel de ref: fagnref="<<fagnref
1107 <<" Jy.Hz (a 1.4GHz), magnref="<<magnref<<" Msol"<<endl;
1108
1109 if(powlaw!=0.) {
1110 // 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)
1111 magnref *= pow(1/(1.+redsh_ref_),powlaw);
1112 if(lp_>0) cout<<" powlaw="<<powlaw<<" -> change magnref to "<<magnref<<" Msol"<<endl;
1113 }
1114
1115 // Les infos en fonction de l'indice "l" selon Oz
1116 vector<double> correction;
1117 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
1118 long nzmod = ((Nz_>10)?Nz_/10:1);
1119 for(long l=0;l<Nz_;l++) {
1120 double z = fabs(DZcom(l));
1121 double zred = interpinv(z);
1122 double dtrc = cosmo_->Dtrcom(zred); // pour variation angle solide
1123 double dlum = cosmo_->Dlum(zred); // pour variation conversion mass HI
1124 double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
1125 double dnu = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
1126 // on a: Mass ~ DNu * Dlum^2 / Dtrcom^2
1127 double corr = dnu/dnu_ref_*pow(dtrc_ref_/dtrc*dlum/dlum_ref_,2.);
1128 // 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)
1129 if(powlaw!=0.) corr *= pow((1.+redsh_ref_)/(1.+zred),powlaw);
1130 correction.push_back(corr);
1131 if(lp_>0 && (l==0 || l==Nz_-1 || l%nzmod==0)) {
1132 cout<<"l="<<l<<" z="<<z<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
1133 <<" dtrc="<<dtrc<<" dlum="<<dlum
1134 <<" -> cor="<<corr<<endl;
1135 }
1136 }
1137
1138 double sum=0., sum2=0., nsum=0.;
1139 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
1140 double a = lsigma*NorRand();
1141 a = magnref*pow(10.,a);
1142 // On met le meme tirage le long de Oz (indice k)
1143 for(long l=0;l<Nz_;l++) {
1144 int_8 ip = IndexR(i,j,l);
1145 data_[ip] += a*correction[l];
1146 }
1147 sum += a; sum2 += a*a; nsum += 1.;
1148 }
1149
1150 if(lp_>0 && nsum>1.) {
1151 sum /= nsum;
1152 sum2 = sum2/nsum - sum*sum;
1153 cout<<"...Mean mass="<<sum<<" Msol , s^2="<<sum2<<" s="<<sqrt(fabs(sum2))<<endl;
1154 }
1155
1156}
1157
1158void GeneFluct3D::AddNoise2Real(double snoise,bool with_evol)
1159// add noise to every pixels (meme les <=0 !)
1160{
1161 if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<with_evol<<endl;
1162 check_array_alloc();
1163
1164 vector<double> correction;
1165 InterpFunc *intercor = NULL;
1166
1167 if(with_evol) {
1168 // Sigma_Noise(en mass) :
1169 // Slim ~ 1/sqrt(DNu) * sqrt(nlobe) en W/m^2Hz
1170 // Flim ~ sqrt(DNu) * sqrt(nlobe) en W/m^2
1171 // Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe) en Msol
1172 // nlobe ~ 1/Dtrcom^2
1173 // Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
1174 if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
1175 char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
1176 cout<<bla<<endl; throw ParmError(bla);
1177 }
1178 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
1179 long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
1180 for(long i=0;i<nsz;i++) {
1181 double d = interpinv.X(i);
1182 double zred = interpinv(d);
1183 double dtrc = cosmo_->Dtrcom(zred); // pour variation angle solide
1184 double dlum = cosmo_->Dlum(zred); // pour variation conversion mass HI
1185 double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
1186 double dnu = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
1187 double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
1188 if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
1189 cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
1190 <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
1191 correction.push_back(corr);
1192 }
1193 intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
1194 }
1195
1196 double dx2=0., dy2=0., dz2=0.;
1197 for(long i=0;i<Nx_;i++) {
1198 dx2 = DXcom(i); dx2 *= dx2;
1199 for(long j=0;j<Ny_;j++) {
1200 dy2 = DYcom(j); dy2 *= dy2;
1201 for(long l=0;l<Nz_;l++) {
1202 double corr = 1.;
1203 if(with_evol) {
1204 dz2 = DZcom(l); dz2 *= dz2; dz2 = sqrt(dx2+dy2+dz2);
1205 corr = (*intercor)(dz2);
1206 }
1207 int_8 ip = IndexR(i,j,l);
1208 data_[ip] += snoise*corr*NorRand();
1209 }
1210 }
1211 }
1212
1213 if(intercor!=NULL) delete intercor;
1214}
Note: See TracBrowser for help on using the repository browser.