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

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

Modif programme pknoise.cc et classes Four3DPk, PkNoiseCalculator pour calcul spectre de bruit en faisant varier la reponse instrumentale avec la frequence, Reza 10/10/2011
File size: 18.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}
130// Constructor
131Four3DPk::Four3DPk(RandomGeneratorInterface& rg, sa_size_t szx, sa_size_t szy, sa_size_t szz)
132 : rg_(rg), fourAmp(szx, szy, szz)
133{
134 SetPrtLevel();
135 SetCellSize();
136}
137
138
139// Generate mass field Fourier Coefficient
140void Four3DPk::ComputeFourierAmp(SpectralShape& pk)
141{
142 // We generate a random gaussian real field
143 // fourAmp represent 3-D fourier transform of a real input array.
144 // The second half of the array along Y and Z contain negative frequencies
145 // double fnorm = 1./sqrt(2.*fourAmp.Size());
146 double fnorm = 1.;
147 double kxx, kyy, kzz;
148 // sa_size_t is large integer type
149 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
150 kzz = (kz>fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
151 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
152 kyy = (ky>fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
153 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
154 double kxx=(double)kx*dkx_;
155 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
156 double amp = sqrt(pk(wk)*fnorm/2.);
157 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp)); // renormalize fourier coeff usin
158 }
159 }
160 }
161 if (prtlev_>2)
162 cout << " Four3DPk::ComputeFourierAmp() done ..." << endl;
163}
164
165
166// Generate mass field Fourier Coefficient
167void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale, bool crmask)
168// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
169// typically = ComovRadialDistance
170{
171 TMatrix<r_4> mask(fourAmp.SizeY(), fourAmp.SizeX());
172 // fourAmp represent 3-D fourier transform of a real input array.
173 // The second half of the array along Y and Z contain negative frequencies
174 double kxx, kyy, kzz, rep, amp;
175 // sa_size_t is large integer type
176 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
177 kzz = (kz>fourAmp.SizeZ()/2) ? -(double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
178 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
179 kyy = (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
180 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
181 kxx=(double)kx*dkx_;
182 rep = resp(kxx*angscale, kyy*angscale);
183 if (crmask&&(kz==0)) mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
184 if (rep<1.e-8) fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
185 else {
186 amp = 1./sqrt(rep)/sqrt(2.);
187 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));
188 }
189 }
190 }
191 }
192 if (prtlev_>2) fourAmp.Show();
193 if (crmask) {
194 POutPersist po("mask.ppf");
195 po << mask;
196 }
197 if (prtlev_>0)
198 cout << " Four3DPk::ComputeNoiseFourierAmp() done ..." << endl;
199}
200
201// Generate mass field Fourier Coefficient
202void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, double angscale)
203// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
204// typically = ComovRadialDistance
205{
206 H21Conversions conv;
207 // fourAmp represent 3-D fourier transform of a real input array.
208 // The second half of the array along Y and Z contain negative frequencies
209 double kxx, kyy, kzz, rep, amp;
210 // sa_size_t is large integer type
211 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
212 conv.setFrequency(f0+kz*df);
213 resp.setLambda(conv.getLambda());
214 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
215 kyy = (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
216 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
217 kxx=(double)kx*dkx_;
218 rep = resp(kxx*angscale, kyy*angscale);
219 if (rep<1.e-8) fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
220 else {
221 amp = 1./sqrt(rep)/sqrt(2.);
222 fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));
223 }
224 }
225 }
226 }
227
228 if (prtlev_>1)
229 cout << " Four3DPk::ComputeNoiseFourierAmp(...) Computing FFT along frequency ..." << endl;
230 TVector< complex<TF> > veczin(fourAmp.SizeZ()), veczout(fourAmp.SizeZ());
231 FFTWServer ffts(true);
232 ffts.setNormalize(true);
233
234 for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
235 for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
236 // veczin=fourAmp(Range(kx), Range(ky), Range::all());
237 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) veczin(kz)=fourAmp(kx,ky,kz);
238 ffts.FFTBackward(veczin,veczout);
239 veczout /= (TF)sqrt((double)fourAmp.SizeZ());
240 // fourAmp(Range(kx), Range(ky), Range::all())=veczout;
241 for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) fourAmp(kx,ky,kz)=veczout(kz);
242 }
243 }
244
245 if (prtlev_>2) fourAmp.Show();
246 if (prtlev_>0)
247 cout << " Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df) done ..." << endl;
248}
249
250// Compute mass field from its Fourier Coefficient
251TArray<TF> Four3DPk::ComputeMassDens()
252{
253 TArray<TF> massdens;
254// Backward fourier transform of the fourierAmp array
255 FFTWServer ffts(true);
256 ffts.setNormalize(true);
257 ffts.FFTBackward(fourAmp, massdens, true);
258 // cout << " Four3DPk::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" << massDens.Size() << endl;
259 cout << " Four3DPk::ComputeMassDens() done NPix=" << massdens.Size() << endl;
260 return massdens;
261}
262
263// Compute power spectrum as a function of wave number k
264// cells with amp^2=re^2+im^2>s2cut are ignored
265// Output : power spectrum (profile histogram)
266HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax)
267{
268 // The second half of the array along Y (matrix rows) contain
269 // negative frequencies
270 // int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
271 // The profile histogram will contain the mean value of FFT amplitude
272 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
273 // if (nbin < 1) nbin = nbh/2;
274 if ((kmax<0.)||(kmax<kmin)) {
275 kmin=0.;
276 double maxx=fourAmp.SizeX()*dkx_;
277 double maxy=fourAmp.SizeY()*dky_/2;
278 double maxz=fourAmp.SizeZ()*dkz_/2;
279 kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
280 }
281 if (nbin<2) nbin=128;
282 HProf hp(kmin, kmax, nbin);
283 hp.SetErrOpt(false);
284 ComputePkCumul(hp, s2cut);
285 return hp;
286}
287
288// Compute power spectrum as a function of wave number k
289// Cumul dans hp - cells with amp^2=re^2+im^2>s2cut are ignored
290void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
291{
292 uint_8 nmodeok=0;
293 // fourAmp represent 3-D fourier transform of a real input array.
294 // The second half of the array along Y and Z contain negative frequencies
295 double kxx, kyy, kzz;
296 // sa_size_t is large integer type
297 // We ignore 0th term in all frequency directions ...
298 for(sa_size_t kz=1; kz<fourAmp.SizeZ(); kz++) {
299 kzz = (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
300 for(sa_size_t ky=1; ky<fourAmp.SizeY(); ky++) {
301 kyy = (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
302 for(sa_size_t kx=1; kx<fourAmp.SizeX(); kx++) { // ignore the 0th coefficient (constant term)
303 double kxx=(double)kx*dkx_;
304 complex<TF> za = fourAmp(kx, ky, kz);
305 if (za.real()>8.e9) continue;
306 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
307 double amp2 = za.real()*za.real()+za.imag()*za.imag();
308 if ((s2cut>1.e-9)&&(amp2>s2cut)) continue;
309 hp.Add(wk, amp2);
310 nmodeok++;
311 }
312 }
313 }
314 if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
315 cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
316 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
317 }
318 return;
319}
320
321// Compute noise power spectrum as a function of wave number k
322// No random generation, computes profile of noise sigma^2
323// cells with amp^2=re^2+im^2>s2cut are ignored
324// Output : noise power spectrum (profile histogram)
325// angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers
326// typically = ComovRadialDistance
327HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt,
328 double angscale, double s2cut, int nbin, double kmin, double kmax)
329{
330 // The second half of the array along Y (matrix rows) contain
331 // negative frequencies
332 // int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
333 // The profile histogram will contain the mean value of noise sigma^2
334 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
335 // if (nbin < 1) nbin = nbh/2;
336 if ((kmax<0.)||(kmax<kmin)) {
337 kmin=0.;
338 double maxx=fourAmp.SizeX()*dkx_;
339 double maxy=fourAmp.SizeY()*dky_/2;
340 double maxz=fourAmp.SizeZ()*dkz_/2;
341 kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
342 }
343 if (nbin<2) nbin=128;
344 HProf hp(kmin, kmax, nbin);
345 hp.SetErrOpt(false);
346 Histo hmcnt(kmin, kmax, nbin);
347 Histo hmcntok(kmin, kmax, nbin);
348
349 uint_8 nmodeok=0;
350 // fourAmp represent 3-D fourier transform of a real input array.
351 // The second half of the array along Y and Z contain negative frequencies
352 double kxx, kyy, kzz;
353 // sa_size_t is large integer type
354 // We ignore 0th term in all frequency directions ...
355 for(sa_size_t kz=1; kz<fourAmp.SizeZ(); kz++) {
356 kzz = (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
357 for(sa_size_t ky=1; ky<fourAmp.SizeY(); ky++) {
358 kyy = (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
359 for(sa_size_t kx=1; kx<fourAmp.SizeX(); kx++) { // ignore the 0th coefficient (constant term)
360 double kxx=(double)kx*dkx_;
361 double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
362 double rep=resp(kxx*angscale, kyy*angscale);
363 double amp2 = (rep>1.e-19)?1./rep:1.e19;
364 hmcnt.Add(wk);
365 if ((s2cut>1.e-9)&&(amp2>s2cut)) continue;
366 hmcntok.Add(wk);
367 hp.Add(wk, amp2);
368 nmodeok++;
369 }
370 }
371 }
372 if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
373 cout << " Four3DPk::ComputeNoisePk/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
374 << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
375 }
376
377 fracmodok=hmcntok/hmcnt;
378
379 char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
380 dt.Clear();
381 dt.AddDoubleColumn(nomcol[0]);
382 dt.AddDoubleColumn(nomcol[1]);
383 dt.AddIntegerColumn(nomcol[2]);
384 dt.AddIntegerColumn(nomcol[3]);
385 dt.AddFloatColumn(nomcol[4]);
386 DataTableRow dtr = dt.EmptyRow();
387 for(int_4 ib=0; ib<hp.NBins(); ib++) {
388 dtr[0]=hp.BinCenter(ib);
389 dtr[1]=hp(ib);
390 dtr[2]=hmcnt(ib);
391 dtr[3]=hmcntok(ib);
392 dtr[4]=fracmodok(ib);
393 dt.AddRow(dtr);
394 }
395 return hp;
396}
397
398//-----------------------------------------------------
399// -- MassDist2D class : 2D mass distribution
400// --- PkNoiseCalculator : Classe de calcul du spectre de bruit PNoise(k)
401// determine par une reponse 2D de l'instrument
402//-----------------------------------------------------
403PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen,
404 const char* tit)
405 : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit)
406{
407 SetFreqRange();
408 SetAngScaleConversion();
409 SetPrtLevel();
410}
411
412HProf PkNoiseCalculator::Compute()
413{
414 Timer tm(title.c_str());
415 tm.Nop();
416 HProf hnd;
417 cout << "PkNoiseCalculator::Compute() " << title << " NGEN=" << NGEN << " S2CUT=" << S2CUT
418 << " Freq0=" << freq0_ << " dFreq=" << dfreq_ << " angscale=" << angscale_ << endl;
419 for(int igen=0; igen<NGEN; igen++) {
420 pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscale_);
421 if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
422 else pkn3d.ComputePkCumul(hnd,S2CUT);
423 if ((prtlev_>0)&&(igen>0)&&(((igen-1)%prtmodulo_)==0))
424 cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
425 }
426 return hnd;
427}
428
429
430//-----------------------------------------------------
431// -- MassDist2D class : 2D mass distribution
432//-----------------------------------------------------
433// Constructor
434MassDist2D::MassDist2D(GenericFunc& pk, int size, double meandens)
435: pkSpec(pk) , sizeA((size>16)?size:16) , massDens(sizeA, sizeA),
436 meanRho(meandens) , fg_fourAmp(false) , fg_massDens(false)
437{
438}
439
440// To the computation job
441void MassDist2D::Compute()
442{
443 ComputeFourierAmp();
444 ComputeMassDens();
445}
446
447// Generate mass field Fourier Coefficient
448void MassDist2D::ComputeFourierAmp()
449{
450 if (fg_fourAmp) return; // job already done
451 // We generate a random gaussian real field
452 double sigma = 1.;
453// The following line fills the array by gaussian random numbers
454//--Replaced-- massDens = RandomSequence(RandomSequence::Gaussian, 0., sigma);
455// Can be replaced by
456 DR48RandGen rg;
457 for(sa_size_t ir=0; ir<massDens.NRows(); ir++) {
458 for(sa_size_t jc=0; jc<massDens.NCols(); jc++) {
459 massDens(ir, jc) = rg.Gaussian(sigma);
460 }
461 }
462// --- End of random filling
463
464 // Compute fourier transform of the random gaussian field -> white noise
465 FFTWServer ffts(true);
466 ffts.setNormalize(true);
467 ffts.FFTForward(massDens, fourAmp);
468
469 // fourAmp represent 2-D fourier transform of a real input array.
470 // The second half of the array along Y (matrix rows) contain
471 // negative frequencies
472// double fnorm = 1./sqrt(2.*fourAmp.Size());
473// PUT smaller value for fnorm and check number of zeros
474 double fnorm = 1.;
475 // sa_size_t is large integer type
476 for(sa_size_t ky=0; ky<fourAmp.NRows(); ky++) {
477 double kyy = ky;
478 if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky; // negative frequencies
479 for(sa_size_t kx=0; kx<fourAmp.NCols(); kx++) {
480 double wk = sqrt((double)(kx*kx+kyy*kyy));
481 double amp = pkSpec(wk)*fnorm;
482 fourAmp(ky, kx) *= amp; // renormalize fourier coeff using
483 }
484 }
485 fg_fourAmp = true;
486 cout << " MassDist2D::ComputeFourierAmp() done ..." << endl;
487}
488
489// Compute mass field from its Fourier Coefficient
490void MassDist2D::ComputeMassDens()
491{
492 if (fg_massDens) return; // job already done
493 if (!fg_fourAmp) ComputeFourierAmp(); // Check fourier amp generation
494
495// Backward fourier transform of the fourierAmp array
496 FFTWServer ffts(true);
497 ffts.setNormalize(true);
498 ffts.FFTBackward(fourAmp, massDens, true);
499// We consider that massDens represents delta rho/rho
500// rho = (delta rho/rho + 1) * MeanDensity
501 massDens += 1.;
502// We remove negative values
503 sa_size_t npbz = 0;
504 for (sa_size_t i=0; i<massDens.NRows(); i++)
505 for (sa_size_t j=0; j<massDens.NCols(); j++)
506 if (massDens(i,j) < 0.) { npbz++; massDens(i,j) = 0.; }
507 massDens *= meanRho;
508 cout << " MassDist2D::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" << massDens.Size() << endl;
509}
510
511// Compute power spectrum as a function of wave number k
512// Output : power spectrum (profile histogram)
513HProf MassDist2D::ReconstructPk(int nbin)
514{
515 // The second half of the array along Y (matrix rows) contain
516 // negative frequencies
517 int nbh = sqrt(2.0)*fourAmp.NCols();
518 // The profile histogram will contain the mean value of FFT amplitude
519 // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
520 if (nbin < 1) nbin = nbh/2;
521 HProf hp(-0.5, nbh-0.5, nbin);
522 hp.SetErrOpt(false);
523
524 for(int ky=0; ky<fourAmp.NRows(); ky++) {
525 double kyy = ky;
526 if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky; // negative frequencies
527 for(int kx=0; kx<fourAmp.NCols(); kx++) {
528 double wk = sqrt((double)(kx*kx+kyy*kyy));
529 complex<r_8> za = fourAmp(ky, kx);
530 double amp = sqrt(za.real()*za.real()+za.imag()*za.imag());
531 hp.Add(wk, amp);
532 }
533 }
534 return hp;
535}
536
Note: See TracBrowser for help on using the repository browser.