diff --git a/work/week12/starting_point/fft.hh b/work/week12/starting_point/fft.hh
index 054c0c3..c51c8da 100644
--- a/work/week12/starting_point/fft.hh
+++ b/work/week12/starting_point/fft.hh
@@ -1,136 +1,147 @@
 #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);
 
 };
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::transform(Matrix<complex>& m_in) {
 
   // Get matrix size
   int N = m_in.size();
 
   // Declare input and output 2D signals
   fftw_complex *in = new fftw_complex[N*N];
   fftw_complex *out = new fftw_complex[N*N];
 
   // Populate input 2D signal
   for (int i = 0; i < N; ++i) {
     for (int j = 0; j < N; ++j) {
       in[i*N + j][0] = m_in(i, j).real();
       in[i*N + j][1] = m_in(i, j).imag();
     }
   }
 
   // Create, execute and destroy FFT plan
   fftw_plan p = fftw_plan_dft_2d(N, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
   fftw_execute(p);
 
   // Initialize and populate transformed matrix from output 2D signal
   Matrix<complex> m_out(N);
   for (int i = 0; i < N; ++i) {
     for (int j = 0; j < N; ++j) {
       m_out(i, j) = complex(out[i*N + j][0], out[i*N + j][1]);
     }
   }
 
   // Destroy plan and fftw signals
   fftw_destroy_plan(p);
   delete [] in;
   delete [] out;
 
   // Return FFT matrix
   return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::itransform(Matrix<complex>& m_in) {
 
   // Get matrix size
   int N = m_in.size();
 
   // Declare input and output 2D signals
   fftw_complex *in = new fftw_complex[N*N];
   fftw_complex *out = new fftw_complex[N*N];
 
   // Populate input 2D signal
   for (int i = 0; i < N; ++i) {
     for (int j = 0; j < N; ++j) {
       in[i*N + j][0] = m_in(i, j).real();
       in[i*N + j][1] = m_in(i, j).imag();
     }
   }
 
   // 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);
 
   // Initialize and populate transformed matrix from output 2D signal
   Matrix<complex> m_out(N);
   for (int i = 0; i < N; ++i) {
     for (int j = 0; j < N; ++j) {
       m_out(i, j) = complex(out[i*N + j][0], out[i*N + j][1]);
       m_out(i, j) /= (N * N);  // dividing by factor N2 needed for 2D arrays
     }
   }
 
   // Destroy plan and fftw signals
   fftw_destroy_plan(p);
   delete [] in;
   delete [] out;
 
   // Return FFT matrix
   return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 
 /* ------------------------------------------------------ */
 
 inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {
 
   // Generate 1D frequency vector
   std::vector<int> vec(size);
   if (size % 2 == 0) {  // if size is even
     for (int i = 0; i < size / 2; ++i) {
       vec[i] = i / size;
     }
     for (int i = size / 2; i < size; ++i) {
       vec[i] = (-size / 2  + (i - size / 2)) / size;
     }
   } else { // if size is odd
     for (int i = 0; i < (size - 1) / 2; ++i) {
       vec[i] = i / size;
     }
     for (int i = (size - 1) / 2; i < size; ++i) {
       vec[i] = (-(size - 1) / 2  + (i - (size - 1) / 2)) / size;
     }
   }
 
   // Populate 2D matrix with frequencies vector
   Matrix<std::complex<int>> m_out(size);
   for (int i = 0; i < size; ++i) {
     for (int j = 0; j < size; ++j) {
       m_out(i, j) = std::complex<int>(vec[i], vec[j]);
     }
   }
 
   return m_out;
 }
 
 
 #endif  // FFT_HH
diff --git a/work/week12/starting_point/material_point.hh b/work/week12/starting_point/material_point.hh
index e8d752e..c6b02a7 100644
--- a/work/week12/starting_point/material_point.hh
+++ b/work/week12/starting_point/material_point.hh
@@ -1,27 +1,36 @@
 #ifndef __MATERIAL_POINT__HH__
 #define __MATERIAL_POINT__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "particle.hh"
 
 //! Class for MaterialPoint
+/**
+This class is meant to be used as the fundamental entity 
+(i.e. a grid point) to compute the spatio-temporal 
+evolution of heat on a 2D grid. To this effect, 
+the class inherits from the Particle class and 
+defines two additional real number attributes:
+ temperature and heat rate, to be used to solve the heat equation.
+
+*/
 class MaterialPoint : public Particle {
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 
 public:
 
   void printself(std::ostream& stream) const override;
   void initself(std::istream& sstr) override;
 
   Real & getTemperature(){return temperature;};
   Real & getHeatRate(){return heat_rate;};
   
 private:
   Real temperature;
   Real heat_rate;
 };
 
 /* -------------------------------------------------------------------------- */
 #endif  //__MATERIAL_POINT__HH__
diff --git a/work/week12/starting_point/material_points_factory.hh b/work/week12/starting_point/material_points_factory.hh
index 85793cd..0db1319 100644
--- a/work/week12/starting_point/material_points_factory.hh
+++ b/work/week12/starting_point/material_points_factory.hh
@@ -1,30 +1,37 @@
 #ifndef __MATERIAL_POINTS_FACTORY__HH__
 #define __MATERIAL_POINTS_FACTORY__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "particles_factory_interface.hh"
 #include "material_point.hh"
 /* -------------------------------------------------------------------------- */
 
+
 //! Factory for material points
+/**
+This class is meant to be used as the interface through which the user can initialize a collection of M
+aterialPoint objects and then run a simulation to compute the spatio-temporal evolution of temperature across these points. 
+To this effect, it inherits from the generic ParticlesFactoryInterface and redefines the createParticle and createSimulation
+ methods in the specific context of material points (e.g. making sure that number of points is square at initialization).
+*/
 class MaterialPointsFactory : public ParticlesFactoryInterface {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 private:
   MaterialPointsFactory() = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   SystemEvolution& createSimulation(const std::string& fname,
 				    Real timestep) override;
 
   std::unique_ptr<Particle> createParticle() override;
 
   static ParticlesFactoryInterface& getInstance();
 };
 
 /* -------------------------------------------------------------------------- */
 #endif  //__MATERIAL_POINTS_FACTORY__HH__
diff --git a/work/week12/starting_point/matrix.hh b/work/week12/starting_point/matrix.hh
index 51776d9..e8a8b17 100644
--- a/work/week12/starting_point/matrix.hh
+++ b/work/week12/starting_point/matrix.hh
@@ -1,111 +1,157 @@
 #ifndef MATRIX_HH
 #define MATRIX_HH
 /* ------------------------------------------------------ */
 #include "my_types.hh"
 #include <algorithm>
 #include <array>
 #include <tuple>
 #include <vector>
 /* ------------------------------------------------------ */
-
+//The *matrix.hh* file defines several template classes to work with matrices of any data type.
+
+
+/**
+ * \Brief Allows iteration over a Matrix object
+ *
+ * This template is an abstract object allowing iteration over Matrix objects.
+ *It contains a pointer to a Matrix object as well as a size and an index field referring to the size of the matrix and the current index.
+ *It also contains several operators:
+ *- the ++ operator handles iteration by just incrementing the index field
+ *- the * operator (i.e. accessor) returns the content of the i-th element of the matrix (i.e. at the current index)
+ *- the != operator checks for inequality of index fields between the current operator and another operator
+ */
 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;
 };
 /* ------------------------------------------------------ */
 
+/**
+
+\Brief Defines a matrix object with its own methods
+
+This template is an abstract data type representing a 2D square (NxN) matrix.
+Its underlying data container is a 1D std::vector.
+It also contains several methods and operators:
+- the size(), rows() and cols() methods all return the size (N) of the matrix, which is the square root of the number of elements in the 1D vector (NxN)
+- the resize() method resizes the matrix by resizing the underlying data vector
+- the accessor operator returns, for a pair of indexes (i, j), the element of the 1D vector at the corresponding serialized index (j * N + i)
+- the /= operator handles element-wise division of the matrix by a scalar
+- the data() method returns the underlying data vector
+- the begin() and end() methods return iterator objects that can point to the first and last elements of the matrix, which can then be used in range for loop to iterate over the matrix elements.
+
+*/
 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 * this->size() + i]; }
 
   Matrix& operator/=(const T& c) {
     std::for_each(storage.begin(), storage.end(), [&c](auto& v) { v /= c; });
     return *this;
   };
 
   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());
   };
 };
 
 /* ------------------------------------------------------ */
 
+/**
+/Brief Iterates over the matrix indexes as tuples (i,j,x): (matrix indexes i, j and x the pointed element )
+
+This template defines a more complex type of abstract operator that inherits from MatrixIterator, 
+but with overwritten accessor element that now returns a 3 elements (i, j, x) tuple, where x 
+is the pointed element, i and j represent the indexes of the element along both matrix dimensions.
+*/
 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]);
   };
 };
 
 /* ------------------------------------------------------ */
+
+
+/**
+/Brief 
+
+This template defines a wrapper around the Matrix template in which the begin and end operators 
+are of type *MatrixIndexIterator*, and can then be used in different types of range for loops.
+The index() method can be used to convert a "simple" Matrix object into an IndexedMatrix object.
+
+Finally, the std::iterator_traits template is used to define the data type used by the MatrixIterator objects when they are then used in other parts of the code.
+
+*/
 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);
 }
 /* ------------------------------------------------------ */
 
 namespace std {
   template <typename T>
   struct iterator_traits<MatrixIterator<T>> {
     using value_type = T;
   };
 }
 
 /* ------------------------------------------------------ */
 
 
 #endif //MATRIX