diff --git a/lectures/pybind_lecture/CMakeLists.txt b/lectures/pybind_lecture/CMakeLists.txt
index b90b6ee..7c426f5 100644
--- a/lectures/pybind_lecture/CMakeLists.txt
+++ b/lectures/pybind_lecture/CMakeLists.txt
@@ -1,19 +1,20 @@
 cmake_minimum_required(VERSION 2.8.12)
 project(example)
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS 1)
 set(CMAKE_CXX_STANDARD 14)
 
-find_package(pybind11 REQUIRED)  # or add_subdirectory(pybind11)
+# find_package(pybind11 REQUIRED)
+add_subdirectory(pybind11)
 
 add_library(pyexample MODULE pyexample.cpp)
 target_link_libraries(pyexample PRIVATE pybind11::module)
 set_target_properties(pyexample PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}"
                                            SUFFIX "${PYTHON_MODULE_EXTENSION}")
 
 file(
   COPY ${CMAKE_CURRENT_SOURCE_DIR}/test.py
   DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/
   )
 
 
diff --git a/lectures/pybind_lecture/pybind11 b/lectures/pybind_lecture/pybind11
new file mode 160000
index 0000000..e2b884c
--- /dev/null
+++ b/lectures/pybind_lecture/pybind11
@@ -0,0 +1 @@
+Subproject commit e2b884c33bcde70b2ea562ffa52dd7ebee276d50
diff --git a/lectures/pybind_lecture/test.py b/lectures/pybind_lecture/test.py
index 1ac87a2..91c7e15 100644
--- a/lectures/pybind_lecture/test.py
+++ b/lectures/pybind_lecture/test.py
@@ -1,121 +1,121 @@
 #!/bin/env python3
 
 import pyexample
-help(pyexample)
+# help(pyexample)
 
-help(pyexample.add)
+# help(pyexample.add)
 print(pyexample.add(1, 2))
 
 help(pyexample.add_withparams)
 print(pyexample.add_withparams(j=2, i=1))
 
 help(pyexample.add_withdefaults)
 print(pyexample.add_withdefaults())
 
 print(pyexample.yopla)
 print(pyexample.global_variable)
 pyexample.global_variable = 1
 print(pyexample.global_variable)
 pyexample.print_var()
 # conclusion ? not copied by reference
 print(pyexample.global_variable_rw)
 pyexample.global_variable_rw.fset(2)
 # this do not work !!!
 # pyexample.global_variable_rw = 2
 pyexample.print_var()
 
 help(pyexample.Animal)
 print(pyexample.Animal)
 a = pyexample.Animal('kitty')
 b = pyexample.Animal()
 a.python_extension()
 print(a)
 print(b)
 a.scream()
 print(a.name)
 a.name = "new"
 print(a.name)
 a.new_member = 'toto'
 print(dir(a))
 print(a.__dict__)
 
 dog = pyexample.Dog('poppy')
 dog.scream()
 help(dog.overloadedFoo)
 
 # automatic upcasting to true type (overhead)
 d = pyexample.makeAnimalUnique()
 print(type(d))
 d.scream()
 
 print(pyexample.ParticleType)
 print(pyexample.ParticleType.__members__)
 
 
 # Python and C++ use fundamentally different ways of managing the
 # memory and lifetime of objects managed by them. This can lead to
 # issues when creating bindings for functions that return a
 # non-trivial type. Just by looking at the type information, it is not
 # clear whether Python should take charge of the returned value and
 # eventually free its resources, or if this is handled on the C++
 # side. For this reason, pybind11 provides a several return value
 # policy annotations that can be passed to the module::def() and
 # class_::def() functions. The default policy is
 # return_value_policy::automatic.
 
 # segfault because del on python object calls delete on c++ pointer
 print("segfault")
 a = pyexample.get_animal()
 print(a)
 del a
 a = pyexample.get_animal()
 print(a)
 
 # solution: copy
 a = pyexample.get_animal_copy()
 print(a)
 a.name = "test"
 del a
 a = pyexample.get_animal_copy()
 print(a)
 
 # solution: by reference (no ownership)
 # memory management made by C++ only
 a = pyexample.get_animal_reference()
 print(a)
 a.name = "test"
 del a
 a = pyexample.get_animal_reference()
 print(a)
 
 # Code with invalid return value policies might access uninitialized
 # memory or free data structures multiple times, which can lead to
 # hard-to-debug non-determinism and segmentation faults, hence it is
 # worth spending the time to understand all the different options in
 # the table of return value policy.
 
 # As an alternative to elaborate call policies and lifetime management
 # logic, consider using smart pointers (see the section on Custom
 # smart pointers for details). Smart pointers can tell whether an
 # object is still referenced from C++ or Python, which generally
 # eliminates the kinds of inconsistencies that can lead to crashes or
 # undefined behavior. For functions returning smart pointers, it is
 # not necessary to specify a return value policy.
 
 a = pyexample.get_shared_ptr()
 a.name = "test"
 print(a)
 
 pyexample.register_animal_shared(a)
 
 pyexample.apply(lambda x: print(x*2))
 
 # https://pybind11.readthedocs.io/en/stable/advanced/exceptions.html
 # automatic management of standard exceptions
 
 try:
     pyexample._raise()
 except Exception as e:
     print(e)
 
 # https://pybind11.readthedocs.io/en/stable/advanced/cast/index.html#
diff --git a/work/week13/starting_point/CMakeLists.txt b/work/week13/starting_point/CMakeLists.txt
index 22e3387..f1027e0 100644
--- a/work/week13/starting_point/CMakeLists.txt
+++ b/work/week13/starting_point/CMakeLists.txt
@@ -1,85 +1,84 @@
 cmake_minimum_required (VERSION 3.1)
 project (Particles)
 cmake_policy(VERSION 3.3)
 
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake")
 
-include_directories(/usr/include/eigen3)
-# find_package(Eigen3 3.3 REQUIRED NO_MODULE)
+find_package(Eigen3 3.3 REQUIRED NO_MODULE)
 
 ################################################################
 # libpart
 ################################################################
 
 add_library(part
   compute_boundary.cc
   compute_verlet_integration.cc
   particle.cc
   planet.cc
   compute_gravity.cc
   csv_reader.cc
   particles_factory_interface.cc
   planets_factory.cc
   compute_contact.cc
   compute_kinetic_energy.cc
   csv_writer.cc
   system.cc
   compute_energy.cc
   compute_potential_energy.cc
   ping_pong_ball.cc
   material_point.cc
   system_evolution.cc
   ping_pong_balls_factory.cc
   compute_interaction.cc
   compute_temperature.cc
   material_points_factory.cc
   compute_temperature_finite_differences.cc
   )
 
-# if (EIGEN3_FOUND)
-#   include_directories(${EIGEN3_INCLUDE_DIR})
-# endif()
+if (EIGEN3_FOUND)
+  include_directories(${EIGEN3_INCLUDE_DIR})
+endif()
 include_directories(.)
 add_executable(particles main.cc)
 target_link_libraries(particles part)
 
 ################################################################
 # Google test
 ################################################################
 
 
 include(GoogleTest)
 enable_testing()
 find_package(GTest)
 if (GTEST_FOUND)
 include_directories(${GTEST_INCLUDE_DIRS})
 add_executable(test_kepler test_kepler.cc)
 add_executable(test_fft test_fft.cc)
 target_link_libraries(test_kepler part ${GTEST_BOTH_LIBRARIES} pthread)
 target_link_libraries(test_fft part ${GTEST_BOTH_LIBRARIES} ${FFTW_LIBRARIES} pthread)
 gtest_discover_tests(test_kepler)
 gtest_discover_tests(test_fft)
 endif()
 
 ################################################################
 # Doxygen
 ################################################################
 
 find_package(Doxygen)
 if (DOXYGEN_FOUND)
 
 # to set other options, read: https://cmake.org/cmake/help/v3.9/module/FindDoxygen.html
 doxygen_add_docs(
     doxygen
     ${PROJECT_SOURCE_DIR}
     COMMENT "Generate html pages"
 )
 add_custom_target(doc DEPENDS doxygen)
 endif(DOXYGEN_FOUND)
 
 option (USE_GSL "Use gsl library" OFF)
 
 if (USE_GSL)
   message("AAAAAA ")
 endif(USE_GSL)
diff --git a/work/week14/particle-pybind/starting_point/.gitignore b/work/week14/particle-pybind/starting_point/.gitignore
new file mode 100644
index 0000000..a3f9ed2
--- /dev/null
+++ b/work/week14/particle-pybind/starting_point/.gitignore
@@ -0,0 +1 @@
+pybind11/
diff --git a/work/week14/particle-pybind/starting_point/CMakeLists.txt b/work/week14/particle-pybind/starting_point/CMakeLists.txt
index d10c10d..ea25663 100644
--- a/work/week14/particle-pybind/starting_point/CMakeLists.txt
+++ b/work/week14/particle-pybind/starting_point/CMakeLists.txt
@@ -1,103 +1,98 @@
 cmake_minimum_required (VERSION 3.1)
 project (Particles)
 cmake_policy(VERSION 3.3)
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS 1)
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake")
 
-################################################################
-# pybind11
-################################################################
-
-find_package(pybind11 REQUIRED)  # or add_subdirectory(pybind11)
-
-add_library(pypart MODULE 
-	material_points_factory.cc
-	particles_factory_interface.cc
-	ping_pong_balls_factory.cc
-	planets_factory.cc
-	csv_writer.cc
-	compute_temperature.cc)
-target_link_libraries(pypart PRIVATE pybind11::module)
-set_target_properties(pypart PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}"
-                                           SUFFIX "${PYTHON_MODULE_EXTENSION}")
-
-file(
-  COPY ${CMAKE_CURRENT_SOURCE_DIR}/main.py
-  DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/
-  )
-
-
 ################################################################
 # FFTW
 ################################################################
 
 set(FFTW_LIBRARY_PATH CACHE PATH "library where to search FFTW")
 find_library (FFTW_LIBRARY fftw3 /usr/include/ ${FFTW_LIBRARY_PATH})
 
 ################################################################
 # libpart
 ################################################################
 
 add_library(part SHARED
   compute_boundary.cc
-  compute_verlet_integration.cc 
-  particle.cc 
+  compute_verlet_integration.cc
+  particle.cc
   planet.cc
-  compute_gravity.cc 
-  csv_reader.cc 
-  particles_factory_interface.cc 
-  planets_factory.cc 
-  compute_contact.cc 
-  compute_kinetic_energy.cc 
-  csv_writer.cc 
-  system.cc 
-  compute_energy.cc 
-  compute_potential_energy.cc 
+  compute_gravity.cc
+  csv_reader.cc
+  particles_factory_interface.cc
+  planets_factory.cc
+  compute_contact.cc
+  compute_kinetic_energy.cc
+  csv_writer.cc
+  system.cc
+  compute_energy.cc
+  compute_potential_energy.cc
   ping_pong_ball.cc
-  material_point.cc 
-  system_evolution.cc 
-  ping_pong_balls_factory.cc 
+  material_point.cc
+  system_evolution.cc
+  ping_pong_balls_factory.cc
   compute_interaction.cc
   compute_temperature.cc
-  material_points_factory.cc 
+  material_points_factory.cc
   )
 
 target_link_libraries(part ${FFTW_LIBRARIES})
 add_executable(particles main.cc)
 target_link_libraries(particles part)
 
+################################################################
+# pybind11
+################################################################
+
+# find_package(pybind11 REQUIRED)   #  ---> for Flavio & Tristan
+add_subdirectory(pybind11)    # ---> for Theo
+
+add_library(pypart MODULE pypart_pybind.cpp)
+target_link_libraries(pypart PRIVATE pybind11::module part)
+set_target_properties(pypart PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}"
+                                           SUFFIX "${PYTHON_MODULE_EXTENSION}")
+
+file(
+  COPY ${CMAKE_CURRENT_SOURCE_DIR}/main.py
+  DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/
+  )
+
+
 ################################################################
 # Google test
 ################################################################
 
 
 include(GoogleTest)
 enable_testing()
 find_package(GTest)
 if (GTEST_FOUND)
 include_directories(${GTEST_INCLUDE_DIRS})
 add_executable(test_kepler test_kepler.cc)
 add_executable(test_fft test_fft.cc)
 target_link_libraries(test_kepler part ${GTEST_BOTH_LIBRARIES} pthread)
 target_link_libraries(test_fft part ${GTEST_BOTH_LIBRARIES} ${FFTW_LIBRARIES} pthread)
 gtest_discover_tests(test_kepler)
 gtest_discover_tests(test_fft)
 endif()
 
 ################################################################
 # Doxygen
 ################################################################
 
 find_package(Doxygen)
 if (DOXYGEN_FOUND)
 
 # to set other options, read: https://cmake.org/cmake/help/v3.9/module/FindDoxygen.html
 doxygen_add_docs(
     doxygen
     ${PROJECT_SOURCE_DIR}
     COMMENT "Generate html pages"
 )
 add_custom_target(doc DEPENDS doxygen)
 endif(DOXYGEN_FOUND)
diff --git a/work/week14/particle-pybind/starting_point/compute_temperature.cc b/work/week14/particle-pybind/starting_point/compute_temperature.cc
index c40f81b..02ba3e1 100644
--- a/work/week14/particle-pybind/starting_point/compute_temperature.cc
+++ b/work/week14/particle-pybind/starting_point/compute_temperature.cc
@@ -1,13 +1,13 @@
 #include "compute_temperature.hh"
 #include "fft.hh"
 #include "material_point.hh"
 #include <cmath>
 
 /* -------------------------------------------------------------------------- */
 
 void ComputeTemperature::compute(System& system) {
 }
 
 /* -------------------------------------------------------------------------- */
 
-#include "pypart.cpp"
\ No newline at end of file
+// #include "pypart_pybind.cpp"
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/csv_writer.cc b/work/week14/particle-pybind/starting_point/csv_writer.cc
index 190fe08..c8ce2a6 100644
--- a/work/week14/particle-pybind/starting_point/csv_writer.cc
+++ b/work/week14/particle-pybind/starting_point/csv_writer.cc
@@ -1,32 +1,32 @@
 #include "csv_writer.hh"
 #include <cstdlib>
 #include <fstream>
 /* -------------------------------------------------------------------------- */
 CsvWriter::CsvWriter(const std::string& filename) : filename(filename) {}
 
 /* -------------------------------------------------------------------------- */
 void CsvWriter::write(System& system) { this->compute(system); }
 
 /* -------------------------------------------------------------------------- */
 
 void CsvWriter::compute(System& system) {
 
   std::ofstream os(filename.c_str());
 
   if (os.is_open() == false) {
     std::cerr << "cannot open file " << filename << std::endl
               << "check that the dumps folder exists" << std::endl;
     std::exit(1);
   }
 
   UInt nb_particles = system.getNbParticles();
 
   for (UInt p = 0; p < nb_particles; ++p) {
     os << system.getParticle(p) << std::endl;
   }
 
   os.close();
 
 }
 
-#include "pypart.cpp"
\ No newline at end of file
+// #include "pypart_pybind.cpp"
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/fft.hh b/work/week14/particle-pybind/starting_point/fft.hh
index 0319abe..fba8f3d 100644
--- a/work/week14/particle-pybind/starting_point/fft.hh
+++ b/work/week14/particle-pybind/starting_point/fft.hh
@@ -1,37 +1,124 @@
 #ifndef FFT_HH
 #define FFT_HH
 /* ------------------------------------------------------ */
 #include "matrix.hh"
 #include "my_types.hh"
 #include <fftw3.h>
 /* ------------------------------------------------------ */
 
+/**
+
+This file defines an FFT structure that is a wrapper around the FFTW library.
+In particular, it defines 3 functions to work with signals of complex numbers:
+- the transform() function to compute the forward FFT of a complex signal
+- the itransform() function to compute the inverse FFT of a complex signal
+- the computeFrequencies() function to compute the sample frequencies of the
+    forward FFT of a complex signal (i.e. the coordinates of the signal in the Fourier space).
+    The Laplacian of these frequencies is then used to sovle the heat equation in the Fourier space.
+
+*/
 struct FFT {
 
   static Matrix<complex> transform(Matrix<complex>& m);
   static Matrix<complex> itransform(Matrix<complex>& m);
 
-  static Matrix<std::complex<int>> computeFrequencies(int size);
+  static Matrix<complex> computeFrequencies(int size);
 
 };
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::transform(Matrix<complex>& m_in) {
+
+  // Check memory allocation of input matrix
+  if (m_in.data() == nullptr) {
+    throw std::runtime_error("no memory allocated for m_in");
+  }
+
+  // Get matrix size and initialize output matrix
+  int N = m_in.size();
+  Matrix<complex> m_out(N);
+
+  // Initialize input and output 2D signals
+  auto in = (fftw_complex*) m_in.data();
+  auto out = (fftw_complex*) m_out.data();
+
+  // Create, execute and destroy FFT plan
+  fftw_plan p = fftw_plan_dft_2d(N, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+  fftw_execute(p);
+
+  // Destroy plan
+  fftw_destroy_plan(p);
+  // fftw_free(in); fftw_free(out);
+
+  // Return FFT matrix
+  return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::itransform(Matrix<complex>& m_in) {
+
+  // Check memory allocation of input matrix
+  if (m_in.data() == nullptr) {
+    throw std::runtime_error("no memory allocated for m_in");
+  }
+
+  // Get matrix size and initialize output matrix
+  int N = m_in.size();
+  Matrix<complex> m_out(N);
+
+  // Initialize input and output 2D signals
+  auto in = (fftw_complex*) m_in.data();
+  auto out = (fftw_complex*) m_out.data();
+
+  // Create, execute and destroy inverse FFT plan
+  fftw_plan p = fftw_plan_dft_2d(N, N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
+  fftw_execute(p);
+
+  // Destroy plan
+  fftw_destroy_plan(p);
+  // fftw_free(in); fftw_free(out);
+
+  // Divide output matrix by factor N2 (needed for 2D arrays)
+  m_out /= (N * N);
+
+  // Return inverse FFT matrix
+  return m_out;
 }
 
 /* ------------------------------------------------------ */
 
+inline Matrix<complex> FFT::computeFrequencies(int size) {
 
-/* ------------------------------------------------------ */
+  // Generate 1D frequency vector
+  std::vector<float> vec(size);
+  if (size % 2 == 0) {  // if size is even
+    for (int i = 0; i < size / 2; ++i) {
+      vec[i] = (float)i;
+    }
+    for (int i = size / 2; i < size; ++i) {
+      vec[i] = -size / 2  + ((float)i - size / 2);
+    }
+  } else { // if size is odd
+    for (int i = 0; i <= (size - 1) / 2; ++i) {
+      vec[i] = (float)i;
+    }
+    for (int i = (size - 1) / 2 + 1; i < size; ++i) {
+      vec[i] = -((size - 1) / 2 + 1)  + ((float)i - (size - 1) / 2);
+    }
+  }
 
-inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {
+  // Populate 2D matrix with frequencies vector
+  Matrix<complex> m_out(size);
+  for (int i = 0; i < size; ++i) {
+    for (int j = 0; j < size; ++j) {
+      m_out(i, j) = complex(vec[j], 0);
+    }
+  }
 
+  return m_out;
 }
 
+
 #endif  // FFT_HH
diff --git a/work/week14/particle-pybind/starting_point/material_points_factory.cc b/work/week14/particle-pybind/starting_point/material_points_factory.cc
index 903b135..589e13f 100644
--- a/work/week14/particle-pybind/starting_point/material_points_factory.cc
+++ b/work/week14/particle-pybind/starting_point/material_points_factory.cc
@@ -1,61 +1,61 @@
 #include "material_points_factory.hh"
 #include "compute_temperature.hh"
 #include "csv_reader.hh"
 #include "csv_writer.hh"
 #include "material_point.hh"
 #include <cmath>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 std::unique_ptr<Particle> MaterialPointsFactory::createParticle() {
   return std::make_unique<MaterialPoint>();
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void MaterialPointsFactory::createDefaultComputes(Real timestep) {
 
   auto compute_temp = std::make_shared<ComputeTemperature>();
   compute_temp->getConductivity() = 1.;
   compute_temp->getL() = 2.;
   compute_temp->getCapacity() = 1.;
   compute_temp->getDensity() = 1.;
   compute_temp->getDeltat() = timestep;
   this->system_evolution->addCompute(compute_temp);
 }
 
 /* -------------------------------------------------------------------------- */
 
 SystemEvolution&
 MaterialPointsFactory::createSimulation(const std::string& fname,
                                         Real timestep) {
 
   this->system_evolution =
       std::make_unique<SystemEvolution>(std::make_unique<System>());
 
   CsvReader reader(fname);
   reader.read(this->system_evolution->getSystem());
 
   // check if it is a square number
   auto N = this->system_evolution->getSystem().getNbParticles();
   int side = std::sqrt(N);
   if (side * side != N)
     throw std::runtime_error("number of particles is not square");
 
   this->createComputes(timestep);
   return *system_evolution;
 }
 
 /* -------------------------------------------------------------------------- */
 
 ParticlesFactoryInterface& MaterialPointsFactory::getInstance() {
   if (not ParticlesFactoryInterface::factory)
     ParticlesFactoryInterface::factory = new MaterialPointsFactory;
 
   return *factory;
 }
 
 /* -------------------------------------------------------------------------- */
 
-#include "pypart.cpp"
\ No newline at end of file
+// #include "pypart_pybind.cpp"
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/matrix.hh b/work/week14/particle-pybind/starting_point/matrix.hh
index 81ca38b..8e8ad98 100644
--- a/work/week14/particle-pybind/starting_point/matrix.hh
+++ b/work/week14/particle-pybind/starting_point/matrix.hh
@@ -1,104 +1,106 @@
 #ifndef MATRIX_HH
 #define MATRIX_HH
 /* ------------------------------------------------------ */
 #include "my_types.hh"
 #include <algorithm>
 #include <array>
 #include <tuple>
 #include <vector>
 /* ------------------------------------------------------ */
 
 template <typename T>
 struct MatrixIterator {
   MatrixIterator(UInt index, UInt size, T* ptr)
       : index(index), size(size), ptr(ptr){};
 
   T& operator*() { return ptr[index]; };
 
   void operator++() { index++; };
 
   bool operator!=(MatrixIterator& other) { return this->index != other.index; };
 
   int index, size;
   T* ptr;
 };
 /* ------------------------------------------------------ */
 
 template <typename T>
 struct Matrix {
 
   Matrix(){};
   Matrix(UInt size) { resize(size); };
 
   UInt rows() { return this->size(); };
   UInt cols() { return this->size(); };
   UInt size() { return std::sqrt(storage.size()); };
 
   void resize(UInt size) { storage.resize(size * size); }
   T& operator()(UInt i, UInt j) { return storage[j + i * this->size()]; }
 
   Matrix& operator/=(const T& c) {
     std::for_each(storage.begin(), storage.end(), [&c](auto& v) { v /= c; });
     return *this;
   };
 
   T* data() { return &storage[0]; };
 
   std::vector<T> storage;
 
   MatrixIterator<T> begin() {
     return MatrixIterator<T>(0, this->size(), this->data());
   };
   MatrixIterator<T> end() {
     return MatrixIterator<T>(storage.size(), this->size(), this->data());
   };
 };
 
 /* ------------------------------------------------------ */
 
 template <typename T>
 struct MatrixIndexIterator : public MatrixIterator<T> {
   MatrixIndexIterator(UInt index, UInt size, T* ptr)
       : MatrixIterator<T>(index, size, ptr){};
 
   std::tuple<UInt, UInt, T&> operator*() {
 
     int i = this->index / this->size;
     int j = this->index % this->size;
     return std::tuple<UInt, UInt, T&>(i, j, this->ptr[this->index]);
   };
 };
 
 /* ------------------------------------------------------ */
 template <typename T>
 struct IndexedMatrix {
 
   IndexedMatrix(Matrix<T>& mat) : mat(mat){};
 
   MatrixIndexIterator<T> begin() {
     return MatrixIndexIterator<T>(0, mat.size(), mat.data());
   };
   MatrixIndexIterator<T> end() {
     return MatrixIndexIterator<T>(mat.storage.size(), mat.size(), mat.data());
   };
 
 private:
   Matrix<T>& mat;
 };
 /* ------------------------------------------------------ */
 
 template <typename T>
 IndexedMatrix<T> index(Matrix<T>& mat) {
   return IndexedMatrix<T>(mat);
 }
 /* ------------------------------------------------------ */
 
-template <typename T>
-struct std::iterator_traits<MatrixIterator<T>> {
-  using value_type = T;
-};
+namespace std {
+  template <typename T>
+  struct iterator_traits<MatrixIterator<T>> {
+    using value_type = T;
+  };
+}
 
 /* ------------------------------------------------------ */
 
 
 #endif //MATRIX
diff --git a/work/week14/particle-pybind/starting_point/particles_factory_interface.cc b/work/week14/particle-pybind/starting_point/particles_factory_interface.cc
index 1485b3d..a244c37 100644
--- a/work/week14/particle-pybind/starting_point/particles_factory_interface.cc
+++ b/work/week14/particle-pybind/starting_point/particles_factory_interface.cc
@@ -1,13 +1,13 @@
 #include "particles_factory_interface.hh"
 #include "planets_factory.hh"
 /* -------------------------------------------------------------------------- */
 ParticlesFactoryInterface& ParticlesFactoryInterface::getInstance() {
 
   return *factory;
 
 }
 
 /* -------------------------------------------------------------------------- */
 ParticlesFactoryInterface* ParticlesFactoryInterface::factory = nullptr;
 
-#include "pypart.cpp"
\ No newline at end of file
+// #include "pypart_pybind.cpp"
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/ping_pong_balls_factory.cc b/work/week14/particle-pybind/starting_point/ping_pong_balls_factory.cc
index 5261dfb..7e2f85a 100644
--- a/work/week14/particle-pybind/starting_point/ping_pong_balls_factory.cc
+++ b/work/week14/particle-pybind/starting_point/ping_pong_balls_factory.cc
@@ -1,47 +1,47 @@
 #include "ping_pong_balls_factory.hh"
 #include "compute_contact.hh"
 #include "compute_verlet_integration.hh"
 #include "csv_reader.hh"
 #include "csv_writer.hh"
 #include "ping_pong_ball.hh"
 #include <cmath>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 std::unique_ptr<Particle> PingPongBallsFactory::createParticle() {
   return std::make_unique<PingPongBall>();
 }
 
 /* -------------------------------------------------------------------------- */
 
 SystemEvolution&
 PingPongBallsFactory::createSimulation(const std::string& fname,
                                        Real timestep) {
   this->system_evolution =
       std::make_unique<SystemEvolution>(std::make_unique<System>());
 
   CsvReader reader(fname);
   reader.read(this->system_evolution->getSystem());
 
   auto contact = std::make_shared<ComputeContact>();
   auto verlet = std::make_shared<ComputeVerletIntegration>(timestep);
 
   contact->setPenalty(1.);
   verlet->addInteraction(contact);
   this->system_evolution->addCompute(verlet);
 
   return *system_evolution;
 }
 
 /* -------------------------------------------------------------------------- */
 
 ParticlesFactoryInterface& PingPongBallsFactory::getInstance() {
   if (not ParticlesFactoryInterface::factory)
     ParticlesFactoryInterface::factory = new PingPongBallsFactory;
 
   return *factory;
 }
 
 /* -------------------------------------------------------------------------- */
 
-#include "pypart.cpp"
\ No newline at end of file
+// #include "pypart_pybind.cpp"
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/planets_factory.cc b/work/week14/particle-pybind/starting_point/planets_factory.cc
index 617509f..c6c166e 100644
--- a/work/week14/particle-pybind/starting_point/planets_factory.cc
+++ b/work/week14/particle-pybind/starting_point/planets_factory.cc
@@ -1,49 +1,49 @@
 #include "planets_factory.hh"
 #include "compute_gravity.hh"
 #include "compute_verlet_integration.hh"
 #include "csv_reader.hh"
 #include "csv_writer.hh"
 #include "planet.hh"
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 
 std::unique_ptr<Particle> PlanetsFactory::createParticle() {
 
   return std::make_unique<Planet>();
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 SystemEvolution& PlanetsFactory::createSimulation(const std::string& fname,
                                                   Real timestep) {
 
   this->system_evolution =
       std::make_unique<SystemEvolution>(std::make_unique<System>());
 
   CsvReader reader(fname);
   reader.read(this->system_evolution->getSystem());
 
   auto gravity = std::make_shared<ComputeGravity>();
   auto verlet = std::make_shared<ComputeVerletIntegration>(timestep);
 
   verlet->addInteraction(gravity);
   this->system_evolution->addCompute(verlet);
 
   return *this->system_evolution;
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 ParticlesFactoryInterface& PlanetsFactory::getInstance() {
 
   if (not ParticlesFactoryInterface::factory)
     ParticlesFactoryInterface::factory = new PlanetsFactory;
 
   return *factory;
 
 }
 /* -------------------------------------------------------------------------- */
 
-#include "pypart.cpp"
\ No newline at end of file
+// #include "pypart_pybind.cpp"
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/pypart.cpp b/work/week14/particle-pybind/starting_point/pypart.cpp
deleted file mode 100644
index e456e3d..0000000
--- a/work/week14/particle-pybind/starting_point/pypart.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#include <functional>
-#include <iostream>
-#include <pybind11/pybind11.h>
-
-namespace py = pybind11;
-
-PYBIND11_MODULE(MaterialPointsFactory, m) {
-
-}
-
-PYBIND11_MODULE(ParticlesFactoryInterface, m) {
-
-}
-
-PYBIND11_MODULE(PingPongBallsFactory, m) {
-
-}
-
-PYBIND11_MODULE(PlanetsFactory, m) {
-
-}
-
-PYBIND11_MODULE(CsvWriter, m) {
-
-}
-
-PYBIND11_MODULE(ComputeTemperature, m) {
-
-}
\ No newline at end of file
diff --git a/work/week14/particle-pybind/starting_point/pypart_pybind.cpp b/work/week14/particle-pybind/starting_point/pypart_pybind.cpp
new file mode 100644
index 0000000..979e88c
--- /dev/null
+++ b/work/week14/particle-pybind/starting_point/pypart_pybind.cpp
@@ -0,0 +1,16 @@
+#include <functional>
+#include <iostream>
+#include <pybind11/pybind11.h>
+
+#include "csv_writer.cc"
+#include "particles_factory_interface.cc"
+#include "ping_pong_balls_factory.cc"
+#include "material_points_factory.cc"
+#include "planets_factory.cc"
+#include "compute_temperature.cc"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(pypart, m) {
+
+}
diff --git a/work/week14/particle-pybind/starting_point/test_fft.cc b/work/week14/particle-pybind/starting_point/test_fft.cc
index f50446c..efd7929 100644
--- a/work/week14/particle-pybind/starting_point/test_fft.cc
+++ b/work/week14/particle-pybind/starting_point/test_fft.cc
@@ -1,42 +1,107 @@
-#include "compute_temperature.hh"
-#include "csv_writer.hh"
-#include "fft.hh"
-#include "material_points_factory.hh"
-#include "my_types.hh"
-#include "vector.hh"
 #include <gtest/gtest.h>
+#include <cstdlib>
+#include <fstream>
+#include <sstream>
+#include "my_types.hh"
+#include "fft.hh"
+
+
 /*****************************************************************/
 TEST(FFT, transform) {
   UInt N = 512;
   Matrix<complex> m(N);
 
   Real k = 2 * M_PI / N;
   for (auto&& entry : index(m)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto& val = std::get<2>(entry);
     val = cos(k * i);
   }
 
   Matrix<complex> res = FFT::transform(m);
 
   for (auto&& entry : index(res)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto& val = std::get<2>(entry);
     if (std::abs(val) > 1e-10)
       std::cout << i << "," << j << " = " << val << std::endl;
 
     if (i == 1 && j == 0)
       ASSERT_NEAR(std::abs(val), N * N / 2, 1e-10);
     else if (i == N - 1 && j == 0)
       ASSERT_NEAR(std::abs(val), N * N / 2, 1e-10);
     else
       ASSERT_NEAR(std::abs(val), 0, 1e-10);
   }
 }
 /*****************************************************************/
 
 TEST(FFT, inverse_transform) {
+  UInt N = 512;
+  Matrix<complex> m(N);
+
+  Real k = 2 * M_PI / N;
+  for (auto&& entry : index(m)) {
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
+    auto& val = std::get<2>(entry);
+    if ((i == 1 && j == 0) || (i == N - 1 && j == 0)) {
+      val = N * N / 2;
+    } else {
+      val = 0;
+    }
+  }
+
+  Matrix<complex> res = FFT::itransform(m);
+
+  for (auto&& entry : index(res)) {
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
+    auto& val = std::get<2>(entry);
+    ASSERT_NEAR(std::abs(val), std::abs(cos(k * i)), 1e-10);
+  }
 }
-/*****************************************************************/
+
+
+void testComputeFrequencies(UInt N) {
+
+  std::string filename = "test_computeFrequencies_n_" + std::to_string(N) + ".txt";
+  std::cout << filename << std::endl;
+
+  // Populate reference frequency matrix from numpy generated file
+  Matrix<std::complex<int>> freqs_ref(N);
+  std::ifstream is(filename.c_str());
+
+  if (is.is_open() == false) {
+    std::cerr << "cannot open file" << std::endl;
+    throw;
+  }
+  std::string line;
+  double tmp;
+  for (int i=0; i < N; ++i) {
+    getline(is, line);
+    std::stringstream sstr(line);
+    for (int j = 0; j < N; ++j) {
+        sstr >> tmp;
+        freqs_ref(i, j) = std::complex<int>(tmp, 0);
+    }
+  }
+  is.close();
+
+  // Create frequency matrix with our own implementation
+  Matrix<complex> freqs = FFT::computeFrequencies(N);
+
+  // Compare the two matrices
+  for (auto&& entry : index(freqs)) {
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
+    auto& val = std::get<2>(entry);
+    auto& ref = freqs_ref(i, j);
+    ASSERT_NEAR(std::abs(val), std::abs(ref), 1e-10);
+  }
+}
+
+TEST(FFT, compute_frequencies_even) { testComputeFrequencies(8); }
+TEST(FFT, compute_frequencies_odd) { testComputeFrequencies(9); }
diff --git a/work/week14/particle-pybind/starting_point/vector.hh b/work/week14/particle-pybind/starting_point/vector.hh
index d12cd46..115ce62 100644
--- a/work/week14/particle-pybind/starting_point/vector.hh
+++ b/work/week14/particle-pybind/starting_point/vector.hh
@@ -1,75 +1,75 @@
 #ifndef __VECTOR__HH__
 #define __VECTOR__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "my_types.hh"
 /* -------------------------------------------------------------------------- */
 #include <array>
 
 /**
  * @brief 3D Vector class
- * 
+ *
  * Note that no constructor is needed for this class
  */
 class Vector {
 public:
   //! Dimension of vector
   static constexpr UInt dim = 3;
 
   // Methods
 public:
   //! access given component
   Real& operator[](UInt i);
   //! access given component (const)
   const Real& operator[](UInt i) const;
-  
+
   //! square of euclidian norm
   Real squaredNorm() const;
 
   // Operators that make sense for vectors
   Vector& operator+=(const Vector& vec);
   Vector& operator-=(const Vector& vec);
   Vector& operator*=(Real val);
   Vector& operator/=(Real val);
 
   //! Copy operator
   Vector& operator=(const Vector& vec);
   //! Assign a value
   Vector& operator=(Real val);
 
 
 public:
   //! Output to stream
   friend std::ostream& operator<<(std::ostream& stream, const Vector& _this);
   //! Initialize from stream
   friend std::istream& operator>>(std::istream& stream, Vector& _this);
 
 private:
   //! Vector values (initilized to zero by default)
-  std::array<Real, dim> values = {0};
+  std::array<Real, dim> values = {{0}};
 };
 
 /* -------------------------------------------------------------------------- */
 /* Separate function definitions */
 /* -------------------------------------------------------------------------- */
 
 /// standard output stream operator
 std::ostream& operator<<(std::ostream& stream, const Vector& _this);
 /// standard input stream operator
 std::istream& operator>>(std::istream& stream, Vector& _this);
 
 // We define the operators +, -, *, / with vectors and scalars
 Vector operator+(const Vector & a, const Vector& b);
 Vector operator-(const Vector & a, const Vector& b);
 
 // Symmetry of multiplication
 Vector operator*(const Vector & a, Real val);
 Vector operator*(Real val, const Vector & a);
 
 // For convenience
 Vector operator/(const Vector & a, Real val);
 
 // Include inline implementation
 #include "vector_inline.hh"
 
 #endif  //__VECTOR__HH__