source: Sophya/trunk/Cosmo/SimLSS/cmvobserv3d.cc@ 3196

Last change on this file since 3196 was 3196, checked in by cmv, 18 years ago

les AGN selon C.Jackson, une premiere approche simplifiee, recodage from Jim Rich. cmv 03/04/2007
File size: 18.5 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#include "timing.h"
10#include "ntuple.h"
11#include "matharr.h"
12
13#include "constcosmo.h"
14#include "cosmocalc.h"
15#include "schechter.h"
16#include "geneutils.h"
17#include "genefluct3d.h"
18
19void usage(void);
20void usage(void)
21{
22 cout<<"cmvobserv3d [...options...]"<<endl
23 <<" -a : init auto de l aleatoire"<<endl
24 <<" -0 : use methode ComputeFourier0"<<endl
25 <<" -G : compute Pk(z=0) and apply growth factor in real space"<<endl
26 <<" -x nx,dx : taille en x (npix,Mpc)"<<endl
27 <<" -y ny,dy : taille en y (npix,Mpc)"<<endl
28 <<" -z nz,dz : taille en z (npix,Mpc)"<<endl
29 <<" -Z zref : redshift median"<<endl
30 <<" -s snoise : sigma du bruit en Msol"<<endl
31 <<" -2 : compute 2D spectrum"<<endl
32 <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
33 <<" -A <log10(S_agn en Msol)>,sigma : distribution des AGN par pixel"<<endl
34 <<" -W : write cube in FITS format"<<endl
35 <<" -P : write cube in PPF format"<<endl
36 <<" -V : compute variance from real space"<<endl
37 <<endl;
38}
39
40int main(int narg,char *arg[])
41{
42 InitTim();
43
44 //-----------------------------------------------------------------
45 // *** Survey definition
46 long nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
47 //long nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
48 //long nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
49
50 // *** Cosmography definition (WMAP)
51 unsigned short flat = 0;
52 double ob0 = 0.0444356;
53 double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
54 double zref = 0.5;
55 double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
56
57 // *** Spectrum and variance definition
58 double ns = 1., as = 1.;
59 double R=8./h100, Rg=R/sqrt(5.);
60 double sigmaR = 1.;
61
62 double kmin=1e-5,kmax=1000.;
63 int npt = 10000;
64 double lkmin=log10(kmin), lkmax=log10(kmax);
65 double eps=1.e-3;
66
67 // *** Schechter mass function definition
68 double h75 = h100 / 0.75;
69 double nstar = 0.006*pow(h75,3.);
70 double mstar = pow(10.,9.8/(h75*h75)); // MSol
71 double alpha = -1.37;
72
73 double schmin=1e8, schmax=1e12;
74 int schnpt = 1000;
75
76 // *** Niveau de bruit
77 double snoise= 0.; // en equivalent MSol
78
79 // *** AGN
80 bool do_agn = false;
81 double lmsol_agn=-99., lsigma_agn=0.; // en equivalent MSol
82
83 // *** type de generation
84 bool computefourier0=false;
85 bool use_growth_factor = false;
86 unsigned short nthread=4;
87
88 // *** What to do
89 bool comp2dspec = false;
90 bool wfits = false;
91 bool wppf = false;
92 bool compvarreal = false;
93
94 // --- Decodage arguments
95
96 char c;
97 while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:A:")) != -1) {
98 switch (c) {
99 case 'a' :
100 Auto_Ini_Ranf(5);
101 break;
102 case '0' :
103 computefourier0 = true;
104 break;
105 case 'G' :
106 use_growth_factor = true;
107 break;
108 case 'x' :
109 sscanf(optarg,"%ld,%lf",&nx,&dx);
110 break;
111 case 'y' :
112 sscanf(optarg,"%ld,%lf",&ny,&dy);
113 break;
114 case 'z' :
115 sscanf(optarg,"%ld,%lf",&nz,&dz);
116 break;
117 case 's' :
118 sscanf(optarg,"%lf",&snoise);
119 break;
120 case 'Z' :
121 sscanf(optarg,"%lf",&zref);
122 break;
123 case '2' :
124 comp2dspec = true;
125 break;
126 case 'M' :
127 sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
128 break;
129 case 'A' :
130 do_agn = true;
131 sscanf(optarg,"%lf,%lf",&lmsol_agn,&lsigma_agn);
132 break;
133 case 'V' :
134 compvarreal = true;
135 break;
136 case 'W' :
137 wfits = true;
138 break;
139 case 'P' :
140 wppf = true;
141 break;
142 case 'h' :
143 default :
144 usage(); return -1;
145 }
146 }
147
148 double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
149
150 string tagobs = "cmvobserv3d.ppf";
151 POutPersist posobs(tagobs);
152
153 cout<<"zref="<<zref<<endl;
154 cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
155 cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"
156 <<", npt="<<npt<<endl;
157 cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
158 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpa="<<alpha<<endl;
159 cout<<"schmin="<<schmin<<" ("<<lschmin
160 <<"), schmax="<<schmax<<" ("<<lschmax<<") Msol"
161 <<", schnpt="<<schnpt<<endl;
162 cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
163 if(do_agn) cout<<"AGN: <log10(Msol)>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
164
165 //-----------------------------------------------------------------
166 cout<<endl<<"\n--- Create Cosmology"<<endl;
167
168 CosmoCalc univ(flat,true,zref+1.);
169 univ.SetInteg(perc,dzinc,dzmax,glorder);
170 univ.SetDynParam(h100,om0,or0,ol0,w0);
171 univ.Print();
172 double loscomref = univ.Dloscom(zref);
173 cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
174 univ.Print(zref);
175
176 //-----------------------------------------------------------------
177 cout<<endl<<"\n--- Create Spectrum"<<endl;
178
179 InitialSpectrum pkini(ns,as);
180
181 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
182 //tf.SetNoOscEnv(2);
183
184 GrowthFactor growth(om0,ol0);
185 // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
186
187 PkSpectrum0 pk0(pkini,tf);
188
189 PkSpectrumZ pkz(pk0,growth,zref);
190
191 //-----------------------------------------------------------------
192 cout<<endl<<"\n--- Create mass function"<<endl;
193
194 Schechter sch(nstar,mstar,alpha);
195 sch.Print();
196
197 //-----------------------------------------------------------------
198 pkz.SetZ(0.);
199 cout<<endl<<"\n--- Compute variance for top-hat R="<<R
200 <<" at z="<<pkz.GetZ()<<endl;
201 VarianceSpectrum varpk_th(pkz,0);
202 double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
203 double pkmax_th = varpk_th(kfind_th);
204 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
205 double k1=kmin, k2=kmax;
206 int rc = varpk_th.FindLimits(R,pkmax_th/1.e4,k1,k2,eps);
207 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
208 <<k2<<" ("<<log10(k2)<<")"<<endl;
209
210 double ldlk = (log10(k2)-log10(k1))/npt;
211 varpk_th.SetInteg(0.01,ldlk,-1.,4);
212 double sr2 = varpk_th.Variance(R,k1,k2);
213 cout<<"varpk_th="<<sr2<<" -> sigma="<<sqrt(sr2)<<endl;
214
215 double normpkz = sigmaR*sigmaR/sr2;
216 pkz.SetScale(normpkz);
217 cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
218
219 pkz.SetZ(zref);
220
221 Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
222 FuncToHisto(pkz,hpkz,true);
223 {
224 tagobs = "hpkz"; posobs.PutObject(hpkz,tagobs);
225 }
226
227 //-----------------------------------------------------------------
228 cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
229 VarianceSpectrum varpk_int(pkz,2);
230
231 double kfind_int = varpk_int.FindMaximum(R,kmin,kmax,eps);
232 double pkmax_int = varpk_int(kfind_int);
233 cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
234 double k1int=kmin, k2int=kmax;
235 int rcint = varpk_int.FindLimits(R,pkmax_int/1.e4,k1int,k2int,eps);
236 cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
237 <<k2int<<" ("<<log10(k2int)<<")"<<endl;
238
239 double ldlkint = (log10(k2int)-log10(k1int))/npt;
240 varpk_int.SetInteg(0.01,ldlkint,-1.,4);
241 double sr2int = varpk_int.Variance(R,k1int,k2int);
242 cout<<"varpk_int="<<sr2int<<" -> sigma="<<sqrt(sr2int)<<endl;
243
244 //-----------------------------------------------------------------
245 cout<<endl<<"\n--- Compute galaxy density number"<<endl;
246
247 sch.SetOutValue(0);
248 cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
249 double ngal_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
250 cout<<"Galaxy density number = "<<ngal_by_mpc3<<" /Mpc^3 between limits"<<endl;
251
252 sch.SetOutValue(1);
253 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
254 double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
255 cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
256 cout<<"Omega_HI at z=0 is "<<mass_by_mpc3/(univ.Rhoc(0.)*GCm3toMsolMpc3_Cst)<<endl
257 <<" at z="<<zref<<" is "<<mass_by_mpc3/(univ.Rhoc(zref)*GCm3toMsolMpc3_Cst)<<endl;
258
259 PrtTim(">>>> End of definition");
260
261 //-----------------------------------------------------------------
262 // FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f with e+f=0 ou 1
263 cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
264
265 TArray< complex<r_8> > pkgen;
266 GeneFluct3D fluct3d(pkgen);
267 fluct3d.SetPrtLevel(2);
268 fluct3d.SetNThread(nthread);
269 fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
270 fluct3d.SetObservator(zref,nz/2.);
271 fluct3d.SetCosmology(univ);
272 fluct3d.SetGrowthFactor(growth);
273 fluct3d.LosComRedshift(0.001);
274 TArray<r_8>& rgen = fluct3d.GetRealArray();
275 cout<<endl; fluct3d.Print();
276
277 double dkmin = fluct3d.GetKincMin();
278 double knyqmax = fluct3d.GetKmax();
279 long nherr = long(knyqmax/dkmin+0.5);
280 cout<<"For HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
281
282 double dktmin = fluct3d.GetKTincMin();
283 double ktnyqmax = fluct3d.GetKTmax();
284 long nherrt = long(ktnyqmax/dktmin+0.5);
285 double dkzmin = fluct3d.GetKinc()[2];
286 double kznyqmax = fluct3d.GetKnyq()[2];
287 long nherrz = long(kznyqmax/dkzmin+0.5);
288 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
289 <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
290
291 //-----------------------------------------------------------------
292 cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
293 // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
294 // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
295 // Vc/Vs = 0.675 -> keff = kmax * (0.675)^(1/3) = kmax * 0.877
296 double knyqmax_mod = 0.877*knyqmax;
297 ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
298 varpk_int.SetInteg(0.01,ldlkint,-1.,4);
299 double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax_mod);
300 cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<" -> sigma="<<sqrt(sr2int_kmax)<<endl;
301
302 PrtTim(">>>> End Initialisation de GeneFluct3D");
303
304 //-----------------------------------------------------------------
305 cout<<"\n--- Computing a realization in Fourier space"<<endl;
306 if(use_growth_factor) pkz.SetZ(0.); else pkz.SetZ(zref);
307 cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
308 if(computefourier0) fluct3d.ComputeFourier0(pkz);
309 else fluct3d.ComputeFourier(pkz);
310 PrtTim(">>>> End Computing a realization in Fourier space");
311
312 if(1) {
313 cout<<"\n--- Checking realization spectra"<<endl;
314 HistoErr hpkgen(0.,knyqmax,nherr);
315 hpkgen.ReCenterBin(); hpkgen.Zero();
316 hpkgen.Show();
317 fluct3d.ComputeSpectrum(hpkgen);
318 {
319 tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
320 }
321 PrtTim(">>>> End Checking realization spectra");
322 }
323
324 if(comp2dspec) {
325 cout<<"\n--- Checking realization 2D spectra"<<endl;
326 Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
327 hpkgen2.ReCenterBin(); hpkgen2.Zero();
328 hpkgen2.Show();
329 fluct3d.ComputeSpectrum2D(hpkgen2);
330 {
331 tagobs = "hpkgen2"; posobs.PutObject(hpkgen2,tagobs);
332 }
333 PrtTim(">>>> End Checking realization 2D spectra");
334 }
335
336 if(1) {
337 cout<<"\n--- Computing convolution by pixel shape"<<endl;
338 fluct3d.FilterByPixel();
339 PrtTim(">>>> End Computing convolution by pixel shape");
340 }
341
342 if(wfits) {
343 fluct3d.WriteFits("!cmvobserv3d_k0.fits");
344 PrtTim(">>>> End WriteFits");
345 }
346 if(wppf) {
347 fluct3d.WritePPF("cmvobserv3d_k0.ppf",false);
348 PrtTim(">>>> End WritePPF");
349 }
350
351 if(1) {
352 cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
353 HistoErr hpkgenf(0.,knyqmax,nherr);
354 hpkgenf.ReCenterBin(); hpkgenf.Zero();
355 hpkgenf.Show();
356 fluct3d.ComputeSpectrum(hpkgenf);
357 {
358 tagobs = "hpkgenf"; posobs.PutObject(hpkgenf,tagobs);
359 }
360 PrtTim(">>>> End Checking realization spectra");
361 }
362
363 if(comp2dspec) {
364 cout<<"\n--- Checking realization 2D spectra after pixel shape convol."<<endl;
365 Histo2DErr hpkgenf2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
366 hpkgenf2.ReCenterBin(); hpkgenf2.Zero();
367 hpkgenf2.Show();
368 fluct3d.ComputeSpectrum2D(hpkgenf2);
369 {
370 tagobs = "hpkgenf2"; posobs.PutObject(hpkgenf2,tagobs);
371 }
372 PrtTim(">>>> End Checking realization 2D spectra");
373 }
374
375 //-----------------------------------------------------------------
376 cout<<"\n--- Computing a realization in real space"<<endl;
377 fluct3d.ComputeReal();
378 double rmin,rmax; rgen.MinMax(rmin,rmax);
379 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
380 PrtTim(">>>> End Computing a realization in real space");
381
382 if(use_growth_factor) {
383 cout<<"\n--- Apply Growth factor"<<endl;
384 cout<<"...D(z=0)="<<growth(0.)<<" D(z="<<zref<<")="<<growth(zref)<<endl;
385 fluct3d.ApplyGrowthFactor(-1);
386 rmin,rmax; rgen.MinMax(rmin,rmax);
387 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
388 PrtTim(">>>> End Applying growth factor");
389 }
390
391 if(wfits) {
392 fluct3d.WriteFits("!cmvobserv3d_r0.fits");
393 PrtTim(">>>> End WriteFits");
394 }
395 if(wppf) {
396 fluct3d.WritePPF("cmvobserv3d_r0.ppf",true);
397 PrtTim(">>>> End WritePPF");
398 }
399
400 int_8 nm;
401 double rm,rs2;
402 if(1) {
403 cout<<"\n--- Check mean and variance in real space"<<endl;
404 int_8 nlowone = fluct3d.NumberOfBad(-1.,1e+200);
405 nm = fluct3d.MeanSigma2(rm,rs2);
406 cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
407 <<rs2<<" -> "<<sqrt(rs2)<<", n(<-1)="<<nlowone<<endl;
408 PrtTim(">>>> End Check mean and variance in real space");
409 }
410
411 if(compvarreal) {
412 cout<<"\n--- Check variance sigmaR in real space"<<endl;
413 double varr;
414 int_8 nvarr = fluct3d.VarianceFrReal(R,varr);
415 cout<<"R="<<R<<" : sigmaR^2="<<varr<<" -> "<<sqrt(varr)<<", n="<<nvarr<<endl;
416 PrtTim(">>>> End Check variance sigmaR in real space");
417 }
418
419 //-----------------------------------------------------------------
420 cout<<endl<<"\n--- Converting fluctuations into mass"<<endl;
421 fluct3d.TurnFluct2Mass();
422 nm = fluct3d.MeanSigma2(rm,rs2);
423 cout<<"1+rgen: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
424 <<rs2<<" -> "<<sqrt(rs2)<<endl;
425 PrtTim(">>>> End Converting fluctuations into mass");
426
427 cout<<"\n--- Converting mass into galaxy number"<<endl;
428 rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
429 cout<<rm<<" galaxies put into survey"<<endl;
430 nm = fluct3d.MeanSigma2(rm,rs2,0.);
431 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
432 <<rs2<<" -> "<<sqrt(rs2)<<endl;
433 PrtTim(">>>> End Converting mass into galaxy number");
434
435 cout<<"\n--- Set negative pixels to BAD"<<endl;
436 nm = fluct3d.SetToVal(0.,1e+200,-999.);
437 cout<<nm<<" negative in survey set to BAD"<<endl;
438 nm = fluct3d.MeanSigma2(rm,rs2,-998.);
439 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
440 <<rs2<<" -> "<<sqrt(rs2)<<endl;
441 PrtTim(">>>> End Set negative pixels to BAD etc...");
442
443 cout<<"\n--- Apply poisson on galaxy number"<<endl;
444 nm = fluct3d.ApplyPoisson();
445 cout<<nm<<" galaxies into survey after poisson"<<endl;
446 nm = fluct3d.MeanSigma2(rm,rs2,-998.);
447 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
448 <<rs2<<" -> "<<sqrt(rs2)<<endl;
449 PrtTim(">>>> End Apply poisson on galaxy number");
450
451 cout<<"\n--- Convert Galaxy number to HI mass"<<endl;
452 long nhmdndm = long( (lschmax-lschmin+1.)*100. + 0.5);
453 Histo hmdndm(lschmin,lschmax,nhmdndm);
454 sch.SetOutValue(1);
455 FuncToHisto(sch,hmdndm,true);
456 FunRan tirhmdndm(hmdndm,true);
457 {
458 tagobs = "hmdndm"; posobs.PutObject(hmdndm,tagobs);
459 Histo hdum1(tirhmdndm);
460 tagobs = "tirhmdndm"; posobs.PutObject(hdum1,tagobs);
461 }
462 double mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
463 cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
464 nm = fluct3d.MeanSigma2(rm,rs2,0.);
465 cout<<"HI mass: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
466 <<rs2<<" -> "<<sqrt(rs2)<<endl;
467 cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
468 PrtTim(">>>> End Convert Galaxy number to HI mass");
469
470 cout<<"\n--- Set BAD pixels to Zero"<<endl;
471 nm = fluct3d.SetToVal(-998.,1e+200,0.);
472 cout<<nm<<" BAD in survey set to zero"<<endl;
473 nm = fluct3d.MeanSigma2(rm,rs2);
474 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
475 <<rs2<<" -> "<<sqrt(rs2)<<endl;
476 PrtTim(">>>> End Set BAD pixels to Zero etc...");
477
478 if(wfits) {
479 fluct3d.WriteFits("!cmvobserv3d_r.fits");
480 PrtTim(">>>> End WriteFits");
481 }
482 if(wppf) {
483 fluct3d.WritePPF("cmvobserv3d_r.ppf",true);
484 PrtTim(">>>> End WritePPF");
485 }
486
487 if(do_agn) {
488 cout<<"\n--- Add AGN: <mass>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
489 fluct3d.AddAGN(lmsol_agn,lsigma_agn);
490 nm = fluct3d.MeanSigma2(rm,rs2);
491 cout<<"HI mass with AGN: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
492 <<rs2<<" -> "<<sqrt(rs2)<<endl;
493 PrtTim(">>>> End Add AGN");
494 }
495
496 if(snoise>0.) {
497 cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
498 fluct3d.AddNoise2Real(snoise);
499 nm = fluct3d.MeanSigma2(rm,rs2);
500 cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
501 <<rs2<<" -> "<<sqrt(rs2)<<endl;
502 PrtTim(">>>> End Add noise");
503 }
504
505 //-----------------------------------------------------------------
506 // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
507
508 if(1) {
509 cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
510 fluct3d.ReComputeFourier();
511 PrtTim(">>>> End ReComputing spectrum");
512 }
513
514 if(wfits) {
515 fluct3d.WriteFits("!cmvobserv3d_k.fits");
516 PrtTim(">>>> End WriteFits");
517 }
518 if(wppf) {
519 fluct3d.WritePPF("cmvobserv3d_k.ppf",false);
520 PrtTim(">>>> End WritePPF");
521 }
522
523 if(1) {
524 cout<<endl<<"\n--- Computing final spectrum"<<endl;
525 HistoErr hpkrec(0.,knyqmax,nherr);
526 hpkrec.ReCenterBin();
527 hpkrec.Show();
528 fluct3d.ComputeSpectrum(hpkrec);
529 tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
530 PrtTim(">>>> End Computing final spectrum");
531 }
532
533 if(comp2dspec) {
534 cout<<"\n--- Computing final 2D spectrum"<<endl;
535 Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
536 hpkrec2.ReCenterBin(); hpkrec2.Zero();
537 hpkrec2.Show();
538 fluct3d.ComputeSpectrum2D(hpkrec2);
539 {
540 tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
541 }
542 PrtTim(">>>> End Computing final 2D spectrum");
543 }
544
545 PrtTim(">>>> End Of Job");
546 return 0;
547}
548
549/*
550######################################################
551readfits cmvobserv3d_k0.fits
552readfits cmvobserv3d_k.fits
553readfits cmvobserv3d_r0.fits
554readfits cmvobserv3d_r.fits
555
556openppf cmvobserv3d_k0.ppf
557openppf cmvobserv3d_k.ppf
558openppf cmvobserv3d_r0.ppf
559openppf cmvobserv3d_r.ppf
560
561# pour le plot 2D d'une slice en Z du 3D: to2d nom_obj3D num_slice
562defscript to2d
563 objaoper $1 sliceyz $2
564 mv sliceyz_${2} ${1}_$2
565 disp ${1}_$2
566 echo display slice $2 of $1
567endscript
568
569to2d $cobj 0
570
571######################################################
572openppf cmvobserv3d.ppf
573
574zone
575set k pow(10.,x)
576n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
577
578zone
579n/plot hpkz.val%x ! ! "nsta connectpoints"
580n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
581n/plot hpkgenf.val%log10(x) x>0 ! "nsta same orange connectpoints"
582n/plot hpkrec.val%log10(x) x>0 ! "nsta same blue connectpoints"
583
584disp hpkgen "hbincont err"
585disp hpkgenf "hbincont err"
586disp hpkrec "hbincont err"
587
588zone 2 2
589imag hpkgen2
590imag hpkgenf2
591imag hpkrec2
592
593zone 2 1
594disp hmdndm "nsta"
595disp tirhmdndm "nsta"
596addline 0 1 20 1 "red"
597
598 */
Note: See TracBrowser for help on using the repository browser.