diff --git a/work/week12/starting_point/test.csv b/work/week12/starting_point/test.csv
new file mode 100644
index 0000000..0fd5177
--- /dev/null
+++ b/work/week12/starting_point/test.csv
@@ -0,0 +1,10 @@
+4.031346456686666357e+01 4.745202600996351805e+01 1.940308014928729818e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.422780306371343340e+00 5.000000000000000000e+00 1.211312698262228052e-02
+4.602478328749229064e+01 1.557006663802866164e+01 2.755035812493834513e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.063556187124171437e+00 5.000000000000000000e+00 2.184984024359264490e-02
+1.060416371625742471e+01 3.773544907146016669e+01 2.001843350167121116e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.594676087448904322e+00 5.000000000000000000e+00 9.668394989608240786e-02
+4.655706344563853349e+01 1.403158588962928910e+01 2.897127099727141442e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.580986344834677482e+00 5.000000000000000000e+00 5.567080689991558279e-02
+7.112404767636071767e+00 4.312096764558515094e-01 2.223146874223951741e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.770610212244721815e+00 5.000000000000000000e+00 1.458395531703818998e-02
+4.463170923972278814e+01 4.319149310941583764e+01 7.579463321327556180e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.010867180192299042e+00 5.000000000000000000e+00 1.888468035024619754e-02
+1.159862413521104330e+01 1.646744678655822725e+01 1.414772546736760539e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.070170682826680819e+00 5.000000000000000000e+00 9.420937623441533182e-02
+1.492623409185320682e+01 1.378683688264756846e+01 3.910164360203076939e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.209275801699564168e+00 5.000000000000000000e+00 3.573290894387221028e-02
+1.748986175971386103e+01 9.508225964676864095e+00 3.080568868856505915e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.917553565743276955e+00 5.000000000000000000e+00 7.664994391488445802e-02
+3.348103820003989206e+00 6.127887005190779668e+00 3.776180931581274436e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.367498102017441486e+00 5.000000000000000000e+00 8.377202814649979679e-02
diff --git a/work/week12/starting_point/test_fft.cc b/work/week12/starting_point/test_fft.cc
index 32a74eb..b8bff84 100644
--- a/work/week12/starting_point/test_fft.cc
+++ b/work/week12/starting_point/test_fft.cc
@@ -1,63 +1,126 @@
 #include "my_types.hh"
 #include "fft.hh"
+#include "compute_temperature.hh"
 #include <gtest/gtest.h>
-
+#include <cstdlib>
 /*****************************************************************/
 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);
   // Matrix<std::complex<int>> a = FFT::computeFrequencies(10);
 
   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);
   }
 }
 /*****************************************************************/
+
+/*****************************************************************/
+
+TEST(FFT, Homogeneous) {
+  UInt N = 512;
+  Matrix<complex> m(N);
+
+
+  //Create a homogeneous matrix of material points
+  for (auto&& entry : index(m)) {
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
+    auto& val = std::get<2>(entry);
+    val = 5;
+  }
+
+  // Apply FFT transformation
+  Matrix<complex> tempMtx_fft = FFT::transform(m);
+
+  // Get FFT coordinates
+  Matrix<std::complex<int>> fft_coord = FFT::computeFrequencies(N);
+
+  // Heat flux
+  Real hv = 0;
+
+  for (auto&& entry : index(tempMtx_fft)) {
+    
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
+    auto& val = std::get<2>(entry);
+    
+    // Generates random particle heat rate
+    Real hr = rand() % 1;
+
+    Real q_squared;
+    if (i == 0 && j == 0){
+      val = 0;
+      q_squared = 1;
+    }else{
+      q_squared = pow(fft_coord(i,j).real(),2) + pow(fft_coord(j,i).real(),2);
+    }
+    // std::cout << "Before val is: " << val << ", " << hr << ", " << q_squared;
+    val = hv - val * hr * q_squared;
+    // std::cout << " After val is: " << val << std::endl;
+  }
+
+    // Inverse FFT transformation
+  Matrix<complex> tempMtx_ifft = FFT::itransform(tempMtx_fft);
+  tempMtx_ifft *= 0.1;
+
+  for (auto&& entry2 : index(tempMtx_ifft)) {
+    int i = std::get<0>(entry2);
+    int j = std::get<1>(entry2);
+    auto& val = std::get<2>(entry2);
+    // std::cout << " IFFT val is: " << val << std::endl;
+    Real new_m = std::abs(m(i,j)) + val.real();
+/*    std::cout << m(i,j)<< std::endl;
+    Real new_m=5;*/
+    // std::cout << " New temp is: " << mp.getTemperature() << std::endl;
+    ASSERT_NEAR(new_m,std::abs(m(i,j)), 1e-10);
+  }
+}
\ No newline at end of file
diff --git a/work/week9/particles/solution/CMakeLists.txt b/work/week9/particles/solution/CMakeLists.txt
new file mode 100644
index 0000000..4da61b7
--- /dev/null
+++ b/work/week9/particles/solution/CMakeLists.txt
@@ -0,0 +1,57 @@
+
+cmake_minimum_required (VERSION 2.6)
+project (Particles)
+
+set(CMAKE_CXX_STANDARD 14)
+
+# Locate GTest
+find_package(GTest REQUIRED)
+include_directories(${GTEST_INCLUDE_DIRS})
+
+
+add_executable(particles
+  main.cc
+  vector.cc
+  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 
+  system_evolution.cc 
+  ping_pong_balls_factory.cc 
+  compute_interaction.cc)
+
+
+add_executable(runtest
+  test_csv.cc
+  vector.cc
+  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 
+  system_evolution.cc 
+  ping_pong_balls_factory.cc 
+  compute_interaction.cc)
+
+target_link_libraries(runtest ${GTEST_LIBRARIES} pthread)
\ No newline at end of file
diff --git a/work/week9/particles/solution/circle_orbit.csv b/work/week9/particles/solution/circle_orbit.csv
new file mode 100644
index 0000000..1d9dfd1
--- /dev/null
+++ b/work/week9/particles/solution/circle_orbit.csv
@@ -0,0 +1,2 @@
+0 0 0 0 0 0 1e-4 0 0 1 sun
+1 0 0 0 1 0 -1e-4 0 0 1e-4 earth
diff --git a/work/week9/particles/solution/compute.hh b/work/week9/particles/solution/compute.hh
new file mode 100644
index 0000000..7ab83db
--- /dev/null
+++ b/work/week9/particles/solution/compute.hh
@@ -0,0 +1,23 @@
+#ifndef __COMPUTE__HH__
+#define __COMPUTE__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "system.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Base class for all compute
+class Compute {
+public:
+  //! Virtual destructor needed because we have subclasses
+  virtual ~Compute() = default;
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! Compute is pure virtual
+  virtual void compute(System& system) = 0;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE__HH__
diff --git a/work/week9/particles/solution/compute_boundary.cc b/work/week9/particles/solution/compute_boundary.cc
new file mode 100644
index 0000000..19766fe
--- /dev/null
+++ b/work/week9/particles/solution/compute_boundary.cc
@@ -0,0 +1,12 @@
+#include "compute_boundary.hh"
+/* -------------------------------------------------------------------------- */
+ComputeBoundary::ComputeBoundary(const Vector& xmin, const Vector& xmax)
+    : xmin(xmin), xmax(xmax) {
+}
+/* -------------------------------------------------------------------------- */
+
+void ComputeBoundary::compute(System& system) {
+
+}
+
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/compute_boundary.hh b/work/week9/particles/solution/compute_boundary.hh
new file mode 100644
index 0000000..c262cbc
--- /dev/null
+++ b/work/week9/particles/solution/compute_boundary.hh
@@ -0,0 +1,24 @@
+#ifndef __COMPUTE_BOUNDARY__HH__
+#define __COMPUTE_BOUNDARY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Compute interaction with simulation box
+class ComputeBoundary : public Compute {
+  // Constructors/Destructors
+public:
+  ComputeBoundary(const Vector& xmin, const Vector& xmax);
+
+  // Methods
+public:
+  void compute(System& system) override;
+
+  // Members
+protected:
+  Vector xmin, xmax;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_BOUNDARY__HH__
diff --git a/work/week9/particles/solution/compute_contact.cc b/work/week9/particles/solution/compute_contact.cc
new file mode 100644
index 0000000..2d5820f
--- /dev/null
+++ b/work/week9/particles/solution/compute_contact.cc
@@ -0,0 +1,9 @@
+#include "compute_contact.hh"
+#include "ping_pong_ball.hh"
+#include <cmath>
+/* -------------------------------------------------------------------------- */
+void ComputeContact::setPenalty(Real penalty) {
+}
+/* -------------------------------------------------------------------------- */
+void ComputeContact::compute(System& system) {
+}
diff --git a/work/week9/particles/solution/compute_contact.hh b/work/week9/particles/solution/compute_contact.hh
new file mode 100644
index 0000000..03ab721
--- /dev/null
+++ b/work/week9/particles/solution/compute_contact.hh
@@ -0,0 +1,26 @@
+#ifndef __COMPUTE_CONTACT__HH__
+#define __COMPUTE_CONTACT__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute_interaction.hh"
+
+//! Compute contact interaction between ping-pong balls
+class ComputeContact : public ComputeInteraction {
+
+  // Virtual implementation
+public:
+  //! Penalty contact implementation
+  void compute(System& system) override;
+
+  // Accessors
+public:
+  //! Set penalty
+  void setPenalty(Real penalty);
+
+  // Members
+protected:
+  Real penalty;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_CONTACT__HH__
diff --git a/work/week9/particles/solution/compute_energy.cc b/work/week9/particles/solution/compute_energy.cc
new file mode 100644
index 0000000..65290c3
--- /dev/null
+++ b/work/week9/particles/solution/compute_energy.cc
@@ -0,0 +1,2 @@
+#include "compute_energy.hh"
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/compute_energy.hh b/work/week9/particles/solution/compute_energy.hh
new file mode 100644
index 0000000..eca44ea
--- /dev/null
+++ b/work/week9/particles/solution/compute_energy.hh
@@ -0,0 +1,19 @@
+#ifndef __COMPUTE_ENERGY__HH__
+#define __COMPUTE_ENERGY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+
+//! Base class for energy computation
+class ComputeEnergy : public Compute {
+
+  // Methods
+public:
+  Real getValue() { return value; }
+
+protected:
+  Real value;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_ENERGY__HH__
diff --git a/work/week9/particles/solution/compute_gravity.cc b/work/week9/particles/solution/compute_gravity.cc
new file mode 100644
index 0000000..647d5c6
--- /dev/null
+++ b/work/week9/particles/solution/compute_gravity.cc
@@ -0,0 +1,47 @@
+#include "compute_gravity.hh"
+#include <cmath>
+/* -------------------------------------------------------------------------- */
+void ComputeGravity::compute(System& system) {
+
+  UInt size = system.getNbParticles();
+  // first of all reset forces to zero
+  for (auto & particle : system)
+    particle.getForce() = 0;
+
+  // the force vector
+  Vector force;
+  // the distance vector
+  Vector v_r;
+  // the normalized vector
+  Vector v_nr;
+  // the square of the distance
+  Real r2 = 0.;
+  // the distance
+  Real r = 0.;
+
+  for (UInt p1 = 0; p1 < size; ++p1) {
+    Particle& par1 = system.getParticle(p1);
+    for (UInt p2 = p1 + 1; p2 < size; ++p2) {
+      Particle& par2 = system.getParticle(p2);
+
+      // compute the distance vector and the square of distance
+
+      v_r = par2.getPosition() - par1.getPosition();
+      r2 = v_r.squaredNorm();
+
+      if (r2 == 0.)
+        continue;
+      // compute the distance
+      r = sqrt(r2);
+      v_nr = v_r / r;
+
+      // compute the pair force
+      force = par1.getMass() * par2.getMass() * G / r2 * v_nr;
+
+      // add up the force for both concerned particles
+      par2.getForce() -= force;
+      par1.getForce() += force;
+    }
+  }
+
+}
diff --git a/work/week9/particles/solution/compute_gravity.hh b/work/week9/particles/solution/compute_gravity.hh
new file mode 100644
index 0000000..77071bd
--- /dev/null
+++ b/work/week9/particles/solution/compute_gravity.hh
@@ -0,0 +1,27 @@
+#ifndef __COMPUTE_GRAVITY__HH__
+#define __COMPUTE_GRAVITY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute_interaction.hh"
+
+//! Compute Newton gravity interaction
+class ComputeGravity : public ComputeInteraction {
+
+  // Virtual implementation
+public:
+  //! Newton gravity implementation
+  void compute(System& system) override;
+
+  // Accessors
+public:
+  //! set the gravitational constant
+  void setG(Real G);
+
+  // Members
+private:
+  //! newton constant
+  Real G = 1.;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_GRAVITY__HH__
diff --git a/work/week9/particles/solution/compute_interaction.cc b/work/week9/particles/solution/compute_interaction.cc
new file mode 100644
index 0000000..d85c12e
--- /dev/null
+++ b/work/week9/particles/solution/compute_interaction.cc
@@ -0,0 +1,3 @@
+#include "compute_interaction.hh"
+#include <cmath>
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/compute_interaction.hh b/work/week9/particles/solution/compute_interaction.hh
new file mode 100644
index 0000000..4374eb4
--- /dev/null
+++ b/work/week9/particles/solution/compute_interaction.hh
@@ -0,0 +1,12 @@
+#ifndef __COMPUTE_INTERACTION__HH__
+#define __COMPUTE_INTERACTION__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+
+//! Base class for interaction computation
+class ComputeInteraction : public Compute {
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_INTERACTION__HH__
diff --git a/work/week9/particles/solution/compute_kinetic_energy.cc b/work/week9/particles/solution/compute_kinetic_energy.cc
new file mode 100644
index 0000000..ce515f0
--- /dev/null
+++ b/work/week9/particles/solution/compute_kinetic_energy.cc
@@ -0,0 +1,6 @@
+#include "compute_kinetic_energy.hh"
+/* -------------------------------------------------------------------------- */
+
+void ComputeKineticEnergy::compute(System& system) {}
+
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/compute_kinetic_energy.hh b/work/week9/particles/solution/compute_kinetic_energy.hh
new file mode 100644
index 0000000..7fce63b
--- /dev/null
+++ b/work/week9/particles/solution/compute_kinetic_energy.hh
@@ -0,0 +1,14 @@
+#ifndef __COMPUTE_KINETIC_ENERGY__HH__
+#define __COMPUTE_KINETIC_ENERGY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute_energy.hh"
+
+//! Compute kinetic energy of system
+class ComputeKineticEnergy : public ComputeEnergy {
+public:
+  void compute(System& system) override;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_KINETIC_ENERGY__HH__
diff --git a/work/week9/particles/solution/compute_potential_energy.cc b/work/week9/particles/solution/compute_potential_energy.cc
new file mode 100644
index 0000000..1174018
--- /dev/null
+++ b/work/week9/particles/solution/compute_potential_energy.cc
@@ -0,0 +1,9 @@
+#include "compute_potential_energy.hh"
+/* -------------------------------------------------------------------------- */
+
+ComputePotentialEnergy::ComputePotentialEnergy(ComputeInteraction& cForces)
+    : cForces(cForces) {}
+
+/* -------------------------------------------------------------------------- */
+
+void ComputePotentialEnergy::compute(System& system) {}
diff --git a/work/week9/particles/solution/compute_potential_energy.hh b/work/week9/particles/solution/compute_potential_energy.hh
new file mode 100644
index 0000000..aa197c4
--- /dev/null
+++ b/work/week9/particles/solution/compute_potential_energy.hh
@@ -0,0 +1,35 @@
+#ifndef __COMPUTE_POTENTIAL_ENERGY__HH__
+#define __COMPUTE_POTENTIAL_ENERGY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute_energy.hh"
+#include "compute_interaction.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Compute potential energy of system
+class ComputePotentialEnergy : public ComputeEnergy {
+
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+
+public:
+  ComputePotentialEnergy(ComputeInteraction& cForces);
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+
+public:
+  void compute(System& system) override;
+
+  /* ------------------------------------------------------------------------ */
+  /* Members                                                                  */
+  /* ------------------------------------------------------------------------ */
+
+protected:
+  ComputeInteraction& cForces;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_POTENTIAL_ENERGY__HH__
diff --git a/work/week9/particles/solution/compute_verlet_integration.cc b/work/week9/particles/solution/compute_verlet_integration.cc
new file mode 100644
index 0000000..e8ee826
--- /dev/null
+++ b/work/week9/particles/solution/compute_verlet_integration.cc
@@ -0,0 +1,43 @@
+#include "compute_verlet_integration.hh"
+
+ComputeVerletIntegration::ComputeVerletIntegration(Real dt) : dt(dt) {}
+
+/* -------------------------------------------------------------------------- */
+
+void ComputeVerletIntegration::setDeltaT(Real dt) {
+
+  this->dt = dt;
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+void ComputeVerletIntegration::compute(System& system) {
+ 
+  UInt size = system.getNbParticles();
+
+  for (auto& par : system) {
+    par.getVelocity() += .5 * dt * par.getForce() / par.getMass();
+    par.getPosition() += dt * par.getVelocity();
+  }
+
+  auto& sun = system.getParticle(0);
+  sun.getPosition() = 0;
+
+  for (auto& interaction : interactions)
+    interaction->compute(system);
+
+  for (auto& par : system) {
+    par.getVelocity() += .5 * dt * par.getForce() / par.getMass();
+  }
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+void ComputeVerletIntegration::addInteraction(
+    std::shared_ptr<ComputeInteraction> interaction) {
+
+  interactions.push_back(interaction);
+
+}
diff --git a/work/week9/particles/solution/compute_verlet_integration.hh b/work/week9/particles/solution/compute_verlet_integration.hh
new file mode 100644
index 0000000..e64876e
--- /dev/null
+++ b/work/week9/particles/solution/compute_verlet_integration.hh
@@ -0,0 +1,33 @@
+#ifndef __COMPUTE_VERLET_INTEGRATION__HH__
+#define __COMPUTE_VERLET_INTEGRATION__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+#include "compute_interaction.hh"
+/* -------------------------------------------------------------------------- */
+#include <cmath>
+
+//! Integrate equation of motion
+class ComputeVerletIntegration : public Compute {
+  using InteractionList = std::vector<std::shared_ptr<ComputeInteraction>>;
+
+  // Constructors/Destructors
+public:
+  ComputeVerletIntegration(Real timestep);
+
+  // Methods
+public:
+  //! Set time step
+  void setDeltaT(Real dt);
+  //! Evolve positions and velocities
+  void compute(System& system) override;
+  //! Add an interaction to the computation of forces
+  void addInteraction(std::shared_ptr<ComputeInteraction> interaction);
+
+private:
+  Real dt;
+  InteractionList interactions;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__COMPUTE_VERLET_INTEGRATION__HH__
diff --git a/work/week9/particles/solution/csv_reader.cc b/work/week9/particles/solution/csv_reader.cc
new file mode 100644
index 0000000..2569f36
--- /dev/null
+++ b/work/week9/particles/solution/csv_reader.cc
@@ -0,0 +1,41 @@
+#include "csv_reader.hh"
+#include "particles_factory_interface.hh"
+#include <fstream>
+#include <sstream>
+/* -------------------------------------------------------------------------- */
+CsvReader::CsvReader(const std::string& filename) : filename(filename) {}
+/* -------------------------------------------------------------------------- */
+void CsvReader::read(System& system) { this->compute(system); }
+/* -------------------------------------------------------------------------- */
+
+void CsvReader::compute(System& system) {
+
+  std::ifstream is(filename.c_str());
+  std::string line;
+
+  if (is.is_open() == false) {
+    std::cerr << "cannot open file " << filename << std::endl;
+    throw;
+  }
+
+  while (is.good()) {
+    getline(is, line);
+
+    if (line[0] == '#' || line.size() == 0)
+      continue;
+
+    auto p = ParticlesFactoryInterface::getInstance().createParticle();
+    std::stringstream sstr(line);
+    sstr >> *p;
+    std::cout << "sstr given to *p" << std::endl;
+    std::cout << p->getPosition() << std::endl;
+    std::cout << p->getVelocity() << std::endl;
+    system.addParticle(std::move(p));
+   
+  }
+
+  is.close();
+
+}
+
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/csv_reader.hh b/work/week9/particles/solution/csv_reader.hh
new file mode 100644
index 0000000..369c6bf
--- /dev/null
+++ b/work/week9/particles/solution/csv_reader.hh
@@ -0,0 +1,32 @@
+#ifndef __CSV_READER__HH__
+#define __CSV_READER__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+#include <string>
+
+//! Read system from csv input file
+class CsvReader : public Compute {
+
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! Construct from filename
+  CsvReader(const std::string& filename);
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! Left here for convenience (calls compute)
+  void read(System& system);
+  //! Read input file into system
+  void compute(System& system);
+
+protected:
+  std::string filename;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__CSV_READER__HH__
diff --git a/work/week9/particles/solution/csv_writer.cc b/work/week9/particles/solution/csv_writer.cc
new file mode 100644
index 0000000..6bbfd6d
--- /dev/null
+++ b/work/week9/particles/solution/csv_writer.cc
@@ -0,0 +1,30 @@
+#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();
+
+}
diff --git a/work/week9/particles/solution/csv_writer.hh b/work/week9/particles/solution/csv_writer.hh
new file mode 100644
index 0000000..0fabca4
--- /dev/null
+++ b/work/week9/particles/solution/csv_writer.hh
@@ -0,0 +1,27 @@
+#ifndef __CSV_WRITER__HH__
+#define __CSV_WRITER__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+#include <string>
+
+//! Write system state to csv file
+class CsvWriter : public Compute {
+  // Constructors/Destructors
+public:
+  //! Construct from filename
+  CsvWriter(const std::string& filename);
+
+  // Methods
+public:
+  //! Write file (calls compute)
+  void write(System& system);
+  //! Write file
+  void compute(System& system);
+
+protected:
+  std::string filename;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__CSV_WRITER__HH__
diff --git a/work/week9/particles/solution/main.cc b/work/week9/particles/solution/main.cc
new file mode 100644
index 0000000..4d94b79
--- /dev/null
+++ b/work/week9/particles/solution/main.cc
@@ -0,0 +1,66 @@
+#include "compute_gravity.hh"
+#include "compute_verlet_integration.hh"
+#include "csv_reader.hh"
+#include "csv_writer.hh"
+#include "my_types.hh"
+#include "ping_pong_balls_factory.hh"
+#include "planets_factory.hh"
+#include "system.hh"
+/* -------------------------------------------------------------------------- */
+#include <cstdlib>
+#include <iostream>
+#include <sstream>
+/* -------------------------------------------------------------------------- */
+
+int main(int argc, char** argv) {
+  if (argc != 6) {
+    std::cout << "Usage: " << argv[0]
+              << " nsteps dump_freq input.csv particle_type timestep yolo"
+              << std::endl;
+    std::cout << "\tparticle type can be: planet, ping_pong and you are in the right folder" << std::endl;
+    std::exit(EXIT_FAILURE);
+  }
+
+
+  // the number of steps to perform
+  Real nsteps;
+  std::stringstream(argv[1]) >> nsteps;
+  // freq to dump
+  int freq;
+  std::stringstream(argv[2]) >> freq;
+  // init file
+  std::string filename = argv[3];
+  // particle type
+  std::string type = argv[4];
+  // timestep
+  Real timestep;
+  std::stringstream(argv[5]) >> timestep;
+
+  if (type == "planet")
+    PlanetsFactory::getInstance();
+  else if (type == "ping_pong")
+    PingPongBallsFactory::getInstance();
+  else {
+    std::cout << "Unknown particle type: " << type << std::endl;
+    std::exit(EXIT_FAILURE);
+  }
+
+  ParticlesFactoryInterface& factory = ParticlesFactoryInterface::getInstance();
+  std::cout << " Instance obtained"<< std::endl;
+
+  SystemEvolution& evol = factory.createSimulation(filename, timestep);
+
+  std::cout << evol.getSystem().getParticle(0).getVelocity()<< std::endl;
+
+  std::cout << " Simulation created"<< std::endl;
+
+
+
+  evol.setNSteps(nsteps);
+  evol.setDumpFreq(freq);
+
+  evol.evolve();
+
+
+  return EXIT_SUCCESS;
+}
diff --git a/work/week9/particles/solution/my_types.hh b/work/week9/particles/solution/my_types.hh
new file mode 100644
index 0000000..8f7ecf8
--- /dev/null
+++ b/work/week9/particles/solution/my_types.hh
@@ -0,0 +1,27 @@
+#ifndef __MY_TYPES_HH__
+#define __MY_TYPES_HH__
+
+/* -------------------------------------------------------------------------- */
+#include <iostream>
+
+typedef unsigned int UInt;
+typedef double Real;
+
+/* -------------------------------------------------------------------------- */
+#define TO_IMPLEMENT                                                           \
+  {                                                                            \
+    std::cerr << std::endl                                                     \
+              << std::endl                                                     \
+              << "*********************************************************\n" \
+              << "To be filled @ " << __FILE__ << ":" << __LINE__ << std::endl \
+              << "*********************************************************\n" \
+              << std::endl                                                     \
+              << std::endl;                                                    \
+    throw;                                                                     \
+  }
+/* -------------------------------------------------------------------------- */
+
+#define EXERCISE_BEGIN_CORRECTION
+#define EXERCISE_END_CORRECTION
+
+#endif /* __MY_TYPES_HH__ */
diff --git a/work/week9/particles/solution/particle.cc b/work/week9/particles/solution/particle.cc
new file mode 100644
index 0000000..a9095bc
--- /dev/null
+++ b/work/week9/particles/solution/particle.cc
@@ -0,0 +1,21 @@
+#include "particle.hh"
+
+void Particle::printself(std::ostream& stream) const {
+
+  stream << " " << position;
+  stream << " " << velocity;
+  stream << " " << force;
+  stream << " " << mass;
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+void Particle::initself(std::istream& sstr) {
+
+  sstr >> position;
+  sstr >> velocity;
+  sstr >> force;
+  sstr >> mass;
+
+}
diff --git a/work/week9/particles/solution/particle.hh b/work/week9/particles/solution/particle.hh
new file mode 100644
index 0000000..db0ed0c
--- /dev/null
+++ b/work/week9/particles/solution/particle.hh
@@ -0,0 +1,61 @@
+#ifndef __PARTICLE__HH__
+#define __PARTICLE__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "vector.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Particle base class
+class Particle {
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! Default destructor (virtual because of virtual functions in this class)
+  virtual ~Particle() = default;
+
+  // Accessors
+public:
+  //! Get mass
+  Real& getMass() { return mass; }
+  //! Get position
+  Vector& getPosition() { return position; }
+  //! Get force
+  Vector& getForce() { return force; }
+  //! Get velocity
+  Vector& getVelocity() { return velocity; }
+
+  // I/O functions
+public:
+  /// function to print the content of the class
+  virtual void printself(std::ostream& stream) const;
+
+  /// function to initilize values from input files
+  virtual void initself(std::istream& sstr);
+
+  // Members
+protected:
+  Real mass;
+  Vector position, force, velocity;
+};
+
+
+/* -------------------------------------------------------------------------- */
+/* Inline functions */
+/* -------------------------------------------------------------------------- */
+
+inline std::istream& operator>>(std::istream& sstr, Particle& _this) {
+  _this.initself(sstr);
+  return sstr;
+}
+
+/* -------------------------------------------------------------------------- */
+inline std::ostream& operator<<(std::ostream& sstr, Particle& _this) {
+  _this.printself(sstr);
+  return sstr;
+}
+
+
+
+/* -------------------------------------------------------------------------- */
+#endif  //__PARTICLE__HH__
diff --git a/work/week9/particles/solution/particles_factory_interface.cc b/work/week9/particles/solution/particles_factory_interface.cc
new file mode 100644
index 0000000..088deeb
--- /dev/null
+++ b/work/week9/particles/solution/particles_factory_interface.cc
@@ -0,0 +1,11 @@
+#include "particles_factory_interface.hh"
+#include "planets_factory.hh"
+/* -------------------------------------------------------------------------- */
+ParticlesFactoryInterface& ParticlesFactoryInterface::getInstance() {
+
+  return *factory;
+
+}
+
+/* -------------------------------------------------------------------------- */
+ParticlesFactoryInterface* ParticlesFactoryInterface::factory = nullptr;
diff --git a/work/week9/particles/solution/particles_factory_interface.hh b/work/week9/particles/solution/particles_factory_interface.hh
new file mode 100644
index 0000000..b04c38b
--- /dev/null
+++ b/work/week9/particles/solution/particles_factory_interface.hh
@@ -0,0 +1,39 @@
+#ifndef __PARTICLES_FACTORY_INTERFACE__HH__
+#define __PARTICLES_FACTORY_INTERFACE__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "system_evolution.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Abstract factory defining interface
+class ParticlesFactoryInterface {
+  // Constructors/Destructors
+protected:
+  //! Instance constructor (protected)
+  ParticlesFactoryInterface() = default;
+
+public:
+  virtual ~ParticlesFactoryInterface() = default;
+
+  // Methods
+public:
+  //! Create a whole simulation from file
+  virtual SystemEvolution& createSimulation(const std::string& fname,
+                                            Real timestep) = 0;
+  //! Create a new particle
+  virtual std::unique_ptr<Particle> createParticle() = 0;
+
+  //! Get singleton instance
+  static ParticlesFactoryInterface& getInstance();
+
+  // Members
+protected:
+  std::vector<Particle*> list_particles;
+  std::unique_ptr<SystemEvolution> system_evolution = nullptr;
+
+  // Standard pointer because constructor is protected (cannot use make_unique)
+  static ParticlesFactoryInterface* factory;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__PARTICLES_FACTORY_INTERFACE__HH__
diff --git a/work/week9/particles/solution/ping_pong_ball.cc b/work/week9/particles/solution/ping_pong_ball.cc
new file mode 100644
index 0000000..60926ea
--- /dev/null
+++ b/work/week9/particles/solution/ping_pong_ball.cc
@@ -0,0 +1,10 @@
+#include "ping_pong_ball.hh"
+
+/* -------------------------------------------------------------------------- */
+void PingPongBall::printself(std::ostream& stream) const {
+}
+
+/* -------------------------------------------------------------------------- */
+
+void PingPongBall::initself(std::istream& sstr) {
+}
diff --git a/work/week9/particles/solution/ping_pong_ball.hh b/work/week9/particles/solution/ping_pong_ball.hh
new file mode 100644
index 0000000..6c5169f
--- /dev/null
+++ b/work/week9/particles/solution/ping_pong_ball.hh
@@ -0,0 +1,27 @@
+#ifndef __PING_PONG_BALL__HH__
+#define __PING_PONG_BALL__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "particle.hh"
+
+//! Class for ping-pong ball
+class PingPongBall : public Particle {
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+
+public:
+  //! Get contact dissipation
+  Real& getContactDissipation() { return contact_dissipation; }
+  //! Get ball radius
+  Real& getRadius() { return radius; }
+
+  void printself(std::ostream& stream) const override;
+  void initself(std::istream& sstr) override;
+
+private:
+  Real radius, contact_dissipation;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__PING_PONG_BALL__HH__
diff --git a/work/week9/particles/solution/ping_pong_balls_factory.cc b/work/week9/particles/solution/ping_pong_balls_factory.cc
new file mode 100644
index 0000000..9a43b06
--- /dev/null
+++ b/work/week9/particles/solution/ping_pong_balls_factory.cc
@@ -0,0 +1,26 @@
+#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() {
+}
+
+/* -------------------------------------------------------------------------- */
+
+SystemEvolution&
+PingPongBallsFactory::createSimulation(const std::string& fname,
+                                       Real timestep) {
+}
+
+/* -------------------------------------------------------------------------- */
+
+ParticlesFactoryInterface& PingPongBallsFactory::getInstance() {
+}
+
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/ping_pong_balls_factory.hh b/work/week9/particles/solution/ping_pong_balls_factory.hh
new file mode 100644
index 0000000..5a75091
--- /dev/null
+++ b/work/week9/particles/solution/ping_pong_balls_factory.hh
@@ -0,0 +1,30 @@
+#ifndef __PING_PING_BALLS_FACTORY__HH__
+#define __PING_PING_BALLS_FACTORY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "particles_factory_interface.hh"
+#include "ping_pong_ball.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Factory for ping-pong balls
+class PingPongBallsFactory : public ParticlesFactoryInterface {
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+private:
+  PingPongBallsFactory() = default;
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+public:
+  SystemEvolution& createSimulation(const std::string& fname,
+				    Real timestep) override;
+
+  std::unique_ptr<Particle> createParticle() override;
+
+  static ParticlesFactoryInterface& getInstance();
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__PING_PING_BALLS_FACTORY__HH__
diff --git a/work/week9/particles/solution/planet.cc b/work/week9/particles/solution/planet.cc
new file mode 100644
index 0000000..e841202
--- /dev/null
+++ b/work/week9/particles/solution/planet.cc
@@ -0,0 +1,9 @@
+#include "planet.hh"
+
+void Planet::initself(std::istream &stream) {
+}
+
+/* -------------------------------------------------------------------------- */
+
+void Planet::printself(std::ostream &stream) const {
+}
diff --git a/work/week9/particles/solution/planet.hh b/work/week9/particles/solution/planet.hh
new file mode 100644
index 0000000..0ebf4c1
--- /dev/null
+++ b/work/week9/particles/solution/planet.hh
@@ -0,0 +1,27 @@
+#ifndef __PLANET__HH__
+#define __PLANET__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "particle.hh"
+/* -------------------------------------------------------------------------- */
+#include <iostream>
+/* -------------------------------------------------------------------------- */
+
+//! Class for planet
+class Planet : public Particle {
+  // Methods
+public:
+  //! Get name
+  std::string& getName() { return name; }
+
+  // Inherited methods
+public:
+  void initself(std::istream& stream) override;
+  void printself(std::ostream& stream) const override;
+
+private:
+  std::string name;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__PLANET__HH__
diff --git a/work/week9/particles/solution/planets_factory.cc b/work/week9/particles/solution/planets_factory.cc
new file mode 100644
index 0000000..78abecb
--- /dev/null
+++ b/work/week9/particles/solution/planets_factory.cc
@@ -0,0 +1,47 @@
+#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;
+
+}
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/planets_factory.hh b/work/week9/particles/solution/planets_factory.hh
new file mode 100644
index 0000000..861ef60
--- /dev/null
+++ b/work/week9/particles/solution/planets_factory.hh
@@ -0,0 +1,31 @@
+#ifndef __PLANETS_FACTORY__HH__
+#define __PLANETS_FACTORY__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "particles_factory_interface.hh"
+#include "planet.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Factory for planet simulations
+class PlanetsFactory : public ParticlesFactoryInterface {
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+private:
+  PlanetsFactory() = default;
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+
+public:
+  SystemEvolution& createSimulation(const std::string& fname,
+                                    Real timestep) override;
+
+  std::unique_ptr<Particle> createParticle() override;
+
+  static ParticlesFactoryInterface& getInstance();
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__PLANETS_FACTORY__HH__
diff --git a/work/week9/particles/solution/system.cc b/work/week9/particles/solution/system.cc
new file mode 100644
index 0000000..4c07c24
--- /dev/null
+++ b/work/week9/particles/solution/system.cc
@@ -0,0 +1,23 @@
+#include "system.hh"
+
+Particle& System::getParticle(UInt i) {
+
+  return *list_particles[i];
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+void System::addParticle(const std::shared_ptr<Particle>& new_particle) {
+
+  list_particles.push_back(new_particle);
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+UInt System::getNbParticles() {
+
+  return list_particles.size();
+
+}
diff --git a/work/week9/particles/solution/system.hh b/work/week9/particles/solution/system.hh
new file mode 100644
index 0000000..967079d
--- /dev/null
+++ b/work/week9/particles/solution/system.hh
@@ -0,0 +1,54 @@
+#ifndef __SYSTEM__HH__
+#define __SYSTEM__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "my_types.hh"
+#include "particle.hh"
+#include <memory>
+#include <vector>
+/* -------------------------------------------------------------------------- */
+
+//! Container for particles
+class System {
+  /*
+   * No need for constructor/destructor with std::unique_ptr
+   */
+
+  // List of particle pointers
+  using ParticleList = std::vector<std::shared_ptr<Particle>>;
+
+  // Methods
+public:
+  //! Remove particle from vector
+  void removeParticle(UInt particle);
+  //! Get particle for specific id
+  Particle& getParticle(UInt i);
+  //! Add a particle to the system
+  void addParticle(const std::shared_ptr<Particle>& new_particle);
+  //! Get number of particles
+  UInt getNbParticles();
+
+
+  //! Iterator class to erase the unique pointer on Particle
+  struct iterator : ParticleList::iterator {
+    iterator(const ParticleList::iterator& it) : ParticleList::iterator(it) {}
+
+    //! Access the underlying particle
+    Particle& operator*() {
+      return *ParticleList::iterator::operator*();
+    }
+  };
+
+  // Iterators
+public:
+  auto begin() { return iterator(list_particles.begin()); }
+  auto end() { return iterator(list_particles.end()); }
+
+
+
+private:
+  ParticleList list_particles;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__SYSTEM__HH__
diff --git a/work/week9/particles/solution/system_evolution.cc b/work/week9/particles/solution/system_evolution.cc
new file mode 100644
index 0000000..f4d3184
--- /dev/null
+++ b/work/week9/particles/solution/system_evolution.cc
@@ -0,0 +1,42 @@
+#include "system_evolution.hh"
+#include "csv_writer.hh"
+/* -------------------------------------------------------------------------- */
+#include <iomanip>
+#include <sstream>
+/* -------------------------------------------------------------------------- */
+
+SystemEvolution::SystemEvolution(std::unique_ptr<System> system)
+    : system(std::move(system)) {}
+
+/* -------------------------------------------------------------------------- */
+
+void SystemEvolution::evolve() {
+
+  for (UInt i = 0; i < nsteps; ++i) {
+
+    for (auto & compute : computes)
+      compute->compute(*system);
+
+    if (i % freq == 0) {
+      std::stringstream sstr;
+      sstr << "dumps/step-" << std::setfill('0') << std::setw(5) << i << ".csv";
+      CsvWriter dumper(sstr.str());
+      dumper.write(*system);
+    }
+  }
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+void SystemEvolution::addCompute(const std::shared_ptr<Compute>& compute) {
+  computes.push_back(compute);
+}
+
+/* -------------------------------------------------------------------------- */
+void SystemEvolution::setNSteps(UInt nsteps) { this->nsteps = nsteps; }
+/* -------------------------------------------------------------------------- */
+void SystemEvolution::setDumpFreq(UInt freq) { this->freq = freq; }
+/* -------------------------------------------------------------------------- */
+System& SystemEvolution::getSystem() { return *system; }
+/* -------------------------------------------------------------------------- */
diff --git a/work/week9/particles/solution/system_evolution.hh b/work/week9/particles/solution/system_evolution.hh
new file mode 100644
index 0000000..3b33052
--- /dev/null
+++ b/work/week9/particles/solution/system_evolution.hh
@@ -0,0 +1,42 @@
+#ifndef __SYSTEM_EVOLUTION__HH__
+#define __SYSTEM_EVOLUTION__HH__
+
+/* -------------------------------------------------------------------------- */
+#include "compute.hh"
+#include "system.hh"
+/* -------------------------------------------------------------------------- */
+
+//! Manager for system evolution
+class SystemEvolution {
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! Construct using existing system (takes ownership)
+  SystemEvolution(std::unique_ptr<System> system);
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! Evolve all time steps
+  void evolve();
+  //! Add compute to list of computes
+  void addCompute(const std::shared_ptr<Compute>& compute);
+  //! Get the system object
+  System& getSystem();
+
+  void setNSteps(UInt nsteps);
+  void setDumpFreq(UInt freq);
+
+  /* ------------------------------------------------------------------------ */
+  /* Members                                                                  */
+  /* ------------------------------------------------------------------------ */
+protected:
+  std::vector<std::shared_ptr<Compute>> computes;
+  std::unique_ptr<System> system;
+  UInt nsteps, freq;
+};
+
+/* -------------------------------------------------------------------------- */
+#endif  //__SYSTEM_EVOLUTION__HH__
diff --git a/work/week9/particles/solution/test_csv.cc b/work/week9/particles/solution/test_csv.cc
new file mode 100644
index 0000000..6805466
--- /dev/null
+++ b/work/week9/particles/solution/test_csv.cc
@@ -0,0 +1,142 @@
+#include <gtest/gtest.h>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+//#include "my_types.cc"
+/*#include "compute_gravity.cc"
+#include "compute_verlet_integration.cc"
+#include "csv_reader.cc"
+#include "csv_writer.cc"
+#include "my_types.hh"
+#include "planets_factory.cc"
+#include "particles_factory_interface.cc"
+#include "system.cc"*/
+
+
+#include "compute_gravity.hh"
+#include "compute_verlet_integration.hh"
+#include "csv_reader.hh"
+#include "csv_writer.hh"
+#include "my_types.hh"
+#include "planets_factory.hh"
+#include "particles_factory_interface.hh"
+#include "system.hh"
+
+
+/* -------------------------------------------------------------------------- */
+
+/* Run command:
+    cd build
+    g++ ../test_csv.cc ../csv_reader.cc ../csv_writer.cc ../system.cc ../particle.cc ../planet.cc ../vector.cc -o test -lgtest_main -lgtest -lpthread -std=c++14
+    ./test
+*/
+
+struct PlanetStruct {
+    std::array<Real, 3> position;
+    std::array<Real, 3> velocity;
+    std::array<Real, 3> force;
+    Real mass;
+    std::string name;
+};
+
+// Fixture class
+class TestCsv : public ::testing::Test {
+protected:
+
+    void SetUp() override {
+        // Create planet structure with fields corresponding to csv file content
+        planet_test.position = {1, 2, 3};
+        planet_test.velocity = {5, 0, 6};
+        planet_test.force = {-1, -1, -2};
+        planet_test.mass = 10;
+        planet_test.name = "test";
+
+        // // Create new system and read planet from file
+        // system = new System();
+        // CsvReader reader("test_reader.csv");
+        // reader.read(*system);
+        std::cout << "Step0" << std::endl;
+        PlanetsFactory::getInstance();
+/*        ParticlesFactoryInterface& factory = ParticlesFactoryInterface::getInstance();
+        SystemEvolution& evol = factory.createSimulation(filename, timestep);
+*/
+        // Extract explicit planet object
+        // planet = dynamic_cast<Planet*>(&system->getParticle(0));
+    
+        
+        // // Write system into output file
+
+        //std::string outfile = "test_writer.csv";
+        // CsvWriter writer(outfile);
+        // writer.write(*system);
+
+        // Get 1st line of output file
+        // std::ifstream is(outfile.c_str());
+        // getline(is, line);
+        // is.close();
+    }
+
+    void TearDown() override {
+        // Clear and delete system
+        // system->removeParticle(0);
+        // delete system;
+        //evol.getSystem().removeParticle(0);
+       //delete factory;
+    }
+
+    //std::cout << "Step1" << std::endl;
+    std::string filename="test_reader.csv";
+    Real nsteps=1;
+    Real timestep=5e-3;
+    PlanetStruct planet_test;
+
+    // System* system;
+   // std::cout << "Step2" << std::endl;
+    //std::unique_ptr<SystemEvolution> system_evolution = nullptr;
+    ParticlesFactoryInterface& factory = ParticlesFactoryInterface::getInstance();
+    //SystemEvolution::SystemEvolution(std::unique_ptr<System> system)
+    //std::cout << "Step3" << std::endl;
+    SystemEvolution& evol = factory.createSimulation(filename, timestep);
+    //std::cout << "Step4" << std::endl;
+    Particle& planet = evol.getSystem().getParticle(0);
+    //std::cout << "Step5" << std::endl;
+    std::string line;
+};
+
+
+// Test CsvReader
+TEST_F(TestCsv, reading) {
+    EXPECT_EQ(planet_test.mass, planet->getMass());
+    for (UInt i = 0; i < 3; ++i) {
+        EXPECT_EQ(planet_test.position[i], planet->getPosition()[i]);
+        EXPECT_EQ(planet_test.velocity[i], planet->getVelocity()[i]);
+        EXPECT_EQ(planet_test.force[i], planet->getForce()[i]);
+    }
+    // EXPECT_EQ(planet_test.name, planet->getName());
+std::cout << "Step6" << std::endl;
+/*    EXPECT_EQ(planet_test.mass, planet.getMass());
+for (UInt i = 0; i < 3; ++i) {
+    std::cout << "Step7" << std::endl;
+        EXPECT_EQ(planet_test.position[i], planet.getPosition()[i]);
+        EXPECT_EQ(planet_test.velocity[i], planet.getVelocity()[i]);
+        EXPECT_EQ(planet_test.force[i], planet.getForce()[i]);
+    }*/
+    // EXPECT_EQ(planet_test.name, evol.getSystem().getParticle(0).getName());
+
+
+}
+
+
+int main(int argc, char **argv) {
+    std::cout << "StepMAIN" << std::endl;
+    testing::InitGoogleTest(&argc, argv);
+    std::cout << "StepTesting" << std::endl;
+    return RUN_ALL_TESTS();
+}
+
+// // Test CsvWriter
+// TEST_F(TestCsv, writing) {
+//     EXPECT_EQ(line, "1 2 3 5 0 6 -1 -1 -2 10 test ");
+// }
diff --git a/work/week9/particles/solution/test_test.cc b/work/week9/particles/solution/test_test.cc
new file mode 100644
index 0000000..8b1868d
--- /dev/null
+++ b/work/week9/particles/solution/test_test.cc
@@ -0,0 +1,21 @@
+#include <cmath>
+#include <exception>
+#include <gtest/gtest.h>
+
+struct NegativeException : public std::exception {};
+
+double mysqrt(double x) {
+
+  if (x < 0.)
+    throw NegativeException();
+
+  return std::sqrt(x);
+}
+
+// MACRO function "TEST" from Google Test
+// defines a test and a context
+TEST(MysqrtTest, positive_integers) {
+  EXPECT_EQ(2, mysqrt(4));
+  EXPECT_EQ(4, mysqrt(16));
+  EXPECT_EQ(3, mysqrt(9));
+}
\ No newline at end of file
diff --git a/work/week9/particles/solution/vector.cc b/work/week9/particles/solution/vector.cc
new file mode 100644
index 0000000..822dae1
--- /dev/null
+++ b/work/week9/particles/solution/vector.cc
@@ -0,0 +1,131 @@
+#include "vector.hh"
+
+Real& Vector::operator[](UInt i) {
+
+  return values[i];
+
+}
+const Real& Vector::operator[](UInt i) const {
+
+  return values[i];
+
+}
+
+Real Vector::squaredNorm() const {
+
+  Real res = 0;
+  for (auto& val : values) {
+    res += val * val;
+  }
+  return res;
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+Vector& Vector::operator+=(const Vector& vec) {
+
+  for (UInt i = 0; i < dim; ++i)
+    values[i] += vec[i];
+
+  return *this;
+}
+
+Vector& Vector::operator-=(const Vector& vec) {
+
+  for (UInt i = 0; i < dim; ++i)
+    values[i] -= vec[i];
+
+  return *this;
+}
+
+Vector& Vector::operator*=(Real val) {
+
+  for (auto& v : values)
+    v *= val;
+
+  return *this;
+}
+
+Vector& Vector::operator/=(Real val) {
+
+  for (auto& v : values)
+    v /= val;
+
+  return *this;
+}
+
+/* -------------------------------------------------------------------------- */
+
+Vector& Vector::operator=(const Vector& vec) {
+
+  std::copy(vec.values.begin(), vec.values.end(), values.begin());
+
+  return *this;
+}
+
+Vector& Vector::operator=(Real val) {
+
+  std::fill(values.begin(), values.end(), val);
+
+  return *this;
+}
+
+/* -------------------------------------------------------------------------- */
+
+Vector operator+(const Vector& a, const Vector& b) {
+
+  Vector res(a);
+  res += b;
+  return res;
+
+}
+
+Vector operator-(const Vector& a, const Vector& b) {
+
+  Vector res(a);
+  res -= b;
+  return res;
+
+}
+
+Vector operator*(const Vector& a, Real val) {
+
+  Vector res(a);
+  res *= val;
+  return res;
+
+}
+
+Vector operator*(Real val, const Vector& a) {
+  return a * val;
+}
+
+Vector operator/(const Vector& a, Real val) {
+
+  Vector res(a);
+  res /= val;
+  return res;
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+/// standard output stream operator
+std::ostream& operator<<(std::ostream& stream, const Vector& _this) {
+
+  for (auto& v : _this.values)
+    stream << v << " ";
+
+  return stream;
+}
+
+/* -------------------------------------------------------------------------- */
+/// standard input stream operator
+std::istream& operator>>(std::istream& stream, Vector& _this) {
+
+  for (auto& v : _this.values)
+    stream >> v;
+
+  return stream;
+}
diff --git a/work/week9/particles/solution/vector.hh b/work/week9/particles/solution/vector.hh
new file mode 100644
index 0000000..1deb5b9
--- /dev/null
+++ b/work/week9/particles/solution/vector.hh
@@ -0,0 +1,72 @@
+#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 const 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};
+};
+
+/* -------------------------------------------------------------------------- */
+/* 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);
+
+#endif  //__VECTOR__HH__
diff --git a/work/week9/particles/starting_point/test_reader.csv b/work/week9/particles/starting_point/test_reader.csv
new file mode 100644
index 0000000..77c99b8
--- /dev/null
+++ b/work/week9/particles/starting_point/test_reader.csv
@@ -0,0 +1 @@
+1 2 3 5 0 6 -1 -1 -2 10 test
\ No newline at end of file