diff --git a/src/surface_statistics.hh b/src/surface_statistics.hh
index ca53377..03d3bae 100644
--- a/src/surface_statistics.hh
+++ b/src/surface_statistics.hh
@@ -1,721 +1,722 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef __SURFACE_STATISTICS_HH__
 #define __SURFACE_STATISTICS_HH__
 /* -------------------------------------------------------------------------- */
 #include "surface.hh"
 #include "fftransform.hh"
 #include "fft_plan_manager.hh"
 #include <fftw3.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 
 class SurfaceStatistics {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   SurfaceStatistics(){};
   virtual ~SurfaceStatistics(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
 
   inline static Real computeMaximum(const Surface<Real> & s);
   inline static Real computeMinimum(const Surface<Real> & s);
   inline static UInt computeAverage(const Surface<UInt> & s,
 				    UInt epsilon=0){SURFACE_FATAL("TO IMPLEMENT");};
 
   inline static unsigned long computeAverage(const Surface<unsigned long> & s,
 					     unsigned long epsilon=0){SURFACE_FATAL("TO IMPLEMENT");};
 
 
   inline static Real computeAverage(const Surface<Real> & s,
 				    Real epsilon= -1.*std::numeric_limits<Real>::max()){
     return computeAverage<Real>(s,epsilon);
   }
 
   template <typename T>
   inline static T computeAverage(const Surface<T> & s,
 				 T epsilon= -1.*std::numeric_limits<T>::max());
   inline static Real computeStdev(const Surface<Real> & s, Real epsilon=-1. * std::numeric_limits<Real>::max());
 
   inline static Real computeSkewness(const Surface<Real> & surface, Real avg, Real std);
   inline static Real computeKurtosis(const Surface<Real> & surface, Real avg, Real std);
   inline static Real computeRMSSlope(const Surface<Real> & surface);
   inline static Real computeSpectralRMSSlope(Surface<Real> & surface);
   inline static Real computeSpectralMeanCurvature(Surface<Real> & surface);
   inline static Real computeSpectralStdev(Surface<Real> & surface);
   //  inline static void computeStatistics(Surface<Real> & surface);
   inline static std::vector<Real> computeMoments(Surface<Real> & surface);
   inline static Real computeAnalyticFractalMoment(UInt order,UInt k1, UInt k2, Real hurst, Real rms, Real L);
   //! compute perimeter
   inline static Real computePerimeter(const Surface<Real> & surface, Real eps = 1e-10);
   //! get real contact area
   inline static Real computeContactArea(const Surface<Real> & surface);
   //! get real contact area
   inline static Real computeContactAreaRatio(const Surface<Real> & surface);
 
 
   template <bool consider_zeros=false>
   inline static std::vector<int>
   computeDistribution(const Surface<Real> & surface,
 		      unsigned int n_bins,
 		      Real min_value=std::numeric_limits<Real>::max(),
 		      Real max_value=-1*std::numeric_limits<Real>::max());
 
   inline static std::vector<complex>
   computeSpectralDistribution(const Surface<Real> & surface,
 			      unsigned int n_bins,
 			      UInt cutoff=std::numeric_limits<UInt>::max());
 
 
 
   //! compute average amplitude
   inline static Real computeAverageAmplitude(const Surface<Real> & surface);
   inline static Real computeAverageAmplitude(const SurfaceComplex<Real> & surface);
 
 
   //! sum all heights
   inline static Real computeSum(const Surface<Real> & surface);
   //! compute auto correlation
   inline static Surface<Real> computeAutocorrelation(Surface<Real> & surface);
   //! compute powerspectrum
   inline static SurfaceComplex<Real> computePowerSpectrum(Surface<Real> & surface);
 
 
 };
 
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 //#include "surface_statistics_inline_impl.cc"
 
 /* -------------------------------------------------------------------------- */
 
 template <bool consider_zeros>
 std::vector<int>
 SurfaceStatistics::computeDistribution(const Surface<Real> & surface,
 				       unsigned int n_bins,
 				       Real min_value,Real max_value){
 
   if (n_bins == 0) SURFACE_FATAL("badly set number of bins");
 
   Real min_heights;
   Real max_heights;
 
   std::vector<int> bins;
   bins.resize(n_bins);
   bins.assign(n_bins,0);
 
 
   min_heights = min_value;
   max_heights = max_value;
 
   if (max_value == -1.*std::numeric_limits<Real>::max() )
     max_heights = computeMaximum(surface);
   if (min_value == std::numeric_limits<Real>::max() )
     min_heights = computeMinimum(surface);
 
   int counter = 0;
   const int n = surface.size();
   for (int i = 0 ; i < n*n ; ++i)
     {
       Real value = surface(i);
       if (!consider_zeros && value < zero_threshold) continue;
       if (value > max_heights || value < min_heights)
 	SURFACE_FATAL("out of bound value: " << min_heights << " <? " << value
 		      << " <? " << max_heights << " : bounds probably badly chosen");
       Real ind = (value - min_heights) / (max_heights - min_heights);
       unsigned int index = ind*n_bins;
       if (index == n_bins) --index;
       ++bins[index];
       ++counter;
     }
   if (consider_zeros && counter != n*n)
     SURFACE_FATAL("problem in counting points");
 
   //  Real offset = (max_heights - min_heights)/n_bins;
   //Real x = 0;
   //  Real f1 = 0, f2 = 0;
   //  for (unsigned int i = 0 ; i < n_bins ; ++i)
   //  {
   //    x  = i*offset+min_heights;
   //f1 = bins[i]/Real(n*n);
   //f2 = bins[i]/Real(counter);
   //}
   Real min_non_zero = std::numeric_limits<Real>::max();
 #pragma omp parallel for reduction(min: min_non_zero)
   for (int i = 0 ; i < n*n ; ++i) {
     Real value = surface(i);
     if (value < min_non_zero && value > 1e-14)
       min_non_zero = value;
   }
   return bins;
 }
 /* -------------------------------------------------------------------------- */
 std::vector<complex>
 SurfaceStatistics::computeSpectralDistribution(const Surface<Real> & surface,
 					       unsigned int n_bins,
 					       UInt cutoff){
 
   if (n_bins == 0) SURFACE_FATAL("badly set number of bins");
 
 
   const int n = surface.size();
 
 
   Real min_value = SurfaceStatistics::computeMinimum(surface);
   Real max_value = SurfaceStatistics::computeMaximum(surface);
 
   std::vector<complex> distrib_fft(n_bins);
   const int N = n_bins;
   Real L = max_value - min_value;
 
 #pragma omp parallel for
   for (int k = 0 ; k < N ; ++k) distrib_fft[k] = 0.;
 
 #pragma omp parallel for
   for (int k = 1 ; k < N/2 ; ++k) {
     for (int i = 0 ; i < n*n ; ++i) {
       Real value = surface(i);
       Real angle = -2.*M_PI*k*(value-min_value)/L;
       complex phase = std::polar(1.,angle);
       distrib_fft[k]   += phase;
       distrib_fft[N-k] += phase;
     }
   }
 
   for (int i = 0 ; i < n*n ; ++i) {
     Real value = surface(i);
     Real angle = -2.*M_PI*N/2*(value-min_value)/L;
     complex phase = std::polar(1.,angle);
     distrib_fft[N/2]   += phase;
   }
 
 #pragma omp parallel for
   for (int k = 1 ; k < N ; ++k) distrib_fft[k] /= n*n;
 
   distrib_fft[0] = std::polar(1.,0.);
 
 #pragma omp parallel for
   for (int k = 1 ; k < N/2 ; ++k) {
     if (k >= (int)cutoff) {
       distrib_fft[k] = 0.;
       distrib_fft[N-k] = 0.;
     }
   }
 
   //proceed to inverse fourier transform
 
 #pragma omp parallel for
   for (int k = 0 ; k < N ; ++k) {
     distrib_fft[k] = std::conj(distrib_fft[k]);
   }
 
   fftw_complex * in  = (fftw_complex*)&distrib_fft[0];
   fftw_plan p = fftw_plan_dft_1d(N, in, in ,FFTW_FORWARD, FFTW_ESTIMATE);
   fftw_execute(p);
   fftw_destroy_plan(p);
 
   std::vector<complex> distrib;
   distrib.resize(N);
 
 #pragma omp parallel for
   for (int k = 0 ; k < N ; ++k) {
     distrib[k] = std::conj(distrib_fft[k]);
   }
 
   return distrib;
 
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeAverageAmplitude(const Surface<Real> & surface){
   Real average_amplitude = 0;
   UInt n = surface.size();
 #pragma omp parallel for collapse(2), reduction(+: average_amplitude)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j)
       average_amplitude += std::abs(surface(i,j));
 
   average_amplitude *= 1./n/n;
   return average_amplitude;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAverageAmplitude(const SurfaceComplex<Real> & surface){
   Real average_amplitude = 0;
   UInt n = surface.size();
 #pragma omp parallel for collapse(2), reduction(+: average_amplitude)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j)
       average_amplitude += std::abs(surface(i,j));
 
   average_amplitude *= 1./n/n;
   return average_amplitude;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeContactArea(const Surface<Real> & surface) {
   Real L = surface.getL();
   return SurfaceStatistics::computeContactAreaRatio(surface)*L*L;
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeContactAreaRatio(const Surface<Real> & surface) {
   UInt area = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+:area)
   for (UInt i = 0; i<n; i++) {
     for (UInt j = 0; j<n; j++) {
       if (surface (i,j) > 1e-10)
         area++;
     }
   }
   return area/Real(n*n);
 }
 
 
 
 
 /* -------------------------------------------------------------------------- */
 
 Surface<Real> SurfaceStatistics::computeAutocorrelation(Surface<Real> & surface){
   SurfaceComplex<Real> correlation_heights = computePowerSpectrum(surface);
   Surface<Real> acf(correlation_heights.size(), correlation_heights.getL());
   FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(acf, correlation_heights);
   plan.backwardTransform();
   return acf;
 }
 /* -------------------------------------------------------------------------- */
 
 
 SurfaceComplex<Real> SurfaceStatistics::computePowerSpectrum(Surface<Real> & surface){
   SurfaceComplex<Real> power_spectrum(surface.size(), surface.getL());
   FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(surface, power_spectrum);
   plan.forwardTransform();
 
   UInt n = power_spectrum.size();
   Real factor = 1./(n*n);
 
 #pragma omp parallel for collapse(2)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j)
       power_spectrum(i,j) *= factor;
 
   power_spectrum.makeItRealBySquare();
   return power_spectrum;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 T SurfaceStatistics::computeAverage(const Surface<T> & surface, T epsilon){
 
   Real avg = 0.0;
   UInt cpt = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: avg, cpt)
   for (UInt i = 0; i < n*n; i++) {
     Real val = surface(i);
     if (std::abs(val) > epsilon){
       avg += val;
       ++cpt;
     }
   }
   avg /= cpt;
   return avg;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeSkewness(const Surface<Real> & surface,
                                         Real avg, Real std){
   Real skw = 0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: skw)
   for (UInt i = 0; i < n*n; i++) {
     Real tmp = surface(i) - avg;
     skw += (tmp*tmp*tmp);
   }
   skw /= (n*std*std*std);
   return skw;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeKurtosis(const Surface<Real> & surface,
                                         Real avg, Real std){
 
   Real kur=0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: kur)
   for (UInt i = 0; i < n*n; i++) {
     Real tmp = surface.size() - avg;
     kur += (tmp*tmp*tmp*tmp);
   }
   kur /= (n*std*std*std*std);
   return kur-3;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMaximum(const Surface<Real> & surface){
   Real _max = -1. * std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(max: _max)
   for (UInt i = 0; i < n*n; i++) {
     _max = std::max(surface(i),_max);
   }
   return _max;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMinimum(const Surface<Real> & surface){
   Real _min = std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(min: _min)
   for (UInt i = 0; i < n*n; i++) {
     _min = std::min(surface(i),_min);
   }
   return _min;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeRMSSlope(const Surface<Real> & surface){
 
   UInt n = surface.size();
   Surface<Real> grad(n,1.);
 
   for (UInt i = 0, m = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       UInt k = i+n*j;
 
       //      xl = ((i==0)?(k+n-1):(k-1));
       UInt xh = ((i==(n-1))?(k-n+1):(k+1));
       //      yl = ((j==0)?(k+n*n-n):(k-n));
       UInt yh = ((j==(n-1))?(k-n*n+n):(k+n));
 
 
       Real gx = surface(xh) - surface(k);
       Real gy = surface(yh) - surface(k);
       grad(m++) = gx*gx + gy*gy;
     }
   }
 
   Real avg = SurfaceStatistics::computeAverage(grad);
   Real slp = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */
   Real scale_factor = surface.getL()/n;
   return slp/scale_factor;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralRMSSlope(Surface<Real> & surface){
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_slopes = 0;
+  // /!\ commented code was causing race conditions with dummy variable in GridHermitian
 #pragma omp parallel for collapse(2), reduction(+:rms_slopes)
   for (UInt i = 1 ; i < n/2 ; ++i)
     for (UInt j = 1 ; j < n/2 ; ++j){
-      rms_slopes += 1.*((i*i)+(j*j))*power(i,j).real();
-      rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
-      rms_slopes += 1.*((i*i)+(j*j))*power(n-i,j).real();
-      rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real();
+      rms_slopes += 2.*((i*i)+(j*j))*power(i,j).real();
+      //rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
+      rms_slopes += 2.*((i*i)+(j*j))*power(n-i,j).real();
+      //rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real();
     }
 
   // external line (i or j == 0)
 #pragma omp parallel for reduction(+:rms_slopes)
   for (UInt i = 1 ; i < n/2 ; ++i){
     rms_slopes += 1.*(i*i)*power(i,0).real();
     rms_slopes += 1.*(i*i)*power(n-i,0).real();
-    rms_slopes += 1.*(i*i)*power(0,i).real();
-    rms_slopes += 1.*(i*i)*power(0,n-i).real();
+    rms_slopes += 2.*(i*i)*power(0,i).real();
+    //rms_slopes += 1.*(i*i)*power(0,n-i).real();
   }
 
   //internal line (i or j == n/2)
 #pragma omp parallel for reduction(+:rms_slopes)
   for (UInt i = 1; i < n/2; ++i) {
     rms_slopes += 0.5*(i*i+n*n/4)*power(i,n/2).real();
     rms_slopes += 0.5*(i*i+n*n/4)*power(n-i,n/2).real();
-    rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,i).real();
-    rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real();
+    rms_slopes += 1.*(i*i+n*n/4)*power(n/2,i).real();
+    //rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real();
   }
 
   {
     //i == n/2 j == n/2
     rms_slopes += 0.25*(n*n/2)*power(n/2,n/2).real();
     //i == 0 j == n/2
     rms_slopes += 1.*(n*n/4)*power(0,n/2).real();
     //i == n/2 j == 0
     rms_slopes += 1.*(n*n/4)*power(n/2,0).real();
   }
   rms_slopes = 2.*M_PI/n*sqrt(rms_slopes);
   Real scale_factor = surface.getL()/n;
   return rms_slopes/scale_factor;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeStdev(const Surface<Real> & surface, Real epsilon){
   UInt n = surface.size();
   UInt cpt = 0;
   Real std = 0.0;
   Real avg = 0.0;
 
 #pragma omp parallel for reduction(+:avg, std, cpt)
   for (UInt i = 0; i < n*n; i++) {
     Real val = surface(i);
     if (val > epsilon){
       avg += val;
       std += val*val;
       ++cpt;
     }
   }
   avg /= cpt;
   return sqrt(std/cpt - avg*avg);
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeSpectralStdev(Surface<Real> & surface){
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real rms_heights = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+:rms_heights)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j){
       rms_heights += power(i,j).real();
     }
   rms_heights = sqrt(rms_heights);
   return rms_heights;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralMeanCurvature(Surface<Real> & surface){
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_curvatures = 0;
 #pragma omp parallel for collapse(2), reduction(+:rms_curvatures)
   for (UInt i = 1 ; i < n/2 ; ++i)
     for (UInt j = 1 ; j < n/2 ; ++j){
       Real coeff = ((i*i)+(j*j));
       coeff *= coeff;
       rms_curvatures += 1.*coeff*power(i,j).real();
       rms_curvatures += 1.*coeff*power(n-i,n-j).real();
       rms_curvatures += 1.*coeff*power(n-i,j).real();
       rms_curvatures += 1.*coeff*power(i,n-j).real();
     }
 
   // external line (i or j == 0)
 #pragma omp parallel for reduction(+:rms_curvatures)
   for (UInt i = 1 ; i < n/2 ; ++i){
     Real coeff = i*i;
     coeff *= coeff;
     rms_curvatures += 4./3.*coeff*power(i,0).real();
     rms_curvatures += 4./3.*coeff*power(n-i,0).real();
     rms_curvatures += 4./3.*coeff*power(0,i).real();
     rms_curvatures += 4./3.*coeff*power(0,n-i).real();
   }
 
   //internal line (i or j == n/2)
 #pragma omp parallel for reduction(+:rms_curvatures)
   for (UInt i = 1; i < n/2; ++i) {
     Real coeff = (i*i+n*n/4);
     coeff *= coeff;
     rms_curvatures += 0.5*coeff*power(i,n/2).real();
     rms_curvatures += 0.5*coeff*power(n-i,n/2).real();
     rms_curvatures += 0.5*coeff*power(n/2,i).real();
     rms_curvatures += 0.5*coeff*power(n/2,n-i).real();
   }
 
   {
     //i == n/2 j == n/2
     rms_curvatures += 0.25*(n*n/2)*(n*n/2)*power(n/2,n/2).real();
     //i == 0 j == n/2
     rms_curvatures += 2./3.*(n*n/4)*(n*n/4)*power(0,n/2).real();
     //i == n/2 j == 0
     rms_curvatures += 2./3.*(n*n/4)*(n*n/4)*power(n/2,0).real();
   }
 
   rms_curvatures *= 9./4.;
   rms_curvatures = M_PI/n/n * sqrt(rms_curvatures);
 
   return rms_curvatures;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computePerimeter(const Surface<Real> & surface, Real eps) {
   int perimeter = 0;
   UInt n = surface.size();
 
 
 #pragma omp parallel for collapse(2), reduction(+:perimeter)
   for (UInt i = 0 ; i < n ; ++i) {
     for (UInt j = 0 ; j < n ; ++j) {
       if ( j > 0 && ( ( surface(i,(j-1)) < eps && surface(i,j) > eps ) ||
                       ( surface(i,(j-1)) > eps && surface(i,j) < eps ) ) )
         perimeter++;
       if ( i > 0 && ( ( surface(j,(i-1)) < eps && surface(j,i) > eps ) ||
                       ( surface(j,(i-1)) > eps && surface(j,i) < eps )) )
         perimeter++;
     }
   }
 
   return perimeter/Real(n);
 }
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeSum(const Surface<Real> & surface){
   Real sum = 0.;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: sum)
   for (UInt i = 0 ; i < n*n ; ++i)
     sum += surface(i);
 
   return sum;
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<Real> SurfaceStatistics::computeMoments(Surface<Real> & surface){
 
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real L = surface.getL();
   UInt n = surface.size();
   Real k_factor = 2.*M_PI/L;
 
   std::vector<std::vector<Real> > moment;
   moment.resize(6);
   for (int i = 0; i < 5; i++)
     moment[i].resize(5);
 
   for (UInt i = 0; i < n/2; i++) {
     Real i0 = pow(Real(i)*k_factor,0);
     Real i2 = pow(Real(i)*k_factor,2);
     Real i4 = pow(Real(i)*k_factor,4);
     Real im0 = pow(-Real(i)*k_factor,0);
     Real im2 = pow(-Real(i)*k_factor,2);
     Real im4 = pow(-Real(i)*k_factor,4);
 
     for (UInt j = 0; j < n/2; j++) {
       Real Phi = power(i,j).real();
 
       Real Phi_10;
       Real Phi_01;
       Real Phi_11;
       if (i==0 && j==0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = 0;
       }
       else if (i==0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = power(i,(n-j)).real();
       }
       else if (j==0) {
         Phi_10 = power(n-i,j).real();
         Phi_11 = 0;
         Phi_01 = 0;
       }
       else {
         Phi_10 = power(n-i,j).real();
         Phi_11 = power(n-i,(n-j)).real();
         Phi_01 = power(i,(n-j)).real();
       }
       Real j0 = pow(Real(j)*k_factor,0);
       Real j2 = pow(Real(j)*k_factor,2);
       Real j4 = pow(Real(j)*k_factor,4);
       Real jm0 = pow(-Real(j)*k_factor,0);
       Real jm2 = pow(-Real(j)*k_factor,2);
       Real jm4 = pow(-Real(j)*k_factor,4);
 
       moment[0][0] += i0*j0*Phi + im0*j0*Phi_10 + i0*jm0*Phi_01 + im0*jm0*Phi_11;
       moment[0][2] += i0*j2*Phi + im0*j2*Phi_10 + i0*jm2*Phi_01 + im0*jm2*Phi_11;
       moment[2][0] += i2*j0*Phi + im2*j0*Phi_10 + i2*jm0*Phi_01 + im2*jm0*Phi_11;
       moment[2][2] += i2*j2*Phi + im2*j2*Phi_10 + i2*jm2*Phi_01 + im2*jm2*Phi_11;
       moment[0][4] += i0*j4*Phi + im0*j4*Phi_10 + i0*jm4*Phi_01 + im0*jm4*Phi_11;
       moment[4][0] += i4*j0*Phi + im4*j0*Phi_10 + i4*jm0*Phi_01 + im4*jm0*Phi_11;
     }
   }
 
   std::vector<Real> mean_moments;
   //m2
   Real mean_m2 = 0.5*(moment[0][2]+moment[2][0]);
   mean_moments.push_back(mean_m2);
   //m4
   Real mean_m4 = (3*moment[2][2]+moment[0][4]+moment[4][0])/3.;
   mean_moments.push_back(mean_m4);
   return mean_moments;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAnalyticFractalMoment(UInt order,UInt k1, UInt k2, Real hurst,Real A, Real L){
   Real xi = k2/k1;
   Real exp = - 2.*hurst;
   Real _k1 = 2.*M_PI/L*k1;
 
   Real Cq;
   switch (order){
   case 0: Cq = 2.*M_PI;break;
   case 2: Cq = M_PI;break;
   case 4: Cq = 3./4.*M_PI;break;
   default: SURFACE_FATAL("invalid order");
   }
   exp = order - 2.*hurst;
   Real res = A * Cq * pow(_k1,exp) / exp * ( pow(xi, exp) - 1 );
 
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif /* __SURFACE_STATISTICS_HH__ */