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

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

Corrections papiers avec les commentaires du referee (1) - Reza 26/10/2011
File size: 20.2 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, Vector& noisevec)
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())||(noisevec.Size() != fourAmp.SizeZ()))
217 throw SzMismatchError("ComputeNoiseFourierAmp(): angscales/noisevec.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 double noisepow=noisevec(kz);
228 if (prtlev_>2)
229 cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc
230 << " NoisePow=" << noisepow << endl;
231 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
232 kyy = (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
233 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
234 kxx=(double)kx*dkx_;
235 rep = resp(kxx*angsc, kyy*angsc);
236 if (rep<1.e-19) rep=1.e-19;
237 else nmodeok++;
238 amp = sqrt(noisepow/rep/2.);
239 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));
240 }
241 }
242 }
243
244 if (prtlev_>1) {
245 cout << " Four3DPk::ComputeNoiseFourierAmp(...) Computing FFT along frequency ... \n"
246 << " NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
247 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
248 }
249 TVector< complex<TF> > veczin(fourAmp.SizeZ()), veczout(fourAmp.SizeZ());
250 FFTWServer ffts(true);
251 ffts.setNormalize(true);
252
253 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
254 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
255 // veczin=fourAmp(Range(kx), Range(ky), Range::all());
256 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) veczin(kz)=fourAmp(kx,ky,kz);
257 ffts.FFTBackward(veczin,veczout);
258 veczout /= (TF)sqrt((double)fourAmp.SizeZ());
259 // fourAmp(Range(kx), Range(ky), Range::all())=veczout;
260 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) fourAmp(kx,ky,kz)=veczout(kz);
261 }
262 }
263
264 // if (prtlev_>2) fourAmp.Show();
265 if (prtlev_>0)
266 cout << " Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df) done ..." << endl;
267}
268
269// Compute mass field from its Fourier Coefficient
270TArray<TF> Four3DPk::ComputeMassDens()
271{
272 TArray<TF> massdens;
273// Backward fourier transform of the fourierAmp array
274 FFTWServer ffts(true);
275 ffts.setNormalize(true);
276 ffts.FFTBackward(fourAmp, massdens, true);
277 // cout << " Four3DPk::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" << massDens.Size() << endl;
278 cout << " Four3DPk::ComputeMassDens() done NPix=" << massdens.Size() << endl;
279 return massdens;
280}
281
282// Compute power spectrum as a function of wave number k
283// cells with amp^2=re^2+im^2>s2cut are ignored
284// Output : power spectrum (profile histogram)
285HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax, bool fgmodcnt)
286{
287 // The second half of the array along Y (matrix rows) contain
288 // negative frequencies
289 // int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
290 // The profile histogram will contain the mean value of FFT amplitude
291 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
292 // if (nbin < 1) nbin = nbh/2;
293 if ((kmax<0.)||(kmax<kmin)) {
294 kmin=0.;
295 double maxx=fourAmp.SizeX()*dkx_;
296 double maxy=fourAmp.SizeY()*dky_/2;
297 double maxz=fourAmp.SizeZ()*dkz_/2;
298 kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
299 }
300 if (nbin<2) nbin=256;
301 hp_pk_p_ = new HProf(kmin, kmax, nbin);
302 hp_pk_p_->SetErrOpt(false);
303 if (fgmodcnt) {
304 hmcnt_p_ = new Histo(kmin, kmax, nbin);
305 hmcntok_p_ = new Histo(kmin, kmax, nbin);
306 }
307 s2cut_=s2cut;
308 ComputePkCumul();
309 return *hp_pk_p_;
310}
311
312// Compute power spectrum as a function of wave number k
313// Cumul dans hp - cells with amp^2=re^2+im^2>s2cut are ignored
314void Four3DPk::ComputePkCumul()
315{
316 uint_8 nmodeok=0;
317 // fourAmp represent 3-D fourier transform of a real input array.
318 // The second half of the array along Y and Z contain negative frequencies
319 double kxx, kyy, kzz;
320 // sa_size_t is large integer type
321 // We ignore 0th term in all frequency directions ...
322 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
323 kzz = (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
324 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
325 kyy = (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
326 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) { // ignore the 0th coefficient (constant term)
327 kxx=(double)kx*dkx_;
328 complex<TF> za = fourAmp(kx, ky, kz);
329 // if (za.real()>8.e19) continue;
330 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
331 double amp2 = za.real()*za.real()+za.imag()*za.imag();
332 if (hmcnt_p_) hmcnt_p_->Add(wk);
333 if ((s2cut_>1.e-9)&&(amp2>s2cut_)) continue;
334 if (hmcntok_p_) hmcntok_p_->Add(wk);
335 hp_pk_p_->Add(wk, amp2);
336 nmodeok++;
337 }
338 }
339 }
340 if ((prtlev_>1)||((prtlev_>0)&&(s2cut_>1.e-9))) {
341 cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
342 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
343 }
344 return;
345}
346
347// Compute noise power spectrum as a function of wave number k
348// No random generation, computes profile of noise sigma^2
349// cells with amp^2=re^2+im^2>s2cut are ignored
350// Output : noise power spectrum (profile histogram)
351// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
352// typically = ComovRadialDistance
353HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, double angscale, double s2cut, int nbin, double kmin, double kmax)
354{
355 // The second half of the array along Y (matrix rows) contain
356 // negative frequencies
357 // int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
358 // The profile histogram will contain the mean value of noise sigma^2
359 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
360 // if (nbin < 1) nbin = nbh/2;
361 double kmax2d=0.;
362 if ((kmax<0.)||(kmax<kmin)) {
363 kmin=0.;
364 double maxx=fourAmp.SizeX()*dkx_;
365 double maxy=fourAmp.SizeY()*dky_/2;
366 double maxz=fourAmp.SizeZ()*dkz_/2;
367 kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
368 kmax2d=sqrt(maxx*maxx+maxy*maxy);
369 }
370 if (nbin<2) nbin=256;
371 hp_pk_p_ = new HProf(kmin, kmax, nbin);
372 hp_pk_p_->SetErrOpt(false);
373 hmcnt_p_ = new Histo(kmin, kmax, nbin);
374 hmcntok_p_ = new Histo(kmin, kmax, nbin);
375 s2cut_=s2cut;
376
377 uint_8 nmodeok=0;
378 // fourAmp represent 3-D fourier transform of a real input array.
379 // The second half of the array along Y and Z contain negative frequencies
380 double kxx, kyy, kzz;
381 // sa_size_t is large integer type
382 // We ignore 0th term in all frequency directions ...
383 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
384 kzz = (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
385 uint_8 nmodeok_kz=0;
386 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
387 kyy = (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
388 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) { // ignore the 0th coefficient (constant term)
389 kxx=(double)kx*dkx_;
390 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
391 double rep=resp(kxx*angscale, kyy*angscale);
392 double amp2 = (rep>1.e-19)?1./rep:1.e19;
393 hmcnt_p_->Add(wk);
394 if ((s2cut_>1.e-9)&&(amp2>s2cut_)) continue;
395 hmcntok_p_->Add(wk);
396 hp_pk_p_->Add(wk, amp2);
397 nmodeok++; nmodeok_kz++;
398 }
399 }
400 if (prtlev_>2)
401 cout << " Four3DPk::ComputeNoisePk(kz="<<kz<<") ModeOK_Kz=" << nmodeok_kz
402 << " FracOK=" << 100.*(double)nmodeok_kz/(double)fourAmp.SizeX()/(double)fourAmp.SizeY() << endl;
403 }
404 if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
405 cout << " Four3DPk::ComputeNoisePk/Info : angscale=" << angscale
406 << " kmin,kmax=" << kmin << "," << kmax << " kmax2d_ang=" << kmax2d*angscale
407 << " \n ... NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
408 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
409 }
410
411 return *hp_pk_p_;
412}
413
414// Fills a data table from the computed P(k) profile histogram and mode count
415Histo Four3DPk::FillPkDataTable(DataTable& dt)
416{
417 if (hp_pk_p_==NULL) throw ParmError("Four3DPk::FillPkDataTable P(k) NOT computed");
418 if ((hmcnt_p_==NULL)||(hmcntok_p_==NULL)) throw ParmError("Four3DPk::FillPkDataTable Mode count histos NOT filled");
419 HProf& hp=(*hp_pk_p_);
420 Histo& hmcnt=(*hmcnt_p_);
421 Histo& hmcntok=(*hmcntok_p_);
422 Histo fracmodok=hmcntok/hmcnt;
423 char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
424 dt.Clear();
425 dt.AddDoubleColumn(nomcol[0]);
426 dt.AddDoubleColumn(nomcol[1]);
427 dt.AddIntegerColumn(nomcol[2]);
428 dt.AddIntegerColumn(nomcol[3]);
429 dt.AddFloatColumn(nomcol[4]);
430 DataTableRow dtr = dt.EmptyRow();
431 for(int_4 ib=0; ib<hp.NBins(); ib++) {
432 dtr[0]=hp.BinCenter(ib);
433 dtr[1]=hp(ib);
434 dtr[2]=hmcnt(ib);
435 dtr[3]=hmcntok(ib);
436 dtr[4]=fracmodok(ib);
437 dt.AddRow(dtr);
438 }
439 return fracmodok;
440}
441
442//-----------------------------------------------------
443// -- MassDist2D class : 2D mass distribution
444// --- PkNoiseCalculator : Classe de calcul du spectre de bruit PNoise(k)
445// determine par une reponse 2D de l'instrument
446//-----------------------------------------------------
447PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen,
448 const char* tit)
449 : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ()), pnoisefac_(pk3.SizeZ())
450{
451 SetFreqRange();
452 SetAngScaleConversion();
453 SetPNoiseFactor();
454 SetPrtLevel();
455}
456
457HProf PkNoiseCalculator::Compute(int nbin, double kmin, double kmax)
458{
459 Timer tm(title.c_str());
460 tm.Nop();
461 HProf hnd;
462 cout << "PkNoiseCalculator::Compute() " << title << " NGEN=" << NGEN << " S2CUT=" << S2CUT
463 << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscales=" << angscales_(0)
464 << " ... " << angscales_(angscales_.Size()-1) << endl;
465 for(int igen=0; igen<NGEN; igen++) {
466 pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_, pnoisefac_);
467 if (igen==0) hnd = pkn3d.ComputePk(S2CUT,nbin,kmin,kmax,true);
468 else pkn3d.ComputePkCumul();
469 if ((prtlev_>0)&&(igen>0)&&(((igen-1)%prtmodulo_)==0))
470 cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
471 }
472 return pkn3d.GetPk() ;
473}
474
475
476//-----------------------------------------------------
477// -- MassDist2D class : 2D mass distribution
478//-----------------------------------------------------
479// Constructor
480MassDist2D::MassDist2D(GenericFunc& pk, int size, double meandens)
481: pkSpec(pk) , sizeA((size>16)?size:16) , massDens(sizeA, sizeA),
482 meanRho(meandens) , fg_fourAmp(false) , fg_massDens(false)
483{
484}
485
486// To the computation job
487void MassDist2D::Compute()
488{
489 ComputeFourierAmp();
490 ComputeMassDens();
491}
492
493// Generate mass field Fourier Coefficient
494void MassDist2D::ComputeFourierAmp()
495{
496 if (fg_fourAmp) return; // job already done
497 // We generate a random gaussian real field
498 double sigma = 1.;
499// The following line fills the array by gaussian random numbers
500//--Replaced-- massDens = RandomSequence(RandomSequence::Gaussian, 0., sigma);
501// Can be replaced by
502 DR48RandGen rg;
503 for(sa_size_t ir=0; ir<massDens.NRows(); ir++) {
504 for(sa_size_t jc=0; jc<massDens.NCols(); jc++) {
505 massDens(ir, jc) = rg.Gaussian(sigma);
506 }
507 }
508// --- End of random filling
509
510 // Compute fourier transform of the random gaussian field -> white noise
511 FFTWServer ffts(true);
512 ffts.setNormalize(true);
513 ffts.FFTForward(massDens, fourAmp);
514
515 // fourAmp represent 2-D fourier transform of a real input array.
516 // The second half of the array along Y (matrix rows) contain
517 // negative frequencies
518// double fnorm = 1./sqrt(2.*fourAmp.Size());
519// PUT smaller value for fnorm and check number of zeros
520 double fnorm = 1.;
521 // sa_size_t is large integer type
522 for(sa_size_t ky=0; ky<fourAmp.NRows(); ky++) {
523 double kyy = ky;
524 if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky; // negative frequencies
525 for(sa_size_t kx=0; kx<fourAmp.NCols(); kx++) {
526 double wk = sqrt((double)(kx*kx+kyy*kyy));
527 double amp = pkSpec(wk)*fnorm;
528 fourAmp(ky, kx) *= amp; // renormalize fourier coeff using
529 }
530 }
531 fg_fourAmp = true;
532 cout << " MassDist2D::ComputeFourierAmp() done ..." << endl;
533}
534
535// Compute mass field from its Fourier Coefficient
536void MassDist2D::ComputeMassDens()
537{
538 if (fg_massDens) return; // job already done
539 if (!fg_fourAmp) ComputeFourierAmp(); // Check fourier amp generation
540
541// Backward fourier transform of the fourierAmp array
542 FFTWServer ffts(true);
543 ffts.setNormalize(true);
544 ffts.FFTBackward(fourAmp, massDens, true);
545// We consider that massDens represents delta rho/rho
546// rho = (delta rho/rho + 1) * MeanDensity
547 massDens += 1.;
548// We remove negative values
549 sa_size_t npbz = 0;
550 for (sa_size_t i=0; i<massDens.NRows(); i++)
551 for (sa_size_t j=0; j<massDens.NCols(); j++)
552 if (massDens(i,j) < 0.) { npbz++; massDens(i,j) = 0.; }
553 massDens *= meanRho;
554 cout << " MassDist2D::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" << massDens.Size() << endl;
555}
556
557// Compute power spectrum as a function of wave number k
558// Output : power spectrum (profile histogram)
559HProf MassDist2D::ReconstructPk(int nbin)
560{
561 // The second half of the array along Y (matrix rows) contain
562 // negative frequencies
563 int nbh = sqrt(2.0)*fourAmp.NCols();
564 // The profile histogram will contain the mean value of FFT amplitude
565 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
566 if (nbin < 1) nbin = nbh/2;
567 HProf hp(-0.5, nbh-0.5, nbin);
568 hp.SetErrOpt(false);
569
570 for(int ky=0; ky<fourAmp.NRows(); ky++) {
571 double kyy = ky;
572 if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky; // negative frequencies
573 for(int kx=0; kx<fourAmp.NCols(); kx++) {
574 double wk = sqrt((double)(kx*kx+kyy*kyy));
575 complex<r_8> za = fourAmp(ky, kx);
576 double amp = sqrt(za.real()*za.real()+za.imag()*za.imag());
577 hp.Add(wk, amp);
578 }
579 }
580 return hp;
581}
582
Note: See TracBrowser for help on using the repository browser.