Changeset 3193 in Sophya for trunk/Cosmo/SimLSS/cmvobserv3d.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/cmvobserv3d.cc (modified) (11 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -31 +31 @@
      <<" -s snoise : sigma du bruit en Msol"<<endl
      <<" -2 : compute 2D spectrum"<<endl
-     <<" -F : write cube in FITS format"<<endl
+     <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
+     <<" -W : write cube in FITS format"<<endl
      <<" -P : write cube in PPF format"<<endl
      <<" -V : compute variance from real space"<<endl
@@ -65 +66 @@
 
  // *** Schechter mass function definition
- double h75 = 0.71 / 0.75;
- double nstar = 0.006*pow(h75,3.); 
+ double h75 = h100 / 0.75;
+ double nstar = 0.006*pow(h75,3.);
  double mstar = pow(10.,9.8/(h75*h75));  // MSol
  double alpha = -1.37;
@@ -72 +73 @@
  double schmin=1e8, schmax=1e12;
  int schnpt = 1000;
- double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
 
  // *** Niveau de bruit
@@ -91 +91 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:")) != -1) {
   switch (c) {
   case 'a' :
@@ -119 +119 @@
   case '2' :
     comp2dspec = true;
+    break;
+  case 'M' :
+    sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
     break;
   case 'V' :
@@ -134 +137 @@
   }
  }
+
+ double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
 
  string tagobs = "cmvobserv3d.ppf";
@@ -157 +162 @@
  univ.Print();
  double loscomref = univ.Dloscom(zref);
- cout<<"zref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
-
- //-----------------------------------------------------------------
-  cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
+ cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
+ univ.Print(zref);
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Spectrum"<<endl;
 
  InitialSpectrum pkini(ns,as);
@@ -168 +174 @@
 
  GrowthFactor growth(om0,ol0);
+ // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
 
  PkSpectrum0 pk0(pkini,tf);
@@ -173 +180 @@
  PkSpectrumZ pkz(pk0,growth,zref);
  
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create mass function"<<endl;
+
  Schechter sch(nstar,mstar,alpha);
+ sch.Print();
 
  //-----------------------------------------------------------------
@@ -234 +245 @@
  double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
  cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
+ cout<<"Omega_HI at z=0 is "<<mass_by_mpc3/(univ.Rhoc(0.)*GCm3toMsolMpc3_Cst)<<endl
+     <<"         at z="<<zref<<" is "<<mass_by_mpc3/(univ.Rhoc(zref)*GCm3toMsolMpc3_Cst)<<endl;
 
  PrtTim(">>>> End of definition");
@@ -463 +476 @@
  }
 
- cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
- fluct3d.AddNoise2Real(snoise);
- nm = fluct3d.MeanSigma2(rm,rs2);
- cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- PrtTim(">>>> End Add noise");
+ if(snoise>0.) {
+   cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
+   fluct3d.AddNoise2Real(snoise);
+   nm = fluct3d.MeanSigma2(rm,rs2);
+   cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+       <<rs2<<" -> "<<sqrt(rs2)<<endl;
+   PrtTim(">>>> End Add noise");
+ }
 
  //-----------------------------------------------------------------