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

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

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