diff --git a/CMakeLists.txt b/CMakeLists.txt
index ecec384db808031e69991f55307f26575f2c07e6..7b4c7e3e7e520b7f478367294ac86c176b5c3a32 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,6 +15,10 @@ add_executable(${PROJECT_NAME}
         src/scene.h
         src/rendering/camera.cpp
         src/rendering/camera.h
+        src/sph/sph_particle.h
+        src/sph/sph_particle.cpp
+        src/sph/sph_system.h
+        src/sph/sph_system.cpp
         )
 
 
diff --git a/src/application.cpp b/src/application.cpp
index a2e6d4792aa158a62f31a5726ac16910cf284a51..b29ba1d4d26230c7e0ce5bb647641a66c8bcf233 100644
--- a/src/application.cpp
+++ b/src/application.cpp
@@ -24,7 +24,7 @@ void Application::init() {
     glLoadIdentity();
     glClearColor(255.0f / 255.0f, 25.0f / 255.0f, 25.0f / 255.0f, 1.0f);
 
-    camera.glSetupCamera();
+    camera->glSetupCamera();
 }
 
 void Application::error_callback(int error, const char* description) {
@@ -106,8 +106,12 @@ void Application::run() {
     // TODO guarantee delta_time to be infinitesimal
     scene->update(delta_time);
     //clear color and depth buffer
+    glClearColor(0,0,0,1);
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
     glLoadIdentity(); //load identity matrix
     camera->glSetupCamera();
     scene->render();
 }
+
+
+
diff --git a/src/application.h b/src/application.h
index 45b638f29390b888b25b1b32eb924a14f39f4d30..ea2e272f3232028b511ceafd87c58cf8ce089a5a 100644
--- a/src/application.h
+++ b/src/application.h
@@ -9,6 +9,7 @@
 
 #include "constants.h"
 #include "scene.h"
+#include "rendering/camera.h"
 
 // TODO make this class a singleton
 class Application {
diff --git a/src/rendering/camera.cpp b/src/rendering/camera.cpp
index fa0b3fd0f6593f626fde8894c8221ad906dfb2da..429007c3d1ecefdb8dea7940817a1ef333fe7bd7 100644
--- a/src/rendering/camera.cpp
+++ b/src/rendering/camera.cpp
@@ -1,6 +1,6 @@
 #include "camera.h"
 
-Camera::Camera(void):
+Camera::Camera():
     aspect_ratio(1),
     pos(0, 25, -50),
     forward(0,0,1),
@@ -9,11 +9,11 @@ Camera::Camera(void):
 {
 }
 
-Camera::~Camera(void) {
+Camera::~Camera() {
 }
 
 
-void Camera::control(float delta_time, bool* inputs) {
+void Camera::control(float delta_time, const bool* inputs) {
     // process camera keys
     for(int i = 0; i < 256; i++) {
         if(!inputs[i]) continue;
@@ -23,18 +23,19 @@ void Camera::control(float delta_time, bool* inputs) {
             case 's': pos -= forward * speed * delta_time; break;
             case 'a': pos -= right * speed * delta_time; break;
             case 'd': pos += right * speed * delta_time; break;
-            case 'q': pos += glm::vec3(0,1,0) * speed * delta_time; break;
+            case 'q': pos += glm::vec3(0,-1,0) * speed * delta_time; break;
             case 'e': pos += glm::vec3(0,1,0) * speed * delta_time; break;
+            default: break;
         }
     }
 }
 
 void Camera::startDrag(int x, int y) {
-   drag_start = glm::vec3(0,0);
+   drag_start = glm::vec2(0,0);
 }
 
 void Camera::drag(int x, int y) {
-    glm::vec2 d(drag_start.x - x, drag_start.y - 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));
@@ -44,7 +45,7 @@ void Camera::drag(int x, int y) {
 }
 
 
-void Camera::glSetupCamera() {
+void Camera::glSetupCamera() const {
     glMatrixMode(GL_PROJECTION);
     glLoadIdentity();
     gluPerspective(45.0f, aspect_ratio, 1.0f, 1000.0f);
diff --git a/src/rendering/camera.h b/src/rendering/camera.h
index 8d2e91e7d1403f4e761064b63cae7fe27b6284c0..8c8ab523559e51343be40d7cc43caa2bc4f5a3eb 100644
--- a/src/rendering/camera.h
+++ b/src/rendering/camera.h
@@ -1,8 +1,8 @@
-#ifndef MADID_SCENE_H
-#define MADID_SCENE_H
+#ifndef MADID_CAMERA_H
+#define MADID_CAMERA_H
 
 #include <glm/glm.hpp>
-#include "utility/gl.h"
+#include "../utility/gl.h"
 
 class Camera {
     float aspect_ratio;
@@ -11,21 +11,21 @@ class Camera {
     glm::vec3 right;
     float speed;
 
-    glm::vec2 drag_start(0,0);
+    glm::vec2 drag_start = glm::vec2(0,0);
 public:
     void setAspectRatio(float ar) {
         aspect_ratio = ar;
     }
 
-    void control(float delta_time, bool* inputs);
+    void control(float delta_time, const bool* inputs);
     void startDrag(int x, int y);
     void drag(int x, int y);
 
-    void glSetupCamera();
+    void glSetupCamera() const;
 
-    Camera(void);
-    ~Camera(void);
+    Camera();
+    ~Camera();
 
-}
+};
 
 #endif
diff --git a/src/scene.cpp b/src/scene.cpp
index 84960c6bbc296b663192d8424aa6a32ceae504d9..befdae600b433669797b35357ebd5dd72ed08697 100644
--- a/src/scene.cpp
+++ b/src/scene.cpp
@@ -9,10 +9,9 @@ Scene::~Scene() {
 }
 
 void Scene::update(double dt) {
-    sph_simulaton
-    return;
+    sph_system->update((float)dt);
 }
 
 void Scene::render() {
-    return;
+    sph_system->render();
 }
diff --git a/src/scene.h b/src/scene.h
index 8841a87474cefe186e6a00c6e3acea464a9670a9..145272dcac61fb42cb4e8012153c16d0ded788ed 100644
--- a/src/scene.h
+++ b/src/scene.h
@@ -2,13 +2,18 @@
 #define MADID_SCENE_H
 
 #include <iostream>
+#include <memory>
+
+#include "sph/sph_system.h"
 
 class Scene {
     bool draw_normals;
     bool draw_surface;
     bool draw_velocity;
-    bool render;
+//    bool render;
     bool use_emitter;
+
+    std::shared_ptr<sph::System> sph_system = std::make_shared<sph::System>();
     
 public:
     Scene();
diff --git a/src/sph/sph_particle.cpp b/src/sph/sph_particle.cpp
index 0a26eb26dd875e24d7e9558c387504f12bf0abaa..c3d3c33708b3dcc3f1064174417ff784c129543b 100644
--- a/src/sph/sph_particle.cpp
+++ b/src/sph/sph_particle.cpp
@@ -1,75 +1,82 @@
-#include "scene.h"
-
-using namespace sph::Particle;
+#include "sph_particle.h"
 
 namespace sph {
-    Particle::Particle():
+
+    Particle::Particle(float x, float y, float z):
+        pos(x, y, z),
+        old_pos(pos),
         radius(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),
         rho(1),
         p(0),
         fixed(false)
-    {}
-
-    Particle::~Particle(){}
-
-    // will need these for hashing
-    Particle::point_type min() const {
-        return Particle::point_type(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; }
+    Particle::~Particle(){}
 
-    // Radius of the particle
-    const float &radius() const { return radius; }
+//    // 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; }
 
 }
diff --git a/src/sph/sph_particle.h b/src/sph/sph_particle.h
index 8e37c3b43a2fcd273a7a9a24f17449a46784a736..252bb811eef39454b3f37d372286c13c093a3416 100644
--- a/src/sph/sph_particle.h
+++ b/src/sph/sph_particle.h
@@ -9,12 +9,15 @@ class Particle {
         typedef glm::vec3 vector;
         typedef vector point_type;
 
-        Particle();
-    private:
-        float radius;       // radius of the particle
-        float radius_sqr;   // squared radius of the particle
+        Particle(float x, float y, float z);
+        ~Particle();
+
+    //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
@@ -25,114 +28,114 @@ class Particle {
         bool fixed;         // Is this particled fixed (locked)?
     public:
         // will need these for hashing
-        point_type min() const;
-        point_type max() const;
+//        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;
+//        bool check(const point_type& candidate) const;
 
         // get position (read only)
-        const point_type& position() const;
-        // get position
-        point_type& position();
+//        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;
+//        const point_type& old_pos() const;
 
         /** Old Position.
         * @return  reference to the old position vector.
         */
-        point_type& old_pos();
+//        point_type& old_pos();
 
         /** Velocity (read only).
         * @return  const reference to the velocity vector.
         */
-        const vector &velocity() const;
+//        const vector &velocity() const;
 
         /** Velocity.
         * @return  reference to the velocity vector.
         */
-        vector &velocity();
+//        vector &velocity();
 
         /** Acceleration (read only).
         * @return  const reference to the acceleration vector.
         */
-        const vector &acceleration() const;
+//        const vector &acceleration() const;
 
         /** Acceleration.
         * @return  reference to the acceleration vector.
         */
-        vector &acceleration();
+//        vector &acceleration();
 
         /** Force (read only).
         * @return  const reference to the external force vector.
         */
-        const vector &force() const;
+//        const vector &force() const;
 
         /** Force.
         * @return  reference to the external force vector.
         */
-        vector &force();
+//        vector &force();
 
         /** Mass (read only).
         * @return  const reference to the mass scalar.
         */
-        const float &mass() const;
+//        const float &mass() const;
 
         /** Mass.
         * @return  reference to the mass scalar.
         */
-        float &mass();
+//        float &mass();
 
         /** Density (read only).
         * @return  const reference to the density value.
         */
-        const float &density() const;
+//        const float &density() const;
 
         /** Density.
         * @return  reference to the density value.
         */
-        float &density();
+//        float &density();
 
         /** Pressure (read only).
         * @return  const reference to the pressure value.
         */
-        const float &pressure() const;
+//        const float &pressure() const;
 
         /** Pressure.
         * @return  reference to the pressure value.
         */
-        float &pressure();
+//        float &pressure();
 
         /** Normal (read only).
         * @return  const reference to the surface normal.
         */
-        const vector &normal() const;
+//        const vector &normal() const;
 
         /** Normal.
         * @return  reference to the surface normal.
         */
-        vector &normal();
+//        vector &normal();
 
         /** Fixed (read only).
         * @return  const reference to the surface normal.
         */
-        const bool &fixed() const;
+//        const bool &fixed() const;
 
         /** Fixed.
         * @return  reference to the surface normal.
         */
-        bool &fixed();
+//        bool &fixed();
 
         /** Fixed.
         * @return  reference to the surface normal.
         */
-        const float &radius() const;
+//        const float &radius() const;
         
-}
+};
 }
 
 #endif
diff --git a/src/sph/sph_system.cpp b/src/sph/sph_system.cpp
index 5b4056ccd1e58b044cdfc8d76a3b8a5f3e437fe7..ce0ef4cae119f121d7416e87f3c42bef0054be4f 100644
--- a/src/sph/sph_system.cpp
+++ b/src/sph/sph_system.cpp
@@ -1 +1,138 @@
 #include "sph_system.h"
+
+
+namespace sph {
+
+    System::System() {
+        std::cout << "Initializing sph particle system with " << N <<
+            " particles." << std::endl;
+
+        for(float y = 0; y < HEIGHT/3; y += H)         
+            for(float z = 0; z < 400; z += H)
+                for (float x = 0; x < WIDTH/4; x += H)
+                    if (particles.size() < N)
+                    {
+                        float x_off = glm::linearRand(.0f, 1.f);
+                        particles.emplace_back(new Particle(x + x_off, y, z));
+                    }
+
+        std::cout << "There are " << particles.size() << "particles." <<
+            std::endl;
+    }
+
+    System::~System() {
+        for(int i = 0; i < particles.size(); i++)
+            delete particles[i];
+    }
+
+    void System::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 r2 = glm::l2Norm(p_j->pos, p_i->pos);
+                
+                // 0 <= r <= d
+                if (r2 <= H2)
+                    p_i->rho += M * POLY6 * glm::pow(H2 - r2, 3.f);
+            }
+            
+            // compute pressure
+            p_i->p = GAS_CONST * (p_i->rho - REST_DENS);
+        }
+
+    }
+
+    void System::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;
+
+                float r_ij = glm::length(p_j->pos - p_i->pos);
+                
+                if(r_ij <= H) {
+                    // compute pressure force contribution
+                    f_pressure += -M * r_ij / (2*p_j->rho) *
+                        POLY6_GRAD_PRESS * glm::pow(H - r_ij, 2.f);
+                    f_viscosity += MU * M * (p_j->v - p_i->v) / p_j->rho *
+                        POLY6_GRAD_VISC * (H - r_ij);
+                }
+
+                glm::vec3 f_gravity = G * p_i->rho;
+
+                p_i->f = f_pressure + f_viscosity + f_gravity;
+            }
+        }
+    }
+
+    void System::integrate() {
+        for (auto &p : particles) {
+            // forward euler
+            p->v += DT * p->f / p->rho;
+            p->pos += DT * p->v;
+
+            //boundary conditions
+            if (p->pos.x - EPS < 0) {
+                p->v.x *= BOUND_DAMPING;
+                p->pos.x = 0.0f;
+                std::cout << "bottom: " << p->pos.x << std::endl;
+            }
+
+            if(p->pos.x + EPS > WIDTH) {
+                p->v.x *= BOUND_DAMPING;
+                p->pos.x = WIDTH - EPS;
+            }
+
+            if(p->pos.y - EPS > HEIGHT) {
+                p->v.y *= BOUND_DAMPING;
+                p->pos.y = HEIGHT - EPS;
+            }
+
+            if(p->pos.y + EPS > HEIGHT) {
+                p->v.y *= BOUND_DAMPING;
+                p->pos.y = HEIGHT - EPS;
+            }
+
+            //TODO DEPTH CONSTANT VALUE
+            if(p->pos.z - EPS < 0) {
+                p->v.z *= BOUND_DAMPING;
+                p->pos.z = EPS;
+            }
+
+            if(p->pos.z + EPS > 400) {
+                p->v.z *= BOUND_DAMPING;
+                p->pos.z = 400 - EPS;
+            }
+        }
+    }
+
+    void System::update(float t) {
+        compute_density_and_pressure();
+        compute_forces();
+
+        if(t > DT)
+            for(float dt = 0.0f; dt < t; dt+=DT)
+                integrate();
+        else
+            integrate();
+    }
+
+    void System::render() {
+
+        // TODO MODERN OPENGL!!!!
+        glEnable(GL_POINT_SMOOTH);
+        glPointSize(H / 2.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();
+        
+    }
+}
diff --git a/src/sph/sph_system.h b/src/sph/sph_system.h
index b11b7303bd3ad09bd2ed27709e5c3fd07df22836..b52b6817def228489443b9f6939f6ba11e581d03 100644
--- a/src/sph/sph_system.h
+++ b/src/sph/sph_system.h
@@ -2,14 +2,73 @@
 #define MADID_SPH_SYSTEM
 
 #include <glm/glm.hpp>
+#include <glm/gtc/random.hpp>
+#include <glm/gtx/norm.hpp>
+#include <iostream>
+#include <vector>
+
+#include "../utility/gl.h"
+
 #include "sph_particle.h"
+#include "../constants.h"
+
+namespace sph {
+    // external (gravitational force)
+    const static glm::vec3 G(0.f, -9.8f, 0.f);
+    // rest density
+    const static float REST_DENS = 1000.0f;
+    // const for equation of state
+    const static float GAS_CONST = 2000.0f;
+    // kernel radius
+    const static float H = 16.f;
+    // (kernel radius)^2 for optimization
+    const static float H2 = H * H;
+    // mass of particles
+    const static float M = 65.f;
+    // viscosity constant
+    const static float MU = 250.f; 
+    // TODO decide the place for this
+    // integration timestep
+    const static float DT = 0.0005f; 
+
+    // epsilon (e.g. at boundary check)
+    const static float EPS = H;
+    // damping when colliding with the boundary
+    const static float BOUND_DAMPING = -0.5f;
+
+    [[maybe_unused]] const static float PI = glm::pi<float>();
+    // smoothing kernel as described in the Müller paper
+    const static float POLY6 = 315.f / (64 * glm::pi<float>() * glm::pow(H, 9));
+    const static float POLY6_GRAD_PRESS = -45.f / (glm::pi<float>()
+            * glm::pow(H, 6));
+    const static float POLY6_GRAD_VISC = 45.f / (glm::pi<float>() * glm::pow(H,
+                6));
+
+    // PARTICLES
+    // number of particles
+    const static int N = 500;
+}
 
 namespace sph {
-    class System {
-        // particle container
-        // collision system
-        // l
-    }
+
+class System {
+
+    std::vector<Particle*> particles;
+    
+    void compute_density_and_pressure();
+    void compute_forces();
+
+    void integrate();
+
+public:
+    System();
+    ~System();
+
+    // step the system by t time
+    void update(float t);
+
+    void render();
+};
 }
 
 #endif