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

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

ajour calcul OmegaHI equiv, cmv 16/10/2007
File size: 25.9 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
23inline double rad2deg(double trad) {return trad/M_PI*180.;}
24inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
25inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
26inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
27inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
28inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
29
30void usage(void);
31void usage(void) {
32 cout<<"cmvdefsurv [options] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
33 <<"----------------"<<endl
34 <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
35 <<" -y adty,atylarg : idem selon Y, si <=0 meme que X"<<endl
36 <<" -z dred,redlarg : resolution et largeur sur la ligne de visee"<<endl
37 <<"-- Unites pour X-Y:"<<endl
38 <<" \'A\' : en angles (pour X-Y) : resolution=ArcMin, largeur=Degre (defaut)"<<endl
39 <<" \'Z\' : en redshift (pour Z) : resolution et largeur en redshift (defaut)"<<endl
40 <<" \'F\' : en frequence (pour Z) : resolution et largeur MHz"<<endl
41 <<" \'M\' : en distance (pour X-Y-Z) : resolution et largeur Mpc"<<endl
42 <<"----------------"<<endl
43 <<" -K k,dk,pk : k(Mpc^-1) dk(Mpc^-1) pk(a z=0 en Mpc^-3) pour estimer la variance cosmique"<<endl
44 <<"----------------"<<endl
45 <<" -O surf,tobs,eta : surface effective (m^2), temps d\'observation (s), efficacite d\'ouverture"<<endl
46 <<" -N Tsys : temperature du system (K)"<<endl
47 <<" -L lobewidth,freqlob : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
48 <<" pour la frequence freqlob en MHz"<<endl
49 <<" Si lobewidth<=0 : l'angle solide du lobe = celui du pixel"<<endl
50 <<" Si freqlob<=0 : la frequence de reference est celle du redshift etudie"<<endl
51 <<" Si freqlob absent : la frequence de reference 1.4 GHz"<<endl
52 <<" -2 : two polarisations measured"<<endl
53 <<" -M : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
54 <<" -F : HI flux factor to be applied for our redshift"<<endl
55 <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
56 <<"----------------"<<endl
57 <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
58 <<" (indnu==0 no evolution with freq.)"<<endl
59 <<"----------------"<<endl
60 <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
61 <<"----------------"<<endl
62 <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
63 <<" redshift : redshift moyen du survey"<<endl
64 <<endl;
65}
66
67int main(int narg,char *arg[])
68{
69 // --- Valeurs fixes
70 // WMAP
71 unsigned short flat = 0;
72 double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
73 // Schechter
74 double h75 = h100 / 0.75;
75 double nstar = 0.006*pow(h75,3.); //
76 double mstar = pow(10.,9.8); // MSol
77 double alpha = -1.37;
78 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
79
80 // --- Arguments
81 double adtx=0., atxlarg=0., dx=0.,txlarg=0.;
82 int nx=0; char unit_x = 'A';
83 double adty=-1., atylarg=-1., dy=0.,tylarg=0.;
84 int ny=0; char unit_y = 'A';
85 double dred=0., redlarg=0., dz=0.,tzlarg=0.;
86 int nz=0; char unit_z = 'Z';
87 double redshift = 0.;
88 double tobs = 6000., surftot = 400000., eta_eff = 1.;
89 // variance cosmique (default = standard SDSSII)
90 double kcosm = 0.05, dkcosm = -1., pkcosm = 40000.;
91 // taille du lobe d'observation en arcmin pour la frequence
92 double lobewidth0 = -1., lobefreq0 = Fr_HyperFin_Par*1.e3;
93 double Tsys=75.;
94 // a 408 MHz (Haslam) + evol index a -2.6
95 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
96 double mhiref = -1.; // reference Mass en HI (def integ schechter)
97 double hifactor = 1.;
98 double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
99 bool ya2polar = false;
100 double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
101 double lflux_agn = -3.;
102
103 // --- Decodage arguments
104 char c;
105 while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:K:")) != -1) {
106 switch (c) {
107 case 'x' :
108 sscanf(optarg,"%lf,%lf,%c",&adtx,&atxlarg,&unit_x);
109 break;
110 case 'y' :
111 sscanf(optarg,"%lf,%lf,%c",&adty,&atylarg,&unit_y);
112 break;
113 case 'z' :
114 sscanf(optarg,"%lf,%lf,%c",&dred,&redlarg,&unit_z);
115 break;
116 case 'O' :
117 sscanf(optarg,"%lf,%lf,%lf",&surftot,&tobs,&eta_eff);
118 break;
119 case 'L' :
120 sscanf(optarg,"%lf,%lf",&lobewidth0,&lobefreq0);
121 break;
122 case 'N' :
123 sscanf(optarg,"%lf",&Tsys);
124 break;
125 case 'S' :
126 sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
127 break;
128 case 'M' :
129 sscanf(optarg,"%lf",&mhiref);
130 break;
131 case 'F' :
132 sscanf(optarg,"%lf",&hifactor);
133 break;
134 case 'V' :
135 sscanf(optarg,"%lf",&vrot);
136 break;
137 case 'U' :
138 sscanf(optarg,"%lf,%lf,%lf,%lf,%hu",&h100,&om0,&ol0,&w0,&flat);
139 break;
140 case '2' :
141 ya2polar = true;
142 facpolar = 1.0;
143 break;
144 case 'A' :
145 sscanf(optarg,"%lf",&lflux_agn);
146 break;
147 case 'K' :
148 sscanf(optarg,"%lf,%lf,%lf",&kcosm,&dkcosm,&pkcosm);
149 break;
150 case 'h' :
151 default :
152 usage(); return -1;
153 }
154 }
155 if(optind>=narg) {usage(); return-1;}
156 sscanf(arg[optind],"%lf",&redshift);
157 if(redshift<=0.) {cout<<"Redshift "<<redshift<<" should be >0"<<endl; return -2;}
158
159 // --- Initialisation de la Cosmologie
160 cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
161 cout<<"h100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
162 <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
163 cout<<"--- Cosmology for z = "<<redshift<<endl;
164 CosmoCalc univ(flat,true,2.*redshift);
165 double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
166 univ.SetInteg(perc,dzinc,dzmax,glorder);
167 univ.SetDynParam(h100,om0,or0,ol0,w0);
168 univ.Print(0.);
169 univ.Print(redshift);
170 GrowthFactor growth(om0,ol0);
171 double growthfac = growth(redshift);
172 cout<<"Facteur de croissance lineaire = "<<growthfac
173 <<" ^2 , ( 1/(1+z)="<<1./(1.+redshift)<<" )"<<endl;
174 double rhoc0 = univ.Rhoc(0.)*GCm3toMsolMpc3_Cst;
175 double rhocz = univ.Rhoc(redshift)*GCm3toMsolMpc3_Cst;
176 cout<<"Rho_c a z=0: "<<rhoc0<<", a z="<<redshift<<": "<<rhocz<<" Msol/Mpc^3"<<endl;
177
178 double dang = univ.Dang(redshift);
179 double dtrcom = univ.Dtrcom(redshift);
180 double dlum = univ.Dlum(redshift);
181 double dloscom = univ.Dloscom(redshift);
182 double dlosdz = univ.Dhubble()/univ.E(redshift);
183 cout<<"dang="<<dang<<" dlum="<<dlum<<" dtrcom="<<dtrcom
184 <<" dloscom="<<dloscom<<" dlosdz="<<dlosdz<<" Mpc"<<endl;
185
186 cout<<"\n1\" -> "<<dang*sec2rad(1.)<<" Mpc = "<<dtrcom*sec2rad(1.)<<" Mpc com"<<endl;
187 cout<<"1\' -> "<<dang*min2rad(1.)<<" Mpc = "<<dtrcom*min2rad(1.)<<" Mpc com"<<endl;
188 cout<<"1d -> "<<dang*deg2rad(1.)<<" Mpc = "<<dtrcom*deg2rad(1.)<<" Mpc com"<<endl;
189
190 cout<<"dz=1 -> "<<dlosdz<<" Mpc com"<<endl;
191
192 cout<<"1 Mpc los com -> dz = "<<1./dlosdz<<endl;
193 cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom)<<"\" = "
194 <<rad2min(1./dtrcom)<<" \' = "<<rad2deg(1./dtrcom)<<" d"<<endl;
195
196 // --- Mise en forme dans les unites appropriees
197 cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
198 if(adty<=0. || atylarg<=0.) {adty=adtx; atylarg=atxlarg; unit_y=unit_x;}
199 cout<<"X values: resolution="<<adtx<<" largeur="<<atxlarg<<" unite="<<unit_x<<endl;
200 if(unit_x == 'A') {
201 nx = int(atxlarg*60./adtx+0.5);
202 adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg);
203 dx = adtx*dtrcom; txlarg = dx*nx;
204 } else if(unit_x == 'M') {
205 nx = int(atxlarg/adtx+0.5);
206 dx = adtx; txlarg = atxlarg;
207 adtx = dx/dtrcom; atxlarg = adtx*nx;
208 } else {
209 cout<<"Unknown unit_x = "<<unit_x<<endl;
210 }
211 cout<<"Y values: resolution="<<adty<<" largeur="<<atylarg<<" unite="<<unit_y<<endl;
212 if(unit_y == 'A') {
213 ny = int(atylarg*60./adty+0.5);
214 adty = min2rad(adty); atylarg = deg2rad(atylarg);
215 dy = adty*dtrcom; tylarg = dy*ny;
216 } else if(unit_y == 'M') {
217 ny = int(atylarg/adty+0.5);
218 dy = adty; tylarg = atylarg;
219 adty = dy/dtrcom; atylarg = adty*ny;
220 } else {
221 cout<<"Unknown unit_y = "<<unit_y<<endl;
222 }
223 cout<<"Z values: resolution="<<dred<<" largeur="<<redlarg<<" unite="<<unit_z<<endl;
224 if(unit_z == 'Z') {
225 nz = int(redlarg/dred+0.5);
226 dz = dred*dlosdz; tzlarg = dz*nz;
227 } else if(unit_z == 'M') {
228 nz = int(redlarg/dred+0.5);
229 dz = dred; tzlarg = redlarg;
230 dred = dz/dlosdz; redlarg = dred*nz;
231 } else if(unit_z == 'F') {
232 nz = int(redlarg/dred+0.5);
233 dred = dred/(Fr_HyperFin_Par*1.e3)*pow(1.+redshift,2.); redlarg = dred*nz;
234 dz = dred*dlosdz; tzlarg = dz*nz;
235 } else {
236 cout<<"Unknown unit_z = "<<unit_z<<endl;
237 }
238
239 double Npix = (double)nx*(double)ny*(double)nz;
240 double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
241 if(redlim[0]<=0.)
242 {cout<<"Lower redshift limit "<<redlim[0]<<" should be >0"<<endl; return -3;}
243 double dtrlim[2] = {univ.Dtrcom(redlim[0]),univ.Dtrcom(redlim[1])};
244 double loslim[2] = {univ.Dloscom(redlim[0]), univ.Dloscom(redlim[1])};
245 double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
246
247 cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
248 <<"dred = "<<dred<<" redlarg = "<<redlarg<<endl
249 <<" dz = "<<dz<<" Mpc redlarg = "<<tzlarg<<" Mpc com, nz = "<<nz<<" pix"<<endl;
250 cout<<"---- Transverse X:"<<endl
251 <<"adtx = "<<rad2min(adtx)<<"\', atxlarg = "<<rad2deg(atxlarg)<<" d"<<endl
252 <<" dx = "<<dx<<" Mpc, txlarg = "<<txlarg<<" Mpc com, nx = "<<nx<<" pix"<<endl;
253 cout<<"---- Transverse Y:"<<endl
254 <<"adty = "<<rad2min(adty)<<"\', atylarg = "<<rad2deg(atylarg)<<" d"<<endl
255 <<" dy = "<<dy<<" Mpc, tylarg = "<<tylarg<<" Mpc com, ny = "<<ny<<" pix"<<endl;
256 cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
257 cout<<" Volume pixel = "<<dx*dy*dz<<" Mpc^3"<<endl;
258 cout<<" Volume total = "<<Npix*dx*dy*dz<<" Mpc^3"<<endl;
259
260 // --- Cosmolographie Transverse
261 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
262 cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
263 <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
264
265 cout<<"... dx = "<<dx<<" Mpc (com), with angle "<<adtx*dtrcom<<endl
266 <<" with angle var_in_z ["<<adtx*dtrlim[0]<<","<<adtx*dtrlim[1]<<"]"<<endl;
267 cout<<"... largx = "<<txlarg<<" Mpc (com), with angle "<<atxlarg*dtrcom<<endl
268 <<" with angle var_in_z ["<<atxlarg*dtrlim[0]<<","<<atxlarg*dtrlim[1]<<"]"<<endl;
269
270 cout<<"... dy = "<<dy<<" Mpc (com), with angle "<<adty*dtrcom<<endl
271 <<" with angle var_in_z ["<<adty*dtrlim[0]<<","<<adty*dtrlim[1]<<"]"<<endl;
272 cout<<"... largy = "<<tylarg<<" Mpc (com), with angle "<<atylarg*dtrcom<<endl
273 <<" with angle var_in_z ["<<atylarg*dtrlim[0]<<","<<atylarg*dtrlim[1]<<"]"<<endl;
274
275 // --- Cosmolographie Line of sight
276 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
277 cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
278 <<loslim[0]<<","<<loslim[1]<<"]"<<endl
279 <<" diff = "
280 <<loslim[1]-loslim[0]<<" Mpc"<<endl;
281
282 cout<<"...dz = "<<dz<<" Mpc (com), with redshift approx "<<dred*dlosdz<<endl;
283 cout<<"...tzlarg = "<<tzlarg<<" Mpc (com), with redshift approx "<<redlarg*dlosdz<<endl;
284
285 // --- Solid Angle & Volume
286 cout<<"\n>>>>\n>>>> Angles solides et Volumes\n>>>>"<<endl;
287 cout<<"--- Solid angle"<<endl;
288 double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
289 cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
290 double angsoltot = AngSol(atxlarg/2.,atylarg/2.,M_PI/2.);
291 cout<<"Total solid angle = "<<angsoltot<<" sr = "<<angsoltot/(4.*M_PI)<<" *4Pi sr"<<endl;
292
293 cout<<"\n--- Volume"<<endl;
294 double dvol = dx*dy*dz;
295 cout<<"Pixel volume comoving = "<<dvol<<" Mpc^3"<<endl;
296 double vol = univ.Vol4Pi(redlim[0],redlim[1])/(4.*M_PI)*angsoltot;
297 cout<<"Volume comoving = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
298 <<"Pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
299
300 // --- Fourier space: k = omega = 2*Pi*Nu
301 cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
302 cout<<"Array size: nx = "<<nx<<", ny = "<<ny<<", nz = "<<nz<<endl;
303 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
304 double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
305 double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
306 cout<<"Resolution: dk_x = "<<dk_x<<" Mpc^-1 (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
307 <<" dk_y = "<<dk_y<<" Mpc^-1 (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl;
308 cout<<"Nyquist: kx = "<<knyq_x<<" Mpc^-1 (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
309 <<" ky = "<<knyq_y<<" Mpc^-1 (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
310 cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1 (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
311 cout<<"Nyquist: kz = "<<knyq_z<<" Mpc^-1 (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
312
313 // --- Variance cosmique
314 // cosmique poisson
315 // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi k^2 dk Vsurvey) * [(1+n*P)/(n*P)]^2
316 // nombre de mode = 1/2 * Vsurvey/(2Pi)^3 * 4Pi*k^2*dk
317 if(kcosm>0. && pkcosm>0.) {
318 double pk = pkcosm*pow(growthfac,2.);
319 cout<<"\n>>>>\n>>>> variance cosmique pour k="<<kcosm<<" Mpc^-1, pk(z=0)="
320 <<pkcosm<<", pk(z="<<redshift<<")="<<pk<<"\n>>>>"<<endl;
321 for(int i=0;i<3;i++) { // la correction de variance du au bruit de poisson
322 double v = 1.1; if(i==1) v=1.5; if(i==2) v=2.0;
323 double ngal = 1./(v-1.)/pk;
324 cout<<" pour "<<ngal<<" gal/Mpc^3 on multiplie par "<<v<<" sigma/P"<<endl;
325 }
326 vector<double> dk; if(dkcosm>0.) dk.push_back(dkcosm);
327 dk.push_back(dk_x); dk.push_back(dk_y); dk.push_back(dk_z);
328 for(int i=0;i<(int)dk.size();i++) { // la variance cosmique pure
329 double vcosm = sqrt( 2.*pow(2.*M_PI,3.)/(4.*M_PI*pow(kcosm,2.)*dk[i]*vol) );
330 double nmode = 0.5*vol/pow(2.*M_PI,3.) * 4.*M_PI*pow(kcosm,2.)*dk[i];
331 cout<<" pour dk = "<<dk[i]<<" Mpc^-1: sigma/P = "<<vcosm
332 <<" , Nmode = "<<nmode<<endl;
333 }
334 }
335
336 // --- Masse de HI
337 cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
338 Schechter sch(nstar,mstar,alpha);
339 sch.SetOutValue(1);
340 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
341 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
342 int npt = 10000;
343 double lnx1=log10(1.e+6), lnx2=log10(1.e+13), dlnx=(lnx2-lnx1)/npt;
344 double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
345 cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
346
347 double masshipix = masshimpc3*dvol;
348 double masshitot = masshimpc3*vol;
349 cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
350 <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
351 cout<<"OmegaHI a z=0: "<<masshimpc3/rhoc0
352 <<", a z="<<redshift<<": "<<masshimpc3/rhocz<<endl;
353 if(mhiref<=0.) mhiref = masshipix;
354
355 sch.SetOutValue(0);
356 cout<<"\nsch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
357 cout<<"Galaxy number density:"<<endl;
358 for(double x=lnx1; x<lnx2-0.5; x+=1.) {
359 double n = IntegrateFuncLog(sch,x,lnx2,0.001,dlnx,10.*dlnx,6);
360 cout<<" m>"<<x<<" Msol: "<<n<<" /Mpc^3, "<<n*dvol<<" /pixel, "
361 <<n*vol<<" in survey"<<endl;
362 }
363 sch.SetOutValue(1);
364
365
366 // --- Survey values
367 cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
368 double surfeff = surftot * eta_eff;
369 double unplusz = 1.+redshift;
370 double nuhiz = Fr_HyperFin_Par / unplusz; // GHz
371 // dnu = NuHi/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
372 // = NuHi*dz/(1.+z0)^2 * 1/[1-(dz/(2*(1+z0)))^2]
373 // ~= NuHi*dz/(1.+z0)^2
374 double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
375 / (1.- pow(dred/.2/unplusz,2.));
376 cout<<" surf="<<surftot<<" m^2, eta="<<eta_eff<<" surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
377 <<" nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
378 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
379
380 double nlobes = 1.;
381 if(lobewidth0>0.) {
382 double lobewidth = lobewidth0; // ArcMin
383 if(lobefreq0<=0.) lobefreq0 = nuhiz*1.e3; // MHz
384 // La taille angulaire du lobe change avec la frequence donc avec le redshift
385 lobewidth *= lobefreq0/(nuhiz*1.e3);
386 cout<<"\n--- Lobe: width="<<lobewidth0<<" pour "<<lobefreq0<<" MHz"<<endl
387 <<" changed to "<<lobewidth<<" pour "<<nuhiz*1.e3<<" MHz"<<endl;
388 double slobe = lobewidth/2.35482; // sigma du lobe en arcmin
389 double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
390 double lobearea = M_PI*lobecyl*lobecyl/4.; // en arcmin^2 (hypothese lobe gaussien)
391 nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
392 cout<<"Beam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
393 <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
394 cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
395 }
396
397 // --- Power emitted by HI
398 cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
399 cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
400
401 double fhi = hifactor*Msol2FluxHI(mhiref,dlum);
402 cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
403 <<" Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
404 <<" in ["<<hifactor*Msol2FluxHI(mhiref,dlumlim[0])
405 <<","<<hifactor*Msol2FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
406 double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
407 cout<<"If spread over pixel depth ("<<dnuhiz<<" GHz), flux density = "<<sfhi*1.e6<<" mu_Jy"<<endl;
408
409 // --- Signal analysis
410 cout<<"\n--- Signal analysis"<<endl;
411 cout<<"Facteur polar = "<<facpolar<<endl;
412
413 PlanckSpectra planck(T_CMB_Par);
414 planck.SetSpectraApprox(PlanckSpectra::RAYLEIGH); // Rayleigh spectra
415 planck.SetSpectraVar(PlanckSpectra::NU); // frequency
416 planck.SetSpectraPower(PlanckSpectra::POWER); // output en W/....
417 planck.SetSpectraUnit(PlanckSpectra::ANGSFLUX); // radiance W/m^2/Sr/Hz
418
419 // Signal
420 double psig_2polar = fhi * surfeff;
421 double tsig_2polar = psig_2polar / k_Boltzman_Cst / (dnuhiz*1e9);
422 double ssig_2polar = psig_2polar / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
423 double psig = facpolar * psig_2polar;
424 double tsig = facpolar * tsig_2polar;
425 double ssig = facpolar * ssig_2polar;
426 cout<<"\nSignal("<<mhiref<<" Msol):"<<endl
427 <<" P="<<psig<<" W"<<endl
428 <<" flux density = "<<ssig*1.e6<<" mu_Jy (for Dnu="<<dnuhiz<<" GHz)"<<endl
429 <<" ("<<ssig_2polar<<" mu_Jy for 2 polars)"<<endl
430 <<" Antenna temperature: tsig="<<tsig<<" K"<<endl;
431
432 // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
433 double doplarge = vrot / SpeedOfLight_Cst * nuhiz;
434 double dzvrot = vrot / SpeedOfLight_Cst * unplusz;
435 cout<<" Doppler width="<<doplarge*1.e3<<" MHz for rotation width of "<<vrot<<" km/s"<<endl
436 <<" dx= "<<dzvrot<<" a z="<<redshift<<endl;
437 if(doplarge>dnuhiz)
438 cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
439
440 // Synchrotron (T en -2.7 -> Flux en -0.7 dans l'approximation Rayleigh)
441 double tsynch = Tsynch408;
442 if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
443 planck.SetTemperature(tsynch);
444 double psynch_2polar = planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
445 double ssynch_2polar = psynch_2polar / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
446 double psynch = facpolar * psynch_2polar;
447 double ssynch = facpolar * ssynch_2polar;
448 cout<<"\nSynchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
449 <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
450 <<" P="<<psynch<<" W for pixel"<<endl
451 <<" flux density = "<<ssynch<<" Jy for pixel solid angle"<<endl;
452
453 // CMB
454 double tcmb = T_CMB_Par;
455 planck.SetTemperature(tcmb);
456 double pcmb_2polar = planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
457 double scmb_2polar = pcmb_2polar / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
458 double pcmb = facpolar * pcmb_2polar;
459 double scmb = facpolar * scmb_2polar;
460 cout<<"\nCMB: T="<<tcmb<<" K"<<endl
461 <<" P="<<pcmb<<" W for pixel"<<endl
462 <<" flux density = "<<scmb<<" Jy for pixel solid angle"<<endl;
463
464 // AGN
465 double flux_agn = pow(10.,lflux_agn);
466 double mass_agn = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum);
467 cout<<"\nAGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="
468 <<flux_agn<<" Jy -> "<<mass_agn<<" equiv. Msol/Hz"<<endl;
469 double flux_agn_pix = flux_agn*(dnuhiz*1e9);
470 double mass_agn_pix = FluxHI2Msol(flux_agn_pix*Jansky2Watt_cst,dlum);
471 double lmass_agn_pix = log10(mass_agn_pix);
472 cout<<"...pixel: f="<<flux_agn_pix<<" 10^-26 W/m^2"
473 <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
474
475 // =====================================================================
476 // ---
477 // --- Noise analysis
478 // ---
479 // --- Puissance du bruit pour un telescope de surface Ae et de BW dNu
480 // Par definition la puissance du bruit est:
481 // Pb = k * Tsys * dNu (W)
482 // Pour une source (non-polarisee) de densite de flux (totale 2 polars)
483 // St (exprimee en Jy = 10^-26 W/m^2/Hz)
484 // Pt = St * Ae * dNu (puissance totale emise en W pour 2 polars)
485 // P1 = 1/2 * St * Ae * dNu (puissance emise en W pour une polar)
486 // la SEFD (system equivalent flux density en Jy) est definie comme
487 // la densite de flux total (2 polars) "St" d'une source (non-polarisee)
488 // dont la puissance P1 mesuree pour une seule polarisation
489 // serait egale a la puissance du bruit. De P1 = Pb on deduit:
490 // SEFD = 2 * k * Tsys / Ae (en Jy)
491 // la puissance du bruit est: Pb = 1/2 * SEFD * Ae * dNu (en W)
492 // la sensibilite Slim tient compte du temps d'integration et de la BW:
493 // le nombre de mesures independantes est "2*dNu*Tobs" donc
494 // Slim = SEFD / sqrt(2*dNu*Tobs) = 2*k*Tsys/[Ae*sqrt(2*dNu*Tobs) (en Jy)
495 // --- Puissance du bruit pour un interferometre
496 // Ae = surface d'un telescope elementaire
497 // N = nombre de telescopes dans l'interferometre (Atot = N*Ae)
498 // La sensibilite Slim en Jy est:
499 // Slim = 2 * k * Tsys / [ Ae * Sqrt(2*N(N-1)/2 *dnu*Tobs) ]
500 // = 2 * k * Tsys / [ Atot/N * Sqrt(2*N(N-1)/2*dnu*Tobs) ]
501 // = 2 * k * Tsys / [ Atot * Sqrt((N-1)/N *dnu*Tobs) ]
502 // - Interferometre a deux antennes:
503 // Slim = 2 * k * Tsys / [ Atot * Sqrt(1/2 *dnu*Tobs) ]
504 // - Interferometre a N antennes (N grand):
505 // Slim -> 2 * k * Tsys / [ Atot * Sqrt(dnu*Tobs) ]
506 // C'est aussi la formule pour un telescope unique de surface Atot
507 // --- On ne mesure qu'une seule polarisation
508 // Ces formules sont valables si on mesure 1 polarisation:
509 // Slim est la densite de flux total "St" (2 polars) d'une source (non-polarisee)
510 // qui donne la meme puissance que le bruit dans un detecteur qui ne
511 // mesure qu'une seule polarisation:
512 // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
513 // S/N = St / Slim
514 // La puissance de bruit est, par definition:
515 // Pb = 1/2 *Slim*Atot*dNu
516 // = k*Tsys*sqrt(2*dNu/Tobs) pour N=2
517 // = k*Tsys*sqrt(dNu/Tobs) pour N>>grand
518 // La densite de flux d'une source a S/N=1 est:
519 // St = Slim
520 // La puissance d'une source a S/N=1 mesuree par un detecteur
521 // qui ne mesure qu'une polar est:
522 // P1_lim = 1/2 *Slim*Atot*dNu
523 // --- On mesure les 2 polarisations avec deux voies d'electronique distinctes
524 // la puissance du signal mesure est multipliee par 2
525 // la puissance du bruit est multipliee par sqrt(2)
526 // on a donc un gain d'un facteur sqrt(2) sur le rapport S/N
527 // (cela revient d'ailleur a doubler le temps de pose: Tobs -> 2*Tobs)
528 // En notant arbitrairement: Slim' = Slim / sqrt(2)
529 // ou Slim est defini par les formules ci-dessus
530 // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
531 // (S/N)_2 = (S/N)_1 * sqrt(2) = (St / Slim) * sqrt(2) = St / Slim'
532 // La densite de flux d'une source a S/N=1 est:
533 // St = Slim' = Slim / sqrt(2)
534 // La puissance d'une source a S/N=1 cumulee par les 2 detecteurs est:
535 // P_lim = St*Atot*dNu = Slim'*Atot*dNu = 1/sqrt(2) *Slim*Atot*dNu
536 // = P1_lim * sqrt(2)
537 // La puissance de bruit cumulee par les 2 detecteurs est, par definition:
538 // Pb = P_lim = Slim'*Atot*dNu = P1_lim * sqrt(2)
539 // = 2*k*Tsys*sqrt(dNu/Tobs) pour N=2
540 // = k*Tsys*sqrt(2*dNu/Tobs) pour N>>grand
541 // =====================================================================
542
543 cout<<"\n---\n--- Noise analysis \n---"<<endl;
544 double psys = k_Boltzman_Cst * Tsys * (dnuhiz*1.e+9);
545 cout<<"Noise: T="<<Tsys<<" K, P="<<psys<<" W (for Dnu="<<dnuhiz<<" GHz)"<<endl;
546
547 cout<<"...Computation assume that noise dominate the signal."<<endl;
548 if(ya2polar)
549 cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<endl;
550
551 double slim,slim_nl,SsN,SsN_nl,smass,smass_nl;
552
553 //---
554 for(unsigned short it=0;it<2;it++) {
555
556 double fac = 1.;
557 if(it==0) { // Interferometre a 2 telescopes
558 fac = 0.5;
559 cout<<"\n...Observation limits for a 2 telescope interferometer (with complex correlator)"<<endl
560 <<" (sensitivity is given for real or complex correlator output)"<<endl;
561 } else if (it==1) { // Interferometre a N>> telescopes
562 fac = 1.;
563 cout<<"\n...Observation limits for a N (large) telescope interferometer (with complex correlator)"<<endl
564 <<" (weak source limit sensitivity in a synthetised image)"<<endl
565 <<" Also valid for a single dish telescope."<<endl;
566 } else continue;
567
568 slim = 2. * k_Boltzman_Cst * Tsys / surfeff
569 / sqrt(fac*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
570 if(ya2polar) slim /= sqrt(2.);
571 SsN = ssig_2polar / slim;
572 smass = mhiref / ssig_2polar * slim;
573 cout<<"for 1 lobe:"<<endl
574 <<" Slim = "<<slim*1.e6<<" mu_Jy"<<endl
575 <<" S/N = "<<SsN<<endl
576 <<" Mass HI = "<<smass<<" Msol"<<endl;
577
578 slim_nl = slim * sqrt(nlobes);
579 SsN_nl = ssig_2polar / slim_nl;
580 smass_nl = mhiref / ssig_2polar * slim_nl;
581 cout<<"for "<<nlobes<<" lobes:"<<endl
582 <<" Flux = "<<slim_nl*1.e6<<" mu_Jy"<<endl
583 <<" S/N = "<<SsN_nl<<endl
584 <<" Mass HI = "<<smass_nl<<" Msol"<<endl;
585
586 }
587
588 return 0;
589}
Note: See TracBrowser for help on using the repository browser.