Changeset 3331 in Sophya

Timestamp:

Oct 2, 2007, 3:11:22 PM (18 years ago)

Author:

cmv

Message:

mise en place evolution seulement sur distance au plan Z , cmv 02/10/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (12 diffs)
• genefluct3d.cc (modified) (6 diffs)
• genefluct3d.h (modified) (2 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -23 +23 @@
      <<" -a : auto init random seed (needed for multiple simul)"<<endl
      <<" -0 : use ComputeFourier0 method (defaut: no, use normal way)"<<endl
-     <<" -G : compute Pk(z=0) and apply growth factor in real space"<<endl
+     <<" -G typevol: compute Pk(z=0) and apply growth factor in real space"<<endl
+     <<"             typevol=1 noise evolved with distance / observateur (def)"<<endl
+     <<"             typevol=2 noise evolved with distance to middle of Z planes"<<endl
      <<"      (default: no, spectrum Pk(z=z_median) for all cube)"<<endl
      <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
@@ -32 +34 @@
      <<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
      <<" -Z zref : redshift for the center of the simulation cube"<<endl
-     <<" -s snoise,evol : gaussian noise sigma in equivalent Msol"<<endl
-     <<"                  if evol>0 noise evolved with distance (def no)"<<endl
+     <<" -s snoise,typevol : gaussian noise sigma in equivalent Msol"<<endl
+     <<"           typevol=0 no evolution (def)"<<endl
+     <<"           typevol=1 noise evolved with distance / observateur"<<endl
+     <<"           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
@@ -95 +99 @@
  // *** Niveau de bruit
  double snoise= 0.;   // en equivalent MSol
- int isnoise_evol = 0;
+ int noise_evol = 0;
 
  // *** AGN
@@ -104 +108 @@
  // *** type de generation
  bool computefourier0=false;
- bool use_growth_factor = false;
+ int use_growth_factor = 0;
  unsigned short nthread=0;
  int filter_by_pixel = 1;
@@ -124 +128 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWSV2GUF:x:y:z:s:Z:M:A:T:N:Q:R:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWSV2UG:F:x:y:z:s:Z:M:A:T:N:Q:R:")) != -1) {
   int nth = 0;
   switch (c) {
@@ -134 +138 @@
     break;
   case 'G' :
-    use_growth_factor = true;
+    sscanf(optarg,"%d",&use_growth_factor);
     break;
   case 'U' :
@@ -152 +156 @@
     break;
   case 's' :
-    sscanf(optarg,"%lf,%d",&snoise,&isnoise_evol);
+    sscanf(optarg,"%lf,%d",&snoise,&noise_evol);
     break;
   case 'Z' :
@@ -212 +216 @@
  cout<<"Filter by pixel = "<<filter_by_pixel<<endl;
  cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
+ cout<<"Use_growth_factor = "<<use_growth_factor<<endl;
  cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
  cout<<"schmin="<<schmin<<" ("<<lschmin
@@ -217 +222 @@
      <<", schnpt="<<schnpt<<endl;
  if(no_poisson) cout<<"No poisson fluctuation, direct conversion to HI mass"<<endl;
- cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<isnoise_evol<<endl;
+ cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
  cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
  if(do_agn)
@@ -385 +390 @@
  //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
- if(use_growth_factor) pkz.SetZ(0.); else pkz.SetZ(zref);
+ if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
  cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
  if(computefourier0) fluct3d.ComputeFourier0(pkz);
@@ -473 +478 @@
  PrtTim(">>>> End Computing a realization in real space");
 
- if(use_growth_factor) {
+ if(use_growth_factor>0) {
    cout<<"\n--- Apply Growth factor"<<endl;
    cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
-   fluct3d.ApplyGrowthFactor();
+   fluct3d.ApplyGrowthFactor(use_growth_factor);
    fluct3d.MinMax(rmin,rmax);
    cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
@@ -607 +612 @@
  double snoisesave = 0.;
  if(snoise>0.) {
-   cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<", evolution="<<isnoise_evol<<endl;
-   fluct3d.AddNoise2Real(snoise,(isnoise_evol>0? true:false));
+   cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<", evolution="<<noise_evol<<endl;
+   fluct3d.AddNoise2Real(snoise,noise_evol);
    snoisesave = snoise;
      nm = fluct3d.MeanSigma2(rm,rs2);

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -744 +744 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ApplyGrowthFactor(void)
+void GeneFluct3D::ApplyGrowthFactor(int type_evol)
 // Apply Growth to real space
 // Using the correspondance between redshift and los comoving distance
 // describe in vector "zred_" "loscom_"
-{
- if(lp_>0) cout<<"--- ApplyGrowthFactor ---"<<endl;
+// type_evol = 1 : evolution de la puissance du bruit avec la distance a l'observateur
+//             2 : evolution de la puissance du bruit avec la distance du plan Z
+//             (tous les pixels d'un plan Z sont mis au meme redshift z que celui du milieu)
+{
+ if(lp_>0) cout<<"--- ApplyGrowthFactor: evol="<<type_evol<<endl;
  check_array_alloc();
 
@@ -756 +759 @@
    cout<<bla<<endl; throw ParmError(bla);
  }
+ if(type_evol<1 || type_evol>2) {
+   char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
@@ -762 +769 @@
  //CHECK: Histo hgr(0.9*zred_[0],1.1*zred_[n-1],1000);
  for(long i=0;i<Nx_;i++) {
-   double dx2 = xobs_[0] - i*Dx_; dx2 *= dx2;
+   double dx2 = DXcom(i); dx2 *= dx2;
    for(long j=0;j<Ny_;j++) {
-     double dy2 = xobs_[1] - j*Dy_; dy2 *= dy2;
+     double dy2 = DYcom(j); dy2 *= dy2;
      for(long l=0;l<Nz_;l++) {
-       double dz2 = xobs_[2] - l*Dz_; dz2 *= dz2;
-       dz2 = sqrt(dx2+dy2+dz2);
-       double z = interpinv(dz2);
+       double dz = DZcom(l);
+       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
+         else dz = fabs(dz); // tous les plans Z au meme redshift
+       double z = interpinv(dz);
        //CHECK: hgr.Add(z);
        double dzgr = (*growth_)(z);   // interpolation par morceau
@@ -1416 +1424 @@
 }
 
-void GeneFluct3D::AddNoise2Real(double snoise,bool with_evol)
+void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
 // add noise to every pixels (meme les <=0 !)
-{
- if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<with_evol<<endl;
- check_array_alloc();
+// type_evol = 0 : pas d'evolution de la puissance du bruit
+//             1 : evolution de la puissance du bruit avec la distance a l'observateur
+//             2 : evolution de la puissance du bruit avec la distance du plan Z
+//                 (tous les plans Z sont mis au meme redshift z de leur milieu)
+{
+ if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
+ check_array_alloc();
+
+ if(type_evol<0) type_evol = 0;
+ if(type_evol>2) {
+   char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
 
  vector<double> correction;
  InterpFunc *intercor = NULL;
 
- if(with_evol) {
+ if(type_evol>0) {
    // Sigma_Noise(en mass) :
    //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
@@ -1454 +1472 @@
  } 
 
- double dx2=0., dy2=0., dz2=0.;
+ double corrlim[2] = {1.,1.};
  for(long i=0;i<Nx_;i++) {
-   dx2 = DXcom(i); dx2 *= dx2;
+   double dx2 = DXcom(i); dx2 *= dx2;
    for(long j=0;j<Ny_;j++) {
-     dy2 = DYcom(j); dy2 *= dy2;
+     double dy2 = DYcom(j); dy2 *= dy2;
      for(long l=0;l<Nz_;l++) {
        double corr = 1.;
-       if(with_evol) {
-         dz2 = DZcom(l); dz2 *= dz2; dz2 = sqrt(dx2+dy2+dz2);
-         corr = (*intercor)(dz2);
+       if(type_evol>0) {
+         double dz = DZcom(l);
+         if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
+           else dz = fabs(dz); // tous les plans Z au meme redshift
+         corr = (*intercor)(dz);
+         if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
        }
        int_8 ip = IndexR(i,j,l);
@@ -1470 +1491 @@
    }
  }
+ if(type_evol>0)
+   cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;
 
  if(intercor!=NULL) delete intercor;

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -93 +93 @@
 
   void ComputeReal(void);
-  void ApplyGrowthFactor(void);
+  void ApplyGrowthFactor(int type_evol=1);
 
   void ReComputeFourier(void);
@@ -118 +118 @@
   double TurnMass2Flux(void);
   void AddAGN(double lfjy,double lsigma,double powlaw=0.);
-  void AddNoise2Real(double snoise,bool with_evol=false);
+  void AddNoise2Real(double snoise,int type_evol=0);
 
   void WriteFits(string cfname,int bitpix=FLOAT_IMG);