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

Timestamp:

Oct 17, 2007, 4:44:53 PM (18 years ago)

Author:

cmv

Message:

modifs sur VarianceFrReal pour calcul; sima8 dans epace reel, cmv 17/10/2007

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -1021 +1021 @@
 int_8 GeneFluct3D::VarianceFrReal(double R,double& var)
 // Recompute MASS variance in spherical top-hat (rayon=R)
+// Par definition: SigmaR^2 = <(M-<M>)^2>/<M>^2
+//                 ou M = masse dans sphere de rayon R
 {
   if(lp_>0) cout<<"--- VarianceFrReal R="<<R<<endl;
  check_array_alloc();
 
- 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;
+ long dnx = long(R/Dx_)+1; if(dnx<=0) dnx = 1;
+ long dny = long(R/Dy_)+1; if(dny<=0) dny = 1;
+ long dnz = long(R/Dz_)+1; if(dnz<=0) dnz = 1;
  if(lp_>0) cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
 
- double sum=0., sum2=0., r2 = R*R; int_8 nsum=0;
-
- 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.; int_8 n=0;
+ double sum=0., sum2=0., sn=0., r2 = R*R;
+ int_8 nsum=0;
+
+ for(long i=dnx;i<Nx_-dnx;i+=2*dnx) {
+   for(long j=dny;j<Ny_-dny;j+=2*dny) {
+     for(long l=dnz;l<Nz_-dnz;l+=2*dnz) {
+       double m=0.; int_8 n=0;
        for(long ii=i-dnx;ii<=i+dnx;ii++) {
          double x = (ii-i)*Dx_; x *= x;
@@ -1044 +1047 @@
              if(x+y+z>r2) continue;
              int_8 ip = IndexR(ii,jj,ll);
-             s += 1.+data_[ip];
+             m += 1.+data_[ip];  // 1+drho/rho
              n++;
            }
          }
        }
-       if(n>0) {sum += s; sum2 += s*s; nsum++;}
-       //cout<<i<<","<<j<<","<<l<<" n="<<n<<" s="<<s<<" sum="<<sum<<" sum2="<<sum2<<endl;
+       if(n>0) {sum += m; sum2 += m*m; nsum++; sn += n;}
+       //cout<<i<<","<<j<<","<<l<<" n="<<n<<" m="<<m<<" sum="<<sum<<" sum2="<<sum2<<endl;
      }
    }
@@ -1056 +1059 @@
 
  if(nsum<=1) {var=0.; return nsum;}
-
  sum /= nsum;
  sum2 = sum2/nsum - sum*sum;
- if(lp_>0) cout<<"...nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
-
- var = sum2/(sum*sum);  // <dM>^2/<M>^2
- if(lp_>0) cout<<"...sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
+ sn /= nsum;
+ if(lp_>0) cout<<"...<n>="<<sn<<", nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
+ var = sum2/(sum*sum);  // <dM^2>/<M>^2
+ if(lp_>0) cout<<"...sigmaR^2 = <(M-<M>)^2>/<M>^2 = "<<var
+               <<" -> sigmaR = "<<sqrt(var)<<endl;
 
  return nsum;