Changeset 3157 in Sophya for trunk/Cosmo/SimLSS/genefluct3d.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/genefluct3d.cc (modified) (5 diffs)

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -32 +32 @@
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
   : T_(T) , Nx_(0) , Ny_(0) , Nz_(0) , array_allocated_(false) , lp_(0)
-  , redshref_(-999.) , kredshref_(-999.)
-{
+  , redshref_(-999.) , kredshref_(0.) , cosmo_(NULL) , growth_(NULL)
+  , loscom_ref_(-999.), loscom_min_(-999.), loscom_max_(-999.)
+
+
+{
+ xobs_[0] = xobs_[1] = xobs_[2] = 0.;
+ zred_.resize(0);
+ loscom_.resize(0);
  SetNThread();
 }
@@ -64 +70 @@
 //     kredshref = indice (en double) correspondant a ce redshift
 //                 Si kredshref<0 alors kredshref=0
+// Exemple: redshref=1.5 kredshref=250.75
+//    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
 {
  if(redshref<0.) redshref = 0.;
  if(kredshref<0.) kredshref = 0.;
- redshref_ = redshref;
+ redshref_  = redshref;
  kredshref_ = kredshref;
+}
+
+void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
+{
+ cosmo_ = &cosmo;
+ if(lp_>1) cosmo_->Print();
+}
+
+void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
+{
+ growth_ = &growth;
 }
 
@@ -164 +183 @@
 }
 
+//-------------------------------------------------------
+long GeneFluct3D::LosComRedshift(double zinc)
+// Given a position of the cube relative to the observer
+// and a cosmology
+// (SetObservator() and SetCosmology() should have been called !)
+// This routine filled:
+//   the vector "zred_" of scanned redshift (by zinc increments)
+//   the vector "loscom_" of corresponding los comoving distance
+// 
+{
+ if(zinc<=0.) zinc = 0.01;
+ if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<endl;
+
+ if(cosmo_ == NULL || redshref_<0.) {
+   cout<<"GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first"<<endl;
+   throw ParmError("");
+ }
+
+ // On calcule les coordonnees de l'observateur
+ // Il est sur un axe centre sur le milieu de la face Oxy
+ double loscom_ref_ = cosmo_->Dloscom(redshref_);
+ xobs_[0] = Nx_/2.*Dx_;
+ xobs_[1] = Ny_/2.*Dy_;
+ xobs_[2] = kredshref_*Dz_ - loscom_ref_;
+
+ // L'observateur est-il dans le cube?
+ bool obs_in_cube = false;
+ if(xobs_[2]>=0. && xobs_[2]<=Nz_*Dz_) obs_in_cube = true;
+
+ // Find MINIMUM los com distance to the observer:
+ // c'est le centre de la face a k=0
+ // (ou zero si l'observateur est dans le cube)
+ loscom_min_ = 0.;
+ if(!obs_in_cube) loscom_min_ = -xobs_[2];
+
+ // Find MAXIMUM los com distance to the observer:
+ // ou que soit positionne l'observateur, la distance
+ // maximal est sur un des coins du cube
+ loscom_max_ = 0.;
+ for(long i=0;i<=1;i++) {
+   double dx2 = xobs_[0] - i*Nx_*Dx_; dx2 *= dx2;
+   for(long j=0;j<=1;j++) {
+     double dy2 = xobs_[1] - j*Ny_*Dy_; dy2 *= dy2;
+     for(long k=0;k<=1;k++) {
+       double dz2 = xobs_[2] - k*Nz_*Dz_; dz2 *= dz2;
+       dz2 = sqrt(dx2+dy2+dz2);
+       if(dz2>loscom_max_) loscom_max_ = dz2;
+     }
+   }
+ }
+ if(lp_>0) {
+   cout<<"...zref="<<redshref_<<" kzref="<<kredshref_<<" losref="<<loscom_ref_<<" Mpc\n"
+       <<"   xobs="<<xobs_[0]<<" , "<<xobs_[1]<<" , "<<xobs_[2]<<" Mpc "
+       <<" in_cube="<<obs_in_cube
+       <<" loscom_min="<<loscom_min_<<" loscom_max="<<loscom_max_<<" Mpc "<<endl;
+ }
+
+ // Fill the corresponding vectors
+ for(double z=0.; ; z+=zinc) {
+   double dlc = cosmo_->Dloscom(z);
+   if(dlc<loscom_min_) {zred_.resize(0); loscom_.resize(0);}
+   zred_.push_back(z);
+   loscom_.push_back(dlc);
+   z += zinc;
+   if(dlc>loscom_max_) break; // on break apres avoir stoque un dlc>dlcmax
+ }
+
+ long n = zred_.size();
+ if(lp_>0) {
+   cout<<"...n="<<n;
+   if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
+   if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
+   cout<<endl;
+ }
+
+ return n;
+}
 
 //-------------------------------------------------------
@@ -323 +419 @@
 
  // --- On remplit avec une realisation
- if(lp_>0) cout<<"...before Fourier realization filling ---"<<endl;
+ if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
  T_(0,0,0) = complex<r_8>(0.);  // on coupe le continue et on l'initialise
  long lmod = Nx_/10; if(lmod<1) lmod=1;
@@ -529 +625 @@
 
 //-------------------------------------------------------------------
+void GeneFluct3D::ApplyGrowthFactor(long npoints)
+// Apply Growth to real space
+// Using the correspondance between redshift and los comoving distance
+// describe in vector "zred_" "loscom_"
+{
+ if(npoints<3) npoints = zred_.size();
+ if(lp_>0) cout<<"--- ApplyGrowthFactor ---  npoints="<<npoints<<endl;
+ check_array_alloc();
+
+ if(growth_ == NULL) {
+   cout<<"GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first"<<endl;
+   throw ParmError("GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first");
+ }
+
+ long n = zred_.size();
+ InverseFunc invfun(zred_,loscom_);
+ vector<double> invlos;
+ invfun.ComputeParab(npoints,invlos);
+
+ InterpFunc interpinv(invfun.YMin(),invfun.YMax(),invlos);
+ unsigned short ok;
+
+ //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;
+   for(long j=0;j<Ny_;j++) {
+     double dy2 = xobs_[1] - j*Dy_; 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);
+       //CHECK: hgr.Add(z);
+       double dzgr = (*growth_)(z);   // interpolation par morceau
+       //double dzgr = growth_->Linear(z,ok);  // interpolation lineaire
+       //double dzgr = growth_->Parab(z,ok);  // interpolation parabolique
+       int_8 ip = IndexR(i,j,l);
+       data_[ip] *= dzgr;
+     }
+   }
+ }
+
+ //CHECK: {POutPersist pos("applygrowth.ppf"); string tag="hgr"; pos.PutObject(hgr,tag);}
+
+}
+
+//-------------------------------------------------------------------
 void GeneFluct3D::ComputeReal(void)
 // Calcule une realisation dans l'espace reel