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

Timestamp:

Oct 3, 2007, 7:22:13 PM (18 years ago)

Author:

cmv

Message:

modifs calcul bruit, mise a niveau cmv 03/10/2007

File:

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

trunk/Cosmo/SimLSS/cmvdefsurv.cc

@@ -14 +14 @@
 #include "planckspectra.h"
 
-/* --- Check Peterson at al. astro-ph/0606104 v1
+/* --- Check Peterson at al. astro-ph/0606104 v1  (pb facteur sqrt(2) sur S/N !)
 cmvdefsurv -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
  --- */
@@ -31 +31 @@
       <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
       <<" -y adty,atylarg : idem selon Y, si <=0 meme que X"<<endl
-      <<" -z dred,redlarg : resolution et largeursur la ligne de visee"<<endl
+      <<" -z dred,redlarg : resolution et largeur sur la ligne de visee"<<endl
       <<"-- Unites pour X-Y:"<<endl
       <<"   \'A\' : en angles (pour X-Y)     : resolution=ArcMin, largeur=Degre (defaut)"<<endl
@@ -44 +44 @@
       <<"            Si lobewidth<=0 : l'angle solide du lobe = celui du pixel"<<endl
       <<"            Si freqlob<=0 : la frequence de reference est celle du redshift etudie"<<endl
+      <<"            Si freqlob absent : la frequence de reference 1.4 GHz"<<endl
       <<" -2 : two polarisations measured"<<endl
       <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
@@ -89 +90 @@
  double hifactor = 1.;
  double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
- double facpolar = 0.5; // si on ne mesure les 2 polars -> 1.0
+ double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
  double lflux_agn = -3.;
 
@@ -319 +320 @@
  cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
 
- double lobewidth = lobewidth0;   // ArcMin
- if(lobefreq0<=0.) lobefreq0 = nuhiz*1.e3; // MHz
- // La taille angulaire du lobe change avec la frequence donc avec le redshift
- lobewidth *= lobefreq0/(nuhiz*1.e3);
- cout<<"\n--- Lobe: width="<<lobewidth0<<" pour "<<lobefreq0<<" MHz"<<endl
-     <<"          changed to "<<lobewidth<<" pour "<<nuhiz*1.e3<<" MHz"<<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)
- double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
- if(lobewidth<=0.) nlobes = 1.;
- cout<<"Beam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
-     <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
- cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
+ double nlobes = 1.;
+ if(lobewidth0>0.) {
+   double lobewidth = lobewidth0;   // ArcMin
+   if(lobefreq0<=0.) lobefreq0 = nuhiz*1.e3; // MHz
+   // La taille angulaire du lobe change avec la frequence donc avec le redshift
+   lobewidth *= lobefreq0/(nuhiz*1.e3);
+   cout<<"\n--- Lobe: width="<<lobewidth0<<" pour "<<lobefreq0<<" MHz"<<endl
+       <<"          changed to "<<lobewidth<<" pour "<<nuhiz*1.e3<<" MHz"<<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)
+   nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
+   cout<<"Beam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
+       <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
+   cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
+ }
 
  // --- Power emitted by HI
@@ -360 +363 @@
  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<<"\nSignal("<<mhiref<<" Msol): P="<<psig<<" W"<<endl;
  cout<<"     flux density = "<<ssig<<" Jy  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
  cout<<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
@@ -378 +381 @@
  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), "
+ cout<<"\nSynchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
      <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
      <<"             P="<<psynch<<" W for pixel"<<endl;
@@ -388 +391 @@
  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 for pixel"<<endl;
+ cout<<"\nCMB: T="<<tcmb<<" K -> P="<<pcmb<<" W for pixel"<<endl;
  cout<<"                       flux density = "<<scmb<<" Jy for pixel solid angle"<<endl;
 
@@ -394 +397 @@
  double flux_agn = pow(10.,lflux_agn);
  double mass_agn = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum);
- cout<<"AGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="
+ cout<<"\nAGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="
      <<flux_agn<<" Jy -> "<<mass_agn<<" equiv. Msol/Hz"<<endl;
  double flux_agn_pix = flux_agn*(dnuhiz*1e9);
@@ -402 +405 @@
      <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
 
+ // ---
  // --- Noise analysis
- cout<<"\n--- Noise analysis"<<endl;
+ // ---
+ // Ae = surface d'un telescope elementaire
+ // N  = nombre de telescopes dans l'interferometre (Atot = N*Ae)
+ // Slim = 2 * k * Tsys / [ Ae * Sqrt(2*N(N-1)/2 *dnu*Tobs) ]
+ //      = 2 * k * Tsys / [ Atot/N * Sqrt(2*N(N-1)/2*dnu*Tobs) ]
+ //      = 2 * k * Tsys / [ Atot * Sqrt((N-1)/N *dnu*Tobs) ]
+ // Interferometre a deux antennes:
+ // Slim =  2 * k * Tsys / [ Atot * Sqrt(1/2 *dnu*Tobs) ]
+ // Interferometre a N antennes (N grand):
+ // Slim ->  2 * k * Tsys / [ Atot * Sqrt(dnu*Tobs) ]
+ //      C'est aussi la formule pour un telescope unique
+ // -
+ // Ces formules sont valables si on mesure 1 polarisation.
+ // Si on mesure 2 polarisations:
+ //   le signal (ssig) est multiplie par 2 (facpolar=1 au lieu de 0.5)
+ //   mais on le lit avec 2 chaines electroniques differentes
+ //   -> le gain n'est que de sqrt(2)
+ cout<<"\n---\n--- Noise analysis \n---"<<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;
-
- double slim_nl = slim * sqrt(nlobes);
- cout<<"Observation flux density limit: "<<slim_nl<<" Jy (in "<<nlobes<<" lobes)"<<endl;
-
- double SsN = ssig/slim;
- cout<<"\nSignal to noise ratio = "<<SsN<<" (1 lobe)"<<endl;
- double SsN_nl = ssig/slim_nl;
- cout<<"\nSignal to noise ratio = "<<SsN_nl<<" ("<<nlobes<<" lobes)"<<endl;
-
- double smass = mhiref/ssig*slim;
- cout<<"\nSigma noise equivalent = "<<smass<<" Msol (1 lobe)"<<endl;
- double smass_nl = mhiref/ssig*slim_nl;
- cout<<"\nSigma noise equivalent = "<<smass_nl<<" Msol ("<<nlobes<<" lobes)"<<endl;
+ double slim,slim_nl,SsN,SsN_nl,smass,smass_nl;
+ cout<<"...Computation assume that noise dominate the signal."<<endl;
+ 
+ double facpolarbruit = 1.;
+ if(fabs(facpolar-1.)<0.01) {
+   facpolarbruit = sqrt(2.);
+   cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<endl;
+ }
+
+ //---
+ cout<<"\n...Observation limits for a 2 telescope interferometer (with complex correlator)"<<endl
+     <<"   (sensitivity is given for real or complex correlator output)"<<endl;
+ slim = 2. * k_Boltzman_Cst * Tsys / surfeff
+               / sqrt(0.5*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
+ slim *= facpolarbruit;
+ SsN = ssig/slim;
+ smass = mhiref/ssig*slim;
+ cout<<"for 1 lobe:"<<endl
+     <<"   Flux    = "<<slim<<" Jy"<<endl
+     <<"   S/N     = "<<SsN<<endl
+     <<"   Mass HI = "<<smass<<" Msol"<<endl;
+ slim_nl = slim * sqrt(nlobes);
+ SsN_nl = ssig/slim_nl;
+ smass_nl = mhiref/ssig*slim_nl;
+ cout<<"for "<<nlobes<<" lobes:"<<endl
+     <<"   Flux    = "<<slim_nl<<" Jy"<<endl
+     <<"   S/N     = "<<SsN_nl<<endl
+     <<"   Mass HI = "<<smass_nl<<" Msol"<<endl;
+
+ //---
+ cout<<"\n...Observation limits for a N (large) telescope interferometer (with complex correlator)"<<endl
+     <<"     (weak source limit sensitivity in a synthetised image)"<<endl
+     <<"     (Also valid for 1 single antenna telescope)"<<endl;
+ slim = 2. * k_Boltzman_Cst * Tsys / surfeff
+               / sqrt((dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
+ slim *= facpolarbruit;
+ SsN = ssig/slim;
+ smass = mhiref/ssig*slim;
+ cout<<"for 1 lobe:"<<endl
+     <<"   Flux    = "<<slim<<" Jy"<<endl
+     <<"   S/N     = "<<SsN<<endl
+     <<"   Mass HI = "<<smass<<" Msol"<<endl;
+ slim_nl = slim * sqrt(nlobes);
+ SsN_nl = ssig/slim_nl;
+ smass_nl = mhiref/ssig*slim_nl;
+ cout<<"for "<<nlobes<<" lobes:"<<endl
+     <<"   Flux    = "<<slim_nl<<" Jy"<<endl
+     <<"   S/N     = "<<SsN_nl<<endl
+     <<"   Mass HI = "<<smass_nl<<" Msol"<<endl;
 
  return 0;