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

Timestamp:

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

Author:

cmv

Message:

cmv 14/06/2007

File:

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

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();
+     }
+   }
+ }
+
+}
+