diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0d85fe72e3d81e15b12c27f0fef3085be9282355..683bc1d702c22478cfe65a64cc2f9ae4ea501381 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,6 +14,7 @@ find_package(glm REQUIRED)
 set(COMMON_SOURCES
     itterative_PCISPH.cu
     spawn_particles.cpp
+    Camera.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/../common/GL_utilities.c
 )
 
diff --git a/src/Camera.cpp b/src/Camera.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6f0035d3b48a24ddeed49c387ce7e30ff75d07ec
--- /dev/null
+++ b/src/Camera.cpp
@@ -0,0 +1,120 @@
+// Camera.cpp
+#include "Camera.h"
+#include <cstdio>
+#include <cstdlib>
+
+// Global pointer to the Camera instance for callback access
+static Camera* g_camera = nullptr;
+
+// Constructor implementation
+Camera::Camera(vec3 position, vec3 target, vec3 up, float fov_, float aspectRatio_, float zNear_, float zFar_)
+    : campos(position), forward(VectorSub(target, position)), upVector(up),
+      fov(fov_), aspectRatio(aspectRatio_), zNear(zNear_), zFar(zFar_),
+      rotationMatrix(IdentityMatrix())
+{
+    // Normalize the forward vector
+    float norm = sqrt(forward.x * forward.x + forward.y * forward.y + forward.z * forward.z);
+    forward = ScalarMult(forward, 1.0f / norm);
+}
+
+
+// Processes keyboard input to control the camera
+void Camera::processKey(unsigned char key)
+{
+    switch (key)
+    {
+        case 'w':
+            campos += ScalarMult(forward, 0.5f);
+            break;
+        case 's':
+            campos -= ScalarMult(forward, 0.5f);
+            break;
+        case 'a':
+            campos -= ScalarMult(CrossProduct(forward, upVector), 0.5f);
+            break;
+        case 'd':
+            campos += ScalarMult(CrossProduct(forward, upVector), 0.5f);
+            break;
+        case 'j':
+            rotateSide(0.03f);
+            break;
+        case 'l':
+            rotateSide(-0.03f);
+            break;
+        case 'i':
+            // Rotate upwards
+            rotateUp(0.05f);
+            break;
+        case 'k':
+            // Rotate downwards
+            rotateUp(-0.05f);
+            break;
+        case 32: // SPACEBAR
+            campos += ScalarMult(upVector, 1.0f);
+            break;
+        case 'c':
+            campos -= ScalarMult(upVector, 1.0f);
+            break;
+        case 'v':
+            // Toggle boundary display
+            toggleDisplayBorder();
+            break;
+        case 0x1b: // ESC
+            exit(0);
+    }
+}
+
+void Camera::reshape(int width, int height)
+{
+    glViewport(0, 0, width, height);
+    aspectRatio = static_cast<float>(width) / static_cast<float>(height);
+}
+
+mat4 Camera::getViewMatrix() const
+{
+    return lookAt(campos, VectorAdd(campos, forward), upVector);
+}
+
+mat4 Camera::getProjectionMatrix() const
+{
+    return perspective(fov, aspectRatio, zNear, zFar);
+}
+
+void Camera::rotateSide(float angle)
+{
+    mat3 rot = mat4tomat3(Ry(angle));
+    forward = MultMat3Vec3(rot, forward);
+}
+
+void Camera::rotateUp(float angle)
+{
+    mat3 rot = mat4tomat3(ArbRotate(CrossProduct(forward, upVector), angle));
+    forward = MultMat3Vec3(rot, forward);
+}
+
+void initCamera(Camera& camera)
+{
+    g_camera = &camera;
+
+    // Register GLUT callbacks
+    glutKeyboardFunc(KeyCallback);
+    glutReshapeFunc(ReshapeCallback);
+}
+
+// GLUT Keyboard callback that forwards input to the Camera instance
+void KeyCallback(unsigned char key, int x, int y)
+{
+    if (g_camera)
+    {
+        g_camera->processKey(key);
+    }
+}
+
+// GLUT Reshape callback that forwards reshape events to the Camera instance
+void ReshapeCallback(int width, int height)
+{
+    if (g_camera)
+    {
+        g_camera->reshape(width, height);
+    }
+}
diff --git a/src/Camera.h b/src/Camera.h
new file mode 100644
index 0000000000000000000000000000000000000000..6a62c34af03af951626561ee72d8f215b62bef3e
--- /dev/null
+++ b/src/Camera.h
@@ -0,0 +1,60 @@
+// Camera.h
+#ifndef CAMERA_H
+#define CAMERA_H
+
+#include "VectorUtils4.h" // Assuming this defines vec3, mat4, etc.
+#include <GL/glut.h>
+
+// Camera class encapsulating all camera-related data and methods
+class Camera {
+private:
+    // Camera attributes
+    vec3 campos;
+    vec3 forward;
+    vec3 upVector;
+
+    // Camera parameters
+    float fov;
+    float aspectRatio;
+    float zNear;
+    float zFar;
+
+    bool displayBorder = false;
+
+    // Rotation matrix
+    mat4 rotationMatrix;
+
+    // Internal methods for camera transformations
+    void rotateSide(float angle);
+    void rotateUp(float angle);
+
+public:
+    // Constructor
+    Camera(vec3 position, vec3 target, vec3 up, float fov, float aspectRatio, float zNear, float zFar);
+
+    // Processes keyboard input
+    void processKey(unsigned char key);
+
+    // Handles window reshape
+    void reshape(int width, int height);
+
+    // Retrieves the view matrix
+    mat4 getViewMatrix() const;
+
+    // Retrieves the projection matrix
+    mat4 getProjectionMatrix() const;
+
+    bool shouldDisplayBorder() const { return displayBorder; }
+    
+    void toggleDisplayBorder() { displayBorder = !displayBorder; }
+
+};
+
+// Initialization function to set up camera and register callbacks
+void initCamera(Camera& camera);
+
+// GLUT callback functions
+void KeyCallback(unsigned char key, int x, int y);
+void ReshapeCallback(int width, int height);
+
+#endif // CAMERA_H
diff --git a/src/PCISPH_sim.h b/src/PCISPH_sim.h
index 28c678e7237fd90e6b2ab3f96925c61aac0eb83d..4ee42e64dfb6500182bc1b52b903338a6cef8975 100644
--- a/src/PCISPH_sim.h
+++ b/src/PCISPH_sim.h
@@ -33,12 +33,15 @@ void neigbourHoodSearch(
         int &gridSizeZ
 ); 
 
-void initializeDeviceArrays(
+void initializeFluidDeviceArrays(
     const float3* h_positions,
     const float3* h_velocities,
+    int numParticles
+);
+
+void initializeBoundryDeviceArrays(
     const float3* h_boundaryPositions,
     const float3* h_boundaryVelocities,
-    int numParticles,
     int numBoundaryParticles
 );
 
diff --git a/src/itterative_PCISPH.cu b/src/itterative_PCISPH.cu
index 4fb4c3644cb0fd4e69dae4cbe4d8c83eef87b2ac..6ab2f69723b038957fd576160b2da81bfd250284 100644
--- a/src/itterative_PCISPH.cu
+++ b/src/itterative_PCISPH.cu
@@ -793,13 +793,10 @@ void unregisterBoundaryVBOFromCUDA()
 }
 
 // Initialize device arrays
-void initializeDeviceArrays(
+void initializeFluidDeviceArrays(
     const float3* h_positions,
     const float3* h_velocities,
-    const float3* h_boundaryPositions,
-    const float3* h_boundaryVelocities,
-    int numParticles,
-    int numBoundaryParticles
+    int numParticles
 )
 {
     // Allocate device memory for velocities
@@ -837,8 +834,21 @@ void initializeDeviceArrays(
     // Allocate device memory for sum_gradW2
     cudaMalloc(&d_sum_gradW2, numParticles * sizeof(float));
 
-    // Allocate bounndary poarticle arrayss
+    // Initialize positions_pred to initial positions
+    cudaMemcpy(d_positions_pred, h_positions, numParticles * sizeof(float3), cudaMemcpyHostToDevice);
+}
 
+void initializeBoundryDeviceArrays(
+    const float3* h_boundaryPositions,
+    const float3* h_boundaryVelocities,
+    int numBoundaryParticles
+)
+{
+    // Allocate cell start and end arrays
+    int gridSizeX = (int)ceil((h_BOUND_X_MAX - h_BOUND_X_MIN) / h_H);
+    int gridSizeY = (int)ceil((h_BOUND_Y_MAX - h_BOUND_Y_MIN) / h_H);
+    int gridSizeZ = (int)ceil((h_BOUND_Z_MAX - h_BOUND_Z_MIN) / h_H);
+    int maxCells = gridSizeX * gridSizeY * gridSizeZ;
 
     // cudaMalloc(&d_boundary_positions, numBoundaryParticles * sizeof(float3));
     cudaMalloc(&d_boundary_velocities, numBoundaryParticles * sizeof(float3));
@@ -851,11 +861,6 @@ void initializeDeviceArrays(
     
     cudaMalloc(&d_boundaryCellStart, maxCells * sizeof(unsigned int));
     cudaMalloc(&d_boundaryCellEnd, maxCells * sizeof(unsigned int));
-
-    
-    
-    // Initialize positions_pred to initial positions
-    cudaMemcpy(d_positions_pred, h_positions, numParticles * sizeof(float3), cudaMemcpyHostToDevice);
 }
 
 
diff --git a/src/main.cpp b/src/main.cpp
index 86addccbfcc62a974255041b55106d06d2a8042e..af5343cbabe0695aa00208722dd5ba92c466a8f7 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -15,12 +15,9 @@
 #include "GL_utilities.h"
 #include "PCISPH_sim.h" // Header for CUDA function declarations
 #include "spawn_particles.h"
+#include "Camera.h"
 #include "constants.h"
 
-//#include <glm/glm.hpp>
-//#include <glm/gtc/matrix_transform.hpp> // For transformation functions
-//#include <glm/gtc/type_ptr.hpp>
-
 #define MAIN
 #include "GL_utilities.h"
 #include "VectorUtils4.h"
@@ -36,12 +33,6 @@ const float pointScale = 10.0f;
 auto previousTime = std::chrono::high_resolution_clock::now();
 float deltaTime = 0.016f;
 
-bool displayBoundary = false;
-
-// Particle data (initialized to 0)
-std::vector<float3> h_positions(NUM_PARTICLES);
-std::vector<float3> h_velocities(NUM_PARTICLES);
-// Acceleration, density, and pressure will be recalculated each time; no need to store them on host
 
 // Shader and VBO handles
 GLuint shaderProgram;
@@ -59,33 +50,7 @@ GLuint pointScaleLoc;
 GLuint uColorLoc;
 
 // CAMEERA VARIABLES
-// Perspective
-float fov = 45.0f; // Field of view in degrees
-float aspectRatio = 800.0f / 600.0f; // Window width / height
-float zNear = 0.1f;
-float zFar = 500.0f;
-
-// Lookat
-vec3 eyePosition(5.0f, 15.0f, 20.0f);
-vec3 targetPosition(0.0f, 4.0f, 0.0f);
-vec3 upVector(0.0f, 1.0f, 0.0f);
-
-
-vec3 campos = {5.0f, 15.0f, 20.0f};
-vec3 forward = {-5.0f, -11.0f, -20.0f};
-mat4 m;
-
-/*
-void printMatrix(const glm::mat4& mat, const char* name)
-{
-    printf("%s:\n", name);
-    const float* m = (const float*)glm::value_ptr(mat);
-    for (int i = 0; i < 4; ++i)
-    {
-        printf("%f %f %f %f\n", m[i], m[i + 4], m[i + 8], m[i + 12]);
-    }
-}
-*/
+Camera* cameraPtr = nullptr;
 
 void printFrameStats() {
     // FPS and deltaTime accumulation
@@ -109,6 +74,60 @@ void printFrameStats() {
     }
 }
 
+void renderBorder()
+{
+    // === Render Boundary Particles ===
+    // Set color for boundary particles (e.g., red)
+    glUniform3f(uColorLoc, 0.29f, 1.0f, 0.57f); // Red color
+
+    // Bind boundary VBO and set up vertex attribute
+    glBindBuffer(GL_ARRAY_BUFFER, vbo_boundary);
+    glEnableVertexAttribArray(positionAttribLoc);
+    glVertexAttribPointer(
+        positionAttribLoc,  // Attribute location
+        3,                  // Number of components (x, y, z)
+        GL_FLOAT,           // Data type
+        GL_FALSE,           // Normalized
+        sizeof(float3),     // Stride (size of one position)
+        (void*)0            // Offset to the position data
+    );
+
+    // Draw boundary particles as points
+    glDrawArrays(GL_POINTS, 0, NUM_BOUNDARY_PARTICLES);
+
+    // Disable vertex attribute
+    glDisableVertexAttribArray(positionAttribLoc);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+
+void renderFluid()
+{
+    // === Render Fluid Particles ===
+    // Set color for fluid particles (e.g., blue)
+    glUniform3f(uColorLoc, 0.29, 0.573, 1.0); // Blue color
+
+    // Bind fluid VBO and set up vertex attribute
+    glBindBuffer(GL_ARRAY_BUFFER, vbo);
+    glEnableVertexAttribArray(positionAttribLoc);
+    glVertexAttribPointer(
+        positionAttribLoc,  // Attribute location
+        3,                  // Number of components (x, y, z)
+        GL_FLOAT,           // Data type
+        GL_FALSE,           // Normalized
+        sizeof(float3),     // Stride (size of one position)
+        (void*)0            // Offset to the position data
+    );
+
+    // Draw fluid particles as points
+    glDrawArrays(GL_POINTS, 0, NUM_PARTICLES);
+
+    // Disable vertex attribute
+    glDisableVertexAttribArray(positionAttribLoc);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+
 // Display callback for OpenGL
 void display()
 {
@@ -134,64 +153,21 @@ void display()
 
     // Use the shader program
     glUseProgram(shaderProgram);
-    
-    mat4 projectionMatrix = perspective(fov, aspectRatio, zNear, zFar);
+
+    // Update the projection and view matrices based on the camera
+    mat4 projectionMatrix = cameraPtr->getProjectionMatrix();
     glUniformMatrix4fv(projectionMatrixLoc, 1, GL_TRUE, projectionMatrix.m);
 
-    mat4 viewMatrix = lookAt(campos, VectorAdd(campos, forward), upVector);
+    mat4 viewMatrix = cameraPtr->getViewMatrix();
     glUniformMatrix4fv(viewMatrixLoc, 1, GL_TRUE, viewMatrix.m);
 
-
-    if (displayBoundary) {
-
-        // === Render Boundary Particles ===
-        // Set color for boundary particles (e.g., red)
-        glUniform3f(uColorLoc, 0.29f, 1.0f, 0.57f); // Red color
-
-        // Bind boundary VBO and set up vertex attribute
-        glBindBuffer(GL_ARRAY_BUFFER, vbo_boundary);
-        glEnableVertexAttribArray(positionAttribLoc);
-        glVertexAttribPointer(
-            positionAttribLoc,  // Attribute location
-            3,                  // Number of components (x, y, z)
-            GL_FLOAT,           // Data type
-            GL_FALSE,           // Normalized
-            sizeof(float3),     // Stride (size of one position)
-            (void*)0            // Offset to the position data
-        );
-
-        // Draw boundary particles as points
-        glDrawArrays(GL_POINTS, 0, NUM_BOUNDARY_PARTICLES);
-
-        // Disable vertex attribute
-        glDisableVertexAttribArray(positionAttribLoc);
-        glBindBuffer(GL_ARRAY_BUFFER, 0);
+    // Render boundary particles if enabled
+    if(cameraPtr->shouldDisplayBorder()) {
+        renderBorder();
     }
-    
-    // === Render Fluid Particles ===
-    // Set color for fluid particles (e.g., blue)
-    glUniform3f(uColorLoc, 0.29, 0.573, 1.0); // Blue color
-
-    // Bind fluid VBO and set up vertex attribute
-    glBindBuffer(GL_ARRAY_BUFFER, vbo);
-    glEnableVertexAttribArray(positionAttribLoc);
-    glVertexAttribPointer(
-        positionAttribLoc,  // Attribute location
-        3,                  // Number of components (x, y, z)
-        GL_FLOAT,           // Data type
-        GL_FALSE,           // Normalized
-        sizeof(float3),     // Stride (size of one position)
-        (void*)0            // Offset to the position data
-    );
-
-    // Draw fluid particles as points
-    glDrawArrays(GL_POINTS, 0, NUM_PARTICLES);
-
-    // Disable vertex attribute
-    glDisableVertexAttribArray(positionAttribLoc);
-    glBindBuffer(GL_ARRAY_BUFFER, 0);
-
 
+    renderFluid();
+    
     // Swap buffers
     glutSwapBuffers();
 
@@ -215,49 +191,63 @@ void initializeBoundaryVBO(const std::vector<float3>& h_boundary_positions) {
     registerBoundaryVBOWithCUDA(vbo_boundary);
 }
 
-void initCamera()
+void initBoundary()
 {
-    modelMatrixLoc = glGetUniformLocation(shaderProgram, "modelMatrix");
-    viewMatrixLoc = glGetUniformLocation(shaderProgram, "viewMatrix");
-    projectionMatrixLoc = glGetUniformLocation(shaderProgram, "projectionMatrix");
+    std::vector<float3> boundary_positions;
+    std::vector<float3> boundary_velocities;
+    int boundaryLayers = 3;
+    float spacing = h_H/1.5;
 
+    createBoundaryParticles(boundary_positions, boundary_velocities, boundaryLayers, spacing);
 
-    // Camera setup
-    mat4 viewMatrix = lookAt(campos,  VectorAdd(campos, forward), upVector);
+    NUM_BOUNDARY_PARTICLES = (int)boundary_positions.size();
+    printf("Num Bounndnary particles: %i \n", NUM_BOUNDARY_PARTICLES);
+    // Initialize device arrays
+    initializeBoundryDeviceArrays(
+        boundary_positions.data(),
+        boundary_velocities.data(),
+        NUM_BOUNDARY_PARTICLES
+    );
 
-    // Create Projection Matrix using glm::perspective
-    mat4 projectionMatrix = perspective(fov, aspectRatio, zNear, zFar);
+    initializeBoundaryVBO(boundary_positions);
 
-    // Model Matrix (Identity matrix)
-    mat4 modelMatrix = IdentityMatrix();
-    
-    // Set uniform matrices
-    glUniformMatrix4fv(modelMatrixLoc, 1, GL_TRUE, modelMatrix.m);
-    glUniformMatrix4fv(viewMatrixLoc, 1, GL_TRUE, viewMatrix.m);
-    glUniformMatrix4fv(projectionMatrixLoc, 1, GL_TRUE, projectionMatrix.m);
+    boundaryNBSearch(NUM_BOUNDARY_PARTICLES);
 }
 
-
-// Initialize OpenGL and CUDA
-void initialize()
+void initParticleSystem()
 {
-    // Initialize constants in CUDA
-    initializeConstants();
-
-    // Set glewExperimental to enable modern OpenGL features
-    glewExperimental = GL_TRUE;
+    // Particle data (initialized to 0)
+    std::vector<float3> h_positions(NUM_PARTICLES);
+    std::vector<float3> h_velocities(NUM_PARTICLES);
+    
+    float girdSpacing = h_H/1.5;
+    // Initialize particle data
+    initParticles(
+        h_positions,
+        h_velocities,
+        NUM_PARTICLES,
+        girdSpacing
+    );
+    
+    initializeFluidDeviceArrays(
+        h_positions.data(),
+        h_velocities.data(),
+        NUM_PARTICLES
+    );
 
-    // Initialize GLEW
-    GLenum err = glewInit();
-    if (err != GLEW_OK) {
-        fprintf(stderr, "Error initializing GLEW: %s\n", glewGetErrorString(err));
-        exit(1);
-    }
+    // Initialize VBO for particle positions
+    glGenBuffers(1, &vbo);
+    glBindBuffer(GL_ARRAY_BUFFER, vbo);
+    // Upload initial positions to the VBO
+    glBufferData(GL_ARRAY_BUFFER, h_positions.size() * sizeof(float3), h_positions.data(), GL_DYNAMIC_DRAW);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
 
-    // Load and compile shaders
-    shaderProgram = loadShaders("../../shaders/particle.vert", "../../shaders/particle.frag");
-    glUseProgram(shaderProgram);
+    // Register VBO with CUDA
+    registerVBOWithCUDA(vbo);
+}
 
+void initShaderUniforms()
+{
     // Query the attribute location for position
     positionAttribLoc = glGetAttribLocation(shaderProgram, "aPosition");
     if (positionAttribLoc == -1) {
@@ -266,7 +256,12 @@ void initialize()
         printf("Attribute 'aPosition' location: %d\n", positionAttribLoc);
     }
 
-    // Query uniform locations
+    // Camera uniform locs
+    modelMatrixLoc = glGetUniformLocation(shaderProgram, "modelMatrix");
+    viewMatrixLoc = glGetUniformLocation(shaderProgram, "viewMatrix");
+    projectionMatrixLoc = glGetUniformLocation(shaderProgram, "projectionMatrix");
+    
+    // Other uniform locs
     pointRadiusLoc = glGetUniformLocation(shaderProgram, "pointRadius");
     pointScaleLoc = glGetUniformLocation(shaderProgram, "pointScale");
     uColorLoc = glGetUniformLocation(shaderProgram, "uColor");
@@ -285,49 +280,48 @@ void initialize()
         fprintf(stderr, "Failed to find 'pointScale' uniform location!\n");
     }
 
-    initCamera();
-
-    float parSpacing = h_H/1.5;
-    // Initialize particle data
-    initParticles(
-        h_positions,
-        h_velocities,
-        NUM_PARTICLES,
-        parSpacing
-    );
+    // Initialize model matrix (to Identity matrix)
+    mat4 modelMatrix = IdentityMatrix();
+    glUniformMatrix4fv(modelMatrixLoc, 1, GL_TRUE, modelMatrix.m);
+}
 
-    // Initialize VBO for particle positions
-    glGenBuffers(1, &vbo);
-    glBindBuffer(GL_ARRAY_BUFFER, vbo);
-    // Upload initial positions to the VBO
-    glBufferData(GL_ARRAY_BUFFER, h_positions.size() * sizeof(float3), h_positions.data(), GL_DYNAMIC_DRAW);
-    glBindBuffer(GL_ARRAY_BUFFER, 0);
+// Initialize OpenGL and CUDA
+void init()
+{
+    // Initialize constants in CUDA
+    initializeConstants();
 
-    // Register VBO with CUDA
-    registerVBOWithCUDA(vbo);
+    // Set glewExperimental to enable modern OpenGL features
+    glewExperimental = GL_TRUE;
 
-    std::vector<float3> boundary_positions;
-    std::vector<float3> boundary_velocities;
-    int boundaryLayers = 3;
-    float spacing = h_H/1.5;
+    // Initialize GLEW
+    GLenum err = glewInit();
+    if (err != GLEW_OK) {
+        fprintf(stderr, "Error initializing GLEW: %s\n", glewGetErrorString(err));
+        exit(1);
+    }
 
-    createBoundaryParticles(boundary_positions, boundary_velocities, boundaryLayers, spacing);
+    // Load and compile shaders
+    shaderProgram = loadShaders("../../shaders/particle.vert", "../../shaders/particle.frag");
+    glUseProgram(shaderProgram);
 
-    NUM_BOUNDARY_PARTICLES = (int)boundary_positions.size();
-    printf("Num Bounndnary particles: %i \n", NUM_BOUNDARY_PARTICLES);
-    // Initialize device arrays
-    initializeDeviceArrays(
-        h_positions.data(),
-        h_velocities.data(),
-        boundary_positions.data(),
-        boundary_velocities.data(),
-        NUM_PARTICLES,
-        NUM_BOUNDARY_PARTICLES
+    initShaderUniforms();
+
+    // Create Camera
+    cameraPtr = new Camera(
+        vec3(5.0f, 15.0f, 20.0f),  // Position
+        vec3(0.0f, 4.0f, 0.0f),    // Target
+        vec3(0.0f, 1.0f, 0.0f),    // Up vector
+        45.0f,                      // FOV
+        800.0f / 600.0f,            // Aspect ratio
+        0.1f,                       // zNear
+        500.0f                      // zFar
     );
 
-    initializeBoundaryVBO(boundary_positions);
-
-    boundaryNBSearch(NUM_BOUNDARY_PARTICLES);
+    initCamera(*cameraPtr);
+    
+    initParticleSystem();
+    initBoundary();
 
     // Enable point size for rendering particles
     glEnable(GL_PROGRAM_POINT_SIZE);
@@ -348,77 +342,6 @@ void cleanup()
 }
 
 
-void Key(unsigned char key,
-         __attribute__((unused)) int x,
-         __attribute__((unused)) int y)
-{
-    
-  switch (key)
-  {
-    case 'w':
-        campos += ScalarMult(forward, 0.01);
-        break;
-    
-    case 's':
-        campos -= ScalarMult(forward, 0.01);
-        break;
-
-    case 'a':
-		campos -= ScalarMult(CrossProduct(forward, upVector), 0.01);
-        break;
-    
-        
-    case 'd':
-		campos += ScalarMult(CrossProduct(forward, upVector), 0.01);
-        break; 
-
-        case 'j':
-        forward = MultMat3Vec3(mat4tomat3(Ry(0.03)), forward);
-        break;
-    
-    case 'l':
-        forward = MultMat3Vec3(mat4tomat3(Ry(-0.03)), forward);
-        break;
-
-    
-    // NOTE: Looking up and down is done by making a side vector and rotation around arbitrary axis!
-	// Note: Press L to rotate camera upwards and press K to look downwards. 
-    case 'i':
-		m = ArbRotate(CrossProduct(forward, upVector), 0.05);
-		forward = MultMat3Vec3(mat4tomat3(m), forward);
-        break;
-        
-    case 'k':
-		m = ArbRotate(CrossProduct(forward, upVector), -0.05);
-		forward = MultMat3Vec3(mat4tomat3(m), forward);
-        break;
-
-    // Move up/down with SPACEBAR and C
-    case 32:
-        campos += ScalarMult(upVector, 1.0);
-        break;
-
-    case 'c':
-        campos -= ScalarMult(upVector, 1.0);
-        break;
-
-    case 'v':
-        displayBoundary = !displayBoundary;
-        break;
-
-    case 0x1b:
-        exit(0);
-  }
-}
-
-void Reshape(int h, int v)
-{
-	glViewport(0, 0, h, v);
-
-    aspectRatio = h / v;
-}
-
-
 int main(int argc, char** argv)
 {
     // Initialize GLUT
@@ -429,13 +352,11 @@ int main(int argc, char** argv)
     glClearColor(0.9f, 0.9f, 0.9f, 1.0f);
 
     // Initialize OpenGL and CUDA
-    initialize();
+    init();
 
     // Register display callback
     glutDisplayFunc(display);
     glutIdleFunc(display);
-    glutKeyboardFunc(Key);
-    glutReshapeFunc(Reshape);
 
     // Cleanup on exit
     atexit(cleanup);