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

Last change on this file since 3892 was 3805, checked in by cmv, 15 years ago

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