Changeset 3129 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Jan 11, 2007, 7:24:49 PM (19 years ago)

Author:

cmv

Message:

passage int en long aux endroits importants cmv 11/01/07

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (4 diffs)
• genefluct3d.cc (modified) (22 diffs)
• genefluct3d.h (modified) (3 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -42 +42 @@
  //-----------------------------------------------------------------
  // *** Survey definition
- int nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
- //int nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
- //int nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
+ long nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
+ //long nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
+ //long nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
 
  // *** Cosmography definition   (WMAP)
@@ -90 +90 @@
     break;
   case 'x' :
-    sscanf(optarg,"%d,%lf",&nx,&dx);
+    sscanf(optarg,"%ld,%lf",&nx,&dx);
     break;
   case 'y' :
-    sscanf(optarg,"%d,%lf",&ny,&dy);
+    sscanf(optarg,"%ld,%lf",&ny,&dy);
     break;
   case 'z' :
-    sscanf(optarg,"%d,%lf",&nz,&dz);
+    sscanf(optarg,"%ld,%lf",&nz,&dz);
     break;
   case 's' :
@@ -230 +230 @@
 
  cout<<"\n--- Checking realization spectra"<<endl;
- int nhprof = int(fluct3d.GetKmax()/dkmin+0.5);
+ long nhprof = long(fluct3d.GetKmax()/dkmin+0.5);
  HProf hpkgen(0.,fluct3d.GetKmax(),nhprof);
  hpkgen.ReCenterBin();
@@ -307 +307 @@
 
  cout<<"\n--- Convert Galaxies number to HI mass"<<endl;
- int nhmdndm = int( (lschmax-lschmin+1.)*100. + 0.5);
+ long nhmdndm = long( (lschmax-lschmin+1.)*100. + 0.5);
  Histo hmdndm(lschmin,lschmax,nhmdndm);
  sch.SetOutValue(1);

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -42 +42 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::SetSize(int nx,int ny,int nz,double dx,double dy,double dz)
+void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
 {
  if(nx<=0 || dx<=0. ) {
@@ -126 +126 @@
  // On tient compte du pb de normalisation de FFTW3
  double sntot = sqrt((double)NRtot_);
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    data[ip] = NorRand()/sntot;
@@ -140 +140 @@
  if(lp>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
  T_(0,0,0) = complex<r_8>(0.);  // on coupe le continue et on l'initialise
- int lmod = Nx_/10; if(lmod<1) lmod=1;
- for(int i=0;i<Nx_;i++) {
-   int ii = (i>Nx_/2) ? Nx_-i : i;
+ long lmod = Nx_/10; if(lmod<1) lmod=1;
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
    if(lp>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
-   for(int j=0;j<Ny_;j++) {
-     int jj = (j>Ny_/2) ? Ny_-j : j;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
      double ky = jj*Dky_;  ky *= ky;
-     for(int l=0;l<NCz_;l++) {
+     for(long l=0;l<NCz_;l++) {
        double kz = l*Dkz_;  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
@@ -209 +209 @@
  // --- On remplit avec une realisation
  if(lp>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
- int lmod = Nx_/10; if(lmod<1) lmod=1;
- for(int i=0;i<Nx_;i++) {
-   int ii = (i>Nx_/2) ? Nx_-i : i;
+ long lmod = Nx_/10; if(lmod<1) lmod=1;
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
    if(lp>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
-   for(int j=0;j<Ny_;j++) {
-     int jj = (j>Ny_/2) ? Ny_-j : j;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
      double ky = jj*Dky_;  ky *= ky;
-     for(int l=0;l<NCz_;l++) {
+     for(long l=0;l<NCz_;l++) {
        double kz = l*Dkz_;  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
@@ -242 +242 @@
 }
 
-int GeneFluct3D::manage_coefficients(void)
+long GeneFluct3D::manage_coefficients(void)
 // Take into account the real and complexe conjugate coefficients
 // because we want a realization of a real data in real space
@@ -249 +249 @@
 
  // 1./ Le Continu et Nyquist sont reels
- int nreal = 0;
- for(int kk=0;kk<2;kk++) {
-   int k=0;  // continu
+ long nreal = 0;
+ for(long kk=0;kk<2;kk++) {
+   long k=0;  // continu
    if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
-   for(int jj=0;jj<2;jj++) {
-     int j=0;
+   for(long jj=0;jj<2;jj++) {
+     long j=0;
      if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
-     for(int ii=0;ii<2;ii++) {
-       int i=0;
+     for(long ii=0;ii<2;ii++) {
+       long i=0;
        if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
        size_t ip = k+NCz_*(j+Ny_*i);
@@ -271 +271 @@
 
  // a./ les lignes et colonnes du continu et de nyquist
- int nconj1 = 0;
- for(int kk=0;kk<2;kk++) {
-   int k=0;  // continu
+ long nconj1 = 0;
+ for(long kk=0;kk<2;kk++) {
+   long k=0;  // continu
    if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
-   for(int jj=0;jj<2;jj++) { // selon j
-     int j=0;
+   for(long jj=0;jj<2;jj++) { // selon j
+     long j=0;
      if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
-     for(int i=1;i<(Nx_+1)/2;i++) {
+     for(long i=1;i<(Nx_+1)/2;i++) {
        size_t ip = k+NCz_*(j+Ny_*i);
        size_t ip1 = k+NCz_*(j+Ny_*(Nx_-i));
@@ -285 +285 @@
      }
    }
-   for(int ii=0;ii<2;ii++) {
-     int i=0;
+   for(long ii=0;ii<2;ii++) {
+     long i=0;
      if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
-     for(int j=1;j<(Ny_+1)/2;j++) {
+     for(long j=1;j<(Ny_+1)/2;j++) {
        size_t ip = k+NCz_*(j+Ny_*i);
        size_t ip1 = k+NCz_*((Ny_-j)+Ny_*i);
@@ -299 +299 @@
 
  // b./ les lignes et colonnes hors continu et de nyquist
- int nconj2 = 0;
- for(int kk=0;kk<2;kk++) {
-   int k=0;  // continu
+ long nconj2 = 0;
+ for(long kk=0;kk<2;kk++) {
+   long k=0;  // continu
    if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
-   for(int j=1;j<(Ny_+1)/2;j++) {
+   for(long j=1;j<(Ny_+1)/2;j++) {
      if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
-     for(int i=1;i<Nx_;i++) {
+     for(long i=1;i<Nx_;i++) {
        if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
        size_t ip = k+NCz_*(j+Ny_*i);
@@ -325 +325 @@
 {
  double s2 = 0.;
- for(int l=0;l<NCz_;l++)
-   for(int j=0;j<Ny_;j++)
-     for(int i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));
+ for(long l=0;l<NCz_;l++)
+   for(long j=0;j<Ny_;j++)
+     for(long i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));
 
  double s20 = 0.;
- for(int j=0;j<Ny_;j++)
-   for(int i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));
+ for(long j=0;j<Ny_;j++)
+   for(long i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));
 
  double s2n = 0.;
  if(Nz_%2==0)
-   for(int j=0;j<Ny_;j++)
-     for(int i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));
+   for(long j=0;j<Ny_;j++)
+     for(long i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));
 
  return 2.*s2 -s20 -s2n;
@@ -353 +353 @@
  if(lp>1) PrtTim("--- FilterByPixel: before ---");
 
- for(int i=0;i<Nx_;i++) {
-   int ii = (i>Nx_/2) ? Nx_-i : i;
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_ *Dx_/2;
    double pkx = pixelfilter(kx); 
-   for(int j=0;j<Ny_;j++) {
-     int jj = (j>Ny_/2) ? Ny_-j : j;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
      double ky = jj*Dky_ *Dy_/2;
      double pky =  pixelfilter(ky); 
-     for(int l=0;l<NCz_;l++) {
+     for(long l=0;l<NCz_;l++) {
        double kz = l*Dkz_ *Dz_/2;
        double pkz =  pixelfilter(kz);
@@ -401 +401 @@
  // On corrige du pb de la normalisation de FFTW3
  double v = (double)NRtot_;
- for(int i=0;i<Nx_;i++)
-   for(int j=0;j<Ny_;j++)
-     for(int l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
+ for(long i=0;i<Nx_;i++)
+   for(long j=0;j<Ny_;j++)
+     for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
 
  Tcontent_ = 3;
@@ -421 +421 @@
 
  // Attention a l'ordre
- for(int i=0;i<Nx_;i++) {
-   int ii = (i>Nx_/2) ? Nx_-i : i;
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
-   for(int j=0;j<Ny_;j++) {
-     int jj = (j>Ny_/2) ? Ny_-j : j;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
      double ky = jj*Dky_;  ky *= ky;
-     for(int l=0;l<NCz_;l++) {
+     for(long l=0;l<NCz_;l++) {
        double kz = l*Dkz_;  kz *= kz;
        double k = sqrt(kx+ky+kz);
@@ -452 +452 @@
 
  double *data = (double *) (&T_(0,0,0));
- int dnx = int(R/Dx_+0.5); if(dnx<=0) dnx = 1;
- int dny = int(R/Dy_+0.5); if(dny<=0) dny = 1;
- int dnz = int(R/Dz_+0.5); if(dnz<=0) dnz = 1;
+ long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
+ long dny = long(R/Dy_+0.5); if(dny<=0) dny = 1;
+ long dnz = long(R/Dz_+0.5); if(dnz<=0) dnz = 1;
  cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
 
  double sum=0., sum2=0., r2 = R*R; size_t nsum=0;
 
- for(int i=dnx;i<Nx_-dnx;i+=dnx) {
-   for(int j=dny;j<Ny_-dny;j+=dny) {
-     for(int l=dnz;l<Nz_-dnz;l+=dnz) {
+ for(long i=dnx;i<Nx_-dnx;i+=dnx) {
+   for(long j=dny;j<Ny_-dny;j+=dny) {
+     for(long l=dnz;l<Nz_-dnz;l+=dnz) {
        double s=0.; size_t n=0;
-       for(int ii=i-dnx;ii<=i+dnx;ii++) {
+       for(long ii=i-dnx;ii<=i+dnx;ii++) {
          double x = (ii-i)*Dx_; x *= x;
-         for(int jj=j-dny;jj<=j+dny;jj++) {
+         for(long jj=j-dny;jj<=j+dny;jj++) {
            double y = (jj-j)*Dy_; y *= y;
-           for(int ll=l-dnz;ll<=l+dnz;ll++) {
+           for(long ll=l-dnz;ll<=l+dnz;ll++) {
              double z = (ll-l)*Dz_; z *= z;
              if(x+y+z>r2) continue;
@@ -503 +503 @@
 
  size_t nbad = 0;
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    double v = data[ip];
@@ -521 +521 @@
  rm = rs2 = 0.;
 
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    double v = data[ip];
@@ -545 +545 @@
 
  size_t nbad = 0;
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    double v = data[ip];
@@ -562 +562 @@
  double *data = (double *) (&T_(0,0,0));
 
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    data[ip] += 1.;
@@ -595 +595 @@
             <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
  double sum = 0.;
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    data[ip] *= dn; // par coherence on multiplie aussi les <=0
@@ -616 +616 @@
 
  double sum = 0.;
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    double v = data[ip];
@@ -646 +646 @@
 
  double sum = 0.;
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    double v = data[ip];
    if(v>0.) {
-     int ngal = int(v+0.1);
+     long ngal = long(v+0.1);
      data[ip] = 0.;
-     for(int i=0;i<ngal;i++) {
+     for(long i=0;i<ngal;i++) {
        double m = massdist.RandomInterp();  // massdist.Random();
        if(axeslog) m = pow(10.,m);
@@ -675 +675 @@
  double *data = (double *) (&T_(0,0,0));
 
- for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    size_t ip = l+NTz_*(j+Ny_*i);
    data[ip] += snoise*NorRand();

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -23 +23 @@
 
   void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
-  void SetSize(int nx,int ny,int nz,double dx,double dy,double dz);  // Mpc
+  void SetSize(long nx,long ny,long nz,double dx,double dy,double dz);  // Mpc
 
   inline void SetZ(double z) {pkz_.SetZ(z);}
@@ -62 +62 @@
 
 protected:
-  int manage_coefficients(void);
+  long manage_coefficients(void);
   double compute_power_carte(void);
   inline double pixelfilter(double x)
@@ -68 +68 @@
 
 
-  int Nx_,Ny_,Nz_,SzK_[3];
-  int NCz_,NTz_;
+  long Nx_,Ny_,Nz_;
+  sa_size_t SzK_[3];
+  long NCz_,NTz_;
   size_t NRtot_;
   double Dx_,Dy_,Dz_;