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

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

modifs d'options de cmvobserv3d.cc
creation de init_fftw() et deplacement a l'initialisation
creation de SetObservator pour positionner l'observateur + ecriture FITS et PPF

(mais le codage des distance/redshift reste a coder)

cmv 19/01/2007
File size: 25.8 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 "timing.h"
11#include "tarray.h"
12#include "pexceptions.h"
13#include "perandom.h"
14#include "srandgen.h"
15
16#include "fabtcolread.h"
17#include "fabtwriter.h"
18#include "fioarr.h"
19
20#include "arrctcast.h"
21
22#include "constcosmo.h"
23#include "integfunc.h"
24#include "geneutils.h"
25
26#include "genefluct3d.h"
27
28//#define FFTW_THREAD
29
30#define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
31
32//-------------------------------------------------------
33GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
34 : T_(T) , Nx_(0) , Ny_(0) , Nz_(0) , array_allocated_(false) , lp_(0)
35 , redshref_(-999.) , kredshref_(-999.)
36{
37 SetNThread();
38}
39
40GeneFluct3D::~GeneFluct3D(void)
41{
42 fftw_destroy_plan(pf_);
43 fftw_destroy_plan(pb_);
44#ifdef FFTW_THREAD
45 if(nthread_>0) fftw_cleanup_threads();
46#endif
47}
48
49//-------------------------------------------------------
50void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
51{
52 setsize(nx,ny,nz,dx,dy,dz);
53 setalloc();
54 setpointers(false);
55 init_fftw();
56}
57
58void GeneFluct3D::SetObservator(double redshref,double kredshref)
59// L'observateur est au redshift z=0
60// est situe sur la "perpendiculaire" a la face x,y
61// issue du centre de cette face
62// Il faut positionner le cube sur l'axe des z cad des redshifts:
63// redshref = redshift de reference
64// Si redshref<0 alors redshref=0
65// kredshref = indice (en double) correspondant a ce redshift
66// Si kredshref<0 alors kredshref=0
67{
68 if(redshref<0.) redshref = 0.;
69 if(kredshref<0.) kredshref = 0.;
70 redshref_ = redshref;
71 kredshref_ = kredshref;
72}
73
74void GeneFluct3D::setsize(long nx,long ny,long nz,double dx,double dy,double dz)
75{
76 if(nx<=0 || dx<=0.) {
77 cout<<"GeneFluct3D::SetSize_Error: bad value(s)"<<endl;
78 throw ParmError("GeneFluct3D::SetSize_Error: bad value(s)");
79 }
80
81 // Les tailles des tableaux
82 Nx_ = nx;
83 Ny_ = ny; if(Ny_ <= 0) Ny_ = Nx_;
84 Nz_ = nz; if(Nz_ <= 0) Nz_ = Nx_;
85 N_.resize(0); N_.push_back(Nx_); N_.push_back(Ny_); N_.push_back(Nz_);
86 NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
87 NCz_ = Nz_/2 +1;
88 NTz_ = 2*NCz_;
89
90 // le pas dans l'espace (Mpc)
91 Dx_ = dx;
92 Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
93 Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
94 D_.resize(0); D_.push_back(Dx_); D_.push_back(Dy_); D_.push_back(Dz_);
95 dVol_ = Dx_*Dy_*Dz_;
96 Vol_ = (Nx_*Dx_)*(Ny_*Dy_)*(Nz_*Dz_);
97
98 // Le pas dans l'espace de Fourier (Mpc^-1)
99 Dkx_ = 2.*M_PI/(Nx_*Dx_);
100 Dky_ = 2.*M_PI/(Ny_*Dy_);
101 Dkz_ = 2.*M_PI/(Nz_*Dz_);
102 Dk_.resize(0); Dk_.push_back(Dkx_); Dk_.push_back(Dky_); Dk_.push_back(Dkz_);
103 Dk3_ = Dkx_*Dky_*Dkz_;
104
105 // La frequence de Nyquist en k (Mpc^-1)
106 Knyqx_ = M_PI/Dx_;
107 Knyqy_ = M_PI/Dy_;
108 Knyqz_ = M_PI/Dz_;
109 Knyq_.resize(0); Knyq_.push_back(Knyqx_); Knyq_.push_back(Knyqy_); Knyq_.push_back(Knyqz_);
110}
111
112void GeneFluct3D::setalloc(void)
113{
114 // Dimensionnement du tableau complex<r_8>
115 // ATTENTION: TArray adresse en memoire a l'envers du C
116 // Tarray(n1,n2,n3) == Carray[n3][n2][n1]
117 sa_size_t SzK_[3] = {NCz_,Ny_,Nx_}; // a l'envers
118 try {
119 T_.ReSize(3,SzK_);
120 array_allocated_ = true;
121 } catch (...) {
122 cout<<"GeneFluct3D::SetSize_Error: Problem allocating T_"<<endl;
123 }
124 T_.SetMemoryMapping(BaseArray::CMemoryMapping);
125}
126
127void GeneFluct3D::setpointers(bool from_real)
128{
129 if(from_real) T_ = ArrCastR2C(R_);
130 else R_ = ArrCastC2R(T_);
131 // On remplit les pointeurs
132 fdata_ = (fftw_complex *) (&T_(0,0,0));
133 data_ = (double *) (&R_(0,0,0));
134}
135
136void GeneFluct3D::check_array_alloc(void)
137// Pour tester si le tableau T_ est alloue
138{
139 if(array_allocated_) return;
140 char bla[90];
141 sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
142 cout<<bla<<endl;
143 throw ParmError(bla);
144}
145
146void GeneFluct3D::init_fftw(void)
147{
148 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
149 if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: before fftw_plan ---");
150#ifdef FFTW_THREAD
151 if(nthread_>0) {
152 cout<<"Computing with "<<nthread_<<" threads"<<endl;
153 fftw_init_threads();
154 fftw_plan_with_nthreads(nthread_);
155 }
156#endif
157 pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
158 pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
159 if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: after fftw_plan ---");
160}
161
162
163//-------------------------------------------------------
164void GeneFluct3D::WriteFits(string cfname,int bitpix)
165{
166 cout<<"GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;
167 try {
168 FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
169 fwrt.WriteKey("NX",Nx_," axe transverse 1");
170 fwrt.WriteKey("NY",Ny_," axe transverse 2");
171 fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");
172 fwrt.WriteKey("DX",Dx_," Mpc");
173 fwrt.WriteKey("DY",Dy_," Mpc");
174 fwrt.WriteKey("DZ",Dz_," Mpc");
175 fwrt.WriteKey("DKX",Dkx_," Mpc^-1");
176 fwrt.WriteKey("DKY",Dky_," Mpc^-1");
177 fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
178 fwrt.WriteKey("ZREF",redshref_," reference redshift");
179 fwrt.WriteKey("KZREF",kredshref_," reference redshift on z axe");
180 fwrt.Write(R_);
181 } catch (PThrowable & exc) {
182 cout<<"Exception : "<<(string)typeid(exc).name()
183 <<" - Msg= "<<exc.Msg()<<endl;
184 return;
185 } catch (...) {
186 cout<<" some other exception was caught !"<<endl;
187 return;
188 }
189}
190
191void GeneFluct3D::ReadFits(string cfname)
192{
193 cout<<"GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;
194 try {
195 FitsImg3DRead fimg(cfname.c_str(),0,5);
196 fimg.Read(R_);
197 long nx = fimg.ReadKeyL("NX");
198 long ny = fimg.ReadKeyL("NY");
199 long nz = fimg.ReadKeyL("NZ");
200 double dx = fimg.ReadKey("DX");
201 double dy = fimg.ReadKey("DY");
202 double dz = fimg.ReadKey("DZ");
203 double zref = fimg.ReadKey("ZREF");
204 double kzref = fimg.ReadKey("KZREF");
205 setsize(nx,ny,nz,dx,dy,dz);
206 setpointers(true);
207 init_fftw();
208 SetObservator(zref,kzref);
209 } catch (PThrowable & exc) {
210 cout<<"Exception : "<<(string)typeid(exc).name()
211 <<" - Msg= "<<exc.Msg()<<endl;
212 return;
213 } catch (...) {
214 cout<<" some other exception was caught !"<<endl;
215 return;
216 }
217}
218
219void GeneFluct3D::WritePPF(string cfname,bool write_real)
220// On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> >
221{
222 cout<<"GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;
223 try {
224 R_.Info()["NX"] = (int_8)Nx_;
225 R_.Info()["NY"] = (int_8)Ny_;
226 R_.Info()["NZ"] = (int_8)Nz_;
227 R_.Info()["DX"] = (r_8)Dx_;
228 R_.Info()["DY"] = (r_8)Dy_;
229 R_.Info()["DZ"] = (r_8)Dz_;
230 R_.Info()["ZREF"] = (r_8)redshref_;
231 R_.Info()["KZREF"] = (r_8)kredshref_;
232 POutPersist pos(cfname.c_str());
233 if(write_real) pos << PPFNameTag("rgen") << R_;
234 else pos << PPFNameTag("pkgen") << T_;
235 } catch (PThrowable & exc) {
236 cout<<"Exception : "<<(string)typeid(exc).name()
237 <<" - Msg= "<<exc.Msg()<<endl;
238 return;
239 } catch (...) {
240 cout<<" some other exception was caught !"<<endl;
241 return;
242 }
243}
244
245void GeneFluct3D::ReadPPF(string cfname)
246{
247 cout<<"GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;
248 try {
249 bool from_real = true;
250 PInPersist pis(cfname.c_str());
251 string name_tag_k = "pkgen";
252 bool found_tag_k = pis.GotoNameTag("pkgen");
253 if(found_tag_k) {
254 cout<<" ...reading spectrun into TArray<complex <r_8> >"<<endl;
255 pis >> PPFNameTag("pkgen") >> T_;
256 from_real = false;
257 } else {
258 cout<<" ...reading space into TArray<r_8>"<<endl;
259 pis >> PPFNameTag("rgen") >> R_;
260 }
261 setpointers(from_real); // a mettre ici pour relire les DVInfo
262 int_8 nx = R_.Info()["NX"];
263 int_8 ny = R_.Info()["NY"];
264 int_8 nz = R_.Info()["NZ"];
265 r_8 dx = R_.Info()["DX"];
266 r_8 dy = R_.Info()["DY"];
267 r_8 dz = R_.Info()["DZ"];
268 r_8 zref = R_.Info()["ZREF"];
269 r_8 kzref = R_.Info()["KZREF"];
270 setsize(nx,ny,nz,dx,dy,dz);
271 init_fftw();
272 SetObservator(zref,kzref);
273 } catch (PThrowable & exc) {
274 cout<<"Exception : "<<(string)typeid(exc).name()
275 <<" - Msg= "<<exc.Msg()<<endl;
276 return;
277 } catch (...) {
278 cout<<" some other exception was caught !"<<endl;
279 return;
280 }
281}
282
283//-------------------------------------------------------
284void GeneFluct3D::Print(void)
285{
286 cout<<"GeneFluct3D(T_alloc="<<array_allocated_<<"):"<<endl;
287 cout<<"Space Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<Nz_<<" ("<<NTz_<<") size="
288 <<NRtot_<<endl;
289 cout<<" Resol: dx="<<Dx_<<" dy="<<Dy_<<" dz="<<Dz_<<" Mpc"
290 <<", dVol="<<dVol_<<", Vol="<<Vol_<<" Mpc^3"<<endl;
291 cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<endl;
292 cout<<" Resol: dkx="<<Dkx_<<" dky="<<Dky_<<" dkz="<<Dkz_<<" Mpc^-1"
293 <<", Dk3="<<Dk3_<<" Mpc^-3"<<endl;
294 cout<<" (2Pi/k: "<<2.*M_PI/Dkx_<<" "<<2.*M_PI/Dky_<<" "<<2.*M_PI/Dkz_<<" Mpc)"<<endl;
295 cout<<" Nyquist: kx="<<Knyqx_<<" ky="<<Knyqy_<<" kz="<<Knyqz_<<" Mpc^-1"
296 <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
297 cout<<" (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
298 cout<<"Redshift "<<redshref_<<" for z axe at k="<<kredshref_<<endl;
299}
300
301//-------------------------------------------------------
302void GeneFluct3D::ComputeFourier0(GenericFunc& pk_at_z)
303// cf ComputeFourier() mais avec autre methode de realisation du spectre
304// (attention on fait une fft pour realiser le spectre)
305{
306
307 // --- realisation d'un tableau de tirage gaussiens
308 if(lp_>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
309 // On tient compte du pb de normalisation de FFTW3
310 double sntot = sqrt((double)NRtot_);
311 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
312 int_8 ip = IndexR(i,j,l);
313 data_[ip] = NorRand()/sntot;
314 }
315 if(lp_>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
316
317 // --- realisation d'un tableau de tirage gaussiens
318 if(lp_>1) PrtTim("--- ComputeFourier0: before fft real ---");
319 fftw_execute(pf_);
320 if(lp_>1) PrtTim("--- ComputeFourier0: after fft real ---");
321
322 // --- On remplit avec une realisation
323 if(lp_>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
324 T_(0,0,0) = complex<r_8>(0.); // on coupe le continue et on l'initialise
325 long lmod = Nx_/10; if(lmod<1) lmod=1;
326 for(long i=0;i<Nx_;i++) {
327 long ii = (i>Nx_/2) ? Nx_-i : i;
328 double kx = ii*Dkx_; kx *= kx;
329 if(lp_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
330 for(long j=0;j<Ny_;j++) {
331 long jj = (j>Ny_/2) ? Ny_-j : j;
332 double ky = jj*Dky_; ky *= ky;
333 for(long l=0;l<NCz_;l++) {
334 double kz = l*Dkz_; kz *= kz;
335 if(i==0 && j==0 && l==0) continue; // Suppression du continu
336 double k = sqrt(kx+ky+kz);
337 // cf normalisation: Peacock, Cosmology, formule 16.38 p504
338 double pk = pk_at_z(k)/Vol_;
339 // ici pas de "/2" a cause de la remarque ci-dessus
340 T_(l,j,i) *= sqrt(pk);
341 }
342 }
343 }
344 if(lp_>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
345
346 double p = compute_power_carte();
347 cout<<"Puissance dans la realisation: "<<p<<endl;
348 if(lp_>1) PrtTim("--- ComputeFourier0: after Computing power ---");
349
350}
351
352//-------------------------------------------------------
353void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
354// Calcule une realisation du spectre "pk_at_z"
355// Attention: dans TArray le premier indice varie le + vite
356// Explication normalisation: see Coles & Lucchin, Cosmology, p264-265
357// FFTW3: on note N=Nx*Ny*Nz
358// f --(FFT)--> F = TF(f) --(FFT^-1)--> fb = TF^-1(F) = TF^-1(TF(f))
359// sum(f(x_i)^2) = S
360// sum(F(nu_i)^2) = S*N
361// sum(fb(x_i)^2) = S*N^2
362{
363 // --- RaZ du tableau
364 T_ = complex<r_8>(0.);
365
366 // --- On remplit avec une realisation
367 if(lp_>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
368 long lmod = Nx_/10; if(lmod<1) lmod=1;
369 for(long i=0;i<Nx_;i++) {
370 long ii = (i>Nx_/2) ? Nx_-i : i;
371 double kx = ii*Dkx_; kx *= kx;
372 if(lp_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
373 for(long j=0;j<Ny_;j++) {
374 long jj = (j>Ny_/2) ? Ny_-j : j;
375 double ky = jj*Dky_; ky *= ky;
376 for(long l=0;l<NCz_;l++) {
377 double kz = l*Dkz_; kz *= kz;
378 if(i==0 && j==0 && l==0) continue; // Suppression du continu
379 double k = sqrt(kx+ky+kz);
380 // cf normalisation: Peacock, Cosmology, formule 16.38 p504
381 double pk = pk_at_z(k)/Vol_;
382 // Explication de la division par 2: voir perandom.cc
383 // ou egalement Coles & Lucchin, Cosmology formula 13.7.2 p279
384 T_(l,j,i) = ComplexGaussRan(sqrt(pk/2.));
385 }
386 }
387 }
388 if(lp_>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");
389
390 manage_coefficients(); // gros effet pour les spectres que l'on utilise !
391 if(lp_>1) PrtTim("--- ComputeFourier: after managing coefficients ---");
392
393 double p = compute_power_carte();
394 cout<<"Puissance dans la realisation: "<<p<<endl;
395 if(lp_>1) PrtTim("--- ComputeFourier: after before Computing power ---");
396
397}
398
399long GeneFluct3D::manage_coefficients(void)
400// Take into account the real and complexe conjugate coefficients
401// because we want a realization of a real data in real space
402{
403 check_array_alloc();
404
405 // 1./ Le Continu et Nyquist sont reels
406 long nreal = 0;
407 for(long kk=0;kk<2;kk++) {
408 long k=0; // continu
409 if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
410 for(long jj=0;jj<2;jj++) {
411 long j=0;
412 if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
413 for(long ii=0;ii<2;ii++) {
414 long i=0;
415 if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
416 int_8 ip = IndexC(i,j,k);
417 //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata_[ip][0]<<","<<fdata_[ip][1]<<")"<<endl;
418 fdata_[ip][1] = 0.; fdata_[ip][0] *= M_SQRT2;
419 nreal++;
420 }
421 }
422 }
423 cout<<"Number of forced real number ="<<nreal<<endl;
424
425 // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist)
426
427 // a./ les lignes et colonnes du continu et de nyquist
428 long nconj1 = 0;
429 for(long kk=0;kk<2;kk++) {
430 long k=0; // continu
431 if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
432 for(long jj=0;jj<2;jj++) { // selon j
433 long j=0;
434 if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
435 for(long i=1;i<(Nx_+1)/2;i++) {
436 int_8 ip = IndexC(i,j,k);
437 int_8 ip1 = IndexC(Nx_-i,j,k);
438 fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
439 nconj1++;
440 }
441 }
442 for(long ii=0;ii<2;ii++) {
443 long i=0;
444 if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
445 for(long j=1;j<(Ny_+1)/2;j++) {
446 int_8 ip = IndexC(i,j,k);
447 int_8 ip1 = IndexC(i,Ny_-j,k);
448 fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
449 nconj1++;
450 }
451 }
452 }
453 cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;
454
455 // b./ les lignes et colonnes hors continu et de nyquist
456 long nconj2 = 0;
457 for(long kk=0;kk<2;kk++) {
458 long k=0; // continu
459 if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;} // Nyquist
460 for(long j=1;j<(Ny_+1)/2;j++) {
461 if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
462 for(long i=1;i<Nx_;i++) {
463 if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
464 int_8 ip = IndexC(i,j,k);
465 int_8 ip1 = IndexC(Nx_-i,Ny_-j,k);
466 fdata_[ip1][0] = fdata_[ip][0]; fdata_[ip1][1] = -fdata_[ip][1];
467 nconj2++;
468 }
469 }
470 }
471 cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
472
473 cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
474
475 return nreal+nconj1+nconj2;
476}
477
478double GeneFluct3D::compute_power_carte(void)
479// Calcul la puissance de la realisation du spectre Pk
480{
481 check_array_alloc();
482
483 double s2 = 0.;
484 for(long l=0;l<NCz_;l++)
485 for(long j=0;j<Ny_;j++)
486 for(long i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));
487
488 double s20 = 0.;
489 for(long j=0;j<Ny_;j++)
490 for(long i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));
491
492 double s2n = 0.;
493 if(Nz_%2==0)
494 for(long j=0;j<Ny_;j++)
495 for(long i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));
496
497 return 2.*s2 -s20 -s2n;
498}
499
500//-------------------------------------------------------------------
501void GeneFluct3D::FilterByPixel(void)
502// Filtrage par la fonction fenetre du pixel (parallelepipede)
503// TF = 1/(dx*dy*dz)*Int[{-dx/2,dx/2},{-dy/2,dy/2},{-dz/2,dz/2}]
504// e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
505// = 2/(k_x*dx) * sin(k_x*dx/2) * (idem y) * (idem z)
506// Gestion divergence en 0: sin(y)/y = 1 - y^2/6*(1-y^2/20)
507// avec y = k_x*dx/2
508{
509 check_array_alloc();
510
511 if(lp_>1) PrtTim("--- FilterByPixel: before ---");
512
513 for(long i=0;i<Nx_;i++) {
514 long ii = (i>Nx_/2) ? Nx_-i : i;
515 double kx = ii*Dkx_ *Dx_/2;
516 double pk_x = pixelfilter(kx);
517 for(long j=0;j<Ny_;j++) {
518 long jj = (j>Ny_/2) ? Ny_-j : j;
519 double ky = jj*Dky_ *Dy_/2;
520 double pk_y = pixelfilter(ky);
521 for(long l=0;l<NCz_;l++) {
522 double kz = l*Dkz_ *Dz_/2;
523 double pk_z = pixelfilter(kz);
524 T_(l,j,i) *= pk_x*pk_y*pk_z;
525 }
526 }
527 }
528
529 if(lp_>1) PrtTim("--- FilterByPixel: after ---");
530}
531
532//-------------------------------------------------------------------
533void GeneFluct3D::ComputeReal(void)
534// Calcule une realisation dans l'espace reel
535{
536 check_array_alloc();
537
538 // On fait la FFT
539 if(lp_>1) PrtTim("--- ComputeReal: before fftw backward ---");
540 fftw_execute(pb_);
541 if(lp_>1) PrtTim("--- ComputeReal: after fftw backward ---");
542}
543
544//-------------------------------------------------------------------
545void GeneFluct3D::ReComputeFourier(void)
546{
547 check_array_alloc();
548
549 // On fait la FFT
550 if(lp_>1) PrtTim("--- ComputeReal: before fftw forward ---");
551 fftw_execute(pf_);
552 // On corrige du pb de la normalisation de FFTW3
553 double v = (double)NRtot_;
554 for(long i=0;i<Nx_;i++)
555 for(long j=0;j<Ny_;j++)
556 for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
557
558 if(lp_>1) PrtTim("--- ComputeReal: after fftw forward ---");
559}
560
561//-------------------------------------------------------------------
562int GeneFluct3D::ComputeSpectrum(HistoErr& herr)
563// Compute spectrum from "T" and fill HistoErr "herr"
564// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
565// cad T(kz,ky,kx) avec 0<kz<kz_nyq -ky_nyq<ky<ky_nyq -kx_nyq<kx<kx_nyq
566{
567 check_array_alloc();
568 if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");
569
570 if(herr.NBins()<0) return -1;
571 herr.Zero();
572
573 // Attention a l'ordre
574 for(long i=0;i<Nx_;i++) {
575 long ii = (i>Nx_/2) ? Nx_-i : i;
576 double kx = ii*Dkx_; kx *= kx;
577 for(long j=0;j<Ny_;j++) {
578 long jj = (j>Ny_/2) ? Ny_-j : j;
579 double ky = jj*Dky_; ky *= ky;
580 for(long l=0;l<NCz_;l++) {
581 double kz = l*Dkz_; kz *= kz;
582 double k = sqrt(kx+ky+kz);
583 double pk = MODULE2(T_(l,j,i));
584 herr.Add(k,pk);
585 }
586 }
587 }
588 herr.ToVariance();
589
590 // renormalize to directly compare to original spectrum
591 double norm = Vol_;
592 herr *= norm;
593
594 if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");
595
596 return 0;
597}
598
599int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr)
600{
601 check_array_alloc();
602 if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");
603
604 if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
605 herr.Zero();
606
607 // Attention a l'ordre
608 for(long i=0;i<Nx_;i++) {
609 long ii = (i>Nx_/2) ? Nx_-i : i;
610 double kx = ii*Dkx_; kx *= kx;
611 for(long j=0;j<Ny_;j++) {
612 long jj = (j>Ny_/2) ? Ny_-j : j;
613 double ky = jj*Dky_; ky *= ky;
614 double kt = sqrt(kx+ky);
615 for(long l=0;l<NCz_;l++) {
616 double kz = l*Dkz_;
617 double pk = MODULE2(T_(l,j,i));
618 herr.Add(kt,kz,pk);
619 }
620 }
621 }
622 herr.ToVariance();
623
624 // renormalize to directly compare to original spectrum
625 double norm = Vol_;
626 herr *= norm;
627
628 if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");
629
630 return 0;
631}
632
633//-------------------------------------------------------
634int_8 GeneFluct3D::VarianceFrReal(double R,double& var)
635// Recompute MASS variance in spherical top-hat (rayon=R)
636{
637 check_array_alloc();
638
639 if(lp_>1) PrtTim("--- VarianceFrReal: before computation ---");
640
641 long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
642 long dny = long(R/Dy_+0.5); if(dny<=0) dny = 1;
643 long dnz = long(R/Dz_+0.5); if(dnz<=0) dnz = 1;
644 cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
645
646 double sum=0., sum2=0., r2 = R*R; int_8 nsum=0;
647
648 for(long i=dnx;i<Nx_-dnx;i+=dnx) {
649 for(long j=dny;j<Ny_-dny;j+=dny) {
650 for(long l=dnz;l<Nz_-dnz;l+=dnz) {
651 double s=0.; int_8 n=0;
652 for(long ii=i-dnx;ii<=i+dnx;ii++) {
653 double x = (ii-i)*Dx_; x *= x;
654 for(long jj=j-dny;jj<=j+dny;jj++) {
655 double y = (jj-j)*Dy_; y *= y;
656 for(long ll=l-dnz;ll<=l+dnz;ll++) {
657 double z = (ll-l)*Dz_; z *= z;
658 if(x+y+z>r2) continue;
659 int_8 ip = IndexR(ii,jj,ll);
660 s += 1.+data_[ip];
661 n++;
662 }
663 }
664 }
665 if(n>0) {sum += s; sum2 += s*s; nsum++;}
666 //cout<<i<<","<<j<<","<<l<<" n="<<n<<" s="<<s<<" sum="<<sum<<" sum2="<<sum2<<endl;
667 }
668 }
669 }
670
671 if(nsum<=1) {var=0.; return nsum;}
672
673 sum /= nsum;
674 sum2 = sum2/nsum - sum*sum;
675 if(lp_) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
676
677 var = sum2/(sum*sum); // <dM>^2/<M>^2
678 if(lp_) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
679
680 if(lp_>1) PrtTim("--- VarianceFrReal: after computation ---");
681 return nsum;
682}
683
684//-------------------------------------------------------
685int_8 GeneFluct3D::NumberOfBad(double vmin,double vmax)
686// number of pixels outside of ]vmin,vmax[ extremites exclues
687// -> vmin and vmax are considered as bad
688{
689 check_array_alloc();
690
691 int_8 nbad = 0;
692 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
693 int_8 ip = IndexR(i,j,l);
694 double v = data_[ip];
695 if(v<=vmin || v>=vmax) nbad++;
696 }
697
698 return nbad;
699}
700
701int_8 GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax)
702// mean,sigma^2 pour pixels avec valeurs ]vmin,vmax[ extremites exclues
703// -> mean and sigma^2 are NOT computed with pixels values vmin and vmax
704{
705 check_array_alloc();
706
707 int_8 n = 0;
708 rm = rs2 = 0.;
709
710 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
711 int_8 ip = IndexR(i,j,l);
712 double v = data_[ip];
713 if(v<=vmin || v>=vmax) continue;
714 rm += v;
715 rs2 += v*v;
716 n++;
717 }
718
719 if(n>1) {
720 rm /= (double)n;
721 rs2 = rs2/(double)n - rm*rm;
722 }
723
724 return n;
725}
726
727int_8 GeneFluct3D::SetToVal(double vmin, double vmax,double val0)
728// set to "val0" if out of range ]vmin,vmax[ extremites exclues
729// -> vmin and vmax are set to val0
730{
731 check_array_alloc();
732
733 int_8 nbad = 0;
734 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
735 int_8 ip = IndexR(i,j,l);
736 double v = data_[ip];
737 if(v<=vmin || v>=vmax) {data_[ip] = val0; nbad++;}
738 }
739
740 return nbad;
741}
742
743//-------------------------------------------------------
744void GeneFluct3D::TurnFluct2Mass(void)
745// d_rho/rho -> Mass (add one!)
746{
747 check_array_alloc();
748
749 if(lp_>1) PrtTim("--- TurnFluct2Mass: before computation ---");
750
751 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
752 int_8 ip = IndexR(i,j,l);
753 data_[ip] += 1.;
754 }
755
756 if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");
757}
758
759double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
760// do NOT treate negative or nul values
761{
762 if(lp_>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
763
764 double m,s2;
765 int_8 ngood = MeanSigma2(m,s2,0.,1e+200);
766 if(lp_) cout<<"TurnMass2MeanNumber: ngood="<<ngood
767 <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;
768
769 // On doit mettre n*Vol galaxies dans notre survey
770 // On en met uniquement dans les pixels de masse >0.
771 // On NE met PAS a zero les pixels <0
772 // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies
773 // comme on ne prend que les pixels >0, on doit normaliser
774 // a la moyenne de <1+d_rho/rho> sur ces pixels
775 // (rappel sur tout les pixels <1+d_rho/rho>=1)
776 double dn = n_by_mpc3*Vol_/m /(double)ngood; // nb de gal a mettre ds 1 px
777 if(lp_) cout<<"galaxy density move from "
778 <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
779 double sum = 0.;
780 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
781 int_8 ip = IndexR(i,j,l);
782 data_[ip] *= dn; // par coherence on multiplie aussi les <=0
783 if(data_[ip]>0.) sum += data_[ip]; // mais on ne les compte pas
784 }
785 if(lp_) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
786
787 if(lp_>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
788 return sum;
789}
790
791double GeneFluct3D::ApplyPoisson(void)
792// do NOT treate negative or nul mass -> let it as it is
793{
794 check_array_alloc();
795
796 if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
797
798 double sum = 0.;
799 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
800 int_8 ip = IndexR(i,j,l);
801 double v = data_[ip];
802 if(v>0.) {
803 unsigned long dn = PoissRandLimit(v,10.);
804 data_[ip] = (double)dn;
805 sum += (double)dn;
806 }
807 }
808 if(lp_) cout<<sum<<" galaxies put into survey"<<endl;
809
810 if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
811 return sum;
812}
813
814double GeneFluct3D::TurnNGal2Mass(FunRan& massdist,bool axeslog)
815// do NOT treate negative or nul mass -> let it as it is
816// INPUT:
817// massdist : distribution de masse (m*dn/dm)
818// axeslog = false : retourne la masse
819// = true : retourne le log10(mass)
820// RETURN la masse totale
821{
822 check_array_alloc();
823
824 if(lp_>1) PrtTim("--- TurnNGal2Mass: before computation ---");
825
826 double sum = 0.;
827 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
828 int_8 ip = IndexR(i,j,l);
829 double v = data_[ip];
830 if(v>0.) {
831 long ngal = long(v+0.1);
832 data_[ip] = 0.;
833 for(long i=0;i<ngal;i++) {
834 double m = massdist.RandomInterp(); // massdist.Random();
835 if(axeslog) m = pow(10.,m);
836 data_[ip] += m;
837 }
838 sum += data_[ip];
839 }
840 }
841 if(lp_) cout<<sum<<" MSol HI mass put into survey"<<endl;
842
843 if(lp_>1) PrtTim("--- TurnNGal2Mass: after computation ---");
844 return sum;
845}
846
847void GeneFluct3D::AddNoise2Real(double snoise)
848// add noise to every pixels (meme les <=0 !)
849{
850 check_array_alloc();
851
852 if(lp_>1) PrtTim("--- AddNoise2Real: before computation ---");
853
854 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
855 int_8 ip = IndexR(i,j,l);
856 data_[ip] += snoise*NorRand();
857 }
858
859 if(lp_>1) PrtTim("--- AddNoise2Real: after computation ---");
860}
Note: See TracBrowser for help on using the repository browser.