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

Last change on this file since 3115 was 3115, checked in by ansari, 19 years ago

Creation initiale du groupe Cosmo avec le repertoire SimLSS de
simulation de distribution de masse 3D des galaxies par CMV+Rz
18/12/2006

File size: 13.5 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 "integfunc.h"
14#include "schechter.h"
15#include "planckspectra.h"
16
17inline double rad2deg(double trad) {return trad/M_PI*180.;}
18inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
19inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
20inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
21inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
22inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
23
24void usage(void);
25void usage(void) {
26 cout<<"cmvdefsurv [-r] -x adtx,atxlarg [-y adty,atylarg] -z dred,redlarg redshift"<<endl
27 <<" -x adtx,atxlarg : resolution en Theta_x (arcmin), largeur (degre)"<<endl
28 <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
29 <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
30 <<" -P : on donne -x -y -z en Mpc aulieu d\'angles et de redshift"<<endl
31 <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
32 <<" -A Tbruit : temperature de bruit (K)"<<endl
33 <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
34 <<" (indnu==0 no evolution with freq.)"<<endl
35 <<" -M : masse de HI de reference (MSol), si <=0 mean schechter"<<endl
36 <<" -F : HI flux factor to be applied for our redshift"<<endl
37 <<" -1 : on ne mesure qu\'une polarisation"<<endl
38 <<" redshift : redshift moyen du survey"<<endl
39 <<endl;
40}
41
42int main(int narg,char *arg[])
43{
44 // --- Valeurs fixes
45 // WMAP
46 unsigned short flat = 0;
47 double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
48 cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="<<or0<<" Ol0="<<ol0<<" w0="<<w0<<endl;
49 // Schechter
50 double h75 = h100 / 0.75;
51 double nstar = 0.006*pow(h75,3.); //
52 double mstar = pow(10.,9.8/(h75*h75)); // MSol
53 double alpha = -1.37;
54 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
55
56
57 // --- Arguments
58 bool inmpc = false;
59 double adtx=1., atxlarg=90., adty=-1., atylarg=-1.;
60 double dx=1.,txlarg=1000., dy=-1.,tylarg=1000., dz=1.,tzlarg=100.;
61 int nx,ny,nz;
62 double redshift = 1., dred=0.01, redlarg=0.3;
63 double tobs = 3600., surfeff = 400000.;
64 double Tbruit=75.;
65 // a 408 MHz (Haslam) + evol index a -2.6
66 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
67 double mhiref = -1.; // reference Mass en HI (def integ schechter)
68 double hifactor = 1.;
69 double facpolar = 1.; // si on ne mesure qu'un polarisation 0.5
70
71 // --- Decodage arguments
72 char c;
73 while((c = getopt(narg,arg,"hP1x:y:z:A:S:O:M:F:")) != -1) {
74 switch (c) {
75 case 'P' :
76 inmpc = true;
77 break;
78 case 'x' :
79 sscanf(optarg,"%lf,%lf",&adtx,&atxlarg);
80 break;
81 case 'y' :
82 sscanf(optarg,"%lf,%lf",&adty,&atylarg);
83 break;
84 case 'z' :
85 sscanf(optarg,"%lf,%lf",&dred,&redlarg);
86 break;
87 case 'O' :
88 sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
89 break;
90 case 'A' :
91 sscanf(optarg,"%lf",&Tbruit);
92 break;
93 case 'S' :
94 sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
95 break;
96 case 'M' :
97 sscanf(optarg,"%lf",&mhiref);
98 break;
99 case 'F' :
100 sscanf(optarg,"%lf",&hifactor);
101 break;
102 case '1' :
103 facpolar = 0.5;
104 break;
105 case 'h' :
106 default :
107 usage(); return -1;
108 }
109 }
110 if(optind>=narg) {usage(); return-1;}
111 sscanf(arg[optind],"%lf",&redshift);
112 if(redshift<=0.) {cout<<"Redshift "<<redshift<<" should be >0"<<endl; return -2;}
113
114 // --- Initialisation de la Cosmologie
115 cout<<"\n--- Cosmology for z = "<<redshift<<endl;
116 CosmoCalc univ(flat,true,2.*redshift);
117 double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
118 univ.SetInteg(perc,dzinc,dzmax,glorder);
119 univ.SetDynParam(h100,om0,or0,ol0,w0);
120 univ.Print(redshift);
121
122 double dang = univ.Dang(redshift);
123 double dtrcom = univ.Dtrcom(redshift);
124 double dlum = univ.Dlum(redshift);
125 double dloscom = univ.Dloscom(redshift);
126 double dlosdz = univ.Dhubble()/univ.E(redshift);
127 cout<<"dang="<<dang<<" dlum="<<dlum<<" dtrcom="<<dtrcom
128 <<" dloscom="<<dloscom<<" dlosdz="<<dlosdz<<" Mpc"<<endl;
129
130 cout<<"\n1\" -> "<<dang*sec2rad(1.)<<" Mpc = "<<dtrcom*sec2rad(1.)<<" Mpc com"<<endl;
131 cout<<"1\' -> "<<dang*min2rad(1.)<<" Mpc = "<<dtrcom*min2rad(1.)<<" Mpc com"<<endl;
132 cout<<"1d -> "<<dang*deg2rad(1.)<<" Mpc = "<<dtrcom*deg2rad(1.)<<" Mpc com"<<endl;
133
134 cout<<"dz=1 -> "<<dlosdz<<" Mpc com"<<endl;
135
136 cout<<"1 Mpc los com -> dz = "<<1./dlosdz<<endl;
137 cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom)<<"\" = "
138 <<rad2min(1./dtrcom)<<" \' = "<<rad2deg(1./dtrcom)<<" d"<<endl;
139
140 // --- Mise en forme dans les unites appropriees
141 if(adty<=0.) adty=adtx;
142 if(atylarg<=0.) atylarg=atxlarg;
143 if(inmpc) {
144 dx = adtx; txlarg = atxlarg; nx = int(txlarg/dx+0.5);
145 dy = adty; txlarg = atxlarg; ny = int(tylarg/dy+0.5);
146 dz = dred; tzlarg = redlarg; nz = int(tzlarg/dz+0.5);
147 adtx = dx/dtrcom; atxlarg = adtx*nx;
148 adty = dy/dtrcom; atylarg = adty*ny;
149 dred = dz/dlosdz; redlarg = dred*nz;
150 } else {
151 adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg); nx = int(atxlarg/adtx+0.5);
152 adty = min2rad(adty); atylarg = deg2rad(atylarg); ny = int(atylarg/adty+0.5);
153 nz = int(redlarg/dred+0.5);
154 dx = adtx*dtrcom; txlarg = dx*nx;
155 dy = adty*dtrcom; tylarg = dy*ny;
156 dz = dred*dlosdz; tzlarg = dz*nz;
157 }
158 double Npix = (double)nx*(double)ny*(double)nz;
159
160 double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
161 if(redlim[0]<=0.)
162 {cout<<"Lower redshift limit "<<redlim[0]<<" should be >0"<<endl; return -3;}
163 double dtrlim[2] = {univ.Dtrcom(redlim[0]),univ.Dtrcom(redlim[1])};
164 double loslim[2] = {univ.Dloscom(redlim[0]), univ.Dloscom(redlim[1])};
165 double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
166
167 cout<<"\n---- Type de donnees: inmpc = "<<inmpc<<endl;
168 cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
169 <<"dred = "<<dred<<" redlarg = "<<redlarg<<endl
170 <<" dz = "<<dz<<" Mpc redlarg = "<<tzlarg<<" Mpc, nz = "<<nz<<" pix"<<endl;
171 cout<<"---- Transverse X:"<<endl
172 <<"adtx = "<<rad2min(adtx)<<"\', atxlarg = "<<rad2deg(atxlarg)<<" d"<<endl
173 <<" dx = "<<dx<<" Mpc, txlarg = "<<txlarg<<" Mpc, nx = "<<nx<<" pix"<<endl;
174 cout<<"---- Transverse Y:"<<endl
175 <<"adty = "<<rad2min(adty)<<"\', atylarg = "<<rad2deg(atylarg)<<" d"<<endl
176 <<" dy = "<<dy<<" Mpc, tylarg = "<<tylarg<<" Mpc, ny = "<<ny<<" pix"<<endl;
177 cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
178
179 // --- Cosmolographie Transverse
180 cout<<"\n--- Transverse"<<endl;
181 cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
182 <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
183
184 cout<<"... dx = "<<dx<<" Mpc (com), with angle "<<adtx*dtrcom<<endl
185 <<" with angle var_in_z ["<<adtx*dtrlim[0]<<","<<adtx*dtrlim[1]<<"]"<<endl;
186 cout<<"... largx = "<<txlarg<<" Mpc (com), with angle "<<atxlarg*dtrcom<<endl
187 <<" with angle var_in_z ["<<atxlarg*dtrlim[0]<<","<<atxlarg*dtrlim[1]<<"]"<<endl;
188
189 cout<<"... dy = "<<dy<<" Mpc (com), with angle "<<adty*dtrcom<<endl
190 <<" with angle var_in_z ["<<adty*dtrlim[0]<<","<<adty*dtrlim[1]<<"]"<<endl;
191 cout<<"... largy = "<<tylarg<<" Mpc (com), with angle "<<atylarg*dtrcom<<endl
192 <<" with angle var_in_z ["<<atylarg*dtrlim[0]<<","<<atylarg*dtrlim[1]<<"]"<<endl;
193
194 // --- Cosmolographie Line of sight
195 cout<<"\n--- Line of Sight"<<endl;
196 cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
197 <<loslim[0]<<","<<loslim[1]<<"]"<<endl
198 <<" diff = "
199 <<loslim[1]-loslim[0]<<" Mpc"<<endl;
200
201 cout<<"...dz = "<<dz<<" Mpc (com), with redshift approx "<<dred*dlosdz<<endl;
202 cout<<"...tzlarg = "<<tzlarg<<" Mpc (com), with redshift approx "<<redlarg*dlosdz<<endl;
203
204 // --- Solid Angle & Volume
205 cout<<"\n--- Solid angle"<<endl;
206 double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
207 cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
208 double angsoltot = AngSol(atxlarg/2.,atylarg/2.,M_PI/2.);
209 cout<<"Total solid angle = "<<angsoltot<<" sr = "<<angsoltot/(4.*M_PI)<<" *4Pi sr"<<endl;
210
211 cout<<"\n--- Volume"<<endl;
212 double dvol = dx*dy*dz;
213 cout<<"Pixel volume comoving = "<<dvol<<" Mpc^3"<<endl;
214 double vol = univ.Vol4Pi(redlim[0],redlim[1])/(4.*M_PI)*angsoltot;
215 cout<<"Volume comoving = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
216 <<"Pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
217
218 // --- Fourier space: k = omega = 2*Pi*Nu
219 cout<<"\n--- Fourier space"<<endl;
220 cout<<"Array size: nx = "<<nx<<", ny = "<<ny<<", nz = "<<nz<<endl;
221 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
222 double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
223 double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
224 cout<<"Resolution: dk_x = "<<dk_x<<" Mpc^-1 (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
225 <<" dk_y = "<<dk_y<<" Mpc^-1 (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl;
226 cout<<"Nyquist: kx = "<<knyq_x<<" Mpc^-1 (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
227 <<" ky = "<<knyq_y<<" Mpc^-1 (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
228 cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1 (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
229 cout<<"Nyquist: kz = "<<knyq_z<<" Mpc^-1 (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
230
231 // --- Masse de HI
232 cout<<"\n--- Mass HI"<<endl;
233 Schechter sch(nstar,mstar,alpha);
234 sch.SetOutValue(1);
235 cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
236 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
237 int npt = 10000;
238 double lnx1=log10(1e-6), lnx2=log10(1e+14), dlnx=(lnx2-lnx1)/npt;
239 double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
240 cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
241
242 double masshipix = masshimpc3*dvol;
243 double masshitot = masshimpc3*vol;
244 cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
245 <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
246 if(mhiref<=0.) mhiref = masshipix;
247
248 // --- Survey values
249 double unplusz = 1.+redshift;
250 double nuhiz = Fr_HyperFin_Par / unplusz; // GHz
251 // dnu = NuHi/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
252 // = NuHi*dz/(1.+z0)^2 * 1/[1-(dz/(2*(1+z0)))^2]
253 double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
254 / (1.- (dred/.2/unplusz)*(dred/.2/unplusz));
255 cout<<"\n--- Observation:"<<endl
256 <<" surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
257 <<" nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
258 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
259
260 // --- Power emitted by HI
261 cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
262 cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
263
264 double fhi = hifactor*FluxHI(mhiref,dlum);
265 cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
266 <<" Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
267 <<" in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
268 <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
269 double sfhi = fhi / (dnuhiz*1.e+9) / Jansky2Watt_cst;
270 cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
271
272 // --- Signal analysis
273 cout<<"\n--- Signal analysis"<<endl;
274 PlanckSpectra planck(T_CMB_Par);
275 planck.SetApprox(1); // Rayleigh spectra
276 planck.SetVar(0); // frequency
277 planck.SetUnitOut(0); // output en W/....
278 planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz
279
280 double hnu = h_Planck_Cst*(nuhiz*1.e+9);
281 cout<<"Facteur polar = "<<facpolar<<", h.nu="<<hnu<<" J"<<endl;
282
283 double psig = facpolar * fhi * surfeff;
284 double esig = psig * tobs;
285 double nsig = esig / hnu;
286 double tsig = psig / k_Boltzman_Cst / (dnuhiz*1.e+9);
287 double ssig = psig / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
288 cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W -> E="<<esig<<" J -> "<<nsig<<" ph"<<endl;
289 cout<<" Antenna temperature: tsig="<<tsig<<" K"<<endl;
290 cout<<" flux density = "<<ssig<<" Jy"<<endl;
291
292 double facnoise = facpolar * tobs * surfeff * angsol * (dnuhiz*1.e+9);
293
294 double tsynch = Tsynch408;
295 if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
296 planck.SetTemperature(tsynch);
297 double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
298 double esynch = psynch * tobs;
299 double nsynch = esynch / hnu;
300 double ssynch = psynch / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
301 cout<<"Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
302 <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
303 <<" P="<<psynch<<" W -> E="<<esynch<<" J -> "<<nsynch<<" ph"<<endl;
304 cout<<" flux density = "<<ssynch<<" Jy"<<endl;
305
306 double tcmb = T_CMB_Par;
307 planck.SetTemperature(tcmb);
308 double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
309 double ecmb = pcmb * tobs;
310 double ncmb = ecmb / hnu;
311 double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
312 cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W -> E="<<ecmb<<" J -> "<<ncmb<<" ph"<<endl;
313 cout<<" flux density = "<<scmb<<" Jy"<<endl;
314
315 // --- Noise analysis
316 cout<<"\n--- Noise analysis"<<endl;
317 double pbruit = k_Boltzman_Cst * Tbruit * (dnuhiz*1.e+9);
318 double ebruit = pbruit * tobs;
319 cout<<"Bruit: T="<<Tbruit<<" K, P="<<pbruit<<" W -> E="<<ebruit<<" J"<<endl;
320
321 double seq = (1./facpolar) * k_Boltzman_Cst * Tbruit / surfeff /Jansky2Watt_cst;
322 cout<<"\nSystem equivalent flux density: "<<seq<<" Jy"<<endl;
323
324 double slim = seq / sqrt(2.*(dnuhiz*1.e+9)*tobs);
325 cout<<"Observation flux density limit: "<<slim<<" Jy"<<endl;
326
327 double SsN = ssig/slim;
328 cout<<"\nSignal to noise ratio = "<<SsN<<endl;
329
330 double smass = mhiref/ssig*slim;
331 cout<<"\nSigma noise equivalent = "<<smass<<" Msol"<<endl;
332
333 return 0;
334}
Note: See TracBrowser for help on using the repository browser.