Changeset 3157 in Sophya for trunk/Cosmo/SimLSS/cmvobserv3d.cc

Timestamp:

Jan 25, 2007, 6:04:48 PM (19 years ago)

Author:

cmv

Message:

intro du facteur de croissance dans la simul cmv 25/01/2007

File:

• trunk/Cosmo/SimLSS/cmvobserv3d.cc (modified) (13 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -12 +12 @@
 
 #include "constcosmo.h"
+#include "cosmocalc.h"
 #include "schechter.h"
 #include "geneutils.h"
@@ -23 +24 @@
      <<" -a : init auto de l aleatoire"<<endl
      <<" -0 : use methode ComputeFourier0"<<endl
+     <<" -G : compute Pk(z=0) and apply growth factor in real space"<<endl
      <<" -x nx,dx : taille en x (npix,Mpc)"<<endl
      <<" -y ny,dy : taille en y (npix,Mpc)"<<endl
@@ -46 +48 @@
 
  // *** Cosmography definition   (WMAP)
+ unsigned short flat = 0;
  double ob0 = 0.0444356;
  double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
  double zref = 0.5;
+ double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
 
  // *** Spectrum and variance definition
@@ -75 +79 @@
  // *** type de generation
  bool computefourier0=false;
+ bool use_growth_factor = false;
  unsigned short nthread=4;
 
@@ -86 +91 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2x:y:z:s:Z:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:")) != -1) {
   switch (c) {
   case 'a' :
@@ -93 +98 @@
   case '0' :
     computefourier0 = true;
+    break;
+  case 'G' :
+    use_growth_factor = true;
     break;
   case 'x' :
@@ -142 +150 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
+ cout<<endl<<"\n--- Create Cosmology"<<endl;
+
+ CosmoCalc univ(flat,true,zref+1.);
+ univ.SetInteg(perc,dzinc,dzmax,glorder);
+ univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.Print();
+ double loscomref = univ.Dloscom(zref);
+ cout<<"zref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
+
+ //-----------------------------------------------------------------
+  cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
 
  InitialSpectrum pkini(ns,as);
@@ -149 +167 @@
  //tf.SetNoOscEnv(2);
 
- GrowthFactor d1(om0,ol0);
+ GrowthFactor growth(om0,ol0);
 
  PkSpectrum0 pk0(pkini,tf);
 
- PkSpectrumZ pkz(pk0,d1,zref);
+ PkSpectrumZ pkz(pk0,growth,zref);
  
  Schechter sch(nstar,mstar,alpha);
@@ -223 +241 @@
  cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
 
- pkz.SetZ(zref);
  TArray< complex<r_8> > pkgen;
  GeneFluct3D fluct3d(pkgen);
@@ -229 +246 @@
  fluct3d.SetNThread(nthread);
  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ fluct3d.SetObservator(zref,nz/2.);
+ fluct3d.SetCosmology(univ);
+ fluct3d.SetGrowthFactor(growth);
+ fluct3d.LosComRedshift(0.001);
  TArray<r_8>& rgen = fluct3d.GetRealArray();
- fluct3d.Print();
+ cout<<endl; fluct3d.Print();
 
  double dkmin = fluct3d.GetKincMin();
@@ -246 +267 @@
      <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
 
+ //-----------------------------------------------------------------
  cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
  // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
@@ -258 +280 @@
  PrtTim(">>>> End Initialisation de GeneFluct3D");
 
+ //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
+ if(use_growth_factor) pkz.SetZ(0.); else pkz.SetZ(zref);
+ cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
  if(computefourier0) fluct3d.ComputeFourier0(pkz);
    else              fluct3d.ComputeFourier(pkz);
@@ -326 +351 @@
  }
 
+ //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in real space"<<endl;
  fluct3d.ComputeReal();
  double rmin,rmax; rgen.MinMax(rmin,rmax);
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
-   PrtTim(">>>> End Computing a realization in real space");
+ PrtTim(">>>> End Computing a realization in real space");
+
+ if(use_growth_factor) {
+   cout<<"\n--- Apply Growth factor"<<endl;
+   cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
+   fluct3d.ApplyGrowthFactor(-1);
+   rmin,rmax; rgen.MinMax(rmin,rmax);
+   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   PrtTim(">>>> End Applying growth factor");
+ }
 
  if(wfits) {