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

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

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