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

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

modifs calcul bruit, mise a niveau cmv 03/10/2007
File size: 20.4 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 "constcosmo.h"
11#include "cosmocalc.h"
12#include "geneutils.h"
13#include "schechter.h"
14#include "planckspectra.h"
15
16/* --- Check Peterson at al. astro-ph/0606104 v1 (pb facteur sqrt(2) sur S/N !)
17cmvdefsurv -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
18 --- */
19
20inline double rad2deg(double trad) {return trad/M_PI*180.;}
21inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
22inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
23inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
24inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
25inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
26
27void usage(void);
28void usage(void) {
29 cout<<"cmvdefsurv [-r] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
30 <<"----------------"<<endl
31 <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
32 <<" -y adty,atylarg : idem selon Y, si <=0 meme que X"<<endl
33 <<" -z dred,redlarg : resolution et largeur sur la ligne de visee"<<endl
34 <<"-- Unites pour X-Y:"<<endl
35 <<" \'A\' : en angles (pour X-Y) : resolution=ArcMin, largeur=Degre (defaut)"<<endl
36 <<" \'Z\' : en redshift (pour Z) : resolution et largeur en redshift (defaut)"<<endl
37 <<" \'F\' : en frequence (pour Z) : resolution et largeur MHz"<<endl
38 <<" \'M\' : en distance (pour X-Y-Z) : resolution et largeur Mpc"<<endl
39 <<"----------------"<<endl
40 <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
41 <<" -N Tsys : temperature du system (K)"<<endl
42 <<" -L lobewidth,freqlob : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
43 <<" pour la frequence freqlob en MHz"<<endl
44 <<" Si lobewidth<=0 : l'angle solide du lobe = celui du pixel"<<endl
45 <<" Si freqlob<=0 : la frequence de reference est celle du redshift etudie"<<endl
46 <<" Si freqlob absent : la frequence de reference 1.4 GHz"<<endl
47 <<" -2 : two polarisations measured"<<endl
48 <<" -M : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
49 <<" -F : HI flux factor to be applied for our redshift"<<endl
50 <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
51 <<"----------------"<<endl
52 <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
53 <<" (indnu==0 no evolution with freq.)"<<endl
54 <<"----------------"<<endl
55 <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
56 <<"----------------"<<endl
57 <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
58 <<" redshift : redshift moyen du survey"<<endl
59 <<endl;
60}
61
62int main(int narg,char *arg[])
63{
64 // --- Valeurs fixes
65 // WMAP
66 unsigned short flat = 0;
67 double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
68 // Schechter
69 double h75 = h100 / 0.75;
70 double nstar = 0.006*pow(h75,3.); //
71 double mstar = pow(10.,9.8/(h75*h75)); // MSol
72 double alpha = -1.37;
73 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
74
75 // --- Arguments
76 double adtx=0., atxlarg=0., dx=0.,txlarg=0.;
77 int nx=0; char unit_x = 'A';
78 double adty=-1., atylarg=-1., dy=0.,tylarg=0.;
79 int ny=0; char unit_y = 'A';
80 double dred=0., redlarg=0., dz=0.,tzlarg=0.;
81 int nz=0; char unit_z = 'Z';
82 double redshift = 0.;
83 double tobs = 6000., surfeff = 400000.;
84 // taille du lobe d'observation en arcmin pour la frequence
85 double lobewidth0 = -1., lobefreq0 = Fr_HyperFin_Par*1.e3;
86 double Tsys=75.;
87 // a 408 MHz (Haslam) + evol index a -2.6
88 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
89 double mhiref = -1.; // reference Mass en HI (def integ schechter)
90 double hifactor = 1.;
91 double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
92 double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
93 double lflux_agn = -3.;
94
95 // --- Decodage arguments
96 char c;
97 while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
98 switch (c) {
99 case 'x' :
100 sscanf(optarg,"%lf,%lf,%c",&adtx,&atxlarg,&unit_x);
101 break;
102 case 'y' :
103 sscanf(optarg,"%lf,%lf,%c",&adty,&atylarg,&unit_y);
104 break;
105 case 'z' :
106 sscanf(optarg,"%lf,%lf,%c",&dred,&redlarg,&unit_z);
107 break;
108 case 'O' :
109 sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
110 break;
111 case 'L' :
112 sscanf(optarg,"%lf,%lf",&lobewidth0,&lobefreq0);
113 break;
114 case 'N' :
115 sscanf(optarg,"%lf",&Tsys);
116 break;
117 case 'S' :
118 sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
119 break;
120 case 'M' :
121 sscanf(optarg,"%lf",&mhiref);
122 break;
123 case 'F' :
124 sscanf(optarg,"%lf",&hifactor);
125 break;
126 case 'V' :
127 sscanf(optarg,"%lf",&vrot);
128 break;
129 case 'U' :
130 sscanf(optarg,"%lf,%lf,%lf,%lf,%hu",&h100,&om0,&ol0,&w0,&flat);
131 break;
132 case '2' :
133 facpolar = 1.0;
134 break;
135 case 'A' :
136 sscanf(optarg,"%lf",&lflux_agn);
137 break;
138 case 'h' :
139 default :
140 usage(); return -1;
141 }
142 }
143 if(optind>=narg) {usage(); return-1;}
144 sscanf(arg[optind],"%lf",&redshift);
145 if(redshift<=0.) {cout<<"Redshift "<<redshift<<" should be >0"<<endl; return -2;}
146
147 // --- Initialisation de la Cosmologie
148 cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
149 cout<<"h100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
150 <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
151 cout<<"--- Cosmology for z = "<<redshift<<endl;
152 CosmoCalc univ(flat,true,2.*redshift);
153 double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
154 univ.SetInteg(perc,dzinc,dzmax,glorder);
155 univ.SetDynParam(h100,om0,or0,ol0,w0);
156 univ.Print(0.);
157 univ.Print(redshift);
158
159 double dang = univ.Dang(redshift);
160 double dtrcom = univ.Dtrcom(redshift);
161 double dlum = univ.Dlum(redshift);
162 double dloscom = univ.Dloscom(redshift);
163 double dlosdz = univ.Dhubble()/univ.E(redshift);
164 cout<<"dang="<<dang<<" dlum="<<dlum<<" dtrcom="<<dtrcom
165 <<" dloscom="<<dloscom<<" dlosdz="<<dlosdz<<" Mpc"<<endl;
166
167 cout<<"\n1\" -> "<<dang*sec2rad(1.)<<" Mpc = "<<dtrcom*sec2rad(1.)<<" Mpc com"<<endl;
168 cout<<"1\' -> "<<dang*min2rad(1.)<<" Mpc = "<<dtrcom*min2rad(1.)<<" Mpc com"<<endl;
169 cout<<"1d -> "<<dang*deg2rad(1.)<<" Mpc = "<<dtrcom*deg2rad(1.)<<" Mpc com"<<endl;
170
171 cout<<"dz=1 -> "<<dlosdz<<" Mpc com"<<endl;
172
173 cout<<"1 Mpc los com -> dz = "<<1./dlosdz<<endl;
174 cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom)<<"\" = "
175 <<rad2min(1./dtrcom)<<" \' = "<<rad2deg(1./dtrcom)<<" d"<<endl;
176
177 // --- Mise en forme dans les unites appropriees
178 cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
179 if(adty<=0. || atylarg<=0.) {adty=adtx; atylarg=atxlarg; unit_y=unit_x;}
180 cout<<"X values: resolution="<<adtx<<" largeur="<<atxlarg<<" unite="<<unit_x<<endl;
181 if(unit_x == 'A') {
182 nx = int(atxlarg*60./adtx+0.5);
183 adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg);
184 dx = adtx*dtrcom; txlarg = dx*nx;
185 } else if(unit_x == 'M') {
186 nx = int(atxlarg/adtx+0.5);
187 dx = adtx; txlarg = atxlarg;
188 adtx = dx/dtrcom; atxlarg = adtx*nx;
189 } else {
190 cout<<"Unknown unit_x = "<<unit_x<<endl;
191 }
192 cout<<"Y values: resolution="<<adty<<" largeur="<<atylarg<<" unite="<<unit_y<<endl;
193 if(unit_y == 'A') {
194 ny = int(atylarg*60./adty+0.5);
195 adty = min2rad(adty); atylarg = deg2rad(atylarg);
196 dy = adty*dtrcom; tylarg = dy*ny;
197 } else if(unit_y == 'M') {
198 ny = int(atylarg/adty+0.5);
199 dy = adty; tylarg = atylarg;
200 adty = dy/dtrcom; atylarg = adty*ny;
201 } else {
202 cout<<"Unknown unit_y = "<<unit_y<<endl;
203 }
204 cout<<"Z values: resolution="<<dred<<" largeur="<<redlarg<<" unite="<<unit_z<<endl;
205 if(unit_z == 'Z') {
206 nz = int(redlarg/dred+0.5);
207 dz = dred*dlosdz; tzlarg = dz*nz;
208 } else if(unit_z == 'M') {
209 nz = int(redlarg/dred+0.5);
210 dz = dred; tzlarg = redlarg;
211 dred = dz/dlosdz; redlarg = dred*nz;
212 } else if(unit_z == 'F') {
213 nz = int(redlarg/dred+0.5);
214 dred = dred/(Fr_HyperFin_Par*1.e3)*pow(1.+redshift,2.); redlarg = dred*nz;
215 dz = dred*dlosdz; tzlarg = dz*nz;
216 } else {
217 cout<<"Unknown unit_z = "<<unit_z<<endl;
218 }
219
220 double Npix = (double)nx*(double)ny*(double)nz;
221 double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
222 if(redlim[0]<=0.)
223 {cout<<"Lower redshift limit "<<redlim[0]<<" should be >0"<<endl; return -3;}
224 double dtrlim[2] = {univ.Dtrcom(redlim[0]),univ.Dtrcom(redlim[1])};
225 double loslim[2] = {univ.Dloscom(redlim[0]), univ.Dloscom(redlim[1])};
226 double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
227
228 cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
229 <<"dred = "<<dred<<" redlarg = "<<redlarg<<endl
230 <<" dz = "<<dz<<" Mpc redlarg = "<<tzlarg<<" Mpc com, nz = "<<nz<<" pix"<<endl;
231 cout<<"---- Transverse X:"<<endl
232 <<"adtx = "<<rad2min(adtx)<<"\', atxlarg = "<<rad2deg(atxlarg)<<" d"<<endl
233 <<" dx = "<<dx<<" Mpc, txlarg = "<<txlarg<<" Mpc com, nx = "<<nx<<" pix"<<endl;
234 cout<<"---- Transverse Y:"<<endl
235 <<"adty = "<<rad2min(adty)<<"\', atylarg = "<<rad2deg(atylarg)<<" d"<<endl
236 <<" dy = "<<dy<<" Mpc, tylarg = "<<tylarg<<" Mpc com, ny = "<<ny<<" pix"<<endl;
237 cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
238
239 // --- Cosmolographie Transverse
240 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
241 cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
242 <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
243
244 cout<<"... dx = "<<dx<<" Mpc (com), with angle "<<adtx*dtrcom<<endl
245 <<" with angle var_in_z ["<<adtx*dtrlim[0]<<","<<adtx*dtrlim[1]<<"]"<<endl;
246 cout<<"... largx = "<<txlarg<<" Mpc (com), with angle "<<atxlarg*dtrcom<<endl
247 <<" with angle var_in_z ["<<atxlarg*dtrlim[0]<<","<<atxlarg*dtrlim[1]<<"]"<<endl;
248
249 cout<<"... dy = "<<dy<<" Mpc (com), with angle "<<adty*dtrcom<<endl
250 <<" with angle var_in_z ["<<adty*dtrlim[0]<<","<<adty*dtrlim[1]<<"]"<<endl;
251 cout<<"... largy = "<<tylarg<<" Mpc (com), with angle "<<atylarg*dtrcom<<endl
252 <<" with angle var_in_z ["<<atylarg*dtrlim[0]<<","<<atylarg*dtrlim[1]<<"]"<<endl;
253
254 // --- Cosmolographie Line of sight
255 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
256 cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
257 <<loslim[0]<<","<<loslim[1]<<"]"<<endl
258 <<" diff = "
259 <<loslim[1]-loslim[0]<<" Mpc"<<endl;
260
261 cout<<"...dz = "<<dz<<" Mpc (com), with redshift approx "<<dred*dlosdz<<endl;
262 cout<<"...tzlarg = "<<tzlarg<<" Mpc (com), with redshift approx "<<redlarg*dlosdz<<endl;
263
264 // --- Solid Angle & Volume
265 cout<<"\n>>>>\n>>>> Angles solides et Volumes\n>>>>"<<endl;
266 cout<<"--- Solid angle"<<endl;
267 double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
268 cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
269 double angsoltot = AngSol(atxlarg/2.,atylarg/2.,M_PI/2.);
270 cout<<"Total solid angle = "<<angsoltot<<" sr = "<<angsoltot/(4.*M_PI)<<" *4Pi sr"<<endl;
271
272 cout<<"\n--- Volume"<<endl;
273 double dvol = dx*dy*dz;
274 cout<<"Pixel volume comoving = "<<dvol<<" Mpc^3"<<endl;
275 double vol = univ.Vol4Pi(redlim[0],redlim[1])/(4.*M_PI)*angsoltot;
276 cout<<"Volume comoving = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
277 <<"Pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
278
279 // --- Fourier space: k = omega = 2*Pi*Nu
280 cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
281 cout<<"Array size: nx = "<<nx<<", ny = "<<ny<<", nz = "<<nz<<endl;
282 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
283 double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
284 double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
285 cout<<"Resolution: dk_x = "<<dk_x<<" Mpc^-1 (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
286 <<" dk_y = "<<dk_y<<" Mpc^-1 (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl;
287 cout<<"Nyquist: kx = "<<knyq_x<<" Mpc^-1 (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
288 <<" ky = "<<knyq_y<<" Mpc^-1 (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
289 cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1 (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
290 cout<<"Nyquist: kz = "<<knyq_z<<" Mpc^-1 (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
291
292 // --- Masse de HI
293 cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
294 Schechter sch(nstar,mstar,alpha);
295 sch.SetOutValue(1);
296 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
297 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
298 int npt = 10000;
299 double lnx1=log10(1e-6), lnx2=log10(1e+14), dlnx=(lnx2-lnx1)/npt;
300 double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
301 cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
302
303 double masshipix = masshimpc3*dvol;
304 double masshitot = masshimpc3*vol;
305 cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
306 <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
307 if(mhiref<=0.) mhiref = masshipix;
308
309 // --- Survey values
310 cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
311 double unplusz = 1.+redshift;
312 double nuhiz = Fr_HyperFin_Par / unplusz; // GHz
313 // dnu = NuHi/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
314 // = NuHi*dz/(1.+z0)^2 * 1/[1-(dz/(2*(1+z0)))^2]
315 // ~= NuHi*dz/(1.+z0)^2
316 double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
317 / (1.- pow(dred/.2/unplusz,2.));
318 cout<<" surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
319 <<" nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
320 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
321
322 double nlobes = 1.;
323 if(lobewidth0>0.) {
324 double lobewidth = lobewidth0; // ArcMin
325 if(lobefreq0<=0.) lobefreq0 = nuhiz*1.e3; // MHz
326 // La taille angulaire du lobe change avec la frequence donc avec le redshift
327 lobewidth *= lobefreq0/(nuhiz*1.e3);
328 cout<<"\n--- Lobe: width="<<lobewidth0<<" pour "<<lobefreq0<<" MHz"<<endl
329 <<" changed to "<<lobewidth<<" pour "<<nuhiz*1.e3<<" MHz"<<endl;
330 double slobe = lobewidth/2.35482; // sigma du lobe en arcmin
331 double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
332 double lobearea = M_PI*lobecyl*lobecyl/4.; // en arcmin^2 (hypothese lobe gaussien)
333 nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
334 cout<<"Beam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
335 <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
336 cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
337 }
338
339 // --- Power emitted by HI
340 cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
341 cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
342
343 double fhi = hifactor*Msol2FluxHI(mhiref,dlum);
344 cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
345 <<" Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
346 <<" in ["<<hifactor*Msol2FluxHI(mhiref,dlumlim[0])
347 <<","<<hifactor*Msol2FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
348 double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
349 cout<<"If spread over pixel depth ("<<dnuhiz<<" GHz), flux density = "<<sfhi<<" Jy"<<endl;
350
351 // --- Signal analysis
352 cout<<"\n--- Signal analysis"<<endl;
353 cout<<"Facteur polar = "<<facpolar<<endl;
354
355 PlanckSpectra planck(T_CMB_Par);
356 planck.SetApprox(1); // Rayleigh spectra
357 planck.SetVar(0); // frequency
358 planck.SetUnitOut(0); // output en W/....
359 planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz
360
361 // Signal
362 double psig = facpolar * fhi * surfeff;
363 double tsig = psig / k_Boltzman_Cst / (dnuhiz*1e9);
364 double ssig = psig / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
365 cout<<"\nSignal("<<mhiref<<" Msol): P="<<psig<<" W"<<endl;
366 cout<<" flux density = "<<ssig<<" Jy (for Dnu="<<dnuhiz<<" GHz)"<<endl;
367 cout<<" Antenna temperature: tsig="<<tsig<<" K"<<endl;
368
369 // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
370 double doplarge = vrot / SpeedOfLight_Cst * nuhiz;
371 double dzvrot = vrot / SpeedOfLight_Cst * unplusz;
372 cout<<" Doppler width="<<doplarge*1.e3<<" MHz for rotation width of "<<vrot<<" km/s"<<endl
373 <<" dx= "<<dzvrot<<" a z="<<redshift<<endl;
374 if(doplarge>dnuhiz)
375 cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
376
377 // Synchrotron (T en -2.7 -> Flux en -0.7 dans l'approximation Rayleigh)
378 double tsynch = Tsynch408;
379 if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
380 planck.SetTemperature(tsynch);
381 double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
382 double ssynch = psynch / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
383 cout<<"\nSynchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
384 <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
385 <<" P="<<psynch<<" W for pixel"<<endl;
386 cout<<" flux density = "<<ssynch<<" Jy for pixel solid angle"<<endl;
387
388 // CMB
389 double tcmb = T_CMB_Par;
390 planck.SetTemperature(tcmb);
391 double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
392 double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
393 cout<<"\nCMB: T="<<tcmb<<" K -> P="<<pcmb<<" W for pixel"<<endl;
394 cout<<" flux density = "<<scmb<<" Jy for pixel solid angle"<<endl;
395
396 // AGN
397 double flux_agn = pow(10.,lflux_agn);
398 double mass_agn = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum);
399 cout<<"\nAGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="
400 <<flux_agn<<" Jy -> "<<mass_agn<<" equiv. Msol/Hz"<<endl;
401 double flux_agn_pix = flux_agn*(dnuhiz*1e9);
402 double mass_agn_pix = FluxHI2Msol(flux_agn_pix*Jansky2Watt_cst,dlum);
403 double lmass_agn_pix = log10(mass_agn_pix);
404 cout<<"...pixel: f="<<flux_agn_pix<<" 10^-26 W/m^2"
405 <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
406
407 // ---
408 // --- Noise analysis
409 // ---
410 // Ae = surface d'un telescope elementaire
411 // N = nombre de telescopes dans l'interferometre (Atot = N*Ae)
412 // Slim = 2 * k * Tsys / [ Ae * Sqrt(2*N(N-1)/2 *dnu*Tobs) ]
413 // = 2 * k * Tsys / [ Atot/N * Sqrt(2*N(N-1)/2*dnu*Tobs) ]
414 // = 2 * k * Tsys / [ Atot * Sqrt((N-1)/N *dnu*Tobs) ]
415 // Interferometre a deux antennes:
416 // Slim = 2 * k * Tsys / [ Atot * Sqrt(1/2 *dnu*Tobs) ]
417 // Interferometre a N antennes (N grand):
418 // Slim -> 2 * k * Tsys / [ Atot * Sqrt(dnu*Tobs) ]
419 // C'est aussi la formule pour un telescope unique
420 // -
421 // Ces formules sont valables si on mesure 1 polarisation.
422 // Si on mesure 2 polarisations:
423 // le signal (ssig) est multiplie par 2 (facpolar=1 au lieu de 0.5)
424 // mais on le lit avec 2 chaines electroniques differentes
425 // -> le gain n'est que de sqrt(2)
426 cout<<"\n---\n--- Noise analysis \n---"<<endl;
427 double psys = k_Boltzman_Cst * Tsys * (dnuhiz*1.e+9);
428 cout<<"Noise: T="<<Tsys<<" K, P="<<psys<<" W (for Dnu="<<dnuhiz<<" GHz)"<<endl;
429
430 double slim,slim_nl,SsN,SsN_nl,smass,smass_nl;
431 cout<<"...Computation assume that noise dominate the signal."<<endl;
432
433 double facpolarbruit = 1.;
434 if(fabs(facpolar-1.)<0.01) {
435 facpolarbruit = sqrt(2.);
436 cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<endl;
437 }
438
439 //---
440 cout<<"\n...Observation limits for a 2 telescope interferometer (with complex correlator)"<<endl
441 <<" (sensitivity is given for real or complex correlator output)"<<endl;
442 slim = 2. * k_Boltzman_Cst * Tsys / surfeff
443 / sqrt(0.5*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
444 slim *= facpolarbruit;
445 SsN = ssig/slim;
446 smass = mhiref/ssig*slim;
447 cout<<"for 1 lobe:"<<endl
448 <<" Flux = "<<slim<<" Jy"<<endl
449 <<" S/N = "<<SsN<<endl
450 <<" Mass HI = "<<smass<<" Msol"<<endl;
451 slim_nl = slim * sqrt(nlobes);
452 SsN_nl = ssig/slim_nl;
453 smass_nl = mhiref/ssig*slim_nl;
454 cout<<"for "<<nlobes<<" lobes:"<<endl
455 <<" Flux = "<<slim_nl<<" Jy"<<endl
456 <<" S/N = "<<SsN_nl<<endl
457 <<" Mass HI = "<<smass_nl<<" Msol"<<endl;
458
459 //---
460 cout<<"\n...Observation limits for a N (large) telescope interferometer (with complex correlator)"<<endl
461 <<" (weak source limit sensitivity in a synthetised image)"<<endl
462 <<" (Also valid for 1 single antenna telescope)"<<endl;
463 slim = 2. * k_Boltzman_Cst * Tsys / surfeff
464 / sqrt((dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
465 slim *= facpolarbruit;
466 SsN = ssig/slim;
467 smass = mhiref/ssig*slim;
468 cout<<"for 1 lobe:"<<endl
469 <<" Flux = "<<slim<<" Jy"<<endl
470 <<" S/N = "<<SsN<<endl
471 <<" Mass HI = "<<smass<<" Msol"<<endl;
472 slim_nl = slim * sqrt(nlobes);
473 SsN_nl = ssig/slim_nl;
474 smass_nl = mhiref/ssig*slim_nl;
475 cout<<"for "<<nlobes<<" lobes:"<<endl
476 <<" Flux = "<<slim_nl<<" Jy"<<endl
477 <<" S/N = "<<SsN_nl<<endl
478 <<" Mass HI = "<<smass_nl<<" Msol"<<endl;
479
480 return 0;
481}
Note: See TracBrowser for help on using the repository browser.