diff --git a/CMakeLists.txt b/CMakeLists.txt
index 11cd58154817f1e557fd4e154661a534f38eb16b..2546307a2f5b9749ba41e27e862654e9b0a1226a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,7 +25,7 @@ add_executable(${PROJECT_NAME}
         src/simulation/sph/sph_simulation.cpp
         src/simulation/sph/sph_simulation.h
 
-        )
+        src/rendering/Drawable.cpp src/rendering/Drawable.h)
 
 
 # opengl
diff --git a/src/application.cpp b/src/application.cpp
index 52c1ce1268f8e16d80d8e7617fe2ea3f4442dfea..1b1e60a59451932e04937fb0553e916dc585b9a0 100644
--- a/src/application.cpp
+++ b/src/application.cpp
@@ -1,6 +1,6 @@
 #include "application.h"
 
-InputHandler* Application::inputHandler = InputHandler::GetInstance();
+InputHandler *Application::inputHandler = InputHandler::GetInstance();
 bool Application::dragging = false;
 
 Application::Application() {
@@ -11,13 +11,13 @@ Application::~Application() {
     glfwDestroyWindow(window);
 }
 
-GLFWwindow* Application::get_window() {
+GLFWwindow *Application::get_window() {
     return window;
 }
 
 void Application::init() {
     this->set_callbacks();
-    
+
     glViewport(0, 0, this->width, this->height);
     glMatrixMode(GL_MODELVIEW);
     glEnable(GL_DEPTH_TEST);
@@ -28,7 +28,7 @@ void Application::init() {
     camera->glSetupCamera();
 }
 
-void Application::error_callback(int error, const char* description) {
+void Application::error_callback(int error, const char *description) {
     std::cerr << "Error: " << description << "." << std::endl;
 }
 
@@ -44,15 +44,15 @@ void Application::resize_callback(GLFWwindow *window, int w, int h) {
     glMatrixMode(GL_MODELVIEW);
 }
 
-void Application::key_callback(   GLFWwindow *window,
-                                        int key, 
-                                        int scancode,
-                                        int action, 
-                                        int mods) {
+void Application::key_callback(GLFWwindow *window,
+                               int key,
+                               int scancode,
+                               int action,
+                               int mods) {
 //    std::cout << key << " was pressed" << std::endl;
-    if(key == GLFW_KEY_ESCAPE)
+    if (key == GLFW_KEY_ESCAPE)
         glfwSetWindowShouldClose(window, GLFW_TRUE);
-    if(action == GLFW_RELEASE)
+    if (action == GLFW_RELEASE)
         InputHandler::KeyRelease(key);
     else
         InputHandler::KeyPress(key);
@@ -61,10 +61,10 @@ void Application::key_callback(   GLFWwindow *window,
 }
 
 
-void Application::mouse_click_callback( GLFWwindow *window,
-                                int button, 
-                                int action, 
-                                int mods) {
+void Application::mouse_click_callback(GLFWwindow *window,
+                                       int button,
+                                       int action,
+                                       int mods) {
     switch (button) {
         case GLFW_MOUSE_BUTTON_1:
             dragging = (action == GLFW_PRESS);
@@ -74,21 +74,21 @@ void Application::mouse_click_callback( GLFWwindow *window,
 
 
 void Application::mouse_motion_callback(GLFWwindow *window,
-                                double x, 
-                                double y) {
+                                        double x,
+                                        double y) {
     printf("mouse");
 }
 
 void Application::set_callbacks() {
-    glfwSetErrorCallback(       error_callback);
-    glfwSetWindowSizeCallback(  window, resize_callback);
-    glfwSetKeyCallback(         window, key_callback);
-    glfwSetMouseButtonCallback( window, mouse_click_callback);
-    glfwSetCursorPosCallback(   window, mouse_motion_callback);
+    glfwSetErrorCallback(error_callback);
+    glfwSetWindowSizeCallback(window, resize_callback);
+    glfwSetKeyCallback(window, key_callback);
+    glfwSetMouseButtonCallback(window, mouse_click_callback);
+    glfwSetCursorPosCallback(window, mouse_motion_callback);
 }
 
 void Application::start() {
-    while(!glfwWindowShouldClose(this->window)) {
+    while (!glfwWindowShouldClose(this->window)) {
         this->run();
         glfwSwapBuffers(window);
         glfwPollEvents();
@@ -103,7 +103,7 @@ void Application::run() {
     double delta_time = glfwGetTime() - time_since_last_frame;
     time_since_last_frame = glfwGetTime();
 
-    if(time_since_last_frame >= 1.0 / 24.0) {
+    if (time_since_last_frame >= 1.0 / 24.0) {
         glfwSetTime(0.0f);
         time_since_last_frame = 0.0f;
         tick = true;
@@ -115,7 +115,7 @@ void Application::run() {
     camera->control(delta_time);
 
     //clear color and depth buffer
-    glClearColor(0,0,0,1);
+    glClearColor(0, 0, 0, 1);
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
     glLoadIdentity(); //load identity matrix
     camera->glSetupCamera();
diff --git a/src/application.h b/src/application.h
index b1a3ec339260f3bf6fb2dbf592cf69f05cbab77b..97aa239c79c3ca594fb5173c9df370b4854603cb 100644
--- a/src/application.h
+++ b/src/application.h
@@ -1,5 +1,6 @@
 #ifndef MADID_APPLICATION_H
 #define MADID_APPLICATION_H
+
 #include <GL/glew.h>
 #include <GLFW/glfw3.h>
 #include <glm/glm.hpp>
@@ -20,13 +21,18 @@ class Application {
     int height = HEIGHT;
 
     static bool dragging;
+
     // Callback functions
-    static void error_callback(int error, const char* description);
+    static void error_callback(int error, const char *description);
+
     static void resize_callback(GLFWwindow *window, int w, int h);
+
     static void key_callback(GLFWwindow *window, int key, int scancode,
-            int action, int mods);
+                             int action, int mods);
+
     static void mouse_click_callback(GLFWwindow *window, int button, int action,
-            int mods);
+                                     int mods);
+
     static void mouse_motion_callback(GLFWwindow *window, double x, double y);
 
     std::shared_ptr<Scene> scene = std::make_shared<Scene>();
@@ -34,6 +40,7 @@ class Application {
     static InputHandler *inputHandler;
 
     void run();
+
     double time_since_last_frame;
 
     std::shared_ptr<Camera> camera = std::make_shared<Camera>();
@@ -41,18 +48,21 @@ class Application {
     // TODO
     //boolean flag, indicates whether a screen capture event should be performed
     //on the next display event.
-    /* bool screen_capture; */ 
-    /* bool zoom_mode; */ 
+    /* bool screen_capture; */
+    /* bool zoom_mode; */
     /* bool pan_mode; */
     /* bool trackball_mode; */
 public:
     Application();
+
     ~Application();
-    
-    GLFWwindow *window; 
-    
+
+    GLFWwindow *window;
+
     GLFWwindow *get_window();
+
     void init();
+
     void set_callbacks();
 
     void start();
diff --git a/src/constants.h b/src/constants.h
index a34aed83aad8bc8fa75c5325ea7287fe9eef3e29..21a56937adb51d90514f9978efae1c9e776cd6b7 100644
--- a/src/constants.h
+++ b/src/constants.h
@@ -4,6 +4,6 @@
 const int WIDTH = 500;
 const int HEIGHT = 500;
 
-char * const NAME = "Madid";
+char *const NAME = "Madid";
 
 #endif
diff --git a/src/main.cpp b/src/main.cpp
index 203a60cf2af6f0d09e89ae49ef733419aaf8d3de..26176b7f220f21146947a8220c17939ff582c086 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,14 +1,22 @@
-// Math
-#include <glm/glm.hpp>
+#include <GL/glew.h>
+#include <GLFW/glfw3.h>
 
 #include <iostream>
 #include <memory>
 
 #include "application.h"
-#include "constants.h"
 
 
 int main(int argc, char **argv) {
+    // start GLEW extension handler
+    glewExperimental = GL_TRUE;
+    glewInit();
+
+    // get version info
+    const GLubyte *renderer = glGetString(GL_RENDERER); // get renderer string
+    const GLubyte *version = glGetString(GL_VERSION); // version as a string
+    std::cerr << "Renderer: " << renderer << std::endl;
+    std::cerr << "OpenGL version supported: " << version << std::endl;
 
     // start GL context and O/S window using the GLFW helper library
     if (!glfwInit()) {
@@ -16,30 +24,21 @@ int main(int argc, char **argv) {
         return 1;
     }
 
+    glfwInit();
+
     /* Creating a GLFW application */
     /* Application creates its GLFW window too */
     std::shared_ptr<Application> application = std::make_shared<Application>();
-    
-    if(!application->window) {
+
+
+    if (!application->window) {
         std::cerr << "ERROR: could not open window with GLFW3" << std::endl;
         glfwTerminate();
         return 1;
     }
-    
-    glfwMakeContextCurrent(application->get_window());
-    
-    // start GLEW extension handler
-    glewExperimental = GL_TRUE;
-    glewInit();
-
-    // get version info
-    const GLubyte *renderer = glGetString(GL_RENDERER); // get renderer string
-    const GLubyte *version = glGetString(GL_VERSION); // version as a string
-    std::cerr << "Renderer: " << renderer << std::endl;
-    std::cerr << "OpenGL version supported: " << version << std::endl;
 
+    glfwMakeContextCurrent(application->get_window());
     application->init();
-    
     application->start();
 
 }
diff --git a/src/rendering/Drawable.cpp b/src/rendering/Drawable.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0bdbd11b7ec9e6669c3859c162ca4c193889fd6e
--- /dev/null
+++ b/src/rendering/Drawable.cpp
@@ -0,0 +1,14 @@
+#include <cstdio>
+#include "Drawable.h"
+
+Drawable::Drawable() {
+    glGenVertexArrays(1, &vao);
+    glBindVertexArray(vao);
+    glGenBuffers(1, &vbo);
+    glBindBuffer(GL_ARRAY_BUFFER, vbo);
+}
+
+Drawable::~Drawable() {
+    glDeleteBuffers(1, &vbo);
+    glDeleteVertexArrays(1, &vao);
+}
diff --git a/src/rendering/Drawable.h b/src/rendering/Drawable.h
new file mode 100644
index 0000000000000000000000000000000000000000..b7e14e95e574be60d10e95c2e6b2c2df58ef9b22
--- /dev/null
+++ b/src/rendering/Drawable.h
@@ -0,0 +1,19 @@
+#ifndef MADID_DRAWABLE_H
+#define MADID_DRAWABLE_H
+
+#include "../utility/gl.h"
+
+class Drawable {
+protected:
+    GLuint vao, vbo;
+
+public:
+    Drawable();
+
+    virtual void Draw() = 0;
+
+    ~Drawable();
+};
+
+
+#endif //MADID_DRAWABLE_H
diff --git a/src/rendering/camera.cpp b/src/rendering/camera.cpp
index 1685c5c9396acb96e3b30cfd0b8ca2c7ba350f2a..8b860f868e6fa407d75b3f17bbfcb11978b76ee9 100644
--- a/src/rendering/camera.cpp
+++ b/src/rendering/camera.cpp
@@ -1,13 +1,12 @@
 #include <iostream>
 #include "camera.h"
 
-Camera::Camera():
-    aspect_ratio(1),
-    pos(1, 1, -3),
-    forward(0,0,1),
-    right(-1,0,0),
-    speed(10) 
-{
+Camera::Camera() :
+        aspect_ratio(1),
+        pos(1, 1, -3),
+        forward(0, 0, 1),
+        right(-1, 0, 0),
+        speed(10) {
     inputHandler = InputHandler::GetInstance();
 }
 
@@ -17,37 +16,37 @@ Camera::~Camera() {
 
 void Camera::control(float delta_time) {
     // process camera keys
-    if(InputHandler::IsPressed(GLFW_KEY_W))
+    if (InputHandler::IsPressed(GLFW_KEY_W))
         pos += forward * speed * delta_time;
 
 
-    if(InputHandler::IsPressed(GLFW_KEY_S))
+    if (InputHandler::IsPressed(GLFW_KEY_S))
         pos -= forward * speed * delta_time;
 
-    if(InputHandler::IsPressed(GLFW_KEY_A))
+    if (InputHandler::IsPressed(GLFW_KEY_A))
         pos -= right * speed * delta_time;
 
-    if(InputHandler::IsPressed(GLFW_KEY_D))
+    if (InputHandler::IsPressed(GLFW_KEY_D))
         pos += right * speed * delta_time;
 
-    if(InputHandler::IsPressed(GLFW_KEY_Q))
-        pos += glm::vec3(0,-1,0) * speed * delta_time;
+    if (InputHandler::IsPressed(GLFW_KEY_Q))
+        pos += glm::vec3(0, -1, 0) * speed * delta_time;
 
-    if(InputHandler::IsPressed(GLFW_KEY_E))
-        pos += glm::vec3(0,1,0) * speed * delta_time;
+    if (InputHandler::IsPressed(GLFW_KEY_E))
+        pos += glm::vec3(0, 1, 0) * speed * delta_time;
 
     std::cout << pos.x << pos.y << pos.z << std::endl;
 }
 
 void Camera::startDrag(int x, int y) {
-   drag_start = glm::vec2(0,0);
+    drag_start = glm::vec2(0, 0);
 }
 
 void Camera::drag(int x, int y) {
-    glm::vec2 d(drag_start.x - (float)x, drag_start.y - (float)y);
-    
+    glm::vec2 d(drag_start.x - (float) x, drag_start.y - (float) y);
+
     forward = glm::normalize(forward);
-    right = cross(forward, glm::vec3(0,1,0));
+    right = cross(forward, glm::vec3(0, 1, 0));
     right = glm::normalize(right);
 
     // TODO 
@@ -63,7 +62,7 @@ void Camera::glSetupCamera() const {
     glLoadIdentity();
     gluLookAt(
             pos.x, pos.y, pos.z,
-            pos.x + forward.x, pos.y + forward.y, pos.z + forward.z, 
+            pos.x + forward.x, pos.y + forward.y, pos.z + forward.z,
             0.0f, 1.0f, 0.0f);
 }
 
diff --git a/src/rendering/camera.h b/src/rendering/camera.h
index 185641506f16690eb3b187c67e795c17adb72bcf..d2d189ec0642fd5f04665bc011887c26e7f715a3 100644
--- a/src/rendering/camera.h
+++ b/src/rendering/camera.h
@@ -12,21 +12,24 @@ class Camera {
     glm::vec3 right;
     float speed;
 
-    InputHandler* inputHandler;
+    InputHandler *inputHandler;
 
-    glm::vec2 drag_start = glm::vec2(0,0);
+    glm::vec2 drag_start = glm::vec2(0, 0);
 public:
     void setAspectRatio(float ar) {
         aspect_ratio = ar;
     }
 
     void control(float delta_time);
+
     void startDrag(int x, int y);
+
     void drag(int x, int y);
 
     void glSetupCamera() const;
 
     Camera();
+
     ~Camera();
 
 };
diff --git a/src/scene.cpp b/src/scene.cpp
index befdae600b433669797b35357ebd5dd72ed08697..52d36e474a7c6aae8d99bfa4bbc0aa079dd40c93 100644
--- a/src/scene.cpp
+++ b/src/scene.cpp
@@ -9,7 +9,7 @@ Scene::~Scene() {
 }
 
 void Scene::update(double dt) {
-    sph_system->update((float)dt);
+    sph_system->update((float) dt);
 }
 
 void Scene::render() {
diff --git a/src/scene.h b/src/scene.h
index f06af3ae792078c1dbf2de13b1c9c117ca7f485e..e19b83060c63fe3d201b1e52fdc693073457fcec 100644
--- a/src/scene.h
+++ b/src/scene.h
@@ -11,16 +11,20 @@ class Scene {
     bool draw_normals;
     bool draw_surface;
     bool draw_velocity;
-//    bool render;
     bool use_emitter;
 
-    std::shared_ptr<sph::SPHSimulation> sph_system = std::make_shared<sph::SPHSimulation>();
-    
+    std::shared_ptr<SPHSimulation> sph_system = std::make_shared<SPHSimulation>();
+
+    // TODO set parameters here
+    std::shared_ptr<PBDSimulation> pbd_system = std::make_shared<PBDSimulation>();
+
 public:
     Scene();
+
     ~Scene();
 
     void update(double dt);
+
     void render();
 };
 
diff --git a/src/simulation/particle.cpp b/src/simulation/particle.cpp
index 7da253bf6b88f3b6ba2067a7f87bc89decb1fc0a..fc9acd9ccee9768b23e8260f7b92aef4bb808b81 100644
--- a/src/simulation/particle.cpp
+++ b/src/simulation/particle.cpp
@@ -1,82 +1,25 @@
 #include "particle.h"
 
-namespace sph {
-
-    Particle::Particle(float x, float y, float z):
+Particle::Particle(float x, float y, float z) :
         pos(x, y, z),
         old_pos(pos),
         radius(1), //TODO
-        radius_sqr(1*1), //TODO
+        radius_sqr(1 * 1), //TODO
         v(0, 0, 0),
-        a(0,0,0),
-        f(0,0,0),
-        n(0,0,0),
-        m(1),
+        a(0, 0, 0),
+        f(0, 0, 0),
+        n(0, 0, 0),
+        w(0), // inverse mass, 1/inf is 0
         rho(1),
         p(0),
-        fixed(false)
-    {
-        
-    }
+        fixed(false) {
+}
 
-    Particle::~Particle(){}
+Particle::~Particle() {}
 
-//    // will need these for hashing
-//    Particle::point_type min() const {
-//        return Particle::point_type(pos.x - Particle::point_type(radius));
-//    }
-//
-//    Particle::point_type max() const {
-//        return Particle::point_type(x + Particle::point_type(radius));
-//    }
-//
-//    // checks if two particles (actually two positions) are within an
-//    // allowed radius.
-//    bool check(const point_type& candidate_pos) const{
-//        Particle::point_type diff(pos - candidate_pos);
-//        return radius_sqrt >= diff*diff;
-//    }
-//    // Position
-//    const point_type& position() const;
-//    point_type& position();
-//
-//    //Old Position
-//    const point_type& old_pos() const { return old_pos; }
-//    point_type& old_pos() { return old_pos; }
-//
-//    // Velocity
-//    const vector &velocity() const { return v; }
-//    vector &velocity() { return v; }
-//
-//    // Acceleration
-//    const vector &acceleration() const { return a; }
-//    vector &acceleration() { return a; }
-//
-//    // Force
-//    const vector &force() const { return f; }
-//    vector &force() { return f; }
-//
-//    // Mass
-//    const float &mass() const { return m; }
-//    float &mass() { return m; }
-//
-//    // Density
-//    const float &density() const { return rho; }
-//    float &density() { return rho; }
-//
-//    // Pressure
-//    const float &pressure() const { return p; }
-//    float &pressure() { return p; }
-//
-//    // Normal vector of the surface
-//    const vector &normal() const { return n; }
-//    vector &normal() { return n; }
-//
-//    // Fixed flag of the particle
-//    const bool &fixed() const { return fixed; }
-//    bool &fixed() { return fixed; }
-//
-//    // Radius of the particle
-//    const float &radius() const { return radius; }
+Particle::Particle(Particle::vector pos, float w, Particle::vector color) :
+        pos(pos),
+        w(w),
+        color(color) {
 
 }
diff --git a/src/simulation/particle.h b/src/simulation/particle.h
index b2a26472b61a2116ca00ea54cfd9ecb5b424a1c3..624e3e148bb508670d94478caa35ceea7b2d3599 100644
--- a/src/simulation/particle.h
+++ b/src/simulation/particle.h
@@ -4,136 +4,32 @@
 #include <glm/glm.hpp>
 
 class Particle {
-    public:
-        typedef glm::vec3 vector;
-        typedef vector point_type;
+public:
+    typedef glm::vec3 vector;
+    typedef vector point_type;
 
-        Particle(float x, float y, float z);
-        ~Particle();
+    Particle(float x, float y, float z);
 
-    //TODO proper scope, maybe getter functions
-    public:
-        point_type pos;     // position of the particle
-        point_type old_pos; // old position of the particle
-        float radius;       // radius of the particle
-        float radius_sqr;   // squared radius of the particle
-        vector v;           // current velocity
-        vector a;           // current acceleration
-        vector f;           // current sum of external forces acting on the particle
-        vector n;           // surface normal if close to/on the surface, else 0
-        float m;            // mass (constant)
-        float rho;          // density (varies)
-        float p;            // pressure (varies)
-        bool fixed;         // Is this particled fixed (locked)?
-    public:
-        // will need these for hashing
-//        point_type min() const;
-//        point_type max() const;
-
-        // checks if two particles (actually two positions) are within an
-        // allowed radius.
-//        bool check(const point_type& candidate) const;
-
-        // get position (read only)
-//        const point_type& position() const;
-//        // get position
-//        point_type& position();
-
-        /** Old Position (read only).
-        * @return  const reference to the old position vector.
-        */
-//        const point_type& old_pos() const;
-
-        /** Old Position.
-        * @return  reference to the old position vector.
-        */
-//        point_type& old_pos();
-
-        /** Velocity (read only).
-        * @return  const reference to the velocity vector.
-        */
-//        const vector &velocity() const;
-
-        /** Velocity.
-        * @return  reference to the velocity vector.
-        */
-//        vector &velocity();
-
-        /** Acceleration (read only).
-        * @return  const reference to the acceleration vector.
-        */
-//        const vector &acceleration() const;
-
-        /** Acceleration.
-        * @return  reference to the acceleration vector.
-        */
-//        vector &acceleration();
-
-        /** Force (read only).
-        * @return  const reference to the external force vector.
-        */
-//        const vector &force() const;
-
-        /** Force.
-        * @return  reference to the external force vector.
-        */
-//        vector &force();
+    Particle(vector pos, float w, vector color);
 
-        /** Mass (read only).
-        * @return  const reference to the mass scalar.
-        */
-//        const float &mass() const;
+    ~Particle();
 
-        /** Mass.
-        * @return  reference to the mass scalar.
-        */
-//        float &mass();
-
-        /** Density (read only).
-        * @return  const reference to the density value.
-        */
-//        const float &density() const;
-
-        /** Density.
-        * @return  reference to the density value.
-        */
-//        float &density();
-
-        /** Pressure (read only).
-        * @return  const reference to the pressure value.
-        */
-//        const float &pressure() const;
-
-        /** Pressure.
-        * @return  reference to the pressure value.
-        */
-//        float &pressure();
-
-        /** Normal (read only).
-        * @return  const reference to the surface normal.
-        */
-//        const vector &normal() const;
-
-        /** Normal.
-        * @return  reference to the surface normal.
-        */
-//        vector &normal();
-
-        /** Fixed (read only).
-        * @return  const reference to the surface normal.
-        */
-//        const bool &fixed() const;
-
-        /** Fixed.
-        * @return  reference to the surface normal.
-        */
-//        bool &fixed();
-
-        /** Fixed.
-        * @return  reference to the surface normal.
-        */
-//        const float &radius() const;
-        
+    //TODO proper scope, maybe getter functions
+public:
+    point_type pos;     // position of the particle
+    point_type old_pos; // old position of the particle
+    float radius;       // radius of the particle
+    float radius_sqr;   // squared radius of the particle
+    vector v;           // current velocity
+    vector a;           // current acceleration
+    vector f;           // current sum of external forces acting on the particle
+    vector n;           // surface normal if close to/on the surface, else 0
+    float w;            // inverse mass (constant)
+    vector tmp_pos;     // temporary position (for PBD simulation)
+    vector color;       // (r,g,b) color of the particle
+    float rho;          // density (varies)
+    float p;            // pressure (varies)
+    bool fixed;         // Is this particled fixed (locked)?
 };
 
 #endif
diff --git a/src/simulation/pbd/pbd_simulation.cpp b/src/simulation/pbd/pbd_simulation.cpp
index 3b4b0b44dc50183a3dc8377645ab5abb73fee37a..369b77d2c315f67af6c8c173e48765dcb40ffa38 100644
--- a/src/simulation/pbd/pbd_simulation.cpp
+++ b/src/simulation/pbd/pbd_simulation.cpp
@@ -1,8 +1,16 @@
 #include "pbd_simulation.h"
+
 void PBDSimulation::addForce(vec3 force) {
     externalForces += force;
 }
 
+PBDSimulation::PBDSimulation() {
+    numberOfParticlesOnASide = 10;
+    nrSegments = 10;
+    lSeg = 0.08f;
+    initParticles();
+}
+
 PBDSimulation::PBDSimulation(size_t _nr_sims, size_t _nr_segments, float _l_seg) :
         numberOfParticlesOnASide(_nr_sims),
         nrSegments(_nr_segments),
@@ -18,9 +26,29 @@ PBDSimulation::PBDSimulation(size_t _nr_sims, size_t _nr_segments, float _l_seg)
 }
 
 void PBDSimulation::initParticles() {
+
+    vec3 startPos = vec3(-1, 0, -1);
+
     for (size_t i = 0; i < numberOfParticlesOnASide; i++) {
         vec3 color = vec3(222.0f, 101.0f, 32.0f);
-        strands.emplace_back(CreateStrand(nrSegments, lSeg, currPos.position * vec3(1, -1, 1), color));
+        std::vector<Particle *> currentStrand;
+
+        for (size_t j = 0; j < numberOfParticlesOnASide; j++) {
+            // mass of the current particle
+            float m = .25f;
+            vec3 currPos = startPos + vec3((float) i * lSeg, 0, (float) j * lSeg);
+
+            // infinite mass on the corners
+            if ((i == 0 && j == 0) ||
+                (i == 0 && j == numberOfParticlesOnASide - 1) ||
+                (i == numberOfParticlesOnASide - 1 && j == 0) ||
+                (i == numberOfParticlesOnASide - 1 && j == numberOfParticlesOnASide - 1)
+                    ) {
+                currentStrand.push_back(new Particle(currPos, 0, color));
+            } else currentStrand.push_back(new Particle(currPos, 1 / m, color));
+        }
+
+        strands.emplace_back(currentStrand);
     }
 }
 
@@ -89,8 +117,10 @@ void PBDSimulation::solve_bending_constraint(Particle *p1, Particle *p2, float d
                 (length(p1->tmp_pos - p2->tmp_pos) - dist) *
                 (p1->tmp_pos - p2->tmp_pos) / length(p1->tmp_pos - p2->tmp_pos);
 
-    p1->tmp_pos += 0.6 * d_p1;
-    p2->tmp_pos += 0.6 * d_p2;
+    // TODO define *(float, vec3) operator...
+    float damping = 0.6f;
+    p1->tmp_pos += vec3(d_p1.x * damping, d_p1.y * damping, d_p1.z * damping);
+    p2->tmp_pos += vec3(d_p2.x * damping, d_p2.y * damping, d_p2.z * damping);
 }
 
 void PBDSimulation::solve_collision_constraint(Particle *p, vec3 &q1, vec3 &q2, vec3 &q3) {
@@ -145,27 +175,6 @@ vec3 PBDSimulation::getExternalForces() const {
     return externalForces;
 }
 
-std::vector<Particle *> PBDSimulation::CreateStrand(size_t n, float l, vec3 startPos, vec3 color) {
-    vec3 currPos = startPos;
-    std::vector<Particle *> currentStrand;
-
-    for (size_t i = 0; i < n; i++) {
-        // random mass
-        // A single strand of hair can weigh between .2 – .5 milligrams,
-        // which is 2.0e-7 kg == 10^(-7) kg
-        float m = util::randomOffsetf(.35f, .15f);
-
-        // first particle's position is infinite
-        if (i == 0) currentStrand.push_back(new Particle(currPos, 0, color));
-        else currentStrand.push_back(new Particle(currPos, 1 / m, color));
-
-        // propagate particles downward
-        currPos.y -= l;
-    }
-
-    return currentStrand;
-}
-
 void PBDSimulation::resetExternalForces() {
     externalForces = vec3(0, 0, 0);
 }
diff --git a/src/simulation/pbd/pbd_simulation.h b/src/simulation/pbd/pbd_simulation.h
index c601f5a0dcfff3995f7b5ac55fe81fc8cbd564e5..e23b4a1bfe74d2bdfcc227da3529568569a578ca 100644
--- a/src/simulation/pbd/pbd_simulation.h
+++ b/src/simulation/pbd/pbd_simulation.h
@@ -5,11 +5,14 @@
 #include <vector>
 #include <glm/glm.hpp>
 
+#include "../../utility/gl.h"
+
 #include "../particle.h"
+#include "../../rendering/Drawable.h"
 
 typedef glm::vec3 vec3;
 
-class PBDSimulation {
+class PBDSimulation : Drawable {
     void solve_distance_constraint(Particle *p1, Particle *p2, float dist);
 
     void solve_bending_constraint(Particle *p1, Particle *p2, float dist);
@@ -37,6 +40,8 @@ public:
 
     void addForce(vec3 force);
 
+    PBDSimulation();
+
     PBDSimulation(size_t _nr_sims, size_t _nr_segments, float _l_seg);
 
     void initParticles();
@@ -48,8 +53,6 @@ public:
     vec3 getExternalForces() const;
 
     void resetExternalForces();
-
-    std::vector<Particle *> CreateStrand(size_t segments, float l, vec3 startPos, vec3 color);
 };
 
 
diff --git a/src/simulation/sph/sph_simulation.cpp b/src/simulation/sph/sph_simulation.cpp
index fcf42dcda885636a4332454cc1305f2f6ff7f12d..75fb76959c60c4bffe43d7bdb8fa0eac061aaed4 100644
--- a/src/simulation/sph/sph_simulation.cpp
+++ b/src/simulation/sph/sph_simulation.cpp
@@ -1,153 +1,152 @@
 #include "sph_simulation.h"
 
-    SPHSimulation::SPHSimulation() {
-        std::cout << "Initializing sph particle system with " <<sph::N<<
-            " particles." << std::endl;
+SPHSimulation::SPHSimulation() {
+    std::cout << "Initializing sph particle system with " << sph::N <<
+              " particles." << std::endl;
 
-        for(float y = 1; y < 5; y += 0.5f)
-            for(float z = 0; z < 2; z += 0.2f)
-                for (float x = 0; x < 2; x += 0.2f)
-                    if (particles.size() < sph::N)
-                    {
+    for (float y = 1; y < 5; y += 0.5f)
+        for (float z = 0; z < 2; z += 0.2f)
+            for (float x = 0; x < 2; x += 0.2f)
+                if (particles.size() < sph::N) {
 
-                        float x_off = glm::linearRand(.0f, 1.f);
+                    float x_off = glm::linearRand(.0f, 1.f);
 //                        particles.emplace_back(new Particle(x + x_off, y, z));
-                        particles.emplace_back(
-                                new Particle(
-                                        glm::linearRand(1.f, 1.1f),
-                                        glm::linearRand(1.f, 1.1f),
-                                        glm::linearRand(1.f, 1.1f)
-                                        )
-                                );
-                    }
-
-        std::cout << "There are " << particles.size() << "particles." <<
-            std::endl;
-    }
-
-    SPHSimulation::~SPHSimulation() {
-        for(auto & particle : particles)
-            delete particle;
-
-    }
+                    particles.emplace_back(
+                            new Particle(
+                                    glm::linearRand(1.f, 1.1f),
+                                    glm::linearRand(1.f, 1.1f),
+                                    glm::linearRand(1.f, 1.1f)
+                            )
+                    );
+                }
 
-    void SPHSimulation::compute_density_and_pressure() {
-        for(auto &p_i : particles) {
-            // compute density
-            p_i->rho = 0.f;
-            for(auto &p_j : particles) {
-                glm::vec3 r_ij = p_j->pos - p_i->pos;
-                // euclidean distance between p_i and p_j
-                float r = glm::l2Norm(p_j->pos, p_i->pos);
-                float r2 = glm::pow(r, 2.f);
-                
-                // 0 <= r <= d
-                if (r2 <= sph::H2)
-                    p_i->rho += sph::M * sph::POLY6 * glm::pow(sph::H2 - r2, 3.f);
-            }
-            
-            // compute pressure
-            p_i->p = sph::GAS_CONST * (p_i->rho - sph::REST_DENS);
+    std::cout << "There are " << particles.size() << "particles." <<
+              std::endl;
+}
+
+SPHSimulation::~SPHSimulation() {
+    for (auto &particle : particles)
+        delete particle;
+
+}
+
+void SPHSimulation::compute_density_and_pressure() {
+    for (auto &p_i : particles) {
+        // compute density
+        p_i->rho = 0.f;
+        for (auto &p_j : particles) {
+            glm::vec3 r_ij = p_j->pos - p_i->pos;
+            // euclidean distance between p_i and p_j
+            float r = glm::l2Norm(p_j->pos, p_i->pos);
+            float r2 = glm::pow(r, 2.f);
+
+            // 0 <= r <= d
+            if (r2 <= sph::H2)
+                p_i->rho += sph::M * sph::POLY6 * glm::pow(sph::H2 - r2, 3.f);
         }
 
+        // compute pressure
+        p_i->p = sph::GAS_CONST * (p_i->rho - sph::REST_DENS);
     }
 
-    void SPHSimulation::compute_forces() {
-        for (auto &p_i : particles) {
-            glm::vec3 f_pressure(0.f, 0.f, 0.f);
-            glm::vec3 f_viscosity(0.f, 0.f, 0.f);
+}
 
-            for(auto &p_j : particles) {
-                if(&p_i == &p_j) continue;
+void SPHSimulation::compute_forces() {
+    for (auto &p_i : particles) {
+        glm::vec3 f_pressure(0.f, 0.f, 0.f);
+        glm::vec3 f_viscosity(0.f, 0.f, 0.f);
 
-                glm::vec3 r_ijv = p_i->pos - p_j->pos;
-                float r_ij = glm::length(r_ijv);
-                
-                if(r_ij <= sph::H) {
-                    // compute pressure force contribution
+        for (auto &p_j : particles) {
+            if (&p_i == &p_j) continue;
 
-                    f_pressure += -sph::M * (p_i->p + p_j->p) / (2*p_j->rho) *
-                        sph::POLY6_GRAD_PRESS * r_ijv/r_ij * glm::pow(sph::H-r_ij, 2.0f);
+            glm::vec3 r_ijv = p_i->pos - p_j->pos;
+            float r_ij = glm::length(r_ijv);
 
-                    f_viscosity += sph::MU * sph::M * (p_j->v - p_i->v) / p_j->rho *
-                        sph::POLY6_GRAD_VISC * (sph::H - r_ij);
-                }
+            if (r_ij <= sph::H) {
+                // compute pressure force contribution
 
-                glm::vec3 f_gravity = sph::G * p_i->rho;
+                f_pressure += -sph::M * (p_i->p + p_j->p) / (2 * p_j->rho) *
+                              sph::POLY6_GRAD_PRESS * r_ijv / r_ij * glm::pow(sph::H - r_ij, 2.0f);
 
-                p_i->f = f_pressure + f_viscosity + f_gravity;
-                p_i->f = f_gravity + f_pressure;// * 0.001f;
+                f_viscosity += sph::MU * sph::M * (p_j->v - p_i->v) / p_j->rho *
+                               sph::POLY6_GRAD_VISC * (sph::H - r_ij);
             }
+
+            glm::vec3 f_gravity = sph::G * p_i->rho;
+
+            p_i->f = f_pressure + f_viscosity + f_gravity;
+            p_i->f = f_gravity + f_pressure;// * 0.001f;
         }
     }
+}
 
-    void SPHSimulation::integrate() {
-        for (auto &p : particles) {
-            // forward euler
-            p->v += sph::DT * p->f / p->rho;
-            p->pos += sph::DT * p->v;
-
-            const float max_b = 1.2f;
-            const float min_b = 0.9f;
+void SPHSimulation::integrate() {
+    for (auto &p : particles) {
+        // forward euler
+        p->v += sph::DT * p->f / p->rho;
+        p->pos += sph::DT * p->v;
 
-            //boundary conditions
-            if (p->pos.x - sph::EPS <= min_b) {
-                p->v.x *= sph::BOUND_DAMPING;
-                p->pos.x = min_b + sph::EPS;
-            }
+        const float max_b = 1.2f;
+        const float min_b = 0.9f;
 
-            if(p->pos.x + sph::EPS > max_b) {
-                p->v.x *= sph::BOUND_DAMPING;
-                p->pos.x = max_b - sph::EPS;
-            }
+        //boundary conditions
+        if (p->pos.x - sph::EPS <= min_b) {
+            p->v.x *= sph::BOUND_DAMPING;
+            p->pos.x = min_b + sph::EPS;
+        }
 
-            // ground
-            if(p->pos.y - sph::EPS <= 0.0f) {
-                p->v.y *= sph::BOUND_DAMPING;
-                p->pos.y = sph::EPS;
-            }
+        if (p->pos.x + sph::EPS > max_b) {
+            p->v.x *= sph::BOUND_DAMPING;
+            p->pos.x = max_b - sph::EPS;
+        }
 
-            if(p->pos.y + sph::EPS >= 51) {
-                p->v.y *= sph::BOUND_DAMPING;
-                p->pos.y = 51 - sph::EPS;
-            }
+        // ground
+        if (p->pos.y - sph::EPS <= 0.0f) {
+            p->v.y *= sph::BOUND_DAMPING;
+            p->pos.y = sph::EPS;
+        }
 
-            //TODO DEPTH CONSTANT VALUE
-            if(p->pos.z - sph::EPS <= min_b) {
-                p->v.z *= sph::BOUND_DAMPING;
-                p->pos.z = min_b + sph::EPS;
-            }
+        if (p->pos.y + sph::EPS >= 51) {
+            p->v.y *= sph::BOUND_DAMPING;
+            p->pos.y = 51 - sph::EPS;
+        }
 
-            if(p->pos.z + sph::EPS >= max_b) {
-                p->v.z *= sph::BOUND_DAMPING;
-                p->pos.z = max_b - sph::EPS;
-            }
+        //TODO DEPTH CONSTANT VALUE
+        if (p->pos.z - sph::EPS <= min_b) {
+            p->v.z *= sph::BOUND_DAMPING;
+            p->pos.z = min_b + sph::EPS;
+        }
 
-            p->v *= 0.98;
+        if (p->pos.z + sph::EPS >= max_b) {
+            p->v.z *= sph::BOUND_DAMPING;
+            p->pos.z = max_b - sph::EPS;
         }
+
+        p->v *= 0.98;
     }
+}
 
-    void SPHSimulation::update(float t) {
-        compute_density_and_pressure();
-        compute_forces();
+void SPHSimulation::update(float t) {
+    compute_density_and_pressure();
+    compute_forces();
 
 //        if(t > sph::DT)
 //            for(float dt = 0.0f; dt < t; dt+=sph::DT)
 //                integrate();
 //        else
-            integrate();
-    }
+    integrate();
+}
 
-    void SPHSimulation::render() {
-        // TODO MODERN OPENGL!!!!
-        glEnable(GL_POINT_SMOOTH);
-        glPointSize( 1.f);
+void SPHSimulation::render() {
+    // TODO MODERN OPENGL!!!!
+    glEnable(GL_POINT_SMOOTH);
+    glPointSize(1.f);
 
-        glColor4f(0.5f, 0.5f, 0.8f, 1.f);
+    glColor4f(0.5f, 0.5f, 0.8f, 1.f);
 
-        glBegin(GL_POINTS);
-        for(auto &p : particles)
-            glVertex3f(p->pos.x, p->pos.y, p->pos.z);
-        glEnd();
-        
-    }
+    glBegin(GL_POINTS);
+    for (auto &p : particles)
+        glVertex3f(p->pos.x, p->pos.y, p->pos.z);
+    glEnd();
+
+}
diff --git a/src/simulation/sph/sph_simulation.h b/src/simulation/sph/sph_simulation.h
index 4d194ef399e6bb405ebe60fc5be64aada4279f3f..30deb8c2fe68ccedc7a0aad5652e3438b4eef99b 100644
--- a/src/simulation/sph/sph_simulation.h
+++ b/src/simulation/sph/sph_simulation.h
@@ -40,9 +40,9 @@ namespace sph {
     // smoothing kernel as described in the Müller paper
     const static float POLY6 = 315.f / (64.f * glm::pi<float>() * glm::pow(H, 9));
     const static float POLY6_GRAD_PRESS = -45.f / (glm::pi<float>()
-            * glm::pow(H, 6));
+                                                   * glm::pow(H, 6));
     const static float POLY6_GRAD_VISC = 45.f / (glm::pi<float>() * glm::pow(H,
-                6));
+                                                                             6));
 
     // PARTICLES
     // number of particles
@@ -52,15 +52,17 @@ namespace sph {
 
 class SPHSimulation {
 
-    std::vector<Particle*> particles;
-    
+    std::vector<Particle *> particles;
+
     void compute_density_and_pressure();
+
     void compute_forces();
 
     void integrate();
 
 public:
     SPHSimulation();
+
     ~SPHSimulation();
 
     // step the system by t time
diff --git a/src/utility/InputHandler.cpp b/src/utility/InputHandler.cpp
index 92a2edf01d99b998e2c92b51e9ee2886e97c1489..9566f56b32446984f382bcb3e2108dd62a57375e 100644
--- a/src/utility/InputHandler.cpp
+++ b/src/utility/InputHandler.cpp
@@ -5,7 +5,7 @@ bool InputHandler::keyPressed[348];
 int InputHandler::modifiers = 0;
 
 InputHandler *InputHandler::GetInstance() {
-    if (singleton_ == nullptr) 
+    if (singleton_ == nullptr)
         singleton_ = new InputHandler();
 
     return singleton_;
@@ -31,7 +31,7 @@ bool InputHandler::IsShiftPressed() {
     return modifiers & GLFW_MOD_SHIFT;
 }
 
-bool InputHandler::IsControlPressed()  {
+bool InputHandler::IsControlPressed() {
     return modifiers & GLFW_MOD_CONTROL;
 }
 
diff --git a/src/utility/InputHandler.h b/src/utility/InputHandler.h
index 12e892f2a462d7a563b5388b25466fee9bbdf921..50780de135e5b8efd58c0b55a2ed3c733958f263 100644
--- a/src/utility/InputHandler.h
+++ b/src/utility/InputHandler.h
@@ -50,7 +50,7 @@ public:
 
     static bool IsPressed(int key);
 
-    static bool IsShiftPressed() ;
+    static bool IsShiftPressed();
 
     static bool IsControlPressed();