Changeset 3350 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Oct 11, 2007, 6:07:47 PM (18 years ago)

Author:

cmv

Message:

add cosmic variance and poisson noise estimate for Pk, cmv 11/10/2007

File:

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

trunk/Cosmo/SimLSS/cmvdefsurv.cc

@@ -8 +8 @@
 #include <unistd.h>
 
+#include <vector>
+
 #include "constcosmo.h"
 #include "cosmocalc.h"
 #include "geneutils.h"
 #include "schechter.h"
+#include "pkspectrum.h"
 #include "planckspectra.h"
 
@@ -27 +30 @@
 void usage(void);
 void usage(void) {
-  cout<<"cmvdefsurv [-r] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
+  cout<<"cmvdefsurv [options] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
       <<"----------------"<<endl
       <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
@@ -37 +40 @@
       <<"   \'F\' : en frequence (pour Z)    : resolution et largeur MHz"<<endl
       <<"   \'M\' : en distance (pour X-Y-Z) : resolution et largeur Mpc"<<endl
+      <<"----------------"<<endl
+      <<" -K k,pk,dk : k(Mpc^-1) dk(Mpc^-1) pk(a z=0 en Mpc^-3) pour estimer la variance cosmique"<<endl
       <<"----------------"<<endl
       <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
@@ -82 +87 @@
  double redshift = 0.;
  double tobs = 6000., surfeff = 400000.;
+ // variance cosmique (default = standard SDSSII)
+ double kcosm = 0.05, dkcosm = -1., pkcosm = 40000.;
  // taille du lobe d'observation en arcmin pour la frequence
  double lobewidth0 = -1., lobefreq0 = Fr_HyperFin_Par*1.e3;
@@ -96 +103 @@
  // --- Decodage arguments
  char c;
-  while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
+  while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:K:")) != -1) {
   switch (c) {
   case 'x' :
@@ -137 +144 @@
   case 'A' :
     sscanf(optarg,"%lf",&lflux_agn);
+    break;
+  case 'K' :
+    sscanf(optarg,"%lf,%lf,%lf",&kcosm,&pkcosm,&dkcosm);
     break;
   case 'h' :
@@ -158 +168 @@
  univ.Print(0.);
  univ.Print(redshift);
+ GrowthFactor growth(om0,ol0);
+ double grothfac = growth(redshift);
+ cout<<"Facteur de croissance lineaire = "<<grothfac
+     <<" ^2 , ( 1/(1+z)="<<1./(1.+redshift)<<" )"<<endl;
 
  double dang    = univ.Dang(redshift);
@@ -291 +305 @@
  cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1  (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
  cout<<"Nyquist:    kz = "<<knyq_z<<" Mpc^-1  (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
+
+ // --- Variance cosmique
+ //                      cosmique                    poisson
+ // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi dk Vsurvey) * [(1+n*P)/(n*P)]^2
+ // nombre de mode = 1/2 * Vsurvey/(2Pi)^3 * 4Pi*k^2*dk
+ if(kcosm>0. && pkcosm>0.) {
+   double pk = pkcosm*pow(grothfac,2.);
+   cout<<"\n>>>>\n>>>> variance cosmique pour k="<<kcosm<<" Mpc^-1, pk(z=0)="
+       <<pkcosm<<", pk(z="<<redshift<<")="<<pk<<"\n>>>>"<<endl;
+   for(int i=0;i<3;i++) {  // la correction de variance du au bruit de poisson
+     double v = 1.1; if(i==1) v=1.5; if(i==2) v=2.0;
+     double ngal = 1./(v-1.)/pk;
+     cout<<"  pour "<<ngal<<" gal/Mpc^3 on multiplie par "<<v<<" sigma/P"<<endl;
+   }
+   vector<double> dk; if(dkcosm>0.) dk.push_back(dkcosm);
+   dk.push_back(dk_x); dk.push_back(dk_y); dk.push_back(dk_z);
+   for(int i=0;i<(int)dk.size();i++) { // la variance cosmique pure
+     double vcosm = sqrt( 2.*pow(2.*M_PI,3.)/(4.*M_PI*dk[i]*vol) );
+     double nmode = 0.5*vol/pow(2.*M_PI,3.) * 4.*M_PI*pow(kcosm,2.)*dk[i];
+     cout<<"  pour dk = "<<dk[i]<<" Mpc^-1:  sigma/P = "<<vcosm
+         <<" , Nmode = "<<nmode<<endl;
+   }
+ }
 
  // --- Masse de HI