Changeset 3199 in Sophya for trunk/Cosmo/SimLSS/genefluct3d.cc

Timestamp:

Apr 5, 2007, 7:19:06 PM (18 years ago)

Author:

cmv

Message:

Intro raffinement dans simul AGN, variation dNu Da Dlum selon Oz cmv 05/04/2007

File:

• trunk/Cosmo/SimLSS/genefluct3d.cc (modified) (15 diffs)

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -21 +21 @@
 #include "constcosmo.h"
 #include "geneutils.h"
+#include "schechter.h"
 
 #include "genefluct3d.h"
@@ -33 +34 @@
   , redshref_(-999.) , kredshref_(0.) , cosmo_(NULL) , growth_(NULL)
   , loscom_ref_(-999.), loscom_min_(-999.), loscom_max_(-999.)
-
-
+  , loscom2zred_min_(0.) , loscom2zred_max_(0.)
 {
  xobs_[0] = xobs_[1] = xobs_[2] = 0.;
  zred_.resize(0);
  loscom_.resize(0);
+ loscom2zred_.resize(0);
  SetNThread();
 }
@@ -76 +77 @@
  redshref_  = redshref;
  kredshref_ = kredshref;
+ if(lp_>0)
+   cout<<"--- GeneFluct3D::SetObservator zref="<<redshref_<<" kref="<<kredshref_<<endl;
 }
 
@@ -94 +97 @@
                 <<" D="<<dx<<","<<dy<<","<<dz<<endl;
  if(nx<=0 || dx<=0.) {
-   cout<<"GeneFluct3D::setsize_Error: bad value(s)"<<endl;
-   throw ParmError("GeneFluct3D::setsize_Error: bad value(s)");
+   char *bla = "GeneFluct3D::setsize_Error: bad value(s";
+   cout<<bla<<endl; throw ParmError(bla);
  }
 
@@ -161 +164 @@
  char bla[90];
  sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
- cout<<bla<<endl;
- throw ParmError(bla);
+ cout<<bla<<endl; throw ParmError(bla);
 }
 
@@ -183 +185 @@
 
 //-------------------------------------------------------
-long GeneFluct3D::LosComRedshift(double zinc)
+long GeneFluct3D::LosComRedshift(double zinc,long npoints)
 // Given a position of the cube relative to the observer
 // and a cosmology
@@ -190 +192 @@
 //   the vector "zred_" of scanned redshift (by zinc increments)
 //   the vector "loscom_" of corresponding los comoving distance
+// -- Input:
+// zinc : redshift increment for computation
+// npoints : number of points required for inverting loscom -> zred
 // 
 {
- if(zinc<=0.) zinc = 0.01;
- if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<endl;
+ if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<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
+   char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ // On calcule les coordonnees de l'observateur dans le repere du cube
+ // cad dans le repere ou l'origine est au centre du pixel i=j=l=0.
+ // L'observateur est sur un axe centre sur le milieu de la face Oxy
  double loscom_ref_ = cosmo_->Dloscom(redshref_);
  xobs_[0] = Nx_/2.*Dx_;
@@ -239 +244 @@
  }
 
- // Fill the corresponding vectors
+ // Fill the corresponding vectors for loscom and zred
+ if(zinc<=0.) zinc = 0.01;
  for(double z=0.; ; z+=zinc) {
    double dlc = cosmo_->Dloscom(z);
@@ -246 +252 @@
    loscom_.push_back(dlc);
    z += zinc;
-   if(dlc>loscom_max_) break; // on break apres avoir stoque un dlc>dlcmax
- }
-
- long n = zred_.size();
+   if(dlc>loscom_max_) break; // on sort apres avoir stoque un dlc>dlcmax
+ }
+
  if(lp_>0) {
-   cout<<"...n="<<n;
+   long n = zred_.size();
+   cout<<"...zred/loscom tables[zinc="<<zinc<<"]: n="<<n;
    if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
    if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
@@ -257 +263 @@
  }
 
- return n;
+ // Compute the parameters and tables needed for inversion loscom->zred
+ if(npoints<3) npoints = zred_.size();
+ InverseFunc invfun(zred_,loscom_);
+ invfun.ComputeParab(npoints,loscom2zred_);
+ loscom2zred_min_ = invfun.YMin();
+ loscom2zred_max_ = invfun.YMax();
+
+ if(lp_>0) {
+   long n = loscom2zred_.size();
+   cout<<"...loscom -> zred[npoints="<<npoints<<"]: n="<<n
+       <<" los_min="<<loscom2zred_min_
+       <<" los_max="<<loscom2zred_max_
+       <<" -> zred=[";
+   if(n>0) cout<<loscom2zred_[0];
+   cout<<",";
+   if(n>1) cout<<loscom2zred_[n-1];
+   cout<<"]"<<endl;
+   if(lp_>1 && n>0)
+     for(int i=0;i<n;i++)
+       if(i==0 || abs(i-n/2)<2 || i==n-1)
+         cout<<"    "<<i<<"  "<<loscom2zred_[i]<<endl;
+ }
+
+ return zred_.size();
 }
 
@@ -624 +653 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ApplyGrowthFactor(long npoints)
+void GeneFluct3D::ApplyGrowthFactor(void)
 // 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;
+ if(lp_>0) cout<<"--- ApplyGrowthFactor ---"<<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);
+   char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
  unsigned short ok;
 
@@ -969 +992 @@
 }
 
-void GeneFluct3D::AddAGN(double lmsol,double lsigma)
+void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
 // Add AGN flux into simulation:
 // --- Procedure:
 // 1. lancer "cmvdefsurv" avec les parametres du survey
+//        (au redshift de reference du survey)
 //    et recuperer l'angle solide "angsol sr" du pixel elementaire
 //    au centre du cube.
@@ -978 +1002 @@
 // 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
 //    cad une moyenne <log10(S)> et un sigma "sig"
-//    Attention: la distribution n'etant pas gaussienne
-// 4. re-lancer "cmvdefsurv" en ajoutant l'info sur les AGN
-//    "cmvdefsurv ... -G <log10(S)> ..."
-//    et recuperer le log10(masse solaire equivalente)
+//    Attention: la distribution n'est pas gaussienne les "mean,sigma"
+//               de l'histo ne sont pas vraiment ce que l'on veut
 // --- Limitations actuelle du code:
-// a. l'angle solide du pixel est pris au centre du cube
-//    et ne varie pas avec la distance a l'interieur du cube
-// b. la taille en dNu des pixels (selon z) est supposee constante
-//    et egale a celle calculee pour le centre du cube
-// c. les AGN sont supposes plats en flux
-// d. le flux des AGN est mis dans une colonne Oz (indice k)
-//    et pas sur la ligne de visee
-// e. la distribution est approximee a une gaussienne
-// .. C'est une approximation pour un observateur loin du centre du cube
-//    et pour un cube peu epais / distance observateur 
+// . les AGN sont supposes en loi de puissance IDENTIQUE pour tout theta,phi
+// . le flux des AGN est mis dans une colonne Oz (indice k) et pas sur la ligne de visee
+// . la distribution est approximee a une gaussienne
+// ... C'est une approximation pour un observateur loin du centre du cube
+//     et pour un cube peu epais / distance observateur 
 // --- Parametres de la routine:
-// lmsol : c'est le <log10(Msol)> qui est la conversion en masse solaire
-//         du flux depose dans un pixel par les AGN
+// llfy : c'est le <log10(S)> du flux depose dans un pixel par les AGN
 // lsigma : c'est le sigma de la distribution
+// powlaw : c'est la pente de ls distribution cad que le flux "lmsol"
+//          et considere comme le flux a 1.4GHz et qu'on suppose une loi
+//             F(nu) = (1.4GHz/nu)^powlaw * F(1.4GHz)
 // - Comme on est en echelle log10():
 //   on tire log10(Msol) + X
@@ -1003 +1022 @@
 // - Pas de probleme de pixel negatif car on a une multiplication!
 {
-  if(lp_>0) cout<<"--- AddAGN: lmsol = "<<lmsol<<" , sigma = "<<lsigma<<endl;
+  if(lp_>0) cout<<"--- AddAGN: <log10(S Jy)> = "<<lfjy<<" , sigma = "<<lsigma<<endl;
   check_array_alloc();
 
-  double m = pow(10.,lmsol);
-
-  double sum=0., sum2=0., nsum=0.;
-  for(long l=0;l<Nz_;l++) {
-    double a = lsigma*NorRand();
-    a = m*pow(10.,a);
-    // On met le meme tirage le long de Oz (indice k)
-    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
-      int_8 ip = IndexR(i,j,l);
-      data_[ip] += a;
-    }
-    sum += a; sum2 += a*a; nsum += 1.;
-  }
-
- if(nsum>1.) {
+ if(cosmo_ == NULL || redshref_<0.| loscom2zred_.size()<1) {
+   char *bla = "GeneFluct3D::AddAGN_Error: set Observator and Cosmology first";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ // La distance angulaire/luminosite/Dnu au centre
+ double dangref = cosmo_->Dang(redshref_);
+ double dlumref = cosmo_->Dlum(redshref_);
+ double dredref  = Dz_/(cosmo_->Dhubble()/cosmo_->E(redshref_));
+ double dnuref = Fr_HyperFin_Par *dredref/pow(1.+redshref_,2.); // GHz
+ double fagnref = pow(10.,lfjy)*(dnuref*1.e9); // Jy.Hz
+ double magnref = FluxHI2Msol(fagnref*Jansky2Watt_cst,dlumref); // Msol
+ if(lp_>0) {
+   cout<<"Au centre: z="<<redshref_<<", dredref="<<dredref<<", dnuref="<<dnuref<<" GHz"<<endl
+       <<" dang="<<dangref<<" Mpc, dlum="<<dlumref<<" Mpc,"
+       <<" fagnref="<<fagnref<<" Jy.Hz (a 1.4GHz), magnref="<<magnref<<" Msol"<<endl;
+ }
+
+ if(powlaw!=0.) {
+   // F(nu) = (nu GHz/1.4 Ghz)^p * F(1.4GHz) et nu = 1.4 GHz / (1+z)
+   magnref *= pow(1/(1.+redshref_),powlaw);
+   if(lp_>0) cout<<" powlaw="<<powlaw<<"  -> change magnref to "<<magnref<<" Msol"<<endl;
+ }
+
+ // Les infos en fonction de l'indice "l" selon Oz
+ vector<double> correction;
+ InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
+ for(long l=0;l<Nz_;l++) {
+   double z = fabs(xobs_[2] - l*Dz_);
+   double zred = interpinv(z);
+   double dang = cosmo_->Dang(zred);  // pour variation angle solide
+   double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
+   double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
+   double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
+   double corr = dnu/dnuref*pow(dangref/dang*dlum/dlumref,2.);
+   if(powlaw!=0.) corr *= pow((1.+redshref_)/(1.+zred),powlaw);
+   correction.push_back(corr);
+   if(lp_>0 && (l==0 || abs(l-Nz_/2)<2 || l==Nz_-1)) {
+     cout<<"l="<<l<<" z="<<z<<" red="<<zred
+         <<" da="<<dang<<" dlu="<<dlum<<" dred="<<dred
+         <<" dnu="<<dnu<<" -> corr="<<corr<<endl;
+   }
+ }
+
+ double sum=0., sum2=0., nsum=0.;
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
+   double a = lsigma*NorRand();
+   a = magnref*pow(10.,a);
+   // On met le meme tirage le long de Oz (indice k)
+   for(long l=0;l<Nz_;l++) {
+     int_8 ip = IndexR(i,j,l);
+     data_[ip] += a*correction[l];
+   }
+   sum += a; sum2 += a*a; nsum += 1.;
+ }
+
+ if(lp_>0 && nsum>1.) {
    sum /= nsum;
    sum2 = sum2/nsum - sum*sum;
@@ -1039 +1100 @@
  }
 }
+
+
+
+//-------------------------------------------------------------------
+//-------------------------------------------------------------------
+//--------------------- BRICOLO A NE PAS GARDER ---------------------
+//-------------------------------------------------------------------
+//-------------------------------------------------------------------
+int GeneFluct3D::ComputeSpectrum_bricolo(HistoErr& herr)
+// Compute spectrum from "T" and fill HistoErr "herr"
+// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
+// cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
+{
+ if(lp_>0) cout<<"--- ComputeSpectrum_bricolo ---"<<endl;
+ check_array_alloc();
+
+ if(herr.NBins()<0) return -1;
+ herr.Zero();
+
+ // Attention a l'ordre
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;  kx *= kx;
+   for(long j=0;j<Ny_;j++) {
+     if(i==0 && j==0) continue;
+     long jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;  ky *= ky;
+     for(long l=1;l<NCz_;l++) {
+       double kz = l*Dkz_;  kz *= kz;
+       double k = sqrt(kx+ky+kz);
+       double pk = MODULE2(T_(l,j,i));
+       herr.Add(k,pk);
+     }
+   }
+ }
+ herr.ToVariance();
+
+ // renormalize to directly compare to original spectrum
+ double norm = Vol_;
+ herr *= norm;
+
+ return 0;
+}
+
+int GeneFluct3D::ComputeSpectrum2D_bricolo(Histo2DErr& herr)
+{
+ if(lp_>0) cout<<"--- ComputeSpectrum2D_bricolo ---"<<endl;
+ check_array_alloc();
+
+ if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
+ herr.Zero();
+
+ // Attention a l'ordre
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;  kx *= kx;
+   for(long j=0;j<Ny_;j++) {
+     if(i==0 && j==0) continue;
+     long jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;  ky *= ky;
+     double kt = sqrt(kx+ky);
+     for(long l=1;l<NCz_;l++) {
+       double kz = l*Dkz_;
+       double pk = MODULE2(T_(l,j,i));
+       herr.Add(kt,kz,pk);
+     }
+   }
+ }
+ herr.ToVariance();
+
+ // renormalize to directly compare to original spectrum
+ double norm = Vol_;
+ herr *= norm;
+
+ return 0;
+}