Changeset 3193 in Sophya for trunk/Cosmo/SimLSS/cmvdefsurv.cc

Timestamp:

Mar 21, 2007, 7:18:25 PM (19 years ago)

Author:

cmv

Message:

petis changements cmv 21/03/2007

File:

• trunk/Cosmo/SimLSS/cmvdefsurv.cc (modified) (12 diffs)

trunk/Cosmo/SimLSS/cmvdefsurv.cc

@@ -15 +15 @@
 #include "planckspectra.h"
 
+/* --- Check Peterson at al. astro-ph/0606104 v1
+cmvdefsurv -z 0.0025 -x 1 -U 0.75,0.3,0.7,-1,1 -V 300 -O 400000,6000 -A 75 -M 6.156e9 -F 3 -2 1.5
+ --- */
+
 inline double rad2deg(double trad) {return trad/M_PI*180.;}
 inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
@@ -28 +32 @@
       <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
       <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
-      <<" -P : on donne -x -y -z en Mpc aulieu d\'angles et de redshift"<<endl
+      <<" -P : on donne -x -y -z en Mpc au lieu d\'angles et de redshift"<<endl
+      <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
       <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
-      <<" -A Tbruit : temperature de bruit (K)"<<endl
+      <<" -A Tsys : temperature du system (K)"<<endl
       <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
       <<"                   (indnu==0 no evolution with freq.)"<<endl
-      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter"<<endl
+      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
       <<" -F  : HI flux factor to be applied for our redshift"<<endl
-      <<" -1 : on ne mesure qu\'une polarisation"<<endl
+      <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
+      <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
+      <<" -2 : two polarisations measured"<<endl
       <<" redshift : redshift moyen du survey"<<endl
       <<endl;
@@ -45 +52 @@
  // WMAP
  unsigned short flat = 0;
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
- cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="<<or0<<" Ol0="<<ol0<<" w0="<<w0<<endl;
+ double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
  // Schechter
  double h75 = h100 / 0.75;
@@ -61 +67 @@
  int nx,ny,nz;
  double redshift = 1., dred=0.01, redlarg=0.3;
- double tobs = 3600., surfeff = 400000.;
- double Tbruit=75.;
+ double tobs = 6000., surfeff = 400000.;
+ double lobewidth = 1.;  // taille du lobe d'observation en arcmin
+ double Tsys=75.;
  // a 408 MHz (Haslam) + evol index a -2.6
  double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
  double mhiref = -1.; // reference Mass en HI (def integ schechter)
  double hifactor = 1.;
- double facpolar = 1.; // si on ne mesure qu'un polarisation 0.5
+ double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
+ double facpolar = 0.5; // si on ne mesure les 2 polars -> 1.0
 
  // --- Decodage arguments
  char c;
-  while((c = getopt(narg,arg,"hP1x:y:z:A:S:O:M:F:")) != -1) {
+  while((c = getopt(narg,arg,"hP2x:y:z:A:S:O:M:F:V:U:L:")) != -1) {
   switch (c) {
   case 'P' :
@@ -88 +96 @@
     sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
     break;
+  case 'L' :
+    sscanf(optarg,"%lf",&lobewidth);
+    break;
   case 'A' :
-    sscanf(optarg,"%lf",&Tbruit);
+    sscanf(optarg,"%lf",&Tsys);
     break;
   case 'S' :
@@ -100 +111 @@
     sscanf(optarg,"%lf",&hifactor);
     break;
-  case '1' :
-    facpolar = 0.5;
+  case 'V' :
+    sscanf(optarg,"%lf",&vrot);
+    break;
+  case 'U' :
+    sscanf(optarg,"%lf,%lf,%lf,%lf,%u",&h100,&om0,&ol0,&w0,&flat);
+    break;
+  case '2' :
+    facpolar = 1.0;
     break;
   case 'h' :
@@ -107 +124 @@
     usage(); return -1;
   }
- }
+ } 
  if(optind>=narg) {usage(); return-1;}
  sscanf(arg[optind],"%lf",&redshift);
@@ -113 +130 @@
 
  // --- Initialisation de la Cosmologie
+ cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
+     <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
  cout<<"\n--- Cosmology for z = "<<redshift<<endl;
  CosmoCalc univ(flat,true,2.*redshift);
@@ -118 +137 @@
  univ.SetInteg(perc,dzinc,dzmax,glorder);
  univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.Print(0.);
  univ.Print(redshift);
 
@@ -258 +278 @@
  cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
 
+ double slobe = lobewidth/2.35482; // sigma du lobe en arcmin
+ double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
+ double lobearea = M_PI*lobecyl*lobecyl/4.;  // en arcmin^2 (hypothese lobe gaussien)
+ cout<<"\nBeam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
+     <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
+ double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
+ cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
+
  // --- Power emitted by HI
  cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
@@ -267 +295 @@
      <<"      in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
      <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
- double sfhi = fhi / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
  cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
 
  // --- Signal analysis
  cout<<"\n--- Signal analysis"<<endl;
+ cout<<"Facteur polar = "<<facpolar<<endl;
+
  PlanckSpectra planck(T_CMB_Par);
  planck.SetApprox(1);  // Rayleigh spectra
@@ -278 +308 @@
  planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz
 
- double hnu = h_Planck_Cst*(nuhiz*1.e+9);
- cout<<"Facteur polar = "<<facpolar<<",   h.nu="<<hnu<<" J"<<endl;
-
+ // Signal
  double psig = facpolar * fhi * surfeff;
- double esig = psig * tobs;
- double nsig = esig / hnu;
- double tsig = psig / k_Boltzman_Cst / (dnuhiz*1.e+9);
- double ssig = psig / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
- cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W -> E="<<esig<<" J -> "<<nsig<<" ph"<<endl;
+ double tsig = psig / k_Boltzman_Cst / (dnuhiz*1e9);
+ double ssig = psig / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
+ cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W"<<endl;
+ cout<<"     flux density = "<<ssig<<" Jy  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
  cout<<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
- cout<<"     flux density = "<<ssig<<" Jy"<<endl;
-
- double facnoise = facpolar * tobs * surfeff * angsol * (dnuhiz*1.e+9);
-
+
+ // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
+ double doplarge = vrot / SpeedOfLight_Cst * nuhiz;
+ cout<<"     Doppler width="<<doplarge*1.e3<<" MHz  for rotation width of "<<vrot<<" km/s"<<endl;
+ if(doplarge>dnuhiz) {
+   cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
+ }
+
+ // Synchrotron
  double tsynch = Tsynch408;
         if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
  planck.SetTemperature(tsynch);
- double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
- double esynch = psynch * tobs;
- double nsynch = esynch / hnu;
- double ssynch = psynch / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
+ double ssynch = psynch / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
  cout<<"Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
      <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
-     <<"             P="<<psynch<<" W -> E="<<esynch<<" J -> "<<nsynch<<" ph"<<endl;
- cout<<"     flux density = "<<ssynch<<" Jy"<<endl;
-
+     <<"             P="<<psynch<<" W"<<endl;
+ cout<<"             flux density = "<<ssynch<<" Jy"<<endl;
+
+ // CMB
  double tcmb = T_CMB_Par;
  planck.SetTemperature(tcmb);
- double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
- double ecmb = pcmb * tobs;
- double ncmb = ecmb / hnu;
+ double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
  double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
- cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W -> E="<<ecmb<<" J -> "<<ncmb<<" ph"<<endl;
- cout<<"     flux density = "<<scmb<<" Jy"<<endl;
+ cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W"<<endl;
+ cout<<"                       flux density = "<<scmb<<" Jy"<<endl;
 
  // --- Noise analysis
  cout<<"\n--- Noise analysis"<<endl;
- double pbruit = k_Boltzman_Cst * Tbruit * (dnuhiz*1.e+9);
- double ebruit = pbruit * tobs;
- cout<<"Bruit: T="<<Tbruit<<" K, P="<<pbruit<<" W -> E="<<ebruit<<" J"<<endl;
-
- double seq = (1./facpolar) * k_Boltzman_Cst * Tbruit / surfeff /Jansky2Watt_cst;
- cout<<"\nSystem equivalent flux density: "<<seq<<" Jy"<<endl;
-
- double slim = seq / sqrt(2.*(dnuhiz*1.e+9)*tobs);
- cout<<"Observation flux density limit: "<<slim<<" Jy"<<endl;
+ double psys = k_Boltzman_Cst * Tsys * (dnuhiz*1.e+9);
+ cout<<"Noise: T="<<Tsys<<" K, P="<<psys<<" W  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
+
+ double slim = 2. * k_Boltzman_Cst * Tsys / surfeff
+               / sqrt(2.*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
+ cout<<"Observation flux density limit: "<<slim<<" Jy (in 1 lobe)"<<endl;
+
+ slim = slim * sqrt(nlobes);
+ cout<<"Observation flux density limit: "<<slim<<" Jy (in "<<nlobes<<" lobes)"<<endl;
 
  double SsN = ssig/slim;