diff --git a/src/map_2d.cpp b/src/map_2d.cpp
index 13a280c..2650acb 100644
--- a/src/map_2d.cpp
+++ b/src/map_2d.cpp
@@ -1,248 +1,246 @@
 /**
  *
  * @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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #include "map_2d.hh"
 #include <fstream>
 #include <sstream>
 #include <fftw3.h>
 #include <omp.h>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 
 template <typename T>
 void Map2d<T>::rescaleHeights(Real C){
   /* compute actual RMS */
   for (UInt i = 0 ; i < n[0] ; ++i)
     for (UInt j = 0 ; j < n[1] ; ++j){
       this->at(i,j) *= C;
     }
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::setGridSize(UInt n0, UInt n1){
   n[0] = n0;
   n[1] = n1;
   data.resize(n[0]*n[1]);
   data.assign(n[0]*n[1],T(0.));
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::setGridSize(const UInt n[2]){
   setGridSize(n[0],n[1]);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::operator=(const T & e){
 
   for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
     this->at(i) = e;
   }
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 Map2d<T>::Map2d(UInt n0, UInt n1, Real L0, Real L1) {
   setGridSize(n0,n1);
   this->L[0] = L0;
   this->L[1] = L1;
   this->plan = NULL;
   generateFFTPlan();
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Map2d<T>::Map2d(const UInt * n, const Real *L) {
   setGridSize(n[0],n[1]);
   this->L[0] = L[0];
   this->L[1] = L[1];
   this->plan = NULL;
   generateFFTPlan();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 Map2d<T>::~Map2d() {
   destroyFFTPlan();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::generateFFTPlan() {
   // this should not do anything
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <>
 void Map2d<complex>::generateFFTPlan() {
   fftw_complex * in  = (fftw_complex*) &data[0];
   this->plan = fftw_plan_dft_2d(n[0], n[1], in, in, FFTW_FORWARD, FFTW_ESTIMATE);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<typename T>
 void Map2d<T>::destroyFFTPlan() {
   if (this->plan) {
     fftw_destroy_plan(this->plan);
     this->plan = NULL;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::FFTTransform(int nthreads){
   SURFACE_FATAL("impossible situation");
   throw;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 template <typename T>
 void Map2d<T>::FFTTransform(Map2d<std::complex<Real> > & output, int nthreads)const{
   output = *this;
   output.FFTTransform(nthreads);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::FFTITransform(Map2d<std::complex<Real> > & output, int nthreads)const{
   output = *this;
   output.FFTITransform(nthreads);
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Map2d<T>::FFTITransform(Map2d<Real> & output, int nthreads) const{
   const UInt * n = this->sizes();
   const Real * Ls = this->getLs();
   Map2d<std::complex<Real> > _output(n,Ls);
-  //  _output = *static_cast< >(this);
-  const Map2d<T> & tmp = *this;
-  _output.copy(tmp);
+  _output = *this;
   _output.FFTITransform(nthreads);
   for (UInt i = 0; i < n[0]*n[1]; i++) {
     output(i) = _output(i).real();
   }
 }
 /* -------------------------------------------------------------------------- */
 
 
 template <>
 void Map2d<complex>::FFTTransform(int nthreads){
   STARTTIMER("FFTW");
   fftw_execute(this->plan);
   STOPTIMER("FFTW");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::FFTITransform(int nthreads){
   SURFACE_FATAL("impossible situation");
   throw;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 template <>
 void Map2d<complex>::FFTITransform(int nthreads){
 
 #pragma omp parallel for
   for (UInt i = 0; i < n[0]*n[1]; ++i) {
     data[i] = std::conj(data[i]);
   }
   FFTTransform(nthreads);
 #pragma omp parallel for
   for (UInt i = 0; i < n[0]*n[1]; ++i) {
     data[i] = std::conj(data[i]);
     data[i] *= 1./(n[0]*n[1]);
   }
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::dumpToTextFile(std::string filename) {
 
   UInt n1 = this->size(0);
   UInt n2 = this->size(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 < n1; i++) {
     for (UInt j=0; j< n2; j++) {
       file << i << " " << j << " ";
       this->writeValue(file,this->at(i,j));
       file << std::endl;
     }
     file << std::endl;
   }
   file.close();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Map2d<T>::writeValue(std::ostream & os, const T & val) {
   os << val;
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void Map2d<complex>::writeValue(std::ostream & os, const complex & val) {
 
   os << std::scientific << val.real();
   os << " " << val.imag();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template class Map2d<complex>;
 template class Map2d<Real>;
 template class Map2d<int>;
 template class Map2d<unsigned int>;
 template class Map2d<bool>;
 template class Map2d<unsigned long>;
 
 __END_TAMAAS__
diff --git a/src/map_2d.hh b/src/map_2d.hh
index 93ed787..6dd4270 100644
--- a/src/map_2d.hh
+++ b/src/map_2d.hh
@@ -1,577 +1,591 @@
 /**
  *
  * @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 "my_vector.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T>
 class Map2d {
 
 
   /* ------------------------------------------------------------------------ */
   /* Typedefs                                                                 */
   /* ------------------------------------------------------------------------ */
 
   //  class iterator;
 
   /* ------------------------------------------------------------------------ */
   /* 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();
 
   //! generate fft plan
   void generateFFTPlan();
   //! destroy fft plan
   void destroyFFTPlan();
   //! perform fft transform
   void FFTTransform(Map2d<std::complex<Real> > & output,int nthreads = 1) const;
   //! perform inverse fft transform
   void FFTITransform(Map2d<std::complex<Real> > & output, int nthreads = 1) const;
   //! perform inverse fft transform
   void FFTITransform(Map2d<Real> & output, int nthreads = 1) const;
   //! perform fft transform
   void FFTTransform(int nthreads = 1);
   // //! inverse fft transform
   void FFTITransform(int nthreads = 1);
   //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);
   //! copy method
   template <template<typename> class surf, typename T2>
   void copy(const surf<T2> & s);
   // set all the members of the map to a given value
   void operator=(const T & e);
   // multiply the map by scalar
   template <typename T2>
   void operator*=(const T2 & e);
   // divide the map by scalar
   template <typename T2>
   void operator/=(const T2 & e);
   // multiply the map by scalar
   void operator*=(const T & e);
   // divide the map by scalar
   void operator/=(const T & e);
   // multiply the map with another map, component by component
   template <template<typename> class surf, typename T2>
   void operator*=(const surf<T2> & e);
   // divide the map with another map, component by component
   template <template<typename> class surf, typename T2>
   void operator/=(const surf<T2> & e);
   // 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 void operator+=(const Map2d<T> & s);
   // access to the i,j-th value in the surface
   inline void operator+=(const T & val);
   // access to the i,j-th value in the surface
   inline void operator-=(const T & val);
   // 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 n[dir];};
   //! get both lateral number of discretization points
   inline const UInt* sizes()const {return n;};
   //! get the lateral size of the sample
   inline const Real * getLs() const {return L;};
   //! get the lateral size of the sample
   inline Real getL(UInt dir) const {return 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 data.getPtr();};
 
   //! Write heights array to text file usable with plot programs
   template <bool consider_zeros=true>
   void dumpToTextFile(std::string filename) const;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 
 protected:
 
   //! number of points in each direction
   UInt n[2];
   //! size in each direction
   Real L[2];
   //! internal storage of data
   MyVector<T> data;
   //! fftw plan
   fftw_plan plan;
 };
 
 /* -------------------------------------------------------------------------- */
 //#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){
 
   if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
 					     << i << " " << j << " "
 					     << n[0] << " " << n[1] << " "
 					     << _at(i,j,n));
   return data[_at(i,j,n)];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::at(UInt i,UInt j)const {
 
   if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
 					     << i << " " << j << " "
 					     << n[0] << " " << n[1] << " "
 					     << _at(i,j,n));
   return data[_at(i,j,n)];
 }
 /* -------------------------------------------------------------------------- */
 #undef _at
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T & Map2d<T>::at(UInt i){
   return data[i];
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline const T & Map2d<T>::at(UInt i) const{
   return data[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 data[i];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline const T & Map2d<T>::operator()(int i) const{
   return data[i];
 }
 
 /* -------------------------------------------------------------------------- */
 // template <typename T>
 // inline T & Map2d<T>::operator[](std::pair<UInt,UInt> c){
 //   return this->at(c.first,c.second);
 // }
 
 // /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::operator-=(const T & e){
 
   for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
     this->at(i) -= e;
   }
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Map2d<T>::operator+=(const T & e){
 
   for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
     this->at(i) += e;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Map2d<T>::operator+=(const Map2d<T> & s){
 
   for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
     this->at(i) += s.at(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 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)n[0] && !periodic) continue;
     if (res.first < 0 && periodic)     res.first += n[0];
     if (res.first >= (int)n[0] && periodic) res.first -= 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)n[1] && !periodic) continue;
     if (res.second < 0 && periodic)     res.second += n[1];
     if (res.second >= (int)n[1] && periodic) res.second -= n[1];
 
     if (ii == 0 && jj == 0) continue;
 
     {
       UInt i = res.first;
       UInt j = res.second;
 
       if (i+j*n[0] >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
 					     << i << " " << j << " "
 					     << n[0] << " " << n[1] << " "
 					     << i+j*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 (n[0] != n1[0] ||
       n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices")
 
   for (UInt i = 0 ; i < n[0] ; ++i)
     for (UInt j = 0 ; j < 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);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 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){
+  SURFACE_FATAL("Cannot copy a complex surface into incompatible type");
+}
+template <>
+inline void Map2d<complex>::operator=(const Map2d<complex> & s){
+  this->copy(s);
+}
+
 /* -------------------------------------------------------------------------- */
 template <typename T>
 template <typename T2>
 void Map2d<T>::operator*=(const T2 & e){
 
   UInt sz = this->n[0]*this->n[1];
 #pragma omp parallel for
   for (unsigned int i = 0; i < sz; ++i) {
     this->at(i) *= e;
   }
 
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <typename T2>
 void Map2d<T>::operator/=(const T2 & e){
 
   UInt sz = this->n[0]*this->n[1];
 #pragma omp parallel for
   for (unsigned int i = 0; i < sz; ++i) {
     this->at(i) /= e;
   }
 
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::operator*=(const T & e){
   
   UInt sz = this->n[0]*this->n[1];
 #pragma omp parallel for
   for (unsigned int i = 0; i < sz; ++i) {
     this->at(i) *= e;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Map2d<T>::operator/=(const T & e){
   
   UInt sz = this->n[0]*this->n[1];
 #pragma omp parallel for
   for (unsigned int i = 0; i < sz; ++i) {
     this->at(i) /= e;
   }
   
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <template<typename> class surf, typename T2>
 void Map2d<T>::operator*=(const surf<T2> & e){
 
   UInt sz = this->n[0]*this->n[1];
 #pragma omp parallel for
   for (unsigned int i = 0; i < sz; ++i) {
     this->at(i) *= e.at(i);
   }
 
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <template<typename> class surf, typename T2>
 void Map2d<T>::operator/=(const surf<T2> & e){
 
   UInt sz = this->n[0]*this->n[1];
 #pragma omp parallel for
   for (unsigned int i = 0; i < sz; ++i) {
     this->at(i) /= e.at(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 __END_TAMAAS__
 
 #endif /* __MAP_2D_HH__ */