Changeset 3773 in Sophya for trunk/Cosmo

Timestamp:

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

Author:

cmv

Message:

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

Location:

trunk/Cosmo/SimLSS

Files:

• cmvrvloscor.cc (modified) (8 diffs)
• genefluct3d.cc (modified) (8 diffs)

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -11 +11 @@
 #include "timing.h"
 #include "dvlist.h"
-#include "ntuple.h"
+#include "histos.h"
 #include "fabtcolread.h"
 
@@ -27 +27 @@
 {
  int nthread = 1;
- SophyaInit();
- InitTim();
+ int_8 nhfill = 100000;
 
  if(narg<=2) {usage(); return -1;}
@@ -35 +34 @@
  try {
  //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+
+ SophyaInit();
+ InitTim();
 
  cout<<"> read rho: "<<arg[1]<<endl;
@@ -52 +54 @@
  cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
 
+ double dmin = min(Dx,min(Dy,Dz));
+ double nmax = max(Nx,max(Ny,Nz));
+ cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
+ Histo hmpc(-dmin*nmax,dmin*nmax,4.*nmax);
+ nhfill = (int_8)Nx*(int_8)Ny*(int_8)Nz/nhfill; if(nhfill<=0) nhfill = 1;
+
  POutPersist pos("cmvrvloscor.ppf");
 
@@ -58 +66 @@
  TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
  rgen = 0.;
- TVector<r_8> R(Nz), Rdis(Nz);
+ TVector<GEN3D_TYPE> R(Nz), V(Nz);
+ TVector<r_8> Rdis(Nz);
 
  cout<<"> filling redshift distorted cube"<<endl;
+ int_8 nread = 0;
  for(int i=0;i<Nx;i++) {
    if(i%(Nx/10)==0) cout<<"i="<<i<<endl;
  for(int j=0;j<Ny;j++) {
-   for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
+   //for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
+   //for(int l=0;l<Nz;l++) V(l) = f3dv.Read(l,j,i);
+   f3dr.Read(j,i,R);
+   f3dv.Read(j,i,V);
    Rdis = 0.;
    for(int l=0;l<Nz;l++) {
-     double d = (1.+Zref) / Href * f3dv.Read(l,j,i);
+     double d = (1.+Zref) / Href * V(l);
+     if(nread%nhfill==0) hmpc.Add(d);
      double lpd = (double)l + d/Dz; // valeur du deplacee
      // on repartit proportionnelement au recouvrement sur 2 pixels
@@ -80 +94 @@
  }
  PrtTim(">>>> End filling redshift distorted cube");
+ pos.PutObject(hmpc,"hmpc");
 
  fluct3d.ReComputeFourier();
@@ -131 +146 @@
 openppf cmvrvloscor.ppf
 
+disp hmpc
+
 n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
 
@@ -144 +161 @@
 
 # proj selon kT (black), selon kZ (red)
-n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta cpts"
+n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
  */

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