Changeset 3267 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Jun 14, 2007, 7:06:35 PM (18 years ago)

Author:

cmv

Message:

cmv 14/06/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (7 diffs)
• cosmocalc.cc (modified) (1 diff)
• genefluct3d.cc (modified) (5 diffs)
• genefluct3d.h (modified) (2 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -39 +39 @@
      <<" -P : write cube in PPF format"<<endl
      <<" -V : compute variance from real space (for check, default: no)"<<endl
-     <<" -K : use power spectrum computation adapted for AGN (bidon!)"<<endl
      <<endl;
 }
@@ -86 +85 @@
  double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
  double powlaw_agn = 0.;
- bool killkz = false;
 
  // *** type de generation
@@ -106 +104 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2GKx:y:z:s:Z:M:A:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:A:")) != -1) {
   switch (c) {
   case 'a' :
@@ -141 +139 @@
     do_agn = true;
     sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
-    break;
-  case 'K' :
-    killkz = true;
     break;
   case 'V' :
@@ -284 +279 @@
  fluct3d.SetNThread(nthread);
  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
- fluct3d.SetObservator(zref,nz/2.);
+ fluct3d.SetObservator(zref,-nz/2.);
  fluct3d.SetCosmology(univ);
  fluct3d.SetGrowthFactor(growth);
@@ -530 +525 @@
    hpkrec.ReCenterBin();
    hpkrec.Show();
-   if(killkz) fluct3d.ComputeSpectrum_bricolo(hpkrec);
-     else     fluct3d.ComputeSpectrum(hpkrec);
+   fluct3d.ComputeSpectrum(hpkrec);
    tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
    PrtTim(">>>> End Computing final spectrum");
@@ -541 +535 @@
    hpkrec2.ReCenterBin(); hpkrec2.Zero();
    hpkrec2.Show();
-   if(killkz) fluct3d.ComputeSpectrum2D_bricolo(hpkrec2);
-     else     fluct3d.ComputeSpectrum2D(hpkrec2);
+   fluct3d.ComputeSpectrum2D(hpkrec2);
    {
    tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);

trunk/Cosmo/SimLSS/cosmocalc.cc

@@ -261 +261 @@
 
 //----------------------------------------------------------
-double CosmoCalc::Dloscom(double z)     /* Mpc */
+double CosmoCalc::Dloscom(double z)     /* Mpc comobile */
 {
   return _Dhubble * NInteg(z);
 }
 
-double CosmoCalc::Dtrcom(double z)     /* Mpc */
+double CosmoCalc::Dtrcom(double z)     /* Mpc comobile */
 {
   double v = NInteg(z);  // Zero curvature

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -69 +69 @@
 //                 Si redshref<0 alors redshref=0
 //     kredshref = indice (en double) correspondant a ce redshift
-//                 Si kredshref<0 alors kredshref=0
+//                 Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
 // 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.;
+ if(kredshref<0.) {
+   if(Nz_<=0) {
+     char *bla = "GeneFluct3D::SetObservator_Error: for kredshref<0 SetSize should be called first";
+     cout<<bla<<endl; throw ParmError(bla);
+   }
+   kredshref = Nz_/2.;
+ }
  redshref_  = redshref;
  kredshref_ = kredshref;
@@ -97 +103 @@
                 <<" D="<<dx<<","<<dy<<","<<dz<<endl;
  if(nx<=0 || dx<=0.) {
-   char *bla = "GeneFluct3D::setsize_Error: bad value(s";
+   char *bla = "GeneFluct3D::setsize_Error: bad value(s) for nn/dx";
    cout<<bla<<endl; throw ParmError(bla);
  }
@@ -228 +234 @@
  loscom_max_ = 0.;
  for(long i=0;i<=1;i++) {
-   double dx2 = xobs_[0] - i*Nx_*Dx_; dx2 *= dx2;
+   double dx2 = xobs_[0] - i*(Nx_-1)*Dx_; dx2 *= dx2;
    for(long j=0;j<=1;j++) {
-     double dy2 = xobs_[1] - j*Ny_*Dy_; dy2 *= dy2;
+     double dy2 = xobs_[1] - j*(Ny_-1)*Dy_; dy2 *= dy2;
      for(long k=0;k<=1;k++) {
-       double dz2 = xobs_[2] - k*Nz_*Dz_; dz2 *= dz2;
+       double dz2 = xobs_[2] - k*(Nz_-1)*Dz_; dz2 *= dz2;
        dz2 = sqrt(dx2+dy2+dz2);
        if(dz2>loscom_max_) loscom_max_ = dz2;
@@ -246 +252 @@
 
  // Fill the corresponding vectors for loscom and zred
+ // Be shure to have one dlc <loscom_min and one >loscom_max
  if(zinc<=0.) zinc = 0.01;
  for(double z=0.; ; z+=zinc) {
@@ -1121 +1128 @@
 }
 
-void GeneFluct3D::AddNoise2Real(double snoise)
+void GeneFluct3D::AddNoise2Real(double snoise,bool with_evol)
 // add noise to every pixels (meme les <=0 !)
 {
- if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<endl;
- check_array_alloc();
-
- for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = IndexR(i,j,l);
-   data_[ip] += snoise*NorRand();
- }
-}
-
-
-
-//-------------------------------------------------------------------
-//-------------------------------------------------------------------
-//--------------------- 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
+ if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<with_evol<<endl;
+ check_array_alloc();
+
  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;
-}
+     for(long l=0;l<Nz_;l++) {
+       int_8 ip = IndexR(i,j,l);
+       data_[ip] += snoise*NorRand();
+     }
+   }
+ }
+
+}
+

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -85 +85 @@
   double TurnMass2Flux(void);
   void AddAGN(double lfjy,double lsigma,double powlaw=0.);
-  void AddNoise2Real(double snoise);
+  void AddNoise2Real(double snoise,bool with_evol=false);
 
   void WriteFits(string cfname,int bitpix=FLOAT_IMG);
@@ -96 +96 @@
   void Print(void);
 
-//-------------------------------------------------------------------
-//--------------------- BRICOLO A NE PAS GARDER ---------------------
-//-------------------------------------------------------------------
-  int  ComputeSpectrum_bricolo(HistoErr& herr);
-  int  ComputeSpectrum2D_bricolo(Histo2DErr& herr);
 //-------------------------------------------------------------------