Changeset 3339 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Oct 4, 2007, 4:03:07 PM (18 years ago)

Author:

cmv

Message:

encore mise au point du calcul du bruit cmv 04/10/2007

File:

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

trunk/Cosmo/SimLSS/cmvdefsurv.cc

@@ -90 +90 @@
  double hifactor = 1.;
  double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
+ bool ya2polar = false;
  double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
  double lflux_agn = -3.;
@@ -131 +132 @@
     break;
   case '2' :
+    ya2polar = true;
     facpolar = 1.0;
     break;
@@ -360 +362 @@
 
  // Signal
- double psig = facpolar * fhi * surfeff;
- double tsig = psig / k_Boltzman_Cst / (dnuhiz*1e9);
- double ssig = psig / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
+ double psig_2polar = fhi * surfeff;
+ double tsig_2polar = psig_2polar / k_Boltzman_Cst / (dnuhiz*1e9);
+ double ssig_2polar = psig_2polar / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
+ double psig = facpolar * psig_2polar;
+ double tsig = facpolar * tsig_2polar;
+ double ssig = facpolar * ssig_2polar;
  cout<<"\nSignal("<<mhiref<<" Msol): P="<<psig<<" W"<<endl;
  cout<<"     flux density = "<<ssig<<" Jy  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
@@ -379 +384 @@
         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*1e9);
- double ssynch = psynch / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
+ double psynch_2polar = planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
+ double ssynch_2polar = psynch_2polar / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
+ double psynch = facpolar * psynch_2polar;
+ double ssynch = facpolar * ssynch_2polar;
  cout<<"\nSynchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
      <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
@@ -389 +396 @@
  double tcmb = T_CMB_Par;
  planck.SetTemperature(tcmb);
- double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
- double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ double pcmb_2polar = planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
+ double scmb_2polar = pcmb_2polar / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ double pcmb = facpolar * pcmb_2polar;
+ double scmb = facpolar * scmb_2polar;
  cout<<"\nCMB: T="<<tcmb<<" K -> P="<<pcmb<<" W for pixel"<<endl;
  cout<<"                       flux density = "<<scmb<<" Jy for pixel solid angle"<<endl;
@@ -405 +414 @@
      <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
 
+ // =====================================================================
  // ---
  // --- Noise analysis
  // ---
+ // --- Puissance du bruit pour un telescope de surface Ae et de BW dNu
+ // Par definition la puissance du bruit est:
+ //     Pb = k * Tsys * dNu  (W)
+ // Pour une source (non-polarisee) de densite de flux (totale 2 polars)
+ //     St (exprimee en Jy=W/m^2/Hz)
+ //     Pt = St * Ae * dNu   (puissance totale emise en W pour 2 polars)
+ //     P1 = 1/2 * St * Ae * dNu   (puissance emise en W pour une polar)
+ // la SEFD (system equivalent flux density en Jy) est definie comme
+ //   la densite de flux total (2 polars) "St" d'une source (non-polarisee)
+ //   dont la puissance P1 mesuree pour une seule polarisation
+ //   serait egale a la puissance du bruit. De P1 = Pb on deduit:
+ //   SEFD = 2 * k * Tsys / Ae    (en Jy)
+ //   la puissance du bruit est: Pb = 1/2 * SEFD * Ae * dNu  (en W)
+ // la sensibilite Slim tient compte du temps d'integration et de la BW:
+ //   le nombre de mesures independantes est "2*dNu*Tobs" donc
+ //   Slim = SEFD / sqrt(2*dNu*Tobs) =  2*k*Tsys/[Ae*sqrt(2*dNu*Tobs) (en Jy)
+ // --- Puissance du bruit pour un interferometre
  // Ae = surface d'un telescope elementaire
  // N  = nombre de telescopes dans l'interferometre (Atot = N*Ae)
+ // La sensibilite Slim en Jy est:
  // 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:
+ // - Interferometre a deux antennes:
  // Slim =  2 * k * Tsys / [ Atot * Sqrt(1/2 *dnu*Tobs) ]
- // Interferometre a N antennes (N grand):
+ // - 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)
+ // --- On ne mesure qu'une seule polarisation
+ // Ces formules sont valables si on mesure 1 polarisation:
+ // Slim est la densite de flux total "St" (2 polars) d'une source (non-polarisee)
+ //   qui donne la meme puissance que le bruit dans un detecteur qui ne
+ //   mesure qu'une seule polarisation:
+ // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
+ //    S/N = St / Slim
+ // La puissance de bruit est, par definition:
+ //    Pb = 1/2 *Slim*Atot*dNu
+ //         = k*Tsys*sqrt(2*dNu/Tobs)  pour N=2
+ //         = k*Tsys*sqrt(dNu/Tobs)    pour N>>grand
+ // La densite de flux d'une source a S/N=1 est:
+ //    St = Slim
+ // La puissance d'une source a S/N=1 mesuree par un detecteur
+ //    qui ne mesure qu'une polar est:
+ //    P1_lim = 1/2 *Slim*Atot*dNu
+ // --- On mesure les 2 polarisations avec deux voies d'electronique distinctes
+ // la puissance du signal mesure est multipliee par 2
+ // la puissance du bruit est multipliee par sqrt(2)
+ // on a donc un gain d'un facteur sqrt(2) sur le rapport S/N
+ // (cela revient d'ailleur a doubler le temps de pose: Tobs -> 2*Tobs)
+ // En notant arbitrairement: Slim' = Slim / sqrt(2)
+ //    ou Slim est defini par les formules ci-dessus
+ // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
+ //   (S/N)_2 = (S/N)_1 * sqrt(2) = (St / Slim) * sqrt(2) = St / Slim'
+ // La densite de flux d'une source a S/N=1 est:
+ //    St = Slim' = Slim / sqrt(2)
+ // La puissance d'une source a S/N=1 cumulee par les 2 detecteurs est:
+ //    P_lim = St*Atot*dNu = Slim'*Atot*dNu = 1/sqrt(2) *Slim*Atot*dNu
+ //          = P1_lim * sqrt(2)
+ // La puissance de bruit cumulee par les 2 detecteurs est, par definition:
+ //    Pb = P_lim = Slim'*Atot*dNu = P1_lim * sqrt(2)
+ //               = 2*k*Tsys*sqrt(dNu/Tobs)   pour N=2
+ //               = k*Tsys*sqrt(2*dNu/Tobs)   pour N>>grand
+ // =====================================================================
+
  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;
 
+ cout<<"...Computation assume that noise dominate the signal."<<endl;
+ if(ya2polar)
+   cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<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;
+
+ //---
+ for(unsigned short it=0;it<2;it++) {
+
+   double fac = 1.;
+   if(it==0) { // Interferometre a 2 telescopes
+     fac = 0.5;
+     cout<<"\n...Observation limits for a 2 telescope interferometer (with complex correlator)"<<endl
+         <<"   (sensitivity is given for real or complex correlator output)"<<endl;
+   } else if (it==1) { // Interferometre a N>> telescopes
+     fac = 1.;
+     cout<<"\n...Observation limits for a N (large) telescope interferometer (with complex correlator)"<<endl
+         <<"     (weak source limit sensitivity in a synthetised image)"<<endl;
+   } else continue;
+
+   slim = 2. * k_Boltzman_Cst * Tsys / surfeff
+               / sqrt(fac*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
+   if(ya2polar) slim /= sqrt(2.);
+   SsN = ssig_2polar/slim;
+   smass = mhiref / ssig_2polar * slim;
+   cout<<"for 1 lobe:"<<endl
+       <<"   Slim    = "<<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 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;