source: Sophya/trunk/Cosmo/RadioBeam/specpk.cc@ 4028

Last change on this file since 4028 was 4028, checked in by ansari, 14 years ago

Debug de calcul Pnoise3D(k) avec D_A(z), Reza 18/10/2011
File size: 20.1 KB
Line 
1
2/* ------------------------ Projet BAORadio --------------------
3 Classes to compute 3D power spectrum and noise power spectrum
4 R. Ansari - Nov 2008 ... Dec 2010
5--------------------------------------------------------------- */
6
7#include "specpk.h"
8#include "radutil.h"
9#include "randr48.h"
10#include "ctimer.h"
11
12//------------------------------------
13// Class SpectralShape
14// -----------------------------------
15
16double Pnu1(double nu)
17{
18 return ( sqrt(sqrt(nu)) / ((nu+1.0)/0.2) *
19 (1+0.2*cos(2*M_PI*(nu-2.)*0.15)*exp(-nu/50.)) );
20}
21
22double Pnu2(double nu)
23{
24 if (nu < 1.e-9) return 0.;
25 return ((1.-exp(-nu/0.5))/nu*(1+0.25*cos(2*M_PI*nu*0.1)*exp(-nu/20.)) );
26}
27
28
29double Pnu3(double nu)
30{
31 return ( log(nu/100.+1)*(1+sin(2*M_PI*nu/300))*exp(-nu/4000) );
32}
33
34
35double Pnu4(double nu)
36{
37 double x = (nu-0.5)/0.05;
38 double rc = 2*exp(-x*x);
39 x = (nu-3.1)/0.27;
40 rc += exp(-x*x);
41 x = (nu-7.6)/1.4;
42 rc += 0.5*exp(-x*x);
43 return ( rc+2.*exp(-x*x) );
44}
45
46//--------------------------------------------------
47// -- SpectralShape class : test P(k) class
48//--------------------------------------------------
49// Constructor
50SpectralShape::SpectralShape(int typ)
51{
52 typ_=typ;
53 SetRenormFac();
54}
55
56// Return the spectral power for a given wave number wk
57double SpectralShape::operator() (double wk)
58{
59 wk/=DeuxPI;
60 double retv=1.;
61 switch (typ_)
62 {
63 case 1:
64 retv=Pnu1(wk);
65 break;
66 case 2:
67 retv=Pnu2(wk);
68 break;
69 case 3:
70 retv=Pnu3(wk);
71 break;
72 case 4:
73 retv=Pnu4(wk);
74 break;
75 default :
76 {
77 // global shape
78 double csp = pow( (2*sin(sqrt(sqrt(wk/7.)))),2.);
79 if (csp < 0.) return 0.;
80
81 // Adding some pics
82 double picpos[5] = {75.,150.,225.,300.,375.,};
83
84 for(int k=0; k<5; k++) {
85 double x0 = picpos[k];
86 if ( (wk > x0-25.) && (wk < x0+25.) ) {
87 double x = (wk-x0);
88 csp *= (1.+0.5*exp(-(x*x)/(2.*5*5)));
89 break;
90 }
91 }
92 retv=csp;
93 }
94 break;
95 }
96 return retv*renorm_fac;
97}
98// Return a vector representing the power spectrum (for checking)
99Histo SpectralShape::GetPk(int n)
100{
101 if (n<16) n = 256;
102 Histo h(0.,1024.*DeuxPI,n);
103 for(int k=0; k<h.NBins(); k++) h(k) = Value((k+0.5)*h.BinWidth());
104 return h;
105}
106
107double SpectralShape::Sommek2Pk(double kmax, int n)
108{
109 double dk=kmax/(double)n;
110 double s=0.;
111 for(int i=1; i<=n; i++) {
112 double ck=(double)i*dk;
113 s += Value(ck)*ck*ck;
114 }
115 return s*dk*4.*M_PI;
116}
117//--------------------------------------------------
118// -- Four2DResponse class : test P(k) class
119
120//---------------------------------------------------------------
121// -- Four3DPk class : 3D fourier amplitudes and power spectrum
122//---------------------------------------------------------------
123// Constructeur avec Tableau des coeff. de Fourier en argument
124Four3DPk::Four3DPk(TArray< complex<TF> > & fourcoedd, RandomGeneratorInterface& rg)
125 : rg_(rg), fourAmp(fourcoedd)
126{
127 SetPrtLevel();
128 SetCellSize();
129 hp_pk_p_=NULL; hmcnt_p_=NULL; hmcntok_p_=NULL; s2cut_=0.;
130}
131// Constructor
132Four3DPk::Four3DPk(RandomGeneratorInterface& rg, sa_size_t szx, sa_size_t szy, sa_size_t szz)
133 : rg_(rg), fourAmp(szx, szy, szz)
134{
135 SetPrtLevel();
136 SetCellSize();
137 hp_pk_p_=NULL; hmcnt_p_=NULL; hmcntok_p_=NULL; s2cut_=0.;
138}
139
140// Destructor
141Four3DPk::~Four3DPk()
142{
143 if (hp_pk_p_) delete hp_pk_p_;
144 if (hmcnt_p_) delete hmcnt_p_;
145 if (hmcntok_p_) delete hmcntok_p_;
146}
147
148// Generate mass field Fourier Coefficient
149void Four3DPk::ComputeFourierAmp(SpectralShape& pk)
150{
151 // We generate a random gaussian real field
152 // fourAmp represent 3-D fourier transform of a real input array.
153 // The second half of the array along Y and Z contain negative frequencies
154 // double fnorm = 1./sqrt(2.*fourAmp.Size());
155 double fnorm = 1.;
156 double kxx, kyy, kzz;
157 // sa_size_t is large integer type
158 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
159 kzz = (kz>fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
160 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
161 kyy = (ky>fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
162 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
163 double kxx=(double)kx*dkx_;
164 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
165 double amp = sqrt(pk(wk)*fnorm/2.);
166 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp)); // renormalize fourier coeff usin
167 }
168 }
169 }
170 if (prtlev_>2)
171 cout << " Four3DPk::ComputeFourierAmp() done ..." << endl;
172}
173
174
175// Generate mass field Fourier Coefficient
176void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale, bool crmask)
177// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
178// typically = ComovRadialDistance
179{
180 TMatrix<r_4> mask(fourAmp.SizeY(), fourAmp.SizeX());
181 // fourAmp represent 3-D fourier transform of a real input array.
182 // The second half of the array along Y and Z contain negative frequencies
183 double kxx, kyy, kzz, rep, amp;
184 // sa_size_t is large integer type
185 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
186 kzz = (kz>fourAmp.SizeZ()/2) ? -(double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
187 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
188 kyy = (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
189 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
190 kxx=(double)kx*dkx_;
191 rep = resp(kxx*angscale, kyy*angscale);
192 if (crmask&&(kz==0)) mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
193 if (rep<1.e-8) fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
194 else {
195 amp = 1./sqrt(rep)/sqrt(2.);
196 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));
197 }
198 }
199 }
200 }
201 if (prtlev_>2) fourAmp.Show();
202 if (crmask) {
203 POutPersist po("mask.ppf");
204 po << mask;
205 }
206 if (prtlev_>0)
207 cout << " Four3DPk::ComputeNoiseFourierAmp() done ..." << endl;
208}
209
210// Generate mass field Fourier Coefficient
211void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales)
212// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
213// typically = ComovRadialDistance
214{
215 uint_8 nmodeok=0;
216 if (angscales.Size() != fourAmp.SizeZ())
217 throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
218 H21Conversions conv;
219 // fourAmp represent 3-D fourier transform of a real input array.
220 // The second half of the array along Y and Z contain negative frequencies
221 double kxx, kyy, kzz, rep, amp;
222 // sa_size_t is large integer type
223 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
224 conv.setFrequency(f0+kz*df);
225 resp.setLambda(conv.getLambda());
226 double angsc=angscales(kz);
227 if (prtlev_>2)
228 cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc << endl;
229 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
230 kyy = (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
231 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
232 kxx=(double)kx*dkx_;
233 rep = resp(kxx*angsc, kyy*angsc);
234 if (rep<1.e-19) rep=1.e-19;
235 else nmodeok++;
236 amp = 1./sqrt(rep)/sqrt(2.);
237 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));
238 }
239 }
240 }
241
242 if (prtlev_>1) {
243 cout << " Four3DPk::ComputeNoiseFourierAmp(...) Computing FFT along frequency ... \n"
244 << " NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
245 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
246 }
247 TVector< complex<TF> > veczin(fourAmp.SizeZ()), veczout(fourAmp.SizeZ());
248 FFTWServer ffts(true);
249 ffts.setNormalize(true);
250
251 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
252 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
253 // veczin=fourAmp(Range(kx), Range(ky), Range::all());
254 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) veczin(kz)=fourAmp(kx,ky,kz);
255 ffts.FFTBackward(veczin,veczout);
256 veczout /= (TF)sqrt((double)fourAmp.SizeZ());
257 // fourAmp(Range(kx), Range(ky), Range::all())=veczout;
258 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) fourAmp(kx,ky,kz)=veczout(kz);
259 }
260 }
261
262 // if (prtlev_>2) fourAmp.Show();
263 if (prtlev_>0)
264 cout << " Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df) done ..." << endl;
265}
266
267// Compute mass field from its Fourier Coefficient
268TArray<TF> Four3DPk::ComputeMassDens()
269{
270 TArray<TF> massdens;
271// Backward fourier transform of the fourierAmp array
272 FFTWServer ffts(true);
273 ffts.setNormalize(true);
274 ffts.FFTBackward(fourAmp, massdens, true);
275 // cout << " Four3DPk::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" << massDens.Size() << endl;
276 cout << " Four3DPk::ComputeMassDens() done NPix=" << massdens.Size() << endl;
277 return massdens;
278}
279
280// Compute power spectrum as a function of wave number k
281// cells with amp^2=re^2+im^2>s2cut are ignored
282// Output : power spectrum (profile histogram)
283HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax, bool fgmodcnt)
284{
285 // The second half of the array along Y (matrix rows) contain
286 // negative frequencies
287 // int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
288 // The profile histogram will contain the mean value of FFT amplitude
289 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
290 // if (nbin < 1) nbin = nbh/2;
291 if ((kmax<0.)||(kmax<kmin)) {
292 kmin=0.;
293 double maxx=fourAmp.SizeX()*dkx_;
294 double maxy=fourAmp.SizeY()*dky_/2;
295 double maxz=fourAmp.SizeZ()*dkz_/2;
296 kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
297 }
298 if (nbin<2) nbin=256;
299 hp_pk_p_ = new HProf(kmin, kmax, nbin);
300 hp_pk_p_->SetErrOpt(false);
301 if (fgmodcnt) {
302 hmcnt_p_ = new Histo(kmin, kmax, nbin);
303 hmcntok_p_ = new Histo(kmin, kmax, nbin);
304 }
305 s2cut_=s2cut;
306 ComputePkCumul();
307 return *hp_pk_p_;
308}
309
310// Compute power spectrum as a function of wave number k
311// Cumul dans hp - cells with amp^2=re^2+im^2>s2cut are ignored
312void Four3DPk::ComputePkCumul()
313{
314 uint_8 nmodeok=0;
315 // fourAmp represent 3-D fourier transform of a real input array.
316 // The second half of the array along Y and Z contain negative frequencies
317 double kxx, kyy, kzz;
318 // sa_size_t is large integer type
319 // We ignore 0th term in all frequency directions ...
320 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
321 kzz = (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
322 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
323 kyy = (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
324 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) { // ignore the 0th coefficient (constant term)
325 kxx=(double)kx*dkx_;
326 complex<TF> za = fourAmp(kx, ky, kz);
327 // if (za.real()>8.e19) continue;
328 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
329 double amp2 = za.real()*za.real()+za.imag()*za.imag();
330 if (hmcnt_p_) hmcnt_p_->Add(wk);
331 if ((s2cut_>1.e-9)&&(amp2>s2cut_)) continue;
332 if (hmcntok_p_) hmcntok_p_->Add(wk);
333 hp_pk_p_->Add(wk, amp2);
334 nmodeok++;
335 }
336 }
337 }
338 if ((prtlev_>1)||((prtlev_>0)&&(s2cut_>1.e-9))) {
339 cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
340 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
341 }
342 return;
343}
344
345// Compute noise power spectrum as a function of wave number k
346// No random generation, computes profile of noise sigma^2
347// cells with amp^2=re^2+im^2>s2cut are ignored
348// Output : noise power spectrum (profile histogram)
349// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
350// typically = ComovRadialDistance
351HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, double angscale, double s2cut, int nbin, double kmin, double kmax)
352{
353 // The second half of the array along Y (matrix rows) contain
354 // negative frequencies
355 // int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
356 // The profile histogram will contain the mean value of noise sigma^2
357 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
358 // if (nbin < 1) nbin = nbh/2;
359 double kmax2d=0.;
360 if ((kmax<0.)||(kmax<kmin)) {
361 kmin=0.;
362 double maxx=fourAmp.SizeX()*dkx_;
363 double maxy=fourAmp.SizeY()*dky_/2;
364 double maxz=fourAmp.SizeZ()*dkz_/2;
365 kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
366 kmax2d=sqrt(maxx*maxx+maxy*maxy);
367 }
368 if (nbin<2) nbin=256;
369 hp_pk_p_ = new HProf(kmin, kmax, nbin);
370 hp_pk_p_->SetErrOpt(false);
371 hmcnt_p_ = new Histo(kmin, kmax, nbin);
372 hmcntok_p_ = new Histo(kmin, kmax, nbin);
373 s2cut_=s2cut;
374
375 uint_8 nmodeok=0;
376 // fourAmp represent 3-D fourier transform of a real input array.
377 // The second half of the array along Y and Z contain negative frequencies
378 double kxx, kyy, kzz;
379 // sa_size_t is large integer type
380 // We ignore 0th term in all frequency directions ...
381 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
382 kzz = (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
383 uint_8 nmodeok_kz=0;
384 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
385 kyy = (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
386 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) { // ignore the 0th coefficient (constant term)
387 kxx=(double)kx*dkx_;
388 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
389 double rep=resp(kxx*angscale, kyy*angscale);
390 double amp2 = (rep>1.e-19)?1./rep:1.e19;
391 hmcnt_p_->Add(wk);
392 if ((s2cut_>1.e-9)&&(amp2>s2cut_)) continue;
393 hmcntok_p_->Add(wk);
394 hp_pk_p_->Add(wk, amp2);
395 nmodeok++; nmodeok_kz++;
396 }
397 }
398 if (prtlev_>2)
399 cout << " Four3DPk::ComputeNoisePk(kz="<<kz<<") ModeOK_Kz=" << nmodeok_kz
400 << " FracOK=" << 100.*(double)nmodeok_kz/(double)fourAmp.SizeX()/(double)fourAmp.SizeY() << endl;
401 }
402 if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
403 cout << " Four3DPk::ComputeNoisePk/Info : angscale=" << angscale
404 << " kmin,kmax=" << kmin << "," << kmax << " kmax2d_ang=" << kmax2d*angscale
405 << " \n ... NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
406 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
407 }
408
409 return *hp_pk_p_;
410}
411
412// Fills a data table from the computed P(k) profile histogram and mode count
413Histo Four3DPk::FillPkDataTable(DataTable& dt)
414{
415 if (hp_pk_p_==NULL) throw ParmError("Four3DPk::FillPkDataTable P(k) NOT computed");
416 if ((hmcnt_p_==NULL)||(hmcntok_p_==NULL)) throw ParmError("Four3DPk::FillPkDataTable Mode count histos NOT filled");
417 HProf& hp=(*hp_pk_p_);
418 Histo& hmcnt=(*hmcnt_p_);
419 Histo& hmcntok=(*hmcntok_p_);
420 Histo fracmodok=hmcntok/hmcnt;
421 char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
422 dt.Clear();
423 dt.AddDoubleColumn(nomcol[0]);
424 dt.AddDoubleColumn(nomcol[1]);
425 dt.AddIntegerColumn(nomcol[2]);
426 dt.AddIntegerColumn(nomcol[3]);
427 dt.AddFloatColumn(nomcol[4]);
428 DataTableRow dtr = dt.EmptyRow();
429 for(int_4 ib=0; ib<hp.NBins(); ib++) {
430 dtr[0]=hp.BinCenter(ib);
431 dtr[1]=hp(ib);
432 dtr[2]=hmcnt(ib);
433 dtr[3]=hmcntok(ib);
434 dtr[4]=fracmodok(ib);
435 dt.AddRow(dtr);
436 }
437 return fracmodok;
438}
439
440//-----------------------------------------------------
441// -- MassDist2D class : 2D mass distribution
442// --- PkNoiseCalculator : Classe de calcul du spectre de bruit PNoise(k)
443// determine par une reponse 2D de l'instrument
444//-----------------------------------------------------
445PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen,
446 const char* tit)
447 : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ())
448{
449 SetFreqRange();
450 SetAngScaleConversion();
451 SetPrtLevel();
452}
453
454HProf PkNoiseCalculator::Compute(int nbin, double kmin, double kmax)
455{
456 Timer tm(title.c_str());
457 tm.Nop();
458 HProf hnd;
459 cout << "PkNoiseCalculator::Compute() " << title << " NGEN=" << NGEN << " S2CUT=" << S2CUT
460 << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscales=" << angscales_(0)
461 << " ... " << angscales_(angscales_.Size()-1) << endl;
462 for(int igen=0; igen<NGEN; igen++) {
463 pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_);
464 if (igen==0) hnd = pkn3d.ComputePk(S2CUT,nbin,kmin,kmax,true);
465 else pkn3d.ComputePkCumul();
466 if ((prtlev_>0)&&(igen>0)&&(((igen-1)%prtmodulo_)==0))
467 cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
468 }
469 return pkn3d.GetPk() ;
470}
471
472
473//-----------------------------------------------------
474// -- MassDist2D class : 2D mass distribution
475//-----------------------------------------------------
476// Constructor
477MassDist2D::MassDist2D(GenericFunc& pk, int size, double meandens)
478: pkSpec(pk) , sizeA((size>16)?size:16) , massDens(sizeA, sizeA),
479 meanRho(meandens) , fg_fourAmp(false) , fg_massDens(false)
480{
481}
482
483// To the computation job
484void MassDist2D::Compute()
485{
486 ComputeFourierAmp();
487 ComputeMassDens();
488}
489
490// Generate mass field Fourier Coefficient
491void MassDist2D::ComputeFourierAmp()
492{
493 if (fg_fourAmp) return; // job already done
494 // We generate a random gaussian real field
495 double sigma = 1.;
496// The following line fills the array by gaussian random numbers
497//--Replaced-- massDens = RandomSequence(RandomSequence::Gaussian, 0., sigma);
498// Can be replaced by
499 DR48RandGen rg;
500 for(sa_size_t ir=0; ir<massDens.NRows(); ir++) {
501 for(sa_size_t jc=0; jc<massDens.NCols(); jc++) {
502 massDens(ir, jc) = rg.Gaussian(sigma);
503 }
504 }
505// --- End of random filling
506
507 // Compute fourier transform of the random gaussian field -> white noise
508 FFTWServer ffts(true);
509 ffts.setNormalize(true);
510 ffts.FFTForward(massDens, fourAmp);
511
512 // fourAmp represent 2-D fourier transform of a real input array.
513 // The second half of the array along Y (matrix rows) contain
514 // negative frequencies
515// double fnorm = 1./sqrt(2.*fourAmp.Size());
516// PUT smaller value for fnorm and check number of zeros
517 double fnorm = 1.;
518 // sa_size_t is large integer type
519 for(sa_size_t ky=0; ky<fourAmp.NRows(); ky++) {
520 double kyy = ky;
521 if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky; // negative frequencies
522 for(sa_size_t kx=0; kx<fourAmp.NCols(); kx++) {
523 double wk = sqrt((double)(kx*kx+kyy*kyy));
524 double amp = pkSpec(wk)*fnorm;
525 fourAmp(ky, kx) *= amp; // renormalize fourier coeff using
526 }
527 }
528 fg_fourAmp = true;
529 cout << " MassDist2D::ComputeFourierAmp() done ..." << endl;
530}
531
532// Compute mass field from its Fourier Coefficient
533void MassDist2D::ComputeMassDens()
534{
535 if (fg_massDens) return; // job already done
536 if (!fg_fourAmp) ComputeFourierAmp(); // Check fourier amp generation
537
538// Backward fourier transform of the fourierAmp array
539 FFTWServer ffts(true);
540 ffts.setNormalize(true);
541 ffts.FFTBackward(fourAmp, massDens, true);
542// We consider that massDens represents delta rho/rho
543// rho = (delta rho/rho + 1) * MeanDensity
544 massDens += 1.;
545// We remove negative values
546 sa_size_t npbz = 0;
547 for (sa_size_t i=0; i<massDens.NRows(); i++)
548 for (sa_size_t j=0; j<massDens.NCols(); j++)
549 if (massDens(i,j) < 0.) { npbz++; massDens(i,j) = 0.; }
550 massDens *= meanRho;
551 cout << " MassDist2D::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" << massDens.Size() << endl;
552}
553
554// Compute power spectrum as a function of wave number k
555// Output : power spectrum (profile histogram)
556HProf MassDist2D::ReconstructPk(int nbin)
557{
558 // The second half of the array along Y (matrix rows) contain
559 // negative frequencies
560 int nbh = sqrt(2.0)*fourAmp.NCols();
561 // The profile histogram will contain the mean value of FFT amplitude
562 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
563 if (nbin < 1) nbin = nbh/2;
564 HProf hp(-0.5, nbh-0.5, nbin);
565 hp.SetErrOpt(false);
566
567 for(int ky=0; ky<fourAmp.NRows(); ky++) {
568 double kyy = ky;
569 if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky; // negative frequencies
570 for(int kx=0; kx<fourAmp.NCols(); kx++) {
571 double wk = sqrt((double)(kx*kx+kyy*kyy));
572 complex<r_8> za = fourAmp(ky, kx);
573 double amp = sqrt(za.real()*za.real()+za.imag()*za.imag());
574 hp.Add(wk, amp);
575 }
576 }
577 return hp;
578}
579
Note: See TracBrowser for help on using the repository browser.