diff --git a/Homework3/src/fft.hh b/Homework3/src/fft.hh
index 0dfb04e..6629904 100644
--- a/Homework3/src/fft.hh
+++ b/Homework3/src/fft.hh
@@ -1,129 +1,131 @@
 #ifndef FFT_HH
 #define FFT_HH
 /* ------------------------------------------------------ */
 #include "matrix.hh"
 #include "my_types.hh"
 #include "fftw3.h"
 /* ------------------------------------------------------ */
 
 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) {
-    fftw_complex *C_in;
-    fftw_complex *C_out;
-    fftw_plan forward_plan;
-    UInt length = m_in.size()*m_in.size();
-    UInt rows = m_in.rows();
-    UInt columns = m_in.cols();
-    Matrix<complex> m_out(rows);
+    fftw_complex *C_in;                     // Input physical domain as a rowmajor vector
+    fftw_complex *C_out;                    // Output spectral domain as a rowmajor vector
+    fftw_plan forward_plan;                 // FFT plan
+    UInt length = std::pow(m_in.size(),2);  // Length of vectors
+    Matrix<complex> m_out(m_in.size());     // The output matrix
     
+    // Forming rowmajor vectors
     C_in = new fftw_complex [length];
     C_out = new fftw_complex [length];
     
+    // Filling the input vector
     for (auto&& entry : index(m_in))
     {
         int i = std::get<0>(entry);
         int j = std::get<1>(entry);
         auto& val = std::get<2>(entry);
-        C_in[i*rows+j][0]=val.real();
-        C_in[i*rows+j][1]=0.0;
+        C_in[i*m_in.size()+j][0]=val.real();
+        C_in[i*m_in.size()+j][1]=0.0;
     }
 
-    forward_plan =  fftw_plan_dft_2d(rows, columns, C_in, C_out, FFTW_FORWARD, FFTW_ESTIMATE);
+    // Executing the plan
+    forward_plan =  fftw_plan_dft_2d(m_in.size(), m_in.size(), C_in, C_out, FFTW_FORWARD, FFTW_ESTIMATE);
     fftw_execute(forward_plan);
     
     for (auto&& entry : index(m_out))
     {
         int i = std::get<0>(entry);
         int j = std::get<1>(entry);
         auto& val = std::get<2>(entry);
-        val=complex(C_out[i*rows+j][0], C_out[i*rows+j][1])/double(rows*columns);
+        val=complex(C_out[i*m_in.size()+j][0], C_out[i*m_in.size()+j][1])/double(std::pow(m_in.size(),2));
     }
     
-    
+    // Clean up
     fftw_destroy_plan(forward_plan);
     fftw_cleanup();
     fftw_free(C_in);
     fftw_free(C_out);
     
     return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::itransform(Matrix<complex>& m_in) {
-    fftw_complex *C_in;
-    fftw_complex *C_out;
-    fftw_plan backward_plan;
-    UInt length = m_in.size()*m_in.size();
-    UInt rows = m_in.rows();
-    UInt columns = m_in.cols();
-    Matrix<complex> m_out(rows);
+    fftw_complex *C_in;                     // Input spectral domain as a rowmajor vector
+    fftw_complex *C_out;                    // Output physical domain as a rowmajor vector
+    fftw_plan backward_plan;                // FFT plan
+    UInt length = std::pow(m_in.size(),2);  // Length of vectors
+    Matrix<complex> m_out(m_in.size());     // The output matrix
 
+    // Forming rowmajor vectors
     C_in = new fftw_complex [length];
     C_out = new fftw_complex [length];
 
+    // Filling the input vector
     for (auto&& entry : index(m_in))
     {
         int i = std::get<0>(entry);
         int j = std::get<1>(entry);
         auto& val = std::get<2>(entry);
-        C_in[i*rows+j][0]=val.real();
-        C_in[i*rows+j][1]=val.imag();
+        C_in[i*m_in.size()+j][0]=val.real();
+        C_in[i*m_in.size()+j][1]=val.imag();
     }
 
-    backward_plan =  fftw_plan_dft_2d(rows, columns, C_in, C_out, FFTW_BACKWARD, FFTW_ESTIMATE);
+    // Executing the plan
+    backward_plan = fftw_plan_dft_2d(m_in.size(), m_in.size(), C_in, C_out, FFTW_BACKWARD, FFTW_ESTIMATE);
     fftw_execute(backward_plan);
 
     for (auto&& entry : index(m_out))
     {
         int i = std::get<0>(entry);
         int j = std::get<1>(entry);
         auto& val = std::get<2>(entry);
-        val=complex(C_out[i*rows+j][0], 0.0);
+        val=complex(C_out[i*m_in.size()+j][0], 0.0);    // Output is real valued function.
     }
 
-
+    // Clean up
     fftw_destroy_plan(backward_plan);
     fftw_cleanup();
     fftw_free(C_in);
     fftw_free(C_out);
 
     return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {
     // Note: Be careful that later in calculating wave vectors
     // you should consider the length of the phisical domain.
     // Here, the real part stores wave number for y direction (i)
     // and the imaginary part stores wave nomber for x direction (j)
     // Note: Number of grids in each direction should be a power of 2
     
     Matrix<std::complex<int>> Wave_number(size);
     
     for (auto&& entry : index(Wave_number))
     {
         int i = std::get<0>(entry);
         int j = std::get<1>(entry);
         auto& val = std::get<2>(entry);
         
         i<size/2 ? val.real(i) : val.real(i-size);
         j<size/2 ? val.imag(j) : val.imag(j-size);
     }
     
     return Wave_number;
 }
 /* ------------------------------------------------------ */
 
 #endif  // FFT_HH
diff --git a/Homework3/src/heat_boundary.cc b/Homework3/src/heat_boundary.cc
index edda708..f841139 100644
--- a/Homework3/src/heat_boundary.cc
+++ b/Homework3/src/heat_boundary.cc
@@ -1,39 +1,37 @@
 #include "heat_boundary.hh"
 
 /* -------------------------------------------------------------------------- */
 
 void Heat_Boundary::set_BC(std::string BC)
 {
     BC_Type = BC;
     
     if(BC == "Dirichlet")
         std::cout<<"Dirichlet boundary condition set for all boundaries."<<std::endl;
     else
         std::cout<<"Periodic boundary condition set in both x- and y- direction."<<std::endl;
 }
 
 void Heat_Boundary::compute(System& system) {
-    UInt nb_particles = system.getNbParticles();
-    UInt size = std::sqrt(nb_particles);
-    UInt id;
+    UInt nb_particles = system.getNbParticles();    // Number of particles (grids)
+    UInt size = std::sqrt(nb_particles);            // Discretization in each side
+    UInt id;                                        // Particle id to call
 
     if (BC_Type == "Dirichlet"){
-        
-        
         for(int i=0; i<size; i++){
             // id = i + j*size   &   j=0
             id = i;
             auto& prtcl = dynamic_cast<MaterialPoint&>(system.getParticle(id));
             prtcl.getTemperature() = 0.0;
         }
         
         for(int j=0; j<size; j++){
             // id = i + j*size   &   i=0
             id = j*size;
             auto& prtcl = dynamic_cast<MaterialPoint&>(system.getParticle(id));
             prtcl.getTemperature() = 0.0;
         }
     }
 }
 
 /* -------------------------------------------------------------------------- */
diff --git a/Homework3/src/heat_boundary.hh b/Homework3/src/heat_boundary.hh
index 73b62f5..e3ae761 100644
--- a/Homework3/src/heat_boundary.hh
+++ b/Homework3/src/heat_boundary.hh
@@ -1,20 +1,20 @@
 #ifndef __HEAT_BOUNDARY__HH__
 #define __HEAT_BOUNDARY__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "compute.hh"
 #include "material_point.hh"
 /* -------------------------------------------------------------------------- */
 
 //! Compute interaction with simulation box
 class Heat_Boundary : public Compute{
-public:
-  void set_BC(std::string BC);
-  void compute(System& system);
-    
-private:
-    std::string BC_Type;
+    public:
+      void set_BC(std::string BC);
+      void compute(System& system);
+        
+    private:
+        std::string BC_Type;
 };
 
 /* -------------------------------------------------------------------------- */
 #endif  //__HEAT_BOUNDARY__HH__