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

Timestamp:

Jun 20, 2007, 6:50:48 PM (18 years ago)

Author:

cmv

Message:

intro evolution bruit avec redshift cmv 20/06/2007

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -31 +31 @@
 //-------------------------------------------------------
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
-  : T_(T) , Nx_(0) , Ny_(0) , Nz_(0) , array_allocated_(false) , lp_(0)
-  , redshref_(-999.) , kredshref_(0.) , cosmo_(NULL) , growth_(NULL)
-  , loscom_ref_(-999.), loscom_min_(-999.), loscom_max_(-999.)
+  : Nx_(0) , Ny_(0) , Nz_(0)
+  , lp_(0)
+  , array_allocated_(false) , T_(T)
+  , cosmo_(NULL) , growth_(NULL)
+  , redsh_ref_(-999.), kredsh_ref_(0.), dred_ref_(-999.)
+  , loscom_ref_(-999.), dtrc_ref_(-999.), dlum_ref_(-999.), dang_ref_(-999.)
+  , nu_ref_(-999.), dnu_ref_ (-999.)
+  , loscom_min_(-999.), loscom_max_(-999.)
   , loscom2zred_min_(0.) , loscom2zred_max_(0.)
 {
@@ -76 +81 @@
  if(kredshref<0.) {
    if(Nz_<=0) {
-     char *bla = "GeneFluct3D::SetObservator_Error: for kredshref<0 SetSize should be called first";
+     char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 SetSize should be called first";
      cout<<bla<<endl; throw ParmError(bla);
    }
    kredshref = Nz_/2.;
  }
- redshref_  = redshref;
- kredshref_ = kredshref;
+ redsh_ref_  = redshref;
+ kredsh_ref_ = kredshref;
  if(lp_>0)
-   cout<<"--- GeneFluct3D::SetObservator zref="<<redshref_<<" kref="<<kredshref_<<endl;
+   cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
 }
 
@@ -206 +211 @@
  if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
 
- if(cosmo_ == NULL || redshref_<0.) {
+ if(cosmo_ == NULL || redsh_ref_<0.) {
    char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
    cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ // La distance angulaire/luminosite/Dnu au pixel de reference
+ dred_ref_ = Dz_/(cosmo_->Dhubble()/cosmo_->E(redsh_ref_));
+ loscom_ref_ = cosmo_->Dloscom(redsh_ref_);
+ dtrc_ref_ = cosmo_->Dtrcom(redsh_ref_);
+ dlum_ref_ = cosmo_->Dlum(redsh_ref_);
+ dang_ref_ = cosmo_->Dang(redsh_ref_);
+ nu_ref_   = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
+ dnu_ref_  = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
+ if(lp_>0) {
+   cout<<"...reference pixel redshref="<<redsh_ref_
+       <<", dredref="<<dred_ref_
+       <<", nuref="<<nu_ref_ <<" GHz"
+       <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
+       <<"   dlosc="<<loscom_ref_<<" Mpc com"
+       <<", dtrc="<<dtrc_ref_<<" Mpc com"
+       <<", dlum="<<dlum_ref_<<" Mpc"
+       <<", dang="<<dang_ref_<<" Mpc"<<endl;
  }
 
@@ -214 +238 @@
  // 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_;
  xobs_[1] = Ny_/2.*Dy_;
- xobs_[2] = kredshref_*Dz_ - loscom_ref_;
+ xobs_[2] = kredsh_ref_*Dz_ - loscom_ref_;
 
  // L'observateur est-il dans le cube?
@@ -229 +252 @@
  if(!obs_in_cube) loscom_min_ = -xobs_[2];
 
+ // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
+ if(loscom_min_<=1.e-50)
+   for(int i=0;i<50;i++)
+     cout<<"ATTENTION TOUTES LES PARTIES DU CODE NE MARCHENT PAS POUR UN OBSERVATEUR DANS LE CUBE"<<endl;
+ // TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED TO BE FIXED
+
+
  // Find MAXIMUM los com distance to the observer:
  // ou que soit positionne l'observateur, la distance
@@ -234 +264 @@
  loscom_max_ = 0.;
  for(long i=0;i<=1;i++) {
-   double dx2 = xobs_[0] - i*(Nx_-1)*Dx_; dx2 *= dx2;
+   double dx2 = DXcom(i*(Nx_-1)); dx2 *= dx2;
    for(long j=0;j<=1;j++) {
-     double dy2 = xobs_[1] - j*(Ny_-1)*Dy_; dy2 *= dy2;
+     double dy2 = DYcom(j*(Ny_-1)); dy2 *= dy2;
      for(long k=0;k<=1;k++) {
-       double dz2 = xobs_[2] - k*(Nz_-1)*Dz_; dz2 *= dz2;
+       double dz2 = DZcom(k*(Nz_-1)); dz2 *= dz2;
        dz2 = sqrt(dx2+dy2+dz2);
        if(dz2>loscom_max_) loscom_max_ = dz2;
@@ -245 +275 @@
  }
  if(lp_>0) {
-   cout<<"...zref="<<redshref_<<" kzref="<<kredshref_<<" losref="<<loscom_ref_<<" Mpc\n"
+   cout<<"...zref="<<redsh_ref_<<" kzref="<<kredsh_ref_<<" 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;
+       <<" loscom_min="<<loscom_min_<<" loscom_max="<<loscom_max_<<" Mpc (com)"<<endl;
  }
 
@@ -312 +342 @@
    fwrt.WriteKey("DKY",Dky_," Mpc^-1");
    fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
-   fwrt.WriteKey("ZREF",redshref_," reference redshift");
-   fwrt.WriteKey("KZREF",kredshref_," reference redshift on z axe");
+   fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
+   fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
    fwrt.Write(R_);
  } catch (PThrowable & exc) {
@@ -364 +394 @@
    R_.Info()["DY"] = (r_8)Dy_;
    R_.Info()["DZ"] = (r_8)Dz_;
-   R_.Info()["ZREF"] = (r_8)redshref_;
-   R_.Info()["KZREF"] = (r_8)kredshref_;
+   R_.Info()["ZREF"] = (r_8)redsh_ref_;
+   R_.Info()["KZREF"] = (r_8)kredsh_ref_;
    POutPersist pos(cfname.c_str());
    if(write_real) pos << PPFNameTag("rgen")  << R_;
@@ -432 +462 @@
      <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
  cout<<"          (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
- cout<<"Redshift "<<redshref_<<" for z axe at k="<<kredshref_<<endl;
+ cout<<"Redshift "<<redsh_ref_<<" for z axe at k="<<kredsh_ref_<<endl;
 }
 
@@ -1044 +1074 @@
 //               de l'histo ne sont pas vraiment ce que l'on veut
 // --- Limitations actuelle du code:
-// . les AGN sont supposes en loi de puissance IDENTIQUE pour tout theta,phi
+// . les AGN sont supposes evoluer avec la meme loi de puissance 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
@@ -1050 +1080 @@
 //     et pour un cube peu epais / distance observateur 
 // --- Parametres de la routine:
-// 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"
+// llfy : c'est le <log10(S)> du flux depose par les AGN
+//        dans l'angle solide du pixel elementaire de reference du cube
+// lsigma : c'est le sigma de la distribution des log10(S)
+// powlaw : c'est la pente de la 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)
@@ -1064 +1095 @@
   check_array_alloc();
 
- if(cosmo_ == NULL || redshref_<0.| loscom2zred_.size()<1) {
+ if(cosmo_ == NULL || redsh_ref_<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;
- }
+ // Le flux des AGN en Jy et en mass solaire
+ double fagnref = pow(10.,lfjy)*(dnu_ref_*1.e9); // Jy.Hz = W/m^2
+ double magnref = FluxHI2Msol(fagnref*Jansky2Watt_cst,dlum_ref_); // Msol
+ if(lp_>0)
+   cout<<"Au pixel de ref: 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);
+   // F(nu) = F(1.4GHz)*(nu GHz/1.4 Ghz)^p = F(1.4GHz)*(1/(1+z))^p , car nu = 1.4 GHz/(1+z)
+   magnref *= pow(1/(1.+redsh_ref_),powlaw);
    if(lp_>0) cout<<" powlaw="<<powlaw<<"  -> change magnref to "<<magnref<<" Msol"<<endl;
  }
@@ -1091 +1116 @@
  vector<double> correction;
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
+ long nzmod = ((Nz_>10)?Nz_/10:1);
  for(long l=0;l<Nz_;l++) {
-   double z = fabs(xobs_[2] - l*Dz_);
+   double z = fabs(DZcom(l));
    double zred = interpinv(z);
-   double dang = cosmo_->Dang(zred);  // pour variation angle solide
+   double dtrc = cosmo_->Dtrcom(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);
+   // on a: Mass ~ DNu * Dlum^2 / Dtrcom^2
+   double corr = dnu/dnu_ref_*pow(dtrc_ref_/dtrc*dlum/dlum_ref_,2.);
+   // F(nu) = F(1.4GHz)*(nu GHz/1.4 Ghz)^p = F(1.4GHz)*(1/(1+z))^p , car nu = 1.4 GHz/(1+z)
+   if(powlaw!=0.) corr *= pow((1.+redsh_ref_)/(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;
+   if(lp_>0 && (l==0 || l==Nz_-1 || l%nzmod==0)) {
+     cout<<"l="<<l<<" z="<<z<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
+         <<" dtrc="<<dtrc<<" dlum="<<dlum
+         <<" -> cor="<<corr<<endl;
    }
  }
@@ -1134 +1162 @@
  check_array_alloc();
 
+ vector<double> correction;
+ InterpFunc *intercor = NULL;
+
+ if(with_evol) {
+   // Sigma_Noise(en mass) :
+   //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
+   //      Flim ~ sqrt(DNu) * sqrt(nlobe)   en W/m^2
+   //      Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe)  en Msol
+   //                                    nlobe ~ 1/Dtrcom^2
+   //      Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
+   if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
+     char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
+     cout<<bla<<endl; throw ParmError(bla);
+   }
+   InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
+   long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
+   for(long i=0;i<nsz;i++) {
+     double d = interpinv.X(i);
+     double zred = interpinv(d);
+     double dtrc = cosmo_->Dtrcom(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 = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
+     if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
+       cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
+           <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
+     correction.push_back(corr);
+   }
+   intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
+ } 
+
+ double dx2=0., dy2=0., dz2=0.;
  for(long i=0;i<Nx_;i++) {
+   dx2 = DXcom(i); dx2 *= dx2;
    for(long j=0;j<Ny_;j++) {
+     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);
+       }
        int_8 ip = IndexR(i,j,l);
-       data_[ip] += snoise*NorRand();
-     }
-   }
- }
-
-}
-
+       data_[ip] += snoise*corr*NorRand();
+     }
+   }
+ }
+
+ if(intercor!=NULL) delete intercor;
+}