source: Sophya/trunk/Cosmo/SimLSS/cmvdefsurv.cc@ 3524

Last change on this file since 3524 was 3516, checked in by cmv, 17 years ago

modifs mineures cmv 9/9/2008
File size: 38.2 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
10#include <vector>
11
12#include "constcosmo.h"
13#include "cosmocalc.h"
14#include "geneutils.h"
15#include "schechter.h"
16#include "pkspectrum.h"
17#include "planckspectra.h"
18
19/* --- Check Peterson at al. astro-ph/0606104 v1 (pb facteur sqrt(2) sur S/N !)
20cmvdefsurv -U 0.75,0.3,0.7,-1,1 -V 300 -z 0.0025,0.2,Z -x 1,90,A -O 400000,6000 -N 75 -M 6.156e9 -F 3 -2 1.5
21 --- */
22
23//-----------------------------------------------------------------------------------------------------------
24inline double rad2deg(double trad) {return trad/M_PI*180.;}
25inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
26inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
27inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
28inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
29inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
30
31inline double sr2deg2(double asr) {return asr*pow(rad2deg(1.),2.);}
32inline double sr2min2(double asr) {return asr*pow(rad2min(1.),2.);}
33inline double deg22sr(double adeg) {return adeg*pow(deg2rad(1.),2.);}
34inline double min22sr(double amin) {return amin*pow(min2rad(1.),2.);}
35
36double LargeurDoppler(double v /* km/s */, double nu);
37double DzFrV(double v /* km/s */, double zred);
38double DNuFrDz(double dzred,double nu_at_0,double zred);
39double DzFrDNu(double dnu_at_0,double nu_at_0,double zred);
40double DzFrDNuApprox(double dnu_at_0,double nu_at_0,double zred);
41double ZFrLos(double loscom /* Mpc com */,CosmoCalc& univ);
42double AngsolEqTelescope(double nu /* GHz */,double telsurf /* m^2 */);
43void usage(void);
44
45void usage(void) {
46 cout<<"cmvdefsurv [options] -x dx,txlarg[,unit_x] -y dy,tylarg[,unit_y] -z dz,tzlarg[,unit_z] redshift"<<endl
47 <<"----------------"<<endl
48 <<" -x dx,txlarg[,unit_x] : resolution et largeur dans le plan transverse selon X"<<endl
49 <<" -y dy,tylarg[,unit_y] : idem selon Y, si <=0 meme que X"<<endl
50 <<" -z dz,tzlarg[,unit_z] : resolution et largeur sur la ligne de visee"<<endl
51 <<"-- Unites pour X-Y:"<<endl
52 <<" \'A\' : en angles (pour X-Y) : resolution=ArcMin, largeur=Degre (defaut)"<<endl
53 <<" \'Z\' : en redshift (pour Z) : resolution et largeur en redshift (defaut)"<<endl
54 <<" \'F\' : en frequence (pour Z) : resolution et largeur MHz"<<endl
55 <<" \'M\' : en longeur comobile (pour X-Y-Z) : resolution et largeur Mpc"<<endl
56 <<"----------------"<<endl
57 <<" -K k,dk,pk : k(Mpc^-1) dk(Mpc^-1) pk(a z=0 en Mpc^-3) pour estimer la variance cosmique"<<endl
58 <<"----------------"<<endl
59 <<" -O surf,tobs,eta : surface effective (m^2), temps d\'observation (s), efficacite d\'ouverture"<<endl
60 <<" -N Tsys : temperature du system (K)"<<endl
61 <<" -L lobearea,freqlob : angle solide du lobe d\'observation en arcmin^2 (def= celle du pixel)"<<endl
62 <<" pour la frequence freqlob en MHz"<<endl
63 <<" Si freqlob<=0 : la frequence de reference est celle du redshift etudie (def)"<<endl
64 <<" -2 : two polarisations measured"<<endl
65 <<" -M : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
66 <<" -F : HI flux factor to be applied for our redshift"<<endl
67 <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
68 <<"----------------"<<endl
69 <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
70 <<" (indnu==0 no evolution with freq.)"<<endl
71 <<"----------------"<<endl
72 <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
73 <<"----------------"<<endl
74 <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
75 <<" redshift : redshift moyen du survey"<<endl
76 <<endl;
77}
78
79//-------------------------------------------------------------------------------------------
80int main(int narg,char *arg[])
81{
82 // ---
83 // --- Valeurs fixes ou par defaut
84 // ---
85
86 //-- WMAP
87 unsigned short flat = 0;
88 double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
89
90 //-- Schechter HIMASS
91 double h75 = h100 / 0.75;
92 double nstar = 0.006*pow(h75,3.); //
93 double mstar = pow(10.,9.8); // MSol
94 double alpha = -1.37;
95 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
96
97 double hifactor = 1.;
98 double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
99
100 //-- CMB , Synchrotron et AGN
101 double tcmb = T_CMB_Par;
102 // a 408 MHz (Haslam) + evol index a -2.6
103 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
104 double lflux_agn = -3.;
105
106 //-- Appareillage
107 bool ya2polar = false;
108 double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
109 double Tsys=75.;
110 double tobs = 6000., surftot = 400000., eta_eff = 1.;
111 // taille du lobe d'observation en arcmin pour la frequence
112 double lobearea0 = -1., lobefreq0 = -1.;
113
114 //-- Variance cosmique (default = standard SDSSII)
115 double kcosm = 0.05, dkcosm = -1., pkcosm = 40000.;
116
117 // --- Pour taille du survey
118 double dx=1., dy=-1., dz=0.0007, txlarg=100., tylarg=-1., tzlarg=0.1;
119 int nx=0, ny=0, nz=0;
120 char unit_x = 'A', unit_y = 'A', unit_z = 'Z';
121 double redshift0 = 0.;
122 double mhiref = -1.; // reference Mass en HI (def integ schechter)
123
124 // ---
125 // --- Decodage arguments
126 // ---
127 char c;
128 while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:K:")) != -1) {
129 switch (c) {
130 case 'x' :
131 sscanf(optarg,"%lf,%lf,%c",&dx,&txlarg,&unit_x);
132 break;
133 case 'y' :
134 sscanf(optarg,"%lf,%lf,%c",&dy,&tylarg,&unit_y);
135 break;
136 case 'z' :
137 sscanf(optarg,"%lf,%lf,%c",&dz,&tzlarg,&unit_z);
138 break;
139 case 'O' :
140 sscanf(optarg,"%lf,%lf,%lf",&surftot,&tobs,&eta_eff);
141 break;
142 case 'L' :
143 sscanf(optarg,"%lf,%lf",&lobearea0,&lobefreq0);
144 break;
145 case 'N' :
146 sscanf(optarg,"%lf",&Tsys);
147 break;
148 case 'S' :
149 sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
150 break;
151 case 'M' :
152 sscanf(optarg,"%lf",&mhiref);
153 break;
154 case 'F' :
155 sscanf(optarg,"%lf",&hifactor);
156 break;
157 case 'V' :
158 sscanf(optarg,"%lf",&vrot);
159 break;
160 case 'U' :
161 sscanf(optarg,"%lf,%lf,%lf,%lf,%hu",&h100,&om0,&ol0,&w0,&flat);
162 break;
163 case '2' :
164 ya2polar = true;
165 facpolar = 1.0;
166 break;
167 case 'A' :
168 sscanf(optarg,"%lf",&lflux_agn);
169 break;
170 case 'K' :
171 sscanf(optarg,"%lf,%lf,%lf",&kcosm,&dkcosm,&pkcosm);
172 break;
173 case 'h' :
174 default :
175 usage(); return -1;
176 }
177 }
178 if(optind>=narg) {usage(); return -1;}
179 sscanf(arg[optind],"%lf",&redshift0);
180 if(redshift0<=0.) {cout<<"Redshift "<<redshift0<<" should be >0"<<endl; return -2;}
181
182 // ---
183 // --- Initialisation de la Cosmologie
184 // ---
185 CosmoCalc univ(flat,true,2.*redshift0);
186 double perc=0.01,dzinc=redshift0/100.,dzmax=dzinc*10.; unsigned short glorder=4;
187 univ.SetInteg(perc,dzinc,dzmax,glorder);
188 univ.SetDynParam(h100,om0,or0,ol0,w0);
189
190 GrowthFactor growth(om0,ol0);
191
192 // ---
193 // --- Mise en forme dans les unites appropriees
194 // ---
195 // ATTENTION: le cube de simulation est en Mpc avec des pixels de taille comobile fixe
196 cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
197 if(dy<=0. || tylarg<=0.) {dy=dx; tylarg=txlarg; unit_y=unit_x;}
198 cout<<"X values: resolution="<<dx<<" largeur="<<txlarg<<" unite="<<unit_x<<endl;
199 if(unit_x == 'A') {
200 nx = int(txlarg*60./dx+0.5);
201 dx = min2rad(dx)*univ.Dtrcom(redshift0);
202 txlarg = dx*nx;
203 } else if(unit_x == 'M') {
204 nx = int(txlarg/dx+0.5);
205 txlarg = dx*nx;
206 } else {
207 cout<<"Unknown unit_x = "<<unit_x<<endl; return -2;
208 }
209 cout<<"Y values: resolution="<<dy<<" largeur="<<tylarg<<" unite="<<unit_y<<endl;
210 if(unit_y == 'A') {
211 ny = int(tylarg*60./dy+0.5);
212 dy = min2rad(dy)*univ.Dtrcom(redshift0);
213 tylarg = dy*ny;
214 } else if(unit_y == 'M') {
215 ny = int(tylarg/dy+0.5);
216 tylarg = dy*ny;
217 } else {
218 cout<<"Unknown unit_y = "<<unit_y<<endl; return -2;
219 }
220 cout<<"Z values: resolution="<<dz<<" largeur="<<tzlarg<<" unite="<<unit_z<<endl;
221 if(unit_z == 'Z') {
222 nz = int(tzlarg/dz+0.5);
223 dz = dz*univ.Dhubble()/univ.E(redshift0);
224 tzlarg = dz*nz;
225 } else if(unit_z == 'M') {
226 nz = int(tzlarg/dz+0.5);
227 tzlarg = dz*nz;
228 } else if(unit_z == 'F') { // unite en MHz
229 nz = int(tzlarg/dz+0.5);
230 dz = DzFrDNu(dz,Fr_HyperFin_Par*1.e3,redshift0); //dzred
231 dz = dz*univ.Dhubble()/univ.E(redshift0);
232 tzlarg = dz*nz;
233 } else {
234 cout<<"Unknown unit_z = "<<unit_z<<endl; return -2;
235 }
236
237 // On estime la valeur du redshift le plus grand
238 double bigred = redshift0 + dz/univ.Dhubble()*univ.E(redshift0) * nz/2.;
239 cout<<"biggest redshift estimated to be "<<bigred<<endl;
240 dzinc=bigred/100.; dzmax=dzinc*10.;
241 univ.SetInteg(perc,dzinc,dzmax,glorder);
242
243 // ---
244 // --- Calcul des valeurs au centre du cube et aux faces limites av/ar
245 // ---
246 cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
247 double adtx[3], adty[3], atxlarg[3], atylarg[3];
248 double loscom[3], redshift[3], unplusz[3], dred[3], growthfac[3], rhocz[3];
249 double dang[3], dtrcom[3], dlum[3], dloscom[3], dlosdz[3];
250 double nuhiz[3], lambdahiz[3], dnuhiz[3];
251 double angsol_pix[3], angsol_tot[3];
252
253 redshift[0] = redshift0;
254 loscom[0] = univ.Dloscom(redshift0);
255 loscom[1] = loscom[0]-tzlarg/2.;
256 loscom[2] = loscom[0]+tzlarg/2.;
257 //if(loscom[1]<=0.) {cout<<"Lower distance limit "<<loscom[1]<<" should be >0"<<endl; return -3;}
258 for(int i=0;i<3;i++) {
259 double l = loscom[i]; if(l<=0.) l = dz/2.;
260 redshift[i] = ZFrLos(l,univ);
261 unplusz[i] = 1. + redshift[i];
262 growthfac[i] = growth(redshift[i]);
263 rhocz[i] = univ.Rhoc(redshift[i])*GCm3toMsolMpc3_Cst;;
264 dang[i] = univ.Dang(redshift[i]);
265 dtrcom[i] = univ.Dtrcom(redshift[i]);
266 dlum[i] = univ.Dlum(redshift[i]);
267 dloscom[i] = univ.Dloscom(redshift[i]);
268 dlosdz[i] = univ.Dhubble()/univ.E(redshift[i]);
269 dred[i] = dz/dlosdz[i];
270 adtx[i] = dx/dtrcom[i];
271 atxlarg[i] = adtx[i]*nx;
272 adty[i] = dy/dtrcom[i];
273 atylarg[i] = adty[i]*ny;
274 nuhiz[i] = Fr_HyperFin_Par / unplusz[i]; // GHz
275 lambdahiz[i] = SpeedOfLight_Cst*1000./(nuhiz[i]*1.e9); // m
276 dnuhiz[i] = DNuFrDz(dred[i],Fr_HyperFin_Par,redshift[i]); // GHz
277 angsol_pix[i] = AngSol(adtx[i]/2.,adty[i]/2.,M_PI/2.);
278 angsol_tot[i] = AngSol(atxlarg[i]/2.,atylarg[i]/2.,M_PI/2.);
279 }
280
281
282 cout<<"--- Cosmology for z = "<<0.<<endl;
283 univ.Print(0.);
284 double rhoc0 = univ.Rhoc(0.)*GCm3toMsolMpc3_Cst;
285
286 for(int i=0;i<3;i++) {
287 cout<<"\n--- Cosmology for z = "<<redshift[i]<<endl;
288 univ.Print(redshift[i]);
289 cout<<"Nu(HI) = "<<Fr_HyperFin_Par*1000./(1.+redshift[i])<<" MHz"<<endl;
290 }
291
292 cout<<endl;
293 cout<<"--- Resume"<<endl;
294 cout<<"Facteur de croissance lineaire = "<<growthfac[0]
295 <<" in ["<<growthfac[1]<<","<<growthfac[2]<<"]"<<endl;
296 cout<<" 1/(1+z) = "<<1./(1.+redshift[0])
297 <<" in ["<<1./(1.+redshift[1])<<","<<1./(1.+redshift[2])<<"]"<<endl;
298 cout<<"Facteur de croissance lineaire^2 = "<<growthfac[0]*growthfac[0]
299 <<" in ["<<growthfac[1]*growthfac[1]<<","<<growthfac[2]*growthfac[2]<<"]"<<endl;
300
301 cout<<"Rho_c (z=0) = "<<rhoc0<<", a z="<<0<<": "<<rhoc0<<" Msol/Mpc^3"<<endl;
302 cout<<"Rho_c a z = "<<rhocz[0]
303 <<" Msol/Mpc^3 in ["<<rhocz[1]<<","<<rhocz[2]<<"]"<<endl;
304
305 cout<<endl;
306 cout<<"dang= "<<dang[0]<<" Mpc in ["<<dang[1]<<","<<dang[2]<<"]"<<endl;
307 cout<<"dtrcom= "<<dtrcom[0]<<" Mpc com in ["<<dtrcom[1]<<","<<dtrcom[2]<<"]"<<endl;
308 cout<<"dlum= "<<dlum[0]<<" Mpc in ["<<dlum[1]<<","<<dlum[2]<<"]"<<endl;
309 cout<<"dloscom= "<<dloscom[0]<<" Mpc com in ["<<dloscom[1]<<","<<dloscom[2]<<"]"<<endl;
310 cout<<"dz=1 -> dlosdz= "<<dlosdz[0]<<" Mpc com in ["<<dlosdz[1]<<","<<dlosdz[2]<<"]"<<endl;
311 cout<<"1 Mpc los com -> dz = "<<1./dlosdz[0]<<" in ["<<1./dlosdz[1]<<","<<1./dlosdz[2]<<"]"<<endl;
312
313 cout<<endl;
314 for(int i=0;i<3;i++) {
315 cout<<"...Redshift="<<redshift[i]<<endl;
316 cout<<"1\" -> "<<dang[i]*sec2rad(1.)<<" Mpc = "<<dtrcom[i]*sec2rad(1.)<<" Mpc com"<<endl;
317 cout<<"1\' -> "<<dang[i]*min2rad(1.)<<" Mpc = "<<dtrcom[i]*min2rad(1.)<<" Mpc com"<<endl;
318 cout<<"1d -> "<<dang[i]*deg2rad(1.)<<" Mpc = "<<dtrcom[i]*deg2rad(1.)<<" Mpc com"<<endl;
319 cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom[i])<<"\" = "
320 <<rad2min(1./dtrcom[i])<<" \' = "<<rad2deg(1./dtrcom[i])<<" d"<<endl;
321 cout<<"dz=1 -> dlosdz= "<<dlosdz[i]<<" Mpc com"<<endl;
322 cout<<"1 Mpc los com -> dz = "<<1./dlosdz[0]<<endl;
323 }
324
325 // ---
326 // --- Cosmolographie Transverse
327 // ---
328 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
329
330 cout<<"dx = "<<dx<<", txlarg = "<<txlarg<<" Mpc (com), nx="<<nx<<endl;
331 cout<<"adtx = "<<adtx[0]<<" rd in ["<<adtx[1]<<","<<adtx[2]<<"]"<<endl;
332 cout<<" = "<<rad2min(adtx[0])<<"\' in ["<<rad2min(adtx[1])<<","<<rad2min(adtx[2])<<"]"<<endl;
333 cout<<"atxlarg = "<<atxlarg[0]<<" rd in ["<<atxlarg[1]<<","<<atxlarg[2]<<"]"<<endl;
334 cout<<" "<<rad2deg(atxlarg[0])<<" d in ["<<rad2deg(atxlarg[1])<<","<<rad2deg(atxlarg[2])<<"]"<<endl;
335
336 if(fabs(dx-dy)>1.e-20 && fabs(txlarg-tylarg)>1.e-20) {
337 cout<<"\ndy = "<<dy<<", tylarg = "<<tylarg<<" Mpc (com), ny="<<ny<<endl;
338 cout<<"adty = "<<adty[0]<<" rd in ["<<adty[1]<<","<<adty[2]<<"]"<<endl;
339 cout<<" = "<<rad2min(adty[0])<<"\' in ["<<rad2min(adty[1])<<","<<rad2min(adty[2])<<"]"<<endl;
340 cout<<"atylarg = "<<atylarg[0]<<" rd in ["<<atylarg[1]<<","<<atylarg[2]<<"]"<<endl;
341 cout<<" "<<rad2deg(atylarg[0])<<" d in ["<<rad2deg(atylarg[1])<<","<<rad2deg(atylarg[2])<<"]"<<endl;
342 }
343
344 // ---
345 // --- Cosmolographie Line of sight
346 // ---
347 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
348 cout<<"dz = "<<dz<<", tzlarg = "<<tzlarg<<" Mpc (com), nz="<<nz<<endl;
349 cout<<"Redshift = "<<redshift[0]<<" in ["<<redshift[1]<<","<<redshift[2]<<"]"<<endl;
350 cout<<"dred = "<<dred[0]<<" in ["<<dred[1]<<","<<dred[2]<<"]"<<endl;
351 cout<<"nu HI = "<<nuhiz[0]<<" GHz in ["<<nuhiz[1]<<","<<nuhiz[2]<<"]"<<endl;
352 cout<<"lambda HI = "<<lambdahiz[0]<<" m in ["<<lambdahiz[1]<<","<<lambdahiz[2]<<"]"<<endl;
353 cout<<"dnu HI= "<<dnuhiz[0]*1e3<<" MHz in ["<<dnuhiz[1]*1e3<<","<<dnuhiz[2]*1e3<<"]"<<endl;
354
355 // ---
356 // --- Volume
357 // ---
358 cout<<"\n>>>>\n>>>> Volumes\n>>>>"<<endl;
359 double Npix = (double)nx*(double)ny*(double)nz;
360 double vol_pixel = dx*dy*dz;
361 double vol_survey = vol_pixel*Npix;
362 cout<<"Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
363 cout<<"Volume pixel = "<<vol_pixel<<" Mpc^3 com"<<endl;
364 cout<<"Volume total = "<<vol_survey<<" Mpc^3 com = "<<vol_survey/1.e9<<" Gpc^3"<<endl;
365
366 double vol = univ.Vol4Pi(redshift[1],redshift[2])/(4.*M_PI)*angsol_tot[0];
367 cout<<"Calcul avec angsol et redshift: vol = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
368 <<" -> pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
369
370 // ---
371 // --- Angles solides
372 // ---
373 cout<<"\n>>>>\n>>>> Angles solides\n>>>>"<<endl;
374 cout<<"Pixel solid angle = "<<angsol_pix[0]<<" sr ("<<angsol_pix[0]/(4.*M_PI)<<" *4Pi)"
375 <<" in ["<<angsol_pix[1]<<","<<angsol_pix[2]<<"]"<<endl;
376 cout<<" = "<<sr2min2(angsol_pix[0])<<"\'^2"
377 <<" in ["<<sr2min2(angsol_pix[1])<<","<<sr2min2(angsol_pix[2])<<"]"<<endl;
378 cout<<"Total solid angle = "<<angsol_tot[0]<<" sr ("<<angsol_tot[0]/(4.*M_PI)<<" *4Pi)"
379 <<" in ["<<angsol_tot[1]<<","<<angsol_tot[2]<<"]"<<endl;
380
381 // ---
382 // --- Fourier space: k = omega = 2*Pi*Nu = 2*Pi*C/Lambda
383 // ---
384 cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
385 cout<<"Array size: nx = "<<nx<<", ny = "<<ny<<", nz = "<<nz<<endl;
386 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
387 double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
388 double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
389 cout<<"Transverse:"<<endl
390 <<" Resolution: dk_x = "<<dk_x<<" Mpc^-1 (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
391 <<" dk_y = "<<dk_y<<" Mpc^-1 (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl
392 <<" Nyquist: kx = "<<knyq_x<<" Mpc^-1 (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
393 <<" ky = "<<knyq_y<<" Mpc^-1 (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
394 cout<<"Line of sight:"<<endl
395 <<" Resolution: dk_z = "<<dk_z<<" Mpc^-1 (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl
396 <<" Nyquist: kz = "<<knyq_z<<" Mpc^-1 (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
397
398 // ---
399 // --- Variance cosmique
400 // ---
401 // cosmique poisson
402 // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi k^2 dk Vsurvey) * [(1+n*P)/(n*P)]^2
403 // nombre de mode = 1/2 * Vsurvey/(2Pi)^3 * 4Pi*k^2*dk
404 if(kcosm>0. && pkcosm>0.) {
405 double pk = pkcosm*pow(growthfac[0],2.);
406 cout<<"\n>>>>\n>>>> variance cosmique pour k="<<kcosm<<" Mpc^-1, pk(z=0)="
407 <<pkcosm<<", pk(z="<<redshift[0]<<")="<<pk<<"\n>>>>"<<endl;
408 for(int i=0;i<3;i++) { // la correction de variance du au bruit de poisson
409 double v = 1.1; if(i==1) v=1.5; if(i==2) v=2.0;
410 double ngal = 1./(v-1.)/pk;
411 cout<<" pour "<<ngal<<" gal/Mpc^3 on multiplie par "<<v<<" sigma/P"<<endl;
412 }
413
414 cout<<endl;
415 vector<double> dk; if(dkcosm>0.) dk.push_back(dkcosm);
416 dk.push_back(dk_x); dk.push_back(dk_y); dk.push_back(dk_z);
417 for(int i=0;i<(int)dk.size();i++) { // la variance cosmique pure
418 double vcosm = sqrt( 2.*pow(2.*M_PI,3.)/(4.*M_PI*pow(kcosm,2.)*dk[i]*vol_survey) );
419 double nmode = 0.5*vol_survey/pow(2.*M_PI,3.) * 4.*M_PI*pow(kcosm,2.)*dk[i];
420 cout<<" pour dk = "<<dk[i]<<" Mpc^-1: sigma/P = "<<vcosm
421 <<" , Nmode = "<<nmode<<endl;
422 }
423 }
424
425 // ---
426 // --- Masse de HI
427 // ---
428 cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
429 Schechter sch(nstar,mstar,alpha);
430 sch.SetOutValue(1);
431 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
432 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
433 int npt = 10000;
434 double lnx1=log10(1.e+6), lnx2=log10(1.e+13), dlnx=(lnx2-lnx1)/npt;
435 double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
436 cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
437
438 double masshipix = masshimpc3*vol_pixel;
439 double masshitot = masshimpc3*vol_survey;
440 cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
441 <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
442 cout<<"OmegaHI a z=0: "<<masshimpc3/rhoc0<<endl;
443 for(int i=0;i<3;i++)
444 cout<<" a z="<<redshift[i]<<": "<<masshimpc3/rhocz[i]<<endl;
445 if(mhiref<=0.) mhiref = masshipix;
446
447 sch.SetOutValue(0);
448 cout<<"\nsch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
449 cout<<"Galaxy number density:"<<endl;
450 for(double x=lnx1; x<lnx2-0.5; x+=1.) {
451 double n = IntegrateFuncLog(sch,x,lnx2,0.001,dlnx,10.*dlnx,6);
452 cout<<" m>10^"<<x<<" Msol: "<<n<<" /Mpc^3, "<<n*vol_pixel<<" /pixel, "
453 <<n*vol_survey<<" in survey"<<endl;
454 }
455 sch.SetOutValue(1);
456
457 // ---
458 // --- Survey values
459 // ---
460 cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
461 double surfeff = surftot * eta_eff;
462 cout<<"surf_tot="<<surftot<<" m^2, eta="<<eta_eff<<" surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl;
463
464 // Angles solides pour un telescope plein
465 double angSEq[4], angEq[4];
466 for(int i=0;i<4;i++) {
467 double nu = Fr_HyperFin_Par;
468 if(i<3) nu = nuhiz[i];
469 angSEq[i] = AngsolEqTelescope(nu,surftot);
470 angEq[i] = 2.*FrAngSol(angSEq[i]);
471 }
472 cout<<"\nFor a "<<surftot<<" m^2 telescope:"<<endl
473 <<" equivalent Omega = "<<angSEq[0]<<" sr in ["<<angSEq[1]<<","<<angSEq[2]<<"]"<<endl
474 <<" angular diameter = "<<angEq[0]<<" rd = "<<rad2min(angEq[0])
475 <<"\' in ["<<angEq[2]<<","<<angEq[2]<<"]"<<endl;
476 cout<<"At z=0: equivalent Omega = "<<angSEq[3]<<" sr"<<endl
477 <<" angular diameter = "<<angEq[3]<<" rd = "<<rad2min(angEq[3])<<"\'"<<endl;
478
479 // Pour une observation ou le lobe est + petit ou grand que le pixel de simulation
480 // La taille angulaire du lobe change avec la frequence donc avec le redshift
481 // La taille du lobe "lobearea0" est donnee pour une frequence de reference "lobefreq0" en MHz
482 double nlobes[3] = {1.,1.,1.};
483 // Si "lobefreq0" negatif c'est la frequence du centre du cube
484 if(lobefreq0<=0.) lobefreq0 = nuhiz[0]*1.e3;
485 // Si "lobearea0" negatif c'est la taille du pixel ramenee a lobefreq0
486 if(lobearea0<=0.) lobearea0 = sr2min2(angsol_pix[0])*pow(nuhiz[0]*1.e3/lobefreq0,2.);
487 cout<<"\n--- Lobe: angsol="<<lobearea0<<"\'^2 pour "<<lobefreq0<<" MHz"<<endl;
488 double lobearea[3];
489 for(int i=0;i<3;i++) {
490 lobearea[i] = lobearea0*pow(lobefreq0/(nuhiz[0]*1.e3),2.);
491 nlobes[i] = sr2min2(angsol_pix[i])/lobearea[i];
492 }
493 cout<<"Lobe cylindrique area "<<lobearea[0]<<"\'^2 in ["<<lobearea[1]<<","<<lobearea[2]<<"]"<<endl;
494 cout<<"Number of beams in one transversal pixel "<<nlobes[0]<<" in ["<<nlobes[1]<<","<<nlobes[2]<<"]"<<endl;
495
496 // ---
497 // --- Signal analysis
498 // ---
499 // --- Temperature d'antenne Ta et temperature de brillance Tb
500 // La puissance d'une source de brillance I non polarisee est:
501 // P = I * S * Omega * dNu
502 // La puissance recue pour une seule polarisation est:
503 // P = 1/2 * I * S * Omega * dNu
504 // La puissance recue pour un telescope (plein) est
505 // en remplacant Omega = lambda^2/S
506 // P = 1/2 * I * lambda^2 * dNu
507 // En appliquant la loi de Rayleigh: I = 2*k*Tb/lambda^2
508 // on obtient
509 // P = 1/2 * 2*k*Tb * dNu = k * Tb * dNu
510 // La temperature d'antenne est definie comme
511 // P = k * Ta * dNu
512 // Donc pour un ciel de temperature de brillance Tb
513 // et pour une antenne qui mesure une seule polarisation
514 // on a: Ta = Tb
515
516 cout<<"\n>>>>\n>>>> Signal Analysis\n>>>>"<<endl;
517 cout<<"Facteur polar = "<<facpolar<<endl;
518
519 PlanckSpectra planck(T_CMB_Par);
520 planck.SetSpectraApprox(PlanckSpectra::RAYLEIGH); // Rayleigh spectra
521 planck.SetSpectraVar(PlanckSpectra::NU); // frequency
522 planck.SetSpectraPower(PlanckSpectra::POWER); // output en W/....
523 planck.SetSpectraUnit(PlanckSpectra::ANGSFLUX); // radiance W/m^2/Sr/Hz
524
525 // ---
526 // --- Signal HI
527 // ---
528 // Power emitted by HI
529 cout<<"--- Power from HI for M = "<<mhiref<<" Msol"<<endl;
530 cout<<"Flux factor = "<<hifactor<<" at redshift = "<<redshift[0]<<endl;
531 cout<<"Luminosite = "<<hifactor*Msol2LumiHI(mhiref)<<" W"<<endl;
532
533 double fhi_ref[3], sfhi_ref[3];
534 for(int i=0;i<3;i++) {
535 fhi_ref[i] = hifactor*Msol2FluxHI(mhiref,dlum[i]);
536 sfhi_ref[i] = fhi_ref[i] / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
537 }
538 cout<<"FluxHI("<<dlum[0]<<" Mpc) all polar: "<<fhi_ref[0]<<" W/m^2 = "
539 <<fhi_ref[0]/Jansky2Watt_cst<<" Jy.Hz"
540 <<" in ["<<fhi_ref[1]/Jansky2Watt_cst<<","<<fhi_ref[2]/Jansky2Watt_cst<<"]"<<endl;
541 cout<<"If spread over pixel depth ("<<dnuhiz[0]<<" GHz):"<<endl
542 <<" flux density = "<<sfhi_ref[0]*1.e6<<" mu_Jy"
543 <<" in ["<<sfhi_ref[1]<<","<<sfhi_ref[2]<<"]"<<endl;
544
545 // Signal HI
546 double psig_2polar[3], tasig_2polar[3], ssig_2polar[3], isig_2polar[3], tsig_2polar[3];
547 double psig[3], tasig[3], ssig[3], isig[3], tsig[3];
548 double doplarge[3], dzvrot[3];
549
550 for(int i=0;i<3;i++) {
551 psig_2polar[i] = fhi_ref[i] * surfeff;
552 tasig_2polar[i] = psig_2polar[i] / k_Boltzman_Cst / (dnuhiz[i]*1e9);
553 ssig_2polar[i] = psig_2polar[i] / surfeff / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
554 isig_2polar[i] = ssig_2polar[i] / angsol_pix[i];
555 tsig_2polar[i] = isig_2polar[i]*Jansky2Watt_cst
556 / (2.*pow(nuhiz[i]*1e9/(SpeedOfLight_Cst*1000.),2.)*k_Boltzman_Cst);
557 psig[i] = facpolar * psig_2polar[i];
558 tasig[i] = facpolar * tasig_2polar[i];
559 ssig[i] = facpolar * ssig_2polar[i];
560 isig[i] = facpolar * isig_2polar[i];
561 tsig[i] = facpolar * tsig_2polar[i];
562 doplarge[i] = LargeurDoppler(vrot,nuhiz[i]);
563 dzvrot[i] = DzFrV(vrot,redshift[i]);
564 }
565 cout<<"\n--- Signal HI("<<mhiref<<" Msol) for Dnu="<<dnuhiz[0]<<" GHz :"<<endl
566 <<" Observation:"<<endl
567 <<" Power="<<psig[0]<<" W in ["<<psig[1]<<","<<psig[2]<<"]"<<endl
568 <<" Flux density = "<<ssig[0]*1.e6<<" mu_Jy in ["<<ssig[1]*1.e6<<","<<ssig[2]*1.e6<<"]"<<endl
569 <<" Intensity = "<<isig[0]<<" Jy/sr in ["<<isig[1]<<","<<isig[2]<<"]"<<endl
570 <<" Antenna temperature = "<<tasig[0]<<" K in ["<<tasig[1]<<","<<tasig[2]<<"]"<<endl
571 <<" Brightness temperature = "<<tsig[0]<<" K in ["<<tsig[1]<<","<<tsig[2]<<"]"<<endl
572 <<" 2 polars:"<<endl
573 <<" Power="<<psig_2polar[0]<<" W in ["<<psig_2polar[1]<<","<<psig_2polar[2]<<"]"<<endl
574 <<" Flux density = "<<ssig_2polar[0]*1.e6<<" mu_Jy in ["<<ssig_2polar[1]*1.e6<<","<<ssig_2polar[2]*1.e6<<"]"<<endl
575 <<" Intensity = "<<isig_2polar[0]<<" Jy/sr in ["<<isig_2polar[1]<<","<<isig_2polar[2]<<"]"<<endl
576 <<" Antenna temperature = "<<tasig_2polar[0]<<" K in ["<<tasig_2polar[1]<<","<<tasig_2polar[2]<<"]"<<endl
577 <<" Brightness temperature = "<<tsig_2polar[0]<<" K in ["<<tsig_2polar[1]<<","<<tsig_2polar[2]<<"]"<<endl;
578
579 // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
580 cout<<" Doppler for rotation width of "<<vrot<<" km/s"<<endl
581 <<" width="<<doplarge[0]*1.e3<<" MHz in ["<<doplarge[1]*1.e3<<","<<doplarge[2]*1.e3<<"]"<<endl
582 <<" dzvrot= "<<dzvrot[0]<<" in ["<<dzvrot[1]<<","<<dzvrot[2]<<"]"<<endl;
583 if(doplarge[0]>dnuhiz[0])
584 cout<<"Warning: doppler width "<<doplarge[0]<<" GHz > "<<dnuhiz[0]<<" GHz redshift bin width"<<endl;
585
586 // ---
587 // --- Synchrotron (T en -2.7 -> Flux en -0.7 dans l'approximation Rayleigh)
588 // ---
589 double tsynch[3];
590 double psynch_2polar[3], tasynch_2polar[3], ssynch_2polar[3], isynch_2polar[3];
591 double psynch[3], tasynch[3], ssynch[3], isynch[3];
592
593 for(int i=0;i<3;i++) {
594 tsynch[i] = Tsynch408;
595 if(fabs(indnu)>1.e-50) tsynch[i] *= pow(nuhiz[i]/nuhaslam,indnu);
596 planck.SetTemperature(tsynch[i]);
597 psynch_2polar[i] = planck(nuhiz[i]*1.e+9) * surfeff * angsol_pix[i] * (dnuhiz[i]*1e9);
598 tasynch_2polar[i] = psynch_2polar[i] / k_Boltzman_Cst / (dnuhiz[i]*1e9);
599 ssynch_2polar[i] = psynch_2polar[i] / surfeff / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
600 isynch_2polar[i] = ssynch_2polar[i] / angsol_pix[i];
601 psynch[i] = facpolar * psynch_2polar[i];
602 tasynch[i] = facpolar * tasynch_2polar[i];
603 ssynch[i] = facpolar * ssynch_2polar[i];
604 isynch[i] = ssynch[i] / angsol_pix[i];
605 }
606 cout<<"\n--- Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "<<endl
607 <<" tsynch="<<tsynch[0]<<" K ("<<nuhiz[0]<<" GHz) in ["<<tsynch[1]<<","<<tsynch[2]<<"]"<<endl
608 <<" Observation:"<<endl
609 <<" Power="<<psynch[0]<<" W for pixel in ["<<psynch[1]<<","<<psynch[2]<<"]"<<endl
610 <<" Flux density = "<<ssynch[0]<<" Jy for pixel solid angle in ["<<ssynch[1]<<","<<ssynch[2]<<"]"<<endl
611 <<" Intensity = "<<isynch[0]<<" Jy/sr in ["<<isynch[1]<<","<<isynch[2]<<"]"<<endl
612 <<" Antenna temperature = "<<tasynch[0]<<" K in ["<<tasynch[1]<<","<<tasynch[2]<<"]"<<endl
613 <<" 2 polars:"<<endl
614 <<" Power="<<psynch_2polar[0]<<" W for pixel in ["<<psynch_2polar[1]<<","<<psynch_2polar[2]<<"]"<<endl
615 <<" Flux density = "<<ssynch_2polar[0]<<" Jy for pixel solid angle in ["<<ssynch_2polar[1]<<","<<ssynch_2polar[2]<<"]"<<endl
616 <<" Intensity = "<<isynch_2polar[0]<<" Jy/sr in ["<<isynch_2polar[1]<<","<<isynch_2polar[2]<<"]"<<endl
617 <<" Antenna temperature = "<<tasynch_2polar[0]<<" K in ["<<tasynch_2polar[1]<<","<<tasynch_2polar[2]<<"]"<<endl;
618
619 // ---
620 // --- CMB
621 // ---
622 planck.SetTemperature(tcmb);
623 double pcmb_2polar[3], tacmb_2polar[3], scmb_2polar[3], icmb_2polar[3];
624 double pcmb[3], tacmb[3], scmb[3], icmb[3];
625
626 for(int i=0;i<3;i++) {
627 pcmb_2polar[i] = planck(nuhiz[i]*1.e+9) * surfeff * angsol_pix[i] * (dnuhiz[i]*1e9);
628 tacmb_2polar[i] = pcmb_2polar[i] / k_Boltzman_Cst / (dnuhiz[i]*1e9);
629 scmb_2polar[i] = pcmb_2polar[i] / surfeff / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
630 icmb_2polar[i] = scmb_2polar[i] / angsol_pix[i];
631 pcmb[i] = facpolar * pcmb_2polar[i];
632 tacmb[i] = facpolar * tacmb_2polar[i];
633 scmb[i] = facpolar * scmb_2polar[i];
634 icmb[i] = scmb[i] / angsol_pix[i];
635 }
636 cout<<"\n--- CMB: T="<<tcmb<<" K"<<endl
637 <<" Observation:"<<endl
638 <<" Power="<<pcmb[0]<<" W for pixel in ["<<pcmb[1]<<","<<pcmb[2]<<"]"<<endl
639 <<" Flux density = "<<scmb[0]<<" Jy for pixel solid angle in ["<<scmb[1]<<","<<scmb[2]<<"]"<<endl
640 <<" Intensity = "<<icmb[0]<<" Jy/sr in ["<<icmb[1]<<","<<icmb[2]<<"]"<<endl
641 <<" Antenna temperature = "<<tacmb[0]<<" K in ["<<tacmb[1]<<","<<tacmb[2]<<"]"<<endl
642 <<" 2 polars:"<<endl
643 <<" Power="<<pcmb_2polar[0]<<" W for pixel in ["<<pcmb_2polar[1]<<","<<pcmb_2polar[2]<<"]"<<endl
644 <<" Flux density = "<<scmb_2polar[0]<<" Jy for pixel solid angle in ["<<scmb_2polar[1]<<","<<scmb_2polar[2]<<"]"<<endl
645 <<" Intensity = "<<icmb_2polar[0]<<" Jy/sr in ["<<icmb_2polar[1]<<","<<icmb_2polar[2]<<"]"<<endl
646 <<" Antenna temperature = "<<tacmb_2polar[0]<<" K in ["<<tacmb_2polar[1]<<","<<tacmb_2polar[2]<<"]"<<endl;
647
648 // ---
649 // --- AGN
650 // ---
651 double flux_agn = pow(10.,lflux_agn);
652 double mass_agn[3], flux_agn_pix[3], mass_agn_pix[3], lmass_agn_pix[3];
653 for(int i=0;i<3;i++) {
654 mass_agn[i] = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum[i]);
655 flux_agn_pix[i] = flux_agn*(dnuhiz[i]*1e9);
656 mass_agn_pix[i] = FluxHI2Msol(flux_agn_pix[i]*Jansky2Watt_cst,dlum[i]);
657 lmass_agn_pix[i] = log10(mass_agn_pix[i]);
658 }
659 cout<<"\n--- AGN: log10(S_agn)="<<lflux_agn<<" S_agn="<<flux_agn<<" Jy :"<<endl
660 <<" mass_agn = "<<mass_agn[0]<<" equiv. Msol/Hz in ["<<mass_agn[1]<<","<<mass_agn[2]<<"]"<<endl
661 <<" flux_agn_pix = "<<flux_agn_pix[0]<<" 10^-26 W/m^2 in ["<<flux_agn_pix[1]<<","<<flux_agn_pix[2]<<"]"<<endl
662 <<" mass_agn_pix = "<<mass_agn_pix[0]<<" Msol in ["<<mass_agn_pix[1]<<","<<mass_agn_pix[2]<<"]"<<endl
663 <<" log10(mass_agn_pix) = "<<lmass_agn_pix[0]<<" in ["<<lmass_agn_pix[1]<<","<<lmass_agn_pix[2]<<"]"<<endl;
664
665 // =====================================================================
666 // ---
667 // --- Noise analysis
668 // ---
669 // --- Puissance du bruit pour un telescope de surface Ae et de BW dNu
670 // Par definition la puissance du bruit est:
671 // Pb = k * Tsys * dNu (W)
672 // Pour une source (non-polarisee) de densite de flux (totale 2 polars)
673 // St (exprimee en Jy = 10^-26 W/m^2/Hz)
674 // Pt = St * Ae * dNu (puissance totale emise en W pour 2 polars)
675 // P1 = 1/2 * St * Ae * dNu (puissance emise en W pour une polar)
676 // la SEFD (system equivalent flux density en Jy) est definie comme
677 // la densite de flux total (2 polars) "St" d'une source (non-polarisee)
678 // dont la puissance P1 mesuree pour une seule polarisation
679 // serait egale a la puissance du bruit. De P1 = Pb on deduit:
680 // SEFD = 2 * k * Tsys / Ae (en Jy)
681 // la puissance du bruit est: Pb = 1/2 * SEFD * Ae * dNu (en W)
682 // la sensibilite Slim tient compte du temps d'integration et de la BW:
683 // le nombre de mesures independantes est "2*dNu*Tobs" donc
684 // Slim = SEFD / sqrt(2*dNu*Tobs) = 2*k*Tsys/[Ae*sqrt(2*dNu*Tobs) (en Jy)
685 // --- Puissance du bruit pour un interferometre
686 // Ae = surface d'un telescope elementaire
687 // N = nombre de telescopes dans l'interferometre (Atot = N*Ae)
688 // La sensibilite Slim en Jy est:
689 // Slim = 2 * k * Tsys / [ Ae * Sqrt(2*N(N-1)/2 *dnu*Tobs) ]
690 // = 2 * k * Tsys / [ Atot/N * Sqrt(2*N(N-1)/2*dnu*Tobs) ]
691 // = 2 * k * Tsys / [ Atot * Sqrt((N-1)/N *dnu*Tobs) ]
692 // - Interferometre a deux antennes:
693 // Slim = 2 * k * Tsys / [ Atot * Sqrt(1/2 *dnu*Tobs) ]
694 // - Interferometre a N antennes (N grand):
695 // Slim -> 2 * k * Tsys / [ Atot * Sqrt(dnu*Tobs) ]
696 // C'est aussi la formule pour un telescope unique de surface Atot
697 // --- On ne mesure qu'une seule polarisation
698 // Ces formules sont valables si on mesure 1 polarisation:
699 // Slim est la densite de flux total "St" (2 polars) d'une source (non-polarisee)
700 // qui donne la meme puissance que le bruit dans un detecteur qui ne
701 // mesure qu'une seule polarisation:
702 // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
703 // S/N = St / Slim
704 // La puissance de bruit est, par definition:
705 // Pb = 1/2 *Slim*Atot*dNu
706 // = k*Tsys*sqrt(2*dNu/Tobs) pour N=2
707 // = k*Tsys*sqrt(dNu/Tobs) pour N>>grand
708 // La densite de flux d'une source a S/N=1 est:
709 // St = Slim
710 // La puissance d'une source a S/N=1 mesuree par un detecteur
711 // qui ne mesure qu'une polar est:
712 // P1_lim = 1/2 *Slim*Atot*dNu
713 // --- On mesure les 2 polarisations avec deux voies d'electronique distinctes
714 // la puissance du signal mesure est multipliee par 2
715 // la puissance du bruit est multipliee par sqrt(2)
716 // on a donc un gain d'un facteur sqrt(2) sur le rapport S/N
717 // (cela revient d'ailleur a doubler le temps de pose: Tobs -> 2*Tobs)
718 // En notant arbitrairement: Slim' = Slim / sqrt(2)
719 // ou Slim est defini par les formules ci-dessus
720 // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
721 // (S/N)_2 = (S/N)_1 * sqrt(2) = (St / Slim) * sqrt(2) = St / Slim'
722 // La densite de flux d'une source a S/N=1 est:
723 // St = Slim' = Slim / sqrt(2)
724 // La puissance d'une source a S/N=1 cumulee par les 2 detecteurs est:
725 // P_lim = St*Atot*dNu = Slim'*Atot*dNu = 1/sqrt(2) *Slim*Atot*dNu
726 // = P1_lim * sqrt(2)
727 // La puissance de bruit cumulee par les 2 detecteurs est, par definition:
728 // Pb = P_lim = Slim'*Atot*dNu = P1_lim * sqrt(2)
729 // = 2*k*Tsys*sqrt(dNu/Tobs) pour N=2
730 // = k*Tsys*sqrt(2*dNu/Tobs) pour N>>grand
731 // =====================================================================
732
733 cout<<"\n>>>>\n>>>> Noise analysis\n>>>>"<<endl;
734 double psys[3];
735 for(int i=0;i<3;i++) psys[i] = k_Boltzman_Cst * Tsys * (dnuhiz[i]*1.e+9);
736 cout<<"Noise: T="<<Tsys<<" K"<<endl
737 <<" P="<<psys[0]<<" W in ["<<psys[1]<<","<<psys[2]<<"]"<<endl;
738
739 cout<<"...Computation assume that noise dominate the signal."<<endl;
740 if(ya2polar)
741 cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<endl;
742
743
744 //---
745 for(unsigned short it=0;it<2;it++) {
746
747 double fac = 1.;
748 if(it==0) { // Interferometre a 2 telescopes
749 fac = 0.5;
750 cout<<"\n...Observation limits for a 2 telescope interferometer (with complex correlator)"<<endl
751 <<" (sensitivity is given for real or complex correlator output)"<<endl;
752 } else if (it==1) { // Interferometre a N>> telescopes
753 fac = 1.;
754 cout<<"\n...Observation limits for a N (large) telescope interferometer (with complex correlator)"<<endl
755 <<" (weak source limit sensitivity in a synthetised image)"<<endl
756 <<" Also valid for a single dish telescope"<<endl
757 <<" (add factor sqrt(2) if ON-OFF with T_ON=t_OFF)"<<endl;
758 } else continue;
759
760 double slim[3], SsN[3], smass[3], pkbruit[3];
761 for(int i=0;i<3;i++) {
762 slim[i] = 2. * k_Boltzman_Cst * Tsys / surfeff
763 / sqrt(fac*(dnuhiz[i]*1.e+9)*tobs) /Jansky2Watt_cst;
764 if(ya2polar) slim[i] /= sqrt(2.);
765 SsN[i] = ssig_2polar[i] / slim[i]; // independant de angsol_pix*surfeff
766 smass[i] = mhiref / ssig_2polar[i] * slim[i];
767 pkbruit[i] = pow(smass[i]/mhiref,2.)*vol_pixel;
768 }
769 cout<<"for 1 lobe:"<<endl
770 <<" Slim = "<<slim[0]*1.e6<<" mu_Jy in ["<<slim[1]*1.e6<<","<<slim[2]*1.e6<<"]"<<endl
771 <<" S/N = "<<SsN[0]<<" in ["<<SsN[1]<<","<<SsN[2]<<"]"<<endl
772 <<" Mass HI = "<<smass[0]<<" Msol in ["<<smass[1]<<","<<smass[2]<<"]"<<endl
773 <<" Pk = "<<pkbruit[0]<<" Mpc^3 in ["<<pkbruit[1]<<","<<pkbruit[2]<<"]"<<endl;
774
775 double slim_nl[3], SsN_nl[3], smass_nl[3], pkbruit_nl[3];
776 for(int i=0;i<3;i++) {
777 slim_nl[i] = slim[i] * sqrt(nlobes[i]);
778 SsN_nl[i] = ssig_2polar[i] / slim_nl[i];
779 smass_nl[i] = mhiref / ssig_2polar[i] * slim_nl[i];
780 pkbruit_nl[i] = pow(smass_nl[i]/mhiref,2.)*vol_pixel;
781 }
782 cout<<"for "<<nlobes[0]<<" lobes[i] in ["<<nlobes[1]<<","<<nlobes[2]<<"] :"<<endl
783 <<" Slim = "<<slim_nl[0]*1.e6<<" mu_Jy in ["<<slim_nl[1]*1.e6<<","<<slim_nl[2]*1.e6<<"]"<<endl
784 <<" S/N = "<<SsN_nl[0]<<" in ["<<SsN_nl[1]<<","<<SsN_nl[2]<<"]"<<endl
785 <<" Mass HI = "<<smass_nl[0]<<" Msol in ["<<smass_nl[1]<<","<<smass_nl[2]<<"]"<<endl
786 <<" Pk = "<<pkbruit_nl[0]<<" Mpc^3 in ["<<pkbruit_nl[1]<<","<<pkbruit_nl[2]<<"]"<<endl;
787
788 }
789
790 return 0;
791 }
792
793
794//-------------------------------------------------------------------------------------------
795double LargeurDoppler(double v, double nu)
796// largeur doppler pour une vitesse v en km/s et une frequence nu
797{
798 return v / SpeedOfLight_Cst * nu;
799}
800
801double DzFrV(double v, double zred)
802// largeur en redshift pour une vitesse v en km/s au redshift zred
803{
804 return v / SpeedOfLight_Cst * (1. + zred);
805}
806
807double DNuFrDz(double dzred,double nu_at_0,double zred)
808// Largeur DNu pour une largeur en redshift "dzred" au redshift "zred"
809// pour la frequence "nu_at_0" a z=0
810// nu = NuHi(z=0)/(1.+z0)
811// dnu = NuHi(z=0)/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
812// = NuHi(z=0)*dz/[ (1+z0)^2 - (dz/2)^2 ]
813// = NuHi(z=0)*dz/(1.+z0)^2 / [ 1 - [dz/(1+z0)/2)]^2 ]
814// = NuHi(z=0)*dz/(1.+z0)^2 / [1 - dz/(1+z0)/2] / [1 + dz/(1+z0)/2]
815// ~= NuHi(z=0)*dz/(1.+z0)^2 (approx. pour dz<<z0 a l'ordre (dz/z0)^2)
816{
817 double zp1 = 1.+zred;
818 return nu_at_0*dzred/(zp1*zp1)/(1.-dzred/zp1/2.)/(1.+dzred/zp1/2.);
819}
820
821double DzFrDNu(double dnu_at_0,double nu_at_0,double zred)
822// Largeur en redshift au redshift "zred" pour une largeur
823// en frequence "dnu_at_0" a la frequence "nu_at_0" a z=0
824{
825 if(dnu_at_0<=0.) return 0.;
826 double zp1 = 1.+zred;
827 double dnusnu0 = dnu_at_0/nu_at_0;
828 return 2./dnusnu0 * (sqrt(1.+(dnusnu0*zp1)*(dnusnu0*zp1)) - 1.);
829}
830double DzFrDNuApprox(double dnu_at_0,double nu_at_0,double zred)
831// idem DzFrDNu mais on utilise l'approximation: dnu=NuHi(z=0)*dz/(1.+z0)^2
832{
833 double zp1 = 1.+zred;
834 return dnu_at_0/nu_at_0 *(zp1*zp1);
835}
836
837double ZFrLos(double loscom,CosmoCalc& univ)
838// Recherche du redshift correspondant a une distance comobile
839// le long de la ligne de visee egale a "loscom" Mpc
840// et pour un univers "univ"
841{
842 double dz = univ.ZMax()/10.; if(dz<=0.) dz = 0.1;
843 double zmin=0., zmax=0.;
844 while(univ.Dloscom(zmax)<loscom) zmax += dz;
845 if(zmax==0.) return 0.;
846 for(int i=0; i<6; i++) {
847 zmin=zmax-dz; if(zmin<0.) zmin=0.;
848 dz /= 10.;
849 for(double z=zmin; z<zmax+dz; z+=dz) {
850 double d = univ.Dloscom(z);
851 if(d<loscom) continue;
852 zmax = z;
853 //cout<<"ZFrLos: z="<<zmax<<" d="<<d<<" / "<<loscom<<endl;
854 break;
855 }
856 }
857 return zmax;
858}
859
860double AngsolEqTelescope(double nu /* GHz */,double telsurf /* m^2 */)
861/*
862Calcul de l'angle solide (sr) equivalent pour un telescope
863de surface totale "telsurf" (m^2) a la frequence "nu" (GHz)
864- Soit D(t) la figure de diffraction du telescope telle que D(t=0)=1
865 (t = angle depuis l'axe optique)
866 telescope circulaire de diametre D:
867 D(t) = [2J1(Pi*D*t/l)/(Pi*D*t/l)]
868 telescope rectangulaire axb:
869 D(t) = [sin(Pi*a*t/l)/(Pi*a*t/l)]*[sin(Pi*b*t/l)/(Pi*b*t/l)]
870- On cherche l'angle solide equivalent (par ex d'un cylindre de hauteur 1)
871 Int[ D(t)^2 dOmega ] = Int[ D(t)^2 2Pi t dt ] = Lambda^2/S
872- En conclusion, pour un ciel d'intensite uniforme I, on a:
873 P = I * S * Omega * dNu = I * S * (Lambda^2/S) * dNu
874*/
875{
876 double lambda = SpeedOfLight_Cst*1000./(nu*1.e9);
877 return lambda*lambda / telsurf;
878}
Note: See TracBrowser for help on using the repository browser.