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

Timestamp:

May 9, 2010, 12:47:58 AM (15 years ago)

Author:

cmv

Message:

speed-up FITS file read, cmv 09/05/2010

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -118 +118 @@
  Nz_ = nz;  if(Nz_ <= 0) Nz_ = Nx_;
  N_.resize(0); N_.push_back(Nx_); N_.push_back(Ny_); N_.push_back(Nz_);
- NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
+ NRtot_ = (int_8)Nx_*(int_8)Ny_*(int_8)Nz_; // nombre de pixels dans le survey
  NCz_ =  Nz_/2 +1;
  NTz_ = 2*NCz_;
@@ -766 +766 @@
 // Take into account the real and complexe conjugate coefficients
 // because we want a realization of a real data in real space
-// On ecrit que: P(k_x,k_y,k_z) = P(-k_x,-k_y,-k_z)
+// On ecrit que: conj(P(k_x,k_y,k_z)) = P(-k_x,-k_y,-k_z)
 //               avec k_x = i, -k_x = N_x - i  etc...
 {
@@ -773 +773 @@
 
  // 1./ Le Continu et Nyquist sont reels
- long nreal = 0;
+ int_8 nreal = 0;
  for(long kk=0;kk<2;kk++) {
    long k=0;  // continu
@@ -795 +795 @@
 
  // a./ les lignes et colonnes du continu et de nyquist
- long nconj1 = 0;
+ int_8  nconj1 = 0;
  for(long kk=0;kk<2;kk++) {
    long k=0;  // continu
@@ -823 +823 @@
 
  // b./ les lignes et colonnes hors continu et de nyquist
- long nconj2 = 0;
+ int_8  nconj2 = 0;
  for(long kk=0;kk<2;kk++) {
    long k=0;  // continu
@@ -840 +840 @@
  if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
 
- if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-nreal<<endl;
+ if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(NRtot_-nconj1-nconj2)-nreal<<endl;
 
  return nreal+nconj1+nconj2;
@@ -897 +897 @@
 //-------------------------------------------------------------------
 void GeneFluct3D::ToVelLoS(void)
+// Le spectre Pk doit etre (dRho/Rho)(k)
 {
  double zpk = compute_pk_redsh_ref_;
@@ -912 +913 @@
        double kz = Kz(l);
        double k2 = kt2 + kz*kz;
-       if(k2<=0.) continue;
-       T_(l,j,i) *= complex<double>(0.,dpsd*kz/k2);
-     }
-   }
- }
+       if(l==0) k2 = 0.; else k2 = dpsd*kz/k2;
+       T_(l,j,i) *= complex<double>(0.,k2);
+     }
+   }
+ }
+
 }