Changeset 3351 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Oct 12, 2007, 7:45:12 PM (18 years ago)

Author:

cmv

Message:

suite des modifs importantes , cmv 12/10/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvdefsurv.cc (modified) (5 diffs)
• cmvobserv3d.cc (modified) (9 diffs)
• genefluct3d.cc (modified) (2 diffs)

trunk/Cosmo/SimLSS/cmvdefsurv.cc

@@ -41 +41 @@
       <<"   \'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
+      <<" -K k,dk,pk : 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
@@ -146 +146 @@
     break;
   case 'K' :
-    sscanf(optarg,"%lf,%lf,%lf",&kcosm,&pkcosm,&dkcosm);
+    sscanf(optarg,"%lf,%lf,%lf",&kcosm,&dkcosm,&pkcosm);
     break;
   case 'h' :
@@ -252 +252 @@
      <<"  dy = "<<dy<<" Mpc,  tylarg = "<<tylarg<<" Mpc com, ny = "<<ny<<" pix"<<endl;
  cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
+ cout<<"     Volume pixel = "<<dx*dy*dz<<" Mpc^3"<<endl;
+ cout<<"     Volume total = "<<Npix*dx*dy*dz<<" Mpc^3"<<endl;
 
  // --- Cosmolographie Transverse
@@ -308 +310 @@
  // --- Variance cosmique
  //                      cosmique                    poisson
- // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi dk Vsurvey) * [(1+n*P)/(n*P)]^2
+ // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi k^2 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.) {
@@ -322 +324 @@
    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 vcosm = sqrt( 2.*pow(2.*M_PI,3.)/(4.*M_PI*pow(kcosm,2.)*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

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -39 +39 @@
      <<"           typevol=2 noise evolved with distance to middle of Z planes"<<endl
      <<" -2 : compute also 2D spectrum (default: no)"<<endl
-     <<" -N scalecube,offsetcube: normalize cube before doing final spectrum (default: automatic)"<<endl
+     <<" -N scalemass: facteur d\'unite pour la masse (default: 1)"<<endl
+     <<"               ex: si on veut unites 10^8 Msol -> scalemass=1.e-8"<<endl
+     <<"               si <0 alors facteur=-scalemass*Mpix"<<endl
      <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
      <<"                         If nsch<0 alors no,bre de points par decade"<<endl
@@ -97 +99 @@
  bool no_poisson = false;
 
+ double scalemass = 1.;
+
  // *** Niveau de bruit
  double snoise= 0.;   // en equivalent MSol
@@ -118 +122 @@
  bool wslice = false;
  bool compvarreal = false;
- double scalecube = -1., offsetcube = 0.;
 
  // --- Decodage arguments
@@ -165 +168 @@
     break;
   case 'N' :
-    sscanf(optarg,"%lf,%lf",&scalecube,&offsetcube);
+    sscanf(optarg,"%lf",&scalemass);
+    if(scalemass==0.) scalemass = 1.;
     break;
   case 'M' :
@@ -223 +227 @@
  if(no_poisson) cout<<"No poisson fluctuation, direct conversion to HI mass"<<endl;
  cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
- cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
+ cout<<"scalemass="<<scalemass<<endl;
  ////if(do_agn)
  ////  cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
@@ -364 +368 @@
  cout<<endl; fluct3d.Print();
  cout<<"\nMean number of galaxies per pixel = "<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
- cout<<"Mean mass per pixel = "<<mass_by_mpc3*fluct3d.GetDVol()<<endl;
+ double mass_by_pixel = mass_by_mpc3 * fluct3d.GetDVol();
+ cout<<"Mean mass per pixel = "<<mass_by_pixel<<endl;
 
  double dkmin = fluct3d.GetKincMin();
@@ -536 +541 @@
  if(no_poisson) {
    cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
-       <<mass_by_mpc3*fluct3d.GetDVol()<<endl;
-   rm = fluct3d.TurnMass2MeanNumber(mass_by_mpc3);
+       <<mass_by_pixel<<endl;
+   rm = fluct3d.TurnMass2MeanNumber(mass_by_mpc3); // ici on doit donner Msol/Mpc^3
  } else {
    cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
@@ -581 +586 @@
  }
 
+ //-----------------------------------------------------------------
  if(wfits) {
    fluct3d.WriteFits("!cmvobserv3d_r.fits");
@@ -613 +619 @@
  }
 
-
- //-----------------------------------------------------------------
- if(scalecube!=0.) {  // Si scalecube==0 pas de normalisation
-   cout<<"\n--- Scale cube rs2ref="<<rs2ref<<endl;
-   if(scalecube<0. && rs2ref>0.) {  // si negatif on scale automatiquement
+ //-----------------------------------------------------------------
+ if(scalemass!=0. && scalemass!=1.) {  // Si scalemass==0 pas de normalisation
+   if(scalemass<0.) scalemass = 1. / (-scalemass * mass_by_pixel);
+   cout<<"\n--- Scale cube scale="<<scalemass<<endl;
+   fluct3d.ScaleOffset(scalemass);
+   snoisesave *= scalemass;
      nm = fluct3d.MeanSigma2(rm,rs2);
-     if(rs2>0.) {scalecube=sqrt(rs2ref)/sqrt(rs2); offsetcube=-rm;}
-   }
-   cout<<"...scale="<<scalecube<<" offset="<<offsetcube<<endl;
-   if(scalecube>0.) {
-     fluct3d.ScaleOffset(scalecube,offsetcube);
-     snoisesave *= scalecube;
-   }
    PrtTim(">>>> End Scale cube");
  }
- 
+
+ //----------------------------------------------------------------- 
  if(wfits) {
    fluct3d.WriteFits("!cmvobserv3d_rf.fits");

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -961 +961 @@
  }
  herr.ToVariance();
+ for(int i=0;i<herr.NBins();i++) herr(i) += sigma2;
 
  // renormalize to directly compare to original spectrum
@@ -1007 +1008 @@
  }
  herr.ToVariance();
+ for(int i=0;i<herr.NBinX();i++)
+    for(int j=0;j<herr.NBinY();j++) herr(i,j) += sigma2;
 
  // renormalize to directly compare to original spectrum