[3115] | 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 | #include "srandgen.h"
|
---|
| 13 | #include "perandom.h"
|
---|
[3120] | 14 | #include "fabtwriter.h"
|
---|
[3115] | 15 |
|
---|
| 16 | #include "arrctcast.h"
|
---|
| 17 |
|
---|
| 18 | #include "constcosmo.h"
|
---|
| 19 | #include "schechter.h"
|
---|
| 20 | #include "geneutils.h"
|
---|
| 21 | #include "integfunc.h"
|
---|
| 22 | #include "genefluct3d.h"
|
---|
| 23 |
|
---|
| 24 | void usage(void);
|
---|
| 25 | void usage(void)
|
---|
| 26 | {
|
---|
| 27 | cout<<"cmvobserv3d [-a] [-0]"<<endl
|
---|
| 28 | <<" -a : init auto de l aleatoire"<<endl
|
---|
| 29 | <<" -0 : use methode ComputeFourier0"<<endl
|
---|
| 30 | <<" -x nx,dx : taille en x (npix,Mpc)"<<endl
|
---|
| 31 | <<" -y ny,dy : taille en y (npix,Mpc)"<<endl
|
---|
| 32 | <<" -z nz,dz : taille en z (npix,Mpc)"<<endl
|
---|
| 33 | <<" -Z zref : redshift median"<<endl
|
---|
[3120] | 34 | <<" -s snoise : sigma du bruit en Msol"<<endl
|
---|
[3115] | 35 | <<endl;
|
---|
| 36 | }
|
---|
| 37 |
|
---|
| 38 | int main(int narg,char *arg[])
|
---|
| 39 | {
|
---|
| 40 | InitTim();
|
---|
| 41 |
|
---|
| 42 | //-----------------------------------------------------------------
|
---|
| 43 | // *** Survey definition
|
---|
[3129] | 44 | long nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
|
---|
| 45 | //long nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
|
---|
| 46 | //long nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
|
---|
[3115] | 47 |
|
---|
| 48 | // *** Cosmography definition (WMAP)
|
---|
| 49 | double ob0 = 0.0444356;
|
---|
| 50 | double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
|
---|
| 51 | double zref = 0.5;
|
---|
| 52 |
|
---|
| 53 | // *** Spectrum and variance definition
|
---|
| 54 | double ns = 1., as = 1.;
|
---|
| 55 | double R=8./h100, Rg=R/sqrt(5.);
|
---|
| 56 | double sigmaR = 1.;
|
---|
| 57 |
|
---|
| 58 | double kmin=1e-5,kmax=1000.;
|
---|
| 59 | int npt = 10000;
|
---|
| 60 | double lkmin=log10(kmin), lkmax=log10(kmax);
|
---|
| 61 | double eps=1.e-3;
|
---|
| 62 |
|
---|
| 63 | // *** Schechter mass function definition
|
---|
| 64 | double h75 = 0.71 / 0.75;
|
---|
| 65 | double nstar = 0.006*pow(h75,3.);
|
---|
| 66 | double mstar = pow(10.,9.8/(h75*h75)); // MSol
|
---|
| 67 | double alpha = -1.37;
|
---|
| 68 |
|
---|
| 69 | double schmin=1e8, schmax=1e12;
|
---|
| 70 | int schnpt = 1000;
|
---|
| 71 | double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
|
---|
| 72 |
|
---|
| 73 | // *** Niveau de bruit
|
---|
| 74 | double snoise= 0.; // en equivalent MSol
|
---|
| 75 |
|
---|
| 76 | // *** type de generation
|
---|
| 77 | bool computefourier0=false;
|
---|
| 78 | unsigned short nthread=4;
|
---|
| 79 |
|
---|
| 80 | // --- Decodage arguments
|
---|
| 81 |
|
---|
| 82 | char c;
|
---|
| 83 | while((c = getopt(narg,arg,"ha0x:y:z:s:Z:")) != -1) {
|
---|
| 84 | switch (c) {
|
---|
| 85 | case 'a' :
|
---|
| 86 | Auto_Ini_Ranf(5);
|
---|
| 87 | break;
|
---|
| 88 | case '0' :
|
---|
| 89 | computefourier0 = true;
|
---|
| 90 | break;
|
---|
| 91 | case 'x' :
|
---|
[3129] | 92 | sscanf(optarg,"%ld,%lf",&nx,&dx);
|
---|
[3115] | 93 | break;
|
---|
| 94 | case 'y' :
|
---|
[3129] | 95 | sscanf(optarg,"%ld,%lf",&ny,&dy);
|
---|
[3115] | 96 | break;
|
---|
| 97 | case 'z' :
|
---|
[3129] | 98 | sscanf(optarg,"%ld,%lf",&nz,&dz);
|
---|
[3115] | 99 | break;
|
---|
| 100 | case 's' :
|
---|
| 101 | sscanf(optarg,"%lf",&snoise);
|
---|
| 102 | break;
|
---|
| 103 | case 'Z' :
|
---|
| 104 | sscanf(optarg,"%lf",&zref);
|
---|
| 105 | break;
|
---|
| 106 | case 'h' :
|
---|
| 107 | default :
|
---|
| 108 | usage(); return -1;
|
---|
| 109 | }
|
---|
| 110 | }
|
---|
| 111 |
|
---|
| 112 | string tagobs = "cmvobserv3d.ppf";
|
---|
| 113 | POutPersist posobs(tagobs);
|
---|
| 114 |
|
---|
| 115 | cout<<"zref="<<zref<<endl;
|
---|
| 116 | cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
|
---|
| 117 | cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"
|
---|
| 118 | <<", npt="<<npt<<endl;
|
---|
| 119 | cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
|
---|
| 120 | cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
|
---|
| 121 | cout<<"schmin="<<schmin<<" ("<<lschmin
|
---|
| 122 | <<"), schmax="<<schmax<<" ("<<lschmax<<") Msol"
|
---|
| 123 | <<", schnpt="<<schnpt<<endl;
|
---|
| 124 | cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
|
---|
| 125 |
|
---|
| 126 | //-----------------------------------------------------------------
|
---|
| 127 | cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
|
---|
| 128 |
|
---|
| 129 | InitialSpectrum pkini(ns,as);
|
---|
| 130 |
|
---|
| 131 | TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
|
---|
| 132 | //tf.SetNoOscEnv(2);
|
---|
| 133 |
|
---|
| 134 | GrowthFactor d1(om0,ol0);
|
---|
| 135 |
|
---|
| 136 | PkSpectrum0 pk0(pkini,tf);
|
---|
| 137 |
|
---|
| 138 | PkSpectrumZ pkz(pk0,d1,zref);
|
---|
| 139 |
|
---|
| 140 | Schechter sch(nstar,mstar,alpha);
|
---|
| 141 |
|
---|
| 142 | //-----------------------------------------------------------------
|
---|
| 143 | pkz.SetZ(0.);
|
---|
| 144 | cout<<endl<<"\n--- Compute variance for top-hat R="<<R
|
---|
| 145 | <<" at z="<<pkz.GetZ()<<endl;
|
---|
| 146 | VarianceSpectrum varpk_th(pkz,0);
|
---|
| 147 | double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
|
---|
| 148 | double pkmax_th = varpk_th(kfind_th);
|
---|
| 149 | cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
|
---|
| 150 | double k1=kmin, k2=kmax;
|
---|
| 151 | int rc = varpk_th.FindLimits(R,pkmax_th/1.e4,k1,k2,eps);
|
---|
| 152 | cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
|
---|
| 153 | <<k2<<" ("<<log10(k2)<<")"<<endl;
|
---|
| 154 |
|
---|
| 155 | double ldlk = (log10(k2)-log10(k1))/npt;
|
---|
| 156 | varpk_th.SetInteg(0.01,ldlk,-1.,4);
|
---|
| 157 | double sr2 = varpk_th.Variance(R,k1,k2);
|
---|
| 158 | cout<<"varpk_th="<<sr2<<" -> sigma="<<sqrt(sr2)<<endl;
|
---|
| 159 |
|
---|
| 160 | double normpkz = sigmaR*sigmaR/sr2;
|
---|
| 161 | pkz.SetScale(normpkz);
|
---|
| 162 | cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
|
---|
| 163 |
|
---|
| 164 | pkz.SetZ(zref);
|
---|
| 165 |
|
---|
[3120] | 166 | Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
|
---|
[3115] | 167 | FuncToHisto(pkz,hpkz,true);
|
---|
| 168 | {
|
---|
| 169 | tagobs = "hpkz"; posobs.PutObject(hpkz,tagobs);
|
---|
| 170 | }
|
---|
| 171 |
|
---|
| 172 | //-----------------------------------------------------------------
|
---|
| 173 | cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
|
---|
| 174 | VarianceSpectrum varpk_int(pkz,2);
|
---|
| 175 |
|
---|
| 176 | double kfind_int = varpk_int.FindMaximum(R,kmin,kmax,eps);
|
---|
| 177 | double pkmax_int = varpk_int(kfind_int);
|
---|
| 178 | cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
|
---|
| 179 | double k1int=kmin, k2int=kmax;
|
---|
| 180 | int rcint = varpk_int.FindLimits(R,pkmax_int/1.e4,k1int,k2int,eps);
|
---|
| 181 | cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
|
---|
| 182 | <<k2int<<" ("<<log10(k2int)<<")"<<endl;
|
---|
| 183 |
|
---|
| 184 | double ldlkint = (log10(k2int)-log10(k1int))/npt;
|
---|
| 185 | varpk_int.SetInteg(0.01,ldlkint,-1.,4);
|
---|
| 186 | double sr2int = varpk_int.Variance(R,k1int,k2int);
|
---|
| 187 | cout<<"varpk_int="<<sr2int<<" -> sigma="<<sqrt(sr2int)<<endl;
|
---|
| 188 |
|
---|
| 189 | PrtTim("End of definition");
|
---|
| 190 |
|
---|
| 191 | //-----------------------------------------------------------------
|
---|
| 192 | cout<<endl<<"\n--- Compute galaxies density number"<<endl;
|
---|
| 193 |
|
---|
| 194 | sch.SetOutValue(0);
|
---|
| 195 | cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
|
---|
| 196 | double ngal_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
|
---|
| 197 | cout<<"Galaxy density number = "<<ngal_by_mpc3<<" /Mpc^3 between limits"<<endl;
|
---|
| 198 |
|
---|
| 199 | sch.SetOutValue(1);
|
---|
| 200 | cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
|
---|
| 201 | double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
|
---|
| 202 | cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
|
---|
| 203 |
|
---|
| 204 | //-----------------------------------------------------------------
|
---|
| 205 | // FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f with e+f=0 ou 1
|
---|
| 206 | cout<<endl<<"\n--- Compute Spectrum Realization"<<endl;
|
---|
| 207 |
|
---|
| 208 | pkz.SetZ(zref);
|
---|
| 209 | TArray< complex<r_8> > pkgen;
|
---|
| 210 | GeneFluct3D fluct3d(pkgen,pkz);
|
---|
| 211 | fluct3d.SetNThread(nthread);
|
---|
| 212 | fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
|
---|
| 213 | fluct3d.Print();
|
---|
| 214 | double knyqmax = fluct3d.GetKmax();
|
---|
| 215 | double dkmin = fluct3d.GetKinc()[0];
|
---|
| 216 |
|
---|
| 217 | cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
|
---|
| 218 | // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
|
---|
| 219 | // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
|
---|
| 220 | // Vc/Vs = 0.675 -> keff = kmax * (0.675)^(1/3) = kmax * 0.877
|
---|
| 221 | knyqmax *= 0.877;
|
---|
| 222 | ldlkint = (log10(knyqmax)-log10(k1int))/npt;
|
---|
| 223 | varpk_int.SetInteg(0.01,ldlkint,-1.,4);
|
---|
| 224 | double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax);
|
---|
| 225 | cout<<"varpk_int(<"<<knyqmax<<")="<<sr2int_kmax<<" -> sigma="<<sqrt(sr2int_kmax)<<endl;
|
---|
| 226 |
|
---|
| 227 | cout<<"\n--- Computing a realization in Fourier space"<<endl;
|
---|
| 228 | if(computefourier0) fluct3d.ComputeFourier0();
|
---|
| 229 | else fluct3d.ComputeFourier();
|
---|
| 230 |
|
---|
| 231 | cout<<"\n--- Checking realization spectra"<<endl;
|
---|
[3129] | 232 | long nhprof = long(fluct3d.GetKmax()/dkmin+0.5);
|
---|
[3115] | 233 | HProf hpkgen(0.,fluct3d.GetKmax(),nhprof);
|
---|
| 234 | hpkgen.ReCenterBin();
|
---|
| 235 | fluct3d.ComputeSpectrum(hpkgen);
|
---|
| 236 | {
|
---|
| 237 | tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
|
---|
| 238 | }
|
---|
| 239 |
|
---|
| 240 | if(1) {
|
---|
| 241 | cout<<"\n--- Computing convolution by pixel shape"<<endl;
|
---|
| 242 | fluct3d.FilterByPixel();
|
---|
| 243 |
|
---|
| 244 | cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
|
---|
| 245 | HProf hpkgenf(hpkgen); hpkgenf.Zero();
|
---|
| 246 | fluct3d.ComputeSpectrum(hpkgenf);
|
---|
| 247 | {
|
---|
| 248 | tagobs = "hpkgenf"; posobs.PutObject(hpkgenf,tagobs);
|
---|
| 249 | }
|
---|
| 250 | }
|
---|
| 251 |
|
---|
| 252 | if(0) {
|
---|
| 253 | cout<<"\n--- Writing cmvobserv3dk.ppf"<<endl;
|
---|
| 254 | string tag = "cmvobserv3dk.ppf";
|
---|
| 255 | POutPersist pos(tag);
|
---|
| 256 | tag = "pkgen"; pos.PutObject(pkgen,tag);
|
---|
| 257 | }
|
---|
| 258 |
|
---|
| 259 | cout<<"\n--- Computing a realization in real space"<<endl;
|
---|
| 260 | fluct3d.ComputeReal();
|
---|
| 261 | r_8 undouble=0.;
|
---|
| 262 | TArray<r_8> rgen = ArrayCast(pkgen,undouble);
|
---|
| 263 | double rmin,rmax; rgen.MinMax(rmin,rmax);
|
---|
| 264 | cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
|
---|
| 265 |
|
---|
| 266 | cout<<"\n--- Check mean and variance in real space"<<endl;
|
---|
| 267 | size_t nlowone = fluct3d.NumberOfBad(-1.,1e+200);
|
---|
| 268 | double rm,rs2; size_t nm;
|
---|
| 269 | nm = fluct3d.MeanSigma2(rm,rs2);
|
---|
| 270 | cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 271 | <<rs2<<" -> "<<sqrt(rs2)<<", n(<-1)="<<nlowone<<endl;
|
---|
| 272 |
|
---|
| 273 | if(1) {
|
---|
| 274 | cout<<"\n--- Check variance sigmaR in real space"<<endl;
|
---|
| 275 | double varr;
|
---|
| 276 | size_t nvarr = fluct3d.VarianceFrReal(R,varr);
|
---|
| 277 | cout<<"R="<<R<<" : sigmaR^2="<<varr<<" -> "<<sqrt(varr)<<", n="<<nvarr<<endl;
|
---|
| 278 | }
|
---|
| 279 |
|
---|
| 280 | //-----------------------------------------------------------------
|
---|
| 281 | cout<<endl<<"\n--- Converting fluctuations into mass"<<endl;
|
---|
| 282 | fluct3d.TurnFluct2Mass();
|
---|
| 283 | nm = fluct3d.MeanSigma2(rm,rs2);
|
---|
| 284 | cout<<"1+rgen: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 285 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 286 |
|
---|
| 287 | cout<<"\n--- Converting mass into galaxy number"<<endl;
|
---|
| 288 | rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
|
---|
| 289 | cout<<rm<<" galaxies put into survey"<<endl;
|
---|
| 290 | nm = fluct3d.MeanSigma2(rm,rs2,0.);
|
---|
| 291 | cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 292 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 293 |
|
---|
| 294 | cout<<"\n--- Set negative pixels to BAD"<<endl;
|
---|
| 295 | nm = fluct3d.SetToVal(0.,1e+200,-999.);
|
---|
| 296 | cout<<nm<<" negative in survey set to BAD"<<endl;
|
---|
| 297 | nm = fluct3d.MeanSigma2(rm,rs2,-998.);
|
---|
| 298 | cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 299 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 300 |
|
---|
| 301 | cout<<"\n--- Apply poisson on galaxie number"<<endl;
|
---|
| 302 | nm = fluct3d.ApplyPoisson();
|
---|
| 303 | cout<<nm<<" galaxies into survey after poisson"<<endl;
|
---|
| 304 | nm = fluct3d.MeanSigma2(rm,rs2,-998.);
|
---|
| 305 | cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 306 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 307 |
|
---|
| 308 | cout<<"\n--- Convert Galaxies number to HI mass"<<endl;
|
---|
[3129] | 309 | long nhmdndm = long( (lschmax-lschmin+1.)*100. + 0.5);
|
---|
[3115] | 310 | Histo hmdndm(lschmin,lschmax,nhmdndm);
|
---|
| 311 | sch.SetOutValue(1);
|
---|
| 312 | FuncToHisto(sch,hmdndm,true);
|
---|
| 313 | FunRan tirhmdndm(hmdndm,true);
|
---|
| 314 | {
|
---|
| 315 | tagobs = "hmdndm"; posobs.PutObject(hmdndm,tagobs);
|
---|
| 316 | Histo hdum1(tirhmdndm);
|
---|
| 317 | tagobs = "tirhmdndm"; posobs.PutObject(hdum1,tagobs);
|
---|
| 318 | }
|
---|
| 319 | double mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
|
---|
| 320 | cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
|
---|
| 321 | nm = fluct3d.MeanSigma2(rm,rs2,0.);
|
---|
| 322 | cout<<"HI mass: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 323 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 324 | cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
|
---|
| 325 |
|
---|
| 326 | cout<<"\n--- Set BAD pixels to Zero"<<endl;
|
---|
| 327 | nm = fluct3d.SetToVal(-998.,1e+200,0.);
|
---|
| 328 | cout<<nm<<" BAD in survey set to zero"<<endl;
|
---|
| 329 | nm = fluct3d.MeanSigma2(rm,rs2);
|
---|
| 330 | cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 331 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 332 |
|
---|
[3120] | 333 | if(1) {
|
---|
| 334 | cout<<"\n---Writing Cube to Fits file"<<endl;
|
---|
| 335 | FitsImg3DWriter fwrt("!cmvobserv3dr.fits",FLOAT_IMG,5);
|
---|
| 336 | fwrt.WriteKey("ZREF",zref," redshift");
|
---|
| 337 | vector<long> n = fluct3d.GetNpix();
|
---|
| 338 | fwrt.WriteKey("NX",n[0]," axe transverse 1");
|
---|
| 339 | fwrt.WriteKey("NY",n[1]," axe transverse 2");
|
---|
| 340 | fwrt.WriteKey("NZ",n[2]," axe longitudinal (redshift)");
|
---|
| 341 | vector<r_8> d;
|
---|
| 342 | d = fluct3d.GetDinc();
|
---|
| 343 | fwrt.WriteKey("DX",d[0]," Mpc");
|
---|
| 344 | fwrt.WriteKey("DY",d[1]," Mpc");
|
---|
| 345 | fwrt.WriteKey("DZ",d[2]," Mpc");
|
---|
| 346 | d = fluct3d.GetKinc();
|
---|
| 347 | fwrt.WriteKey("DKX",d[0]," Mpc^-1");
|
---|
| 348 | fwrt.WriteKey("DKY",d[1]," Mpc^-1");
|
---|
| 349 | fwrt.WriteKey("DKZ",d[2]," Mpc^-1");
|
---|
| 350 | fwrt.WriteKey("SNOISE",snoise," Msol");
|
---|
| 351 | fwrt.Write(rgen);
|
---|
| 352 | }
|
---|
| 353 |
|
---|
[3115] | 354 | cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
|
---|
| 355 | fluct3d.AddNoise2Real(snoise);
|
---|
| 356 | nm = fluct3d.MeanSigma2(rm,rs2);
|
---|
| 357 | cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
|
---|
| 358 | <<rs2<<" -> "<<sqrt(rs2)<<endl;
|
---|
| 359 |
|
---|
| 360 | if(0) {
|
---|
| 361 | cout<<"\n--- Writing cmvobserv3dr.ppf"<<endl;
|
---|
| 362 | string tag = "cmvobserv3dr.ppf";
|
---|
| 363 | POutPersist pos(tag);
|
---|
| 364 | tag = "rgen"; pos.PutObject(rgen,tag);
|
---|
| 365 | }
|
---|
| 366 |
|
---|
| 367 | //-----------------------------------------------------------------
|
---|
| 368 | // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
|
---|
| 369 | if(1) {
|
---|
| 370 | cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
|
---|
| 371 | fluct3d.ReComputeFourier();
|
---|
| 372 | HProf hpkrec(0.,fluct3d.GetKmax(),nhprof);
|
---|
| 373 | hpkrec.ReCenterBin();
|
---|
| 374 | fluct3d.ComputeSpectrum(hpkrec);
|
---|
| 375 | tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
|
---|
| 376 | }
|
---|
| 377 |
|
---|
| 378 | return 0;
|
---|
| 379 | }
|
---|
| 380 |
|
---|
| 381 | /*
|
---|
| 382 | openppf cmvobserv3dk.ppf
|
---|
| 383 | openppf cmvobserv3dr.ppf
|
---|
| 384 | openppf cmvobserv3d.ppf
|
---|
| 385 |
|
---|
| 386 | objaoper pkgen sliceyz 0
|
---|
| 387 |
|
---|
[3120] | 388 | n/plot hpkz.val%x ! ! "nsta connectpoints"
|
---|
| 389 | n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
|
---|
| 390 | n/plot hpkgenf.val%log10(x) x>0 ! "nsta same orange connectpoints"
|
---|
| 391 | n/plot hpkrec.val%log10(x) x>0 ! "nsta same blue connectpoints"
|
---|
[3115] | 392 |
|
---|
| 393 | set k pow(10.,x)
|
---|
| 394 | n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
|
---|
| 395 |
|
---|
| 396 | defscript rgensl
|
---|
| 397 | objaoper rgen sliceyz $1
|
---|
| 398 | disp sliceyz_$1
|
---|
| 399 | endscript
|
---|
| 400 |
|
---|
| 401 | rgensl 0
|
---|
| 402 |
|
---|
[3120] | 403 | disp hmdndm
|
---|
| 404 | disp tirhmdndm
|
---|
| 405 | addline 0 1 20 1 "red"
|
---|
| 406 |
|
---|
[3115] | 407 | */
|
---|