diff --git a/src/map_2d.hh b/src/map_2d.hh
index d769dd9..42ec432 100644
--- a/src/map_2d.hh
+++ b/src/map_2d.hh
@@ -1,416 +1,416 @@
 /**
  *
  * @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 __MAP_2D_HH__
 #define __MAP_2D_HH__
 /* -------------------------------------------------------------------------- */
 #include <vector>
 #include <iostream>
 #include <fftw3.h>
 #include "tamaas.hh"
 #include <fstream>
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T>
 class Map2d : public Grid<T, 2> {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   Map2d(UInt n1, UInt n2, Real L0 = 1., Real L1 = 1.);
   Map2d(const UInt * n, const Real * L);
 
   virtual ~Map2d();
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
 
   template <bool periodic>
   std::vector<std::pair<UInt,UInt> > getNeighborIndexes(UInt i, UInt j);
   template <bool periodic>
   std::vector<std::pair<UInt,UInt> > getNeighborIndexes(std::pair<UInt,UInt> i);
 
   // //! return an iterator over the neighbors
   // iterator beginNeighborhood();
   // //! return the end iterator over the neighbors
   // iterator endNeighborhood();
 
   //compute compute component by component product and put result in object
   template <typename T1, typename T2>
   void fullProduct(Map2d<T1> & s1,Map2d<T2> &s2);
   //! set the grid size
   void setGridSize(const UInt n[2]);
   //! set the grid size
   void setGridSize(UInt n0,UInt n1);
   //! rescale height profile
   void rescaleHeights(Real factor);
   //! Write heights array to text file usable with plot programs
   void dumpToTextFile(std::string filename);
   //! set to scalar value
   void operator=(const T & val);
   //! copy method
   template <template<typename> class surf, typename T2>
   void copy(const surf<T2> & s);
   // copy the map passed as parameter
   template <template<typename> class surf, typename T2>
   void operator=(const surf<T2> & s);
   // temporary hack
   void operator=(const Map2d<complex> & s);
   // copy the map passed as const parameter
   template <template<typename> class surf, typename T2>
   void operator=(surf<T2> & s);
 
 /* -------------------------------------------------------------------------- */
 
 private:
 
   //! write an entry to a stringstream
   void writeValue(std::ostream & os, const T & val);
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 
 public:
 
   // access to the i,j-th value in the surface
   inline T & operator()(UInt i,UInt j);
   // access to the i,j-th value in the surface
   inline T & operator()(std::pair<UInt,UInt> i);
   // access to the i,j-th value in the surface
   inline const T & operator()(UInt i,UInt j)const;
   // access to the i,j-th value in the surface
   inline const T & operator()(std::pair<UInt,UInt> i)const;
   // access to the i,j-th value in the surface
   inline T & at(UInt i,UInt j);
   // access to the i-th value in the array
   inline T & at(UInt i);
   // access to the i,j-th value in the surface
   inline T & at(std::pair<UInt,UInt> i);
   // access to the i,j-th value in the surface
   inline const T & at(UInt i,UInt j) const;
   // access to the i,j-th value in the surface
   inline const T & at(std::pair<UInt,UInt> i) const;
   // access to the i,j-th value in the surface
   inline const T & at(UInt i) const;
   // access to the i-th value in the array
   inline T & operator()(int i);
   // access to the i-th value in the array
   //  inline T & operator[](std::pair<UInt,UInt> i);
   // access to the i-th value in the array
   inline const T & operator()(int i) const;
   //! get lateral number of discretization points
   inline UInt size(UInt dir)const {return this->n[dir];};
   //! get both lateral number of discretization points
   inline const UInt* sizes()const {return this->n;};
   //! get the lateral size of the sample
   inline const Real * getLs() const {return this->L;};
   //! get the lateral size of the sample
   inline Real getL(UInt dir) const {return this->L[dir];};
   //! get the lateral size of the sample
   inline void setL(Real L[2]){this->L[0] = L[0];this->L[1] = L[1];};
   //! get the lateral size of the sample
   inline void setL(Real L0,Real L1){this->L[0] = L0;this->L[1] = L1;};
 
   const T * getInternalData() const{return this->data.getPtr();};
 
   //! Write heights array to text file usable with plot programs
   template <bool consider_zeros=true>
   void dumpToTextFile(std::string filename) const;
 };
 
 /* -------------------------------------------------------------------------- */
 //#define _at(I,J,N) I+J*N[1]
 #define _at(I,J,N) (I)*(N[1])+(J)
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline T & Map2d<T>::at(UInt i,UInt j){
   return this->Grid<T, 2>::operator()(i, j, 0);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::at(UInt i,UInt j)const {
   return this->Grid<T, 2>::operator()(i, j, 0);
 }
 /* -------------------------------------------------------------------------- */
 #undef _at
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T & Map2d<T>::at(UInt i){
   return this->Grid<T, 2>::operator()(i);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline const T & Map2d<T>::at(UInt i) const{
   return this->Grid<T, 2>::operator()(i);
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 template <typename T>
 inline T & Map2d<T>::operator()(UInt i,UInt j){
   return this->at(i,j);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T & Map2d<T>::at(std::pair<UInt,UInt> i){
   return this->at(i.first,i.second);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::at(std::pair<UInt,UInt> i)const {
   return this->at(i.first,i.second);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T & Map2d<T>::operator()(std::pair<UInt,UInt> i){
   return this->at(i.first,i.second);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::operator()(UInt i,UInt j)const{
   return this->at(i,j);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::operator()(std::pair<UInt,UInt> i)const {
   return this->at(i.first,i.second);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T & Map2d<T>::operator()(int i){
   return this->Grid<T, 2>::operator()(i);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::operator()(int i) const{
   return this->Grid<T, 2>::operator()(i);
 }
 
 /* -------------------------------------------------------------------------- */
 // template <typename T>
 // inline T & Map2d<T>::operator[](std::pair<UInt,UInt> c){
 //   return this->at(c.first,c.second);
 // }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <bool periodic>
 std::vector<std::pair<UInt,UInt> > Map2d<T>::getNeighborIndexes(UInt i, UInt j){
 
   std::vector<std::pair<UInt,UInt> > vres;
 
 
   std::pair<int,int> res;
 
   for (int ii = -1; ii < 2; ++ii) {
     res.first = ii+static_cast<int>(i);
     if (res.first < 0 && !periodic)    continue;
     if (res.first >= (int)this->n[0] && !periodic) continue;
     if (res.first < 0 && periodic)     res.first += this->n[0];
     if (res.first >= (int)this->n[0] && periodic) res.first -= this->n[0];
 
     for (int jj = -1; jj < 2; ++jj) {
 
     res.second = jj+static_cast<int>(j);
     if (res.second < 0 && !periodic)    continue;
     if (res.second >= (int)this->n[1] && !periodic) continue;
     if (res.second < 0 && periodic)     res.second += this->n[1];
     if (res.second >= (int)this->n[1] && periodic) res.second -= this->n[1];
 
     if (ii == 0 && jj == 0) continue;
 
     {
       UInt i = res.first;
       UInt j = res.second;
 
       if (i+j*this->n[0] >= this->n[0]*this->n[1]) SURFACE_FATAL("out of bound access "
 					     << i << " " << j << " "
 					     << this->n[0] << " " << this->n[1] << " "
 					     << i+j*this->n[0]);
     }
     vres.push_back(res);
     }
   }
   return vres;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <bool periodic>
 std::vector<std::pair<UInt,UInt> > Map2d<T>
 ::getNeighborIndexes(std::pair<UInt,UInt> i){
   return this->getNeighborIndexes<periodic>(i.first,i.second);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 template <typename T1, typename T2>
 void Map2d<T>::fullProduct(Map2d<T1> & s1,Map2d<T2> & s2){
 
   auto n1 = s1.sizes();
   auto n2 = s2.sizes();
 
   if (n1[0] != n2[0] ||
       n1[1] != n2[1]) SURFACE_FATAL("fullproduct of non compatible matrices");
 
   if (this->n[0] != n1[0] ||
       this->n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices")
 
   for (UInt i = 0 ; i < this->n[0] ; ++i)
     for (UInt j = 0 ; j < this->n[1] ; ++j){
       this->at(i,j) = s1(i,j)* s2(i,j);
     }
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 template <bool consider_zeros>
 void Map2d<T>::dumpToTextFile(std::string filename) const {
 
   UInt m = this->n[0];
   UInt n = this->n[1];
   std::cout << "Writing surface file " << filename << std::endl;
   std::ofstream file(filename.c_str());
   if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename);
 
   file.precision(15);
   for (UInt i = 0; i < m; i++) {
     for (UInt j=0; j< n; j++) {
       if (!consider_zeros &&
 	  std::norm(this->at(i,j)) < zero_threshold*zero_threshold)
 	continue;
 
       file << i << " " << j << " ";
       file << std::scientific << this->at(i,j);
       file << std::endl;
     }
     file << std::endl;
   }
   file.close();
 }
 
 
 // template <typename T>
 // class Map2d<T>::iterator {
 
 // public:
 
 //   iterator(std::vector<T> & data,
 // 	   std::queue<std::pair<UInt,UInt> > & indices):
 //     data(data),indices(indices){
 //   }
 
 //   inline void operator++(){
 //     indices.pop();
 //   }
 
 //   inline T & operator *(){
 //     return data[indices.front];
 //   }
 
 //   bool operator != (Map2d<T>::iterator & it){
 //     return (indices.front == it.indices.front);
 //   }
 
 // private:
 
 //   std::vector<T> & data;
 //   std::queue<std::pair<UInt,UInt> >  indices;
 
 //   bool finished;
 // };
 
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 template <template<typename> class surf, typename T2>
 void Map2d<T>::copy(const surf<T2> & s){
   //this->setGridSize(s.sizes());
   //UInt m = this->sizes()[0];
   //UInt n = this->sizes()[1];
   //#pragma omp parallel for
   //for (UInt i = 0; i < m*n; i++) {
   //  this->at(i) = s(i);
   //}
   this->Grid<T, 2>::copy(s);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <template<typename> class surf, typename T2>
 void Map2d<T>::operator=(surf<T2> & s){
   this->copy(s);
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 template <template<typename> class surf, typename T2>
 void Map2d<T>::operator=(const surf<T2> & s){
   this->copy(s);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
-void Map2d<T>::operator=(const Map2d<complex> & s){
+void Map2d<T>::operator=(const Map2d<complex> &){
   SURFACE_FATAL("Cannot copy a complex surface into incompatible type");
 }
 template <>
 inline void Map2d<complex>::operator=(const Map2d<complex> & s){
   this->copy(s);
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif /* __MAP_2D_HH__ */