Changeset 3155 in Sophya for trunk/Cosmo
- Timestamp:
- Jan 19, 2007, 6:18:38 PM (19 years ago)
- Location:
- trunk/Cosmo/SimLSS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cosmo/SimLSS/cmvobserv3d.cc
r3154 r3155 135 135 <<", npt="<<npt<<endl; 136 136 cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl; 137 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alp ha="<<alpha<<endl;137 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpa="<<alpha<<endl; 138 138 cout<<"schmin="<<schmin<<" ("<<lschmin 139 139 <<"), schmax="<<schmax<<" ("<<lschmax<<") Msol" … … 204 204 cout<<"varpk_int="<<sr2int<<" -> sigma="<<sqrt(sr2int)<<endl; 205 205 206 PrtTim("End of definition");207 208 206 //----------------------------------------------------------------- 209 207 cout<<endl<<"\n--- Compute galaxy density number"<<endl; … … 219 217 cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl; 220 218 219 PrtTim(">>>> End of definition"); 220 221 221 //----------------------------------------------------------------- 222 222 // FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f with e+f=0 ou 1 223 cout<<endl<<"\n--- Compute Spectrum Realization"<<endl;223 cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl; 224 224 225 225 pkz.SetZ(zref); 226 226 TArray< complex<r_8> > pkgen; 227 227 GeneFluct3D fluct3d(pkgen); 228 fluct3d.SetPrtLevel( 3);228 fluct3d.SetPrtLevel(2); 229 229 fluct3d.SetNThread(nthread); 230 230 fluct3d.SetSize(nx,ny,nz,dx,dy,dz); … … 256 256 cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<" -> sigma="<<sqrt(sr2int_kmax)<<endl; 257 257 258 PrtTim(">>>> End Initialisation de GeneFluct3D"); 259 258 260 cout<<"\n--- Computing a realization in Fourier space"<<endl; 259 261 if(computefourier0) fluct3d.ComputeFourier0(pkz); 260 262 else fluct3d.ComputeFourier(pkz); 263 PrtTim(">>>> End Computing a realization in Fourier space"); 261 264 262 265 if(1) { … … 269 272 tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs); 270 273 } 274 PrtTim(">>>> End Checking realization spectra"); 271 275 } 272 276 … … 280 284 tagobs = "hpkgen2"; posobs.PutObject(hpkgen2,tagobs); 281 285 } 286 PrtTim(">>>> End Checking realization 2D spectra"); 282 287 } 283 288 … … 285 290 cout<<"\n--- Computing convolution by pixel shape"<<endl; 286 291 fluct3d.FilterByPixel(); 287 } 288 289 if(wfits) fluct3d.WriteFits("!cmvobserv3d_k0.fits"); 290 if(wppf) fluct3d.WritePPF("cmvobserv3d_k0.ppf",false); 292 PrtTim(">>>> End Computing convolution by pixel shape"); 293 } 294 295 if(wfits) { 296 fluct3d.WriteFits("!cmvobserv3d_k0.fits"); 297 PrtTim(">>>> End WriteFits"); 298 } 299 if(wppf) { 300 fluct3d.WritePPF("cmvobserv3d_k0.ppf",false); 301 PrtTim(">>>> End WritePPF"); 302 } 291 303 292 304 if(1) { … … 299 311 tagobs = "hpkgenf"; posobs.PutObject(hpkgenf,tagobs); 300 312 } 313 PrtTim(">>>> End Checking realization spectra"); 301 314 } 302 315 … … 310 323 tagobs = "hpkgenf2"; posobs.PutObject(hpkgenf2,tagobs); 311 324 } 325 PrtTim(">>>> End Checking realization 2D spectra"); 312 326 } 313 327 … … 316 330 double rmin,rmax; rgen.MinMax(rmin,rmax); 317 331 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl; 318 319 if(wfits) fluct3d.WriteFits("!cmvobserv3d_r0.fits"); 320 if(wppf) fluct3d.WritePPF("cmvobserv3d_r0.ppf",true); 332 PrtTim(">>>> End Computing a realization in real space"); 333 334 if(wfits) { 335 fluct3d.WriteFits("!cmvobserv3d_r0.fits"); 336 PrtTim(">>>> End WriteFits"); 337 } 338 if(wppf) { 339 fluct3d.WritePPF("cmvobserv3d_r0.ppf",true); 340 PrtTim(">>>> End WritePPF"); 341 } 321 342 322 343 int_8 nm; … … 328 349 cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 329 350 <<rs2<<" -> "<<sqrt(rs2)<<", n(<-1)="<<nlowone<<endl; 351 PrtTim(">>>> End Check mean and variance in real space"); 330 352 } 331 353 … … 335 357 int_8 nvarr = fluct3d.VarianceFrReal(R,varr); 336 358 cout<<"R="<<R<<" : sigmaR^2="<<varr<<" -> "<<sqrt(varr)<<", n="<<nvarr<<endl; 359 PrtTim(">>>> End Check variance sigmaR in real space"); 337 360 } 338 361 … … 343 366 cout<<"1+rgen: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 344 367 <<rs2<<" -> "<<sqrt(rs2)<<endl; 368 PrtTim(">>>> End Converting fluctuations into mass"); 345 369 346 370 cout<<"\n--- Converting mass into galaxy number"<<endl; … … 350 374 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 351 375 <<rs2<<" -> "<<sqrt(rs2)<<endl; 376 PrtTim(">>>> End Converting mass into galaxy number"); 352 377 353 378 cout<<"\n--- Set negative pixels to BAD"<<endl; … … 357 382 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 358 383 <<rs2<<" -> "<<sqrt(rs2)<<endl; 384 PrtTim(">>>> End Set negative pixels to BAD etc..."); 359 385 360 386 cout<<"\n--- Apply poisson on galaxy number"<<endl; … … 364 390 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 365 391 <<rs2<<" -> "<<sqrt(rs2)<<endl; 392 PrtTim(">>>> End Apply poisson on galaxy number"); 366 393 367 394 cout<<"\n--- Convert Galaxy number to HI mass"<<endl; … … 382 409 <<rs2<<" -> "<<sqrt(rs2)<<endl; 383 410 cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl; 411 PrtTim(">>>> End Convert Galaxy number to HI mass"); 384 412 385 413 cout<<"\n--- Set BAD pixels to Zero"<<endl; … … 389 417 cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 390 418 <<rs2<<" -> "<<sqrt(rs2)<<endl; 391 392 if(wfits) fluct3d.WriteFits("!cmvobserv3d_r.fits"); 393 if(wppf) fluct3d.WritePPF("cmvobserv3d_r.ppf",true); 419 PrtTim(">>>> End Set BAD pixels to Zero etc..."); 420 421 if(wfits) { 422 fluct3d.WriteFits("!cmvobserv3d_r.fits"); 423 PrtTim(">>>> End WriteFits"); 424 } 425 if(wppf) { 426 fluct3d.WritePPF("cmvobserv3d_r.ppf",true); 427 PrtTim(">>>> End WritePPF"); 428 } 394 429 395 430 cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl; … … 398 433 cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = " 399 434 <<rs2<<" -> "<<sqrt(rs2)<<endl; 435 PrtTim(">>>> End Add noise"); 400 436 401 437 //----------------------------------------------------------------- … … 405 441 cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl; 406 442 fluct3d.ReComputeFourier(); 407 } 408 409 if(wfits) fluct3d.WriteFits("!cmvobserv3d_k.fits"); 410 if(wppf) fluct3d.WritePPF("cmvobserv3d_k.ppf",false); 443 PrtTim(">>>> End ReComputing spectrum"); 444 } 445 446 if(wfits) { 447 fluct3d.WriteFits("!cmvobserv3d_k.fits"); 448 PrtTim(">>>> End WriteFits"); 449 } 450 if(wppf) { 451 fluct3d.WritePPF("cmvobserv3d_k.ppf",false); 452 PrtTim(">>>> End WritePPF"); 453 } 411 454 412 455 if(1) { … … 417 460 fluct3d.ComputeSpectrum(hpkrec); 418 461 tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs); 462 PrtTim(">>>> End Computing final spectrum"); 419 463 } 420 464 … … 428 472 tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs); 429 473 } 430 } 431 474 PrtTim(">>>> End Computing final 2D spectrum"); 475 } 476 477 PrtTim(">>>> End Of Job"); 432 478 return 0; 433 479 } -
trunk/Cosmo/SimLSS/genefluct3d.cc
r3154 r3155 8 8 #include <unistd.h> 9 9 10 #include "timing.h"11 10 #include "tarray.h" 12 11 #include "pexceptions.h" … … 74 73 void GeneFluct3D::setsize(long nx,long ny,long nz,double dx,double dy,double dz) 75 74 { 75 if(lp_>1) cout<<"--- GeneFluct3D::setsize: N="<<nx<<","<<ny<<","<<nz 76 <<" D="<<dx<<","<<dy<<","<<dz<<endl; 76 77 if(nx<=0 || dx<=0.) { 77 cout<<"GeneFluct3D:: SetSize_Error: bad value(s)"<<endl;78 throw ParmError("GeneFluct3D:: SetSize_Error: bad value(s)");78 cout<<"GeneFluct3D::setsize_Error: bad value(s)"<<endl; 79 throw ParmError("GeneFluct3D::setsize_Error: bad value(s)"); 79 80 } 80 81 … … 112 113 void GeneFluct3D::setalloc(void) 113 114 { 115 if(lp_>1) cout<<"--- GeneFluct3D::setalloc ---"<<endl; 114 116 // Dimensionnement du tableau complex<r_8> 115 117 // ATTENTION: TArray adresse en memoire a l'envers du C … … 120 122 array_allocated_ = true; 121 123 } catch (...) { 122 cout<<"GeneFluct3D:: SetSize_Error: Problem allocating T_"<<endl;124 cout<<"GeneFluct3D::setalloc_Error: Problem allocating T_"<<endl; 123 125 } 124 126 T_.SetMemoryMapping(BaseArray::CMemoryMapping); … … 127 129 void GeneFluct3D::setpointers(bool from_real) 128 130 { 131 if(lp_>1) cout<<"--- GeneFluct3D::setpointers ---"<<endl; 129 132 if(from_real) T_ = ArrCastR2C(R_); 130 133 else R_ = ArrCastC2R(T_); … … 147 150 { 148 151 // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init) 149 if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: before fftw_plan ---");152 if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl; 150 153 #ifdef FFTW_THREAD 151 154 if(nthread_>0) { 152 cout<<" Computing with "<<nthread_<<" threads"<<endl;155 cout<<"...Computing with "<<nthread_<<" threads"<<endl; 153 156 fftw_init_threads(); 154 157 fftw_plan_with_nthreads(nthread_); 155 158 } 156 159 #endif 160 if(lp_>1) cout<<"...forward plan"<<endl; 157 161 pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE); 162 if(lp_>1) cout<<"...backward plan"<<endl; 158 163 pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE); 159 if(lp_>1) PrtTim("--- GeneFluct3D::init_fftw: after fftw_plan ---");160 164 } 161 165 … … 164 168 void GeneFluct3D::WriteFits(string cfname,int bitpix) 165 169 { 166 cout<<" GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl;170 cout<<"--- GeneFluct3D::WriteFits: Writing Cube to "<<cfname<<endl; 167 171 try { 168 172 FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5); … … 191 195 void GeneFluct3D::ReadFits(string cfname) 192 196 { 193 cout<<" GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl;197 cout<<"--- GeneFluct3D::ReadFits: Reading Cube from "<<cfname<<endl; 194 198 try { 195 199 FitsImg3DRead fimg(cfname.c_str(),0,5); … … 220 224 // On ecrit soit le TArray<r_8> ou le TArray<complex <r_8> > 221 225 { 222 cout<<" GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl;226 cout<<"--- GeneFluct3D::WritePPF: Writing Cube (real="<<write_real<<") to "<<cfname<<endl; 223 227 try { 224 228 R_.Info()["NX"] = (int_8)Nx_; … … 245 249 void GeneFluct3D::ReadPPF(string cfname) 246 250 { 247 cout<<" GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl;251 cout<<"--- GeneFluct3D::ReadPPF: Reading Cube from "<<cfname<<endl; 248 252 try { 249 253 bool from_real = true; … … 306 310 307 311 // --- realisation d'un tableau de tirage gaussiens 308 if(lp_> 1) PrtTim("--- ComputeFourier0: before gaussian filling ---");312 if(lp_>0) cout<<"--- ComputeFourier0: before gaussian filling ---"<<endl; 309 313 // On tient compte du pb de normalisation de FFTW3 310 314 double sntot = sqrt((double)NRtot_); … … 313 317 data_[ip] = NorRand()/sntot; 314 318 } 315 if(lp_>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");316 319 317 320 // --- realisation d'un tableau de tirage gaussiens 318 if(lp_> 1) PrtTim("--- ComputeFourier0: before fft real ---");321 if(lp_>0) cout<<"...before fft real ---"<<endl; 319 322 fftw_execute(pf_); 320 if(lp_>1) PrtTim("--- ComputeFourier0: after fft real ---");321 323 322 324 // --- On remplit avec une realisation 323 if(lp_> 1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");325 if(lp_>0) cout<<"...before Fourier realization filling ---"<<endl; 324 326 T_(0,0,0) = complex<r_8>(0.); // on coupe le continue et on l'initialise 325 327 long lmod = Nx_/10; if(lmod<1) lmod=1; … … 327 329 long ii = (i>Nx_/2) ? Nx_-i : i; 328 330 double kx = ii*Dkx_; kx *= kx; 329 if(lp_> 1&& i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;331 if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl; 330 332 for(long j=0;j<Ny_;j++) { 331 333 long jj = (j>Ny_/2) ? Ny_-j : j; … … 342 344 } 343 345 } 344 if(lp_>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---"); 345 346 347 if(lp_>0) cout<<"...computing power"<<endl; 346 348 double p = compute_power_carte(); 347 cout<<"Puissance dans la realisation: "<<p<<endl; 348 if(lp_>1) PrtTim("--- ComputeFourier0: after Computing power ---"); 349 if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl; 349 350 350 351 } … … 365 366 366 367 // --- On remplit avec une realisation 367 if(lp_> 1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");368 if(lp_>0) cout<<"--- ComputeFourier ---"<<endl; 368 369 long lmod = Nx_/10; if(lmod<1) lmod=1; 369 370 for(long i=0;i<Nx_;i++) { 370 371 long ii = (i>Nx_/2) ? Nx_-i : i; 371 372 double kx = ii*Dkx_; kx *= kx; 372 if(lp_> 1&& i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;373 if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl; 373 374 for(long j=0;j<Ny_;j++) { 374 375 long jj = (j>Ny_/2) ? Ny_-j : j; … … 386 387 } 387 388 } 388 if(lp_>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");389 389 390 390 manage_coefficients(); // gros effet pour les spectres que l'on utilise ! 391 if(lp_>1) PrtTim("--- ComputeFourier: after managing coefficients ---"); 392 391 392 if(lp_>0) cout<<"...computing power"<<endl; 393 393 double p = compute_power_carte(); 394 cout<<"Puissance dans la realisation: "<<p<<endl; 395 if(lp_>1) PrtTim("--- ComputeFourier: after before Computing power ---"); 394 if(lp_>0) cout<<"Puissance dans la realisation: "<<p<<endl; 396 395 397 396 } … … 401 400 // because we want a realization of a real data in real space 402 401 { 402 if(lp_>1) cout<<"...managing coefficients"<<endl; 403 403 check_array_alloc(); 404 404 … … 421 421 } 422 422 } 423 cout<<"Number of forced real number ="<<nreal<<endl;423 if(lp_>1) cout<<"Number of forced real number ="<<nreal<<endl; 424 424 425 425 // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist) … … 451 451 } 452 452 } 453 cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;453 if(lp_>1) cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl; 454 454 455 455 // b./ les lignes et colonnes hors continu et de nyquist … … 469 469 } 470 470 } 471 cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;472 473 cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;471 if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl; 472 473 if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl; 474 474 475 475 return nreal+nconj1+nconj2; … … 507 507 // avec y = k_x*dx/2 508 508 { 509 check_array_alloc(); 510 511 if(lp_>1) PrtTim("--- FilterByPixel: before ---"); 509 if(lp_>0) cout<<"--- FilterByPixel ---"<<endl; 510 check_array_alloc(); 512 511 513 512 for(long i=0;i<Nx_;i++) { … … 527 526 } 528 527 529 if(lp_>1) PrtTim("--- FilterByPixel: after ---");530 528 } 531 529 … … 534 532 // Calcule une realisation dans l'espace reel 535 533 { 534 if(lp_>0) cout<<"--- ComputeReal ---"<<endl; 536 535 check_array_alloc(); 537 536 538 537 // On fait la FFT 539 if(lp_>1) PrtTim("--- ComputeReal: before fftw backward ---");540 538 fftw_execute(pb_); 541 if(lp_>1) PrtTim("--- ComputeReal: after fftw backward ---");542 539 } 543 540 … … 545 542 void GeneFluct3D::ReComputeFourier(void) 546 543 { 544 if(lp_>0) cout<<"--- ReComputeFourier ---"<<endl; 547 545 check_array_alloc(); 548 546 549 547 // On fait la FFT 550 if(lp_>1) PrtTim("--- ComputeReal: before fftw forward ---");551 548 fftw_execute(pf_); 552 549 // On corrige du pb de la normalisation de FFTW3 … … 556 553 for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v); 557 554 558 if(lp_>1) PrtTim("--- ComputeReal: after fftw forward ---");559 555 } 560 556 … … 565 561 // cad T(kz,ky,kx) avec 0<kz<kz_nyq -ky_nyq<ky<ky_nyq -kx_nyq<kx<kx_nyq 566 562 { 567 check_array_alloc();568 if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");563 if(lp_>0) cout<<"--- ComputeSpectrum ---"<<endl; 564 check_array_alloc(); 569 565 570 566 if(herr.NBins()<0) return -1; … … 592 588 herr *= norm; 593 589 594 if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");595 596 590 return 0; 597 591 } … … 599 593 int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr) 600 594 { 601 check_array_alloc();602 if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");595 if(lp_>0) cout<<"--- ComputeSpectrum2D ---"<<endl; 596 check_array_alloc(); 603 597 604 598 if(herr.NBinX()<0 || herr.NBinY()<0) return -1; … … 626 620 herr *= norm; 627 621 628 if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");629 630 622 return 0; 631 623 } … … 635 627 // Recompute MASS variance in spherical top-hat (rayon=R) 636 628 { 637 check_array_alloc(); 638 639 if(lp_>1) PrtTim("--- VarianceFrReal: before computation ---"); 629 if(lp_>0) cout<<"--- VarianceFrReal ---"<<endl; 630 check_array_alloc(); 640 631 641 632 long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1; 642 633 long dny = long(R/Dy_+0.5); if(dny<=0) dny = 1; 643 634 long dnz = long(R/Dz_+0.5); if(dnz<=0) dnz = 1; 644 cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;635 if(lp_>0) cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl; 645 636 646 637 double sum=0., sum2=0., r2 = R*R; int_8 nsum=0; … … 673 664 sum /= nsum; 674 665 sum2 = sum2/nsum - sum*sum; 675 if(lp_ ) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;666 if(lp_>0) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl; 676 667 677 668 var = sum2/(sum*sum); // <dM>^2/<M>^2 678 if(lp_) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl; 679 680 if(lp_>1) PrtTim("--- VarianceFrReal: after computation ---"); 669 if(lp_>0) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl; 670 681 671 return nsum; 682 672 } … … 745 735 // d_rho/rho -> Mass (add one!) 746 736 { 747 check_array_alloc();748 749 if(lp_>1) PrtTim("--- TurnFluct2Mass: before computation ---"); 737 if(lp_>0) cout<<"--- TurnFluct2Mass ---"<<endl; 738 check_array_alloc(); 739 750 740 751 741 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) { … … 753 743 data_[ip] += 1.; 754 744 } 755 756 if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");757 745 } 758 746 … … 760 748 // do NOT treate negative or nul values 761 749 { 762 if(lp_> 1) PrtTim("--- TurnMass2MeanNumber: before computation ---");750 if(lp_>0) cout<<"--- TurnMass2MeanNumber ---"<<endl; 763 751 764 752 double m,s2; 765 753 int_8 ngood = MeanSigma2(m,s2,0.,1e+200); 766 if(lp_ ) cout<<"TurnMass2MeanNumber:ngood="<<ngood767 <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;754 if(lp_>0) cout<<"...ngood="<<ngood 755 <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl; 768 756 769 757 // On doit mettre n*Vol galaxies dans notre survey … … 775 763 // (rappel sur tout les pixels <1+d_rho/rho>=1) 776 764 double dn = n_by_mpc3*Vol_/m /(double)ngood; // nb de gal a mettre ds 1 px 777 if(lp_ ) cout<<"galaxy density move from "778 <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;765 if(lp_>0) cout<<"...galaxy density move from " 766 <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl; 779 767 double sum = 0.; 780 768 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) { … … 783 771 if(data_[ip]>0.) sum += data_[ip]; // mais on ne les compte pas 784 772 } 785 if(lp_) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl; 786 787 if(lp_>1) PrtTim("--- TurnMass2MeanNumber: after computation ---"); 773 if(lp_>0) cout<<sum<<"...galaxies put into survey / "<<n_by_mpc3*Vol_<<endl; 774 788 775 return sum; 789 776 } … … 792 779 // do NOT treate negative or nul mass -> let it as it is 793 780 { 794 check_array_alloc(); 795 796 if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---"); 781 if(lp_>0) cout<<"--- ApplyPoisson ---"<<endl; 782 check_array_alloc(); 797 783 798 784 double sum = 0.; … … 806 792 } 807 793 } 808 if(lp_) cout<<sum<<" galaxies put into survey"<<endl; 809 810 if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---"); 794 if(lp_>0) cout<<sum<<" galaxies put into survey"<<endl; 795 811 796 return sum; 812 797 } … … 820 805 // RETURN la masse totale 821 806 { 822 check_array_alloc(); 823 824 if(lp_>1) PrtTim("--- TurnNGal2Mass: before computation ---"); 807 if(lp_>0) cout<<"--- TurnNGal2Mass ---"<<endl; 808 check_array_alloc(); 825 809 826 810 double sum = 0.; … … 839 823 } 840 824 } 841 if(lp_) cout<<sum<<" MSol HI mass put into survey"<<endl; 842 843 if(lp_>1) PrtTim("--- TurnNGal2Mass: after computation ---"); 825 if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl; 826 844 827 return sum; 845 828 } … … 848 831 // add noise to every pixels (meme les <=0 !) 849 832 { 850 check_array_alloc(); 851 852 if(lp_>1) PrtTim("--- AddNoise2Real: before computation ---"); 833 if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<endl; 834 check_array_alloc(); 853 835 854 836 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) { … … 856 838 data_[ip] += snoise*NorRand(); 857 839 } 858 859 if(lp_>1) PrtTim("--- AddNoise2Real: after computation ---"); 860 } 840 }
Note:
See TracChangeset
for help on using the changeset viewer.