diff --git a/src/MakeFile_Erik.txt b/src/MakeFile_Alt.txt
similarity index 82%
rename from src/MakeFile_Erik.txt
rename to src/MakeFile_Alt.txt
index c27eeef0d9fd8ca53b291d4066c2b89b519e155a..75d9e148140da94e89569734cc9ea1433689d08f 100644
--- a/src/MakeFile_Erik.txt
+++ b/src/MakeFile_Alt.txt
@@ -1,3 +1,3 @@
 nvcc -c kernels_cuda.cu -o kernels_cuda.o -I/usr/local/cuda/include
-nvcc main2.cpp ../common/GL_utilities.c ../common/Linux/MicroGlut.c kernels_cuda.o -I/usr/local/cuda/include -I../common -I../common/Linux -DGL_GLEXT_PROTOTYPES -I/usr/include/GL -L/usr/local/cuda/lib64 -lGL -lGLU -lcudart -lX11 -lXt -lm -o simple_particle_system
+nvcc main3.cpp ../common/GL_utilities.c ../common/Linux/MicroGlut.c kernels_cuda.o -I/usr/local/cuda/include -I../common -I../common/Linux -DGL_GLEXT_PROTOTYPES -I/usr/include/GL -L/usr/local/cuda/lib64 -lGL -lGLU -lcudart -lX11 -lXt -lm -o simple_particle_system
 ./simple_particle_system
\ No newline at end of file
diff --git a/src/main2.cpp b/src/main2.cpp
deleted file mode 100644
index ce874a41fba2c59214e2bb2903913bcc94c9b091..0000000000000000000000000000000000000000
--- a/src/main2.cpp
+++ /dev/null
@@ -1,319 +0,0 @@
-// main.cpp: This example uses FreeGLUT for OpenGL and CUDA interoperation.
-// It creates a simple particle simulation and uses CUDA to update particle positions with gravity.
-
-#include "Particle.h"  // Include the Particle struct definition
-#include <vector>
-#include <cstdio> // For fprintf and stderr
-
-#ifdef _WIN32
-    #include <GL/glew.h>  // Include GLEW before other OpenGL headers
-#endif
-
-//#include <GL/freeglut.h>
-#include "kernels_cuda.h"  // Header for CUDA function declarations
-
-//#include <glm/glm.hpp>
-//#include <glm/gtc/matrix_transform.hpp> // For transformation functions
-//#include <glm/gtc/type_ptr.hpp> 
-
-#define MAIN
-#include "MicroGlut.h"
-#include "VectorUtils4.h"
-#include "GL_utilities.h"
-
-
-
-// Number of particles
-const int NUM_PARTICLES = 12800; // Should match a square grid for simplicity
-const int GRID_SIZE = 16;       // Number of particles per row and column in the grid
-const float GRID_SPACING = 1.2f;  // Spacing between particles in the grid
-const float PARTICLE_SIZE = 1.0f;
-
-const float pointRadius = 20.0f;  // Point size in world space
-const float pointScale = 30.0f;
-
-// Particle data (initialized to 0)
-std::vector<Particle> h_particles(NUM_PARTICLES);
-
-// Shader and VBO handles
-GLuint shaderProgram;
-GLuint positionAttribLoc;
-// Vertex Buffer Object (VBO) for particle positions
-GLuint vbo;
-
-// Uniform locations
-GLuint modelMatrixLoc;
-GLuint viewMatrixLoc;
-GLuint projectionMatrixLoc;
-GLuint pointRadiusLoc;
-GLuint pointScaleLoc;
-
-/*
-void printMatrix(const mat4& mat, const char* name)
-{
-    printf("%s:\n", name);
-    const float* m = static_cast<const float*>(*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]);
-    }
-}
-
-*/
-
-
-// Display callback for OpenGL
-void display()
-{
-
-    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-
-    glEnable(GL_BLEND);
-    //glBlendFunc(GL_SRC_ALPHA, GL_ONE); // Additive blending
-    glDepthMask(GL_FALSE);
-
-    // Update particle positions using CUDA
-    updateParticles();
-
-    // Use the shader program
-    glUseProgram(shaderProgram);
-
-    // Bind VBO and set up vertex attribute for position
-    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(Particle),                   // Stride (size of one Particle struct)
-        (void*)offsetof(Particle, x)        // Offset to the position data in the struct
-    );
-
-    // Set a fixed color via uniform (e.g., white)
-    // Assuming your shader uses a uniform color
-    // If not, you can modify the fragment shader to use a fixed color
-    // Example: glUniform3f(colorUniformLoc, 1.0f, 1.0f, 1.0f);
-
-    // Draw particles as points
-    glDrawArrays(GL_POINTS, 0, NUM_PARTICLES);
-
-    // Disable vertex attribute
-    glDisableVertexAttribArray(positionAttribLoc);
-    glBindBuffer(GL_ARRAY_BUFFER, 0);
-
-    glutSwapBuffers();
-
-    glDepthMask(GL_TRUE);
-    glDisable(GL_BLEND);
-}
-
-void initParticles() {
-     // Calculate particles per axis
-    int particlesPerAxis = static_cast<int>(std::ceil(std::pow(NUM_PARTICLES, 1.0 / 3.0)));
-    int totalGridParticles = particlesPerAxis * particlesPerAxis * particlesPerAxis;
-
-    // Initialize particle data in a grid
-    int index = 0;
-    for (int i = 0; i < particlesPerAxis && index < NUM_PARTICLES; ++i) {
-        for (int j = 0; j < particlesPerAxis && index < NUM_PARTICLES; ++j) {
-            for (int k = 0; k < particlesPerAxis && index < NUM_PARTICLES; ++k) {
-                // Center particles around the origin
-                h_particles[index].x = (i - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING;
-                h_particles[index].y = (j - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + 10;
-                h_particles[index].z = (k - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING;
-
-                // Initialize other particle properties
-                h_particles[index].mass = 1.0f;
-                h_particles[index].vx = 0.0f;
-                h_particles[index].vy = 0.0f;
-                h_particles[index].vz = 0.0f;
-                h_particles[index].ax = 0.0f;
-                h_particles[index].ay = 0.0f;
-                h_particles[index].az = 0.0f;
-
-                ++index;
-                if (index >= NUM_PARTICLES) break;
-            }
-            if (index >= NUM_PARTICLES) break;
-        }
-        if (index >= NUM_PARTICLES) break;
-    }
-    /*
-    int index = 0;
-    for (int i = 0; i < GRID_SIZE; ++i) {
-        for (int j = 0; j < GRID_SIZE; ++j) {
-            h_particles[index].x = i * GRID_SPACING - (GRID_SIZE * GRID_SPACING) / 2.0f;  // X position
-            h_particles[index].y = j * GRID_SPACING - (GRID_SIZE * GRID_SPACING) / 2.0f;  // Y position
-            if (i % 2 == 0){
-                h_particles[index].z = i*5.0f;                                                  // Z position
-            } else {
-                h_particles[index].z = 100.0f;                                                  // Z position
-            }
-            h_particles[index].mass = 1.0f;                                               // Mass
-            h_particles[index].vx = 0.0f;                                                 // Initial velocity
-            h_particles[index].vy = 0.0f;
-            h_particles[index].vz = 0.0f;
-            h_particles[index].ax = 0.0f;                                                 // Initial acceleration
-            h_particles[index].ay = 0.0f;
-            h_particles[index].az = 0.0f;
-            ++index;
-        }
-    }
-    */
-    
-}
-
-// Initialize OpenGL and CUDA
-void initialize()
-{
-    // Init cuda constants
-    initializeConstants();
-    
-    // Set glewExperimental to enable modern OpenGL features
-    //glewExperimental = GL_TRUE;
-
-    // Initialize GLEW
-    /*
-    GLenum err = glewInit();
-    if (err != GLEW_OK) {
-        fprintf(stderr, "Error initializing GLEW: %s\n", glewGetErrorString(err));
-        exit(1);
-    }
-    */
-
-    // Load and compile shaders
-    shaderProgram = loadShaders("../shaders/particle.vert", "../shaders/particle.frag");
-    glUseProgram(shaderProgram);
-
-    // Query the attribute location for position
-    positionAttribLoc = glGetAttribLocation(shaderProgram, "aPosition");
-    if (positionAttribLoc == -1) {
-        fprintf(stderr, "Could not find attribute 'aPosition'\n");
-    } else {
-        printf("Attribute 'aPosition' location: %d\n", positionAttribLoc);
-    }
-
-    // Replace vec3 with glm::vec3
-    vec3 eyePosition(0.0f, 65.0f, 100.0f);
-    vec3 targetPosition(0.0f, 0.0f, 0.0f);
-    vec3 upVector(0.0f, 1.0f, 0.0f);
-
-    // Create View Matrix using glm::lookAt
-    mat4 viewMatrix = lookAt(eyePosition, targetPosition, upVector);
-
-    // Create Projection Matrix using glm::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;
-    mat4 projectionMatrix = perspective(fov * 3.1415/180, aspectRatio, zNear, zFar);
-
-    // Model Matrix (Identity matrix)
-    mat4 modelMatrix = mat4(1.0f);
-
-    // Query uniform locations
-    modelMatrixLoc = glGetUniformLocation(shaderProgram, "modelMatrix");
-    viewMatrixLoc = glGetUniformLocation(shaderProgram, "viewMatrix");
-    projectionMatrixLoc = glGetUniformLocation(shaderProgram, "projectionMatrix");
-    pointRadiusLoc = glGetUniformLocation(shaderProgram, "pointRadius");
-    pointScaleLoc = glGetUniformLocation(shaderProgram, "pointScale");
-
-
-    // Enable programmable point size
-
-    //printMatrix(modelMatrix, "Model Matrix");
-    //printMatrix(viewMatrix, "View Matrix");
-    //printMatrix(projectionMatrix, "Projection Matrix");
-
-    glUniformMatrix4fv(modelMatrixLoc, 1, GL_FALSE, modelMatrix.m);
-    glUniformMatrix4fv(viewMatrixLoc, 1, GL_FALSE, viewMatrix.m);
-    glUniformMatrix4fv(projectionMatrixLoc, 1, GL_FALSE, projectionMatrix.m);
-    
-    // Set point size uniform
-    if (pointRadiusLoc != -1) {
-        glUniform1f(pointRadiusLoc, pointRadius); // Adjust as needed
-    } else {
-        fprintf(stderr, "Failed to find 'pointRadius' uniform location!\n");
-    }
-
-    if (pointScaleLoc != -1) {
-        glUniform1f(pointScaleLoc, pointScale); // Adjust as needed
-    } else {
-        fprintf(stderr, "Failed to find 'pointScale' uniform location!\n");
-    }
-
-    
-
-    initParticles();
-    // Initialize VBO for particle data
-    glGenBuffers(1, &vbo);
-    glBindBuffer(GL_ARRAY_BUFFER, vbo);
-    // Upload initial particle data to the VBO
-    glBufferData(GL_ARRAY_BUFFER, h_particles.size() * sizeof(Particle), h_particles.data(), GL_DYNAMIC_DRAW);
-
-    // Map VBO for reading (optional, for debugging)
-    glBindBuffer(GL_ARRAY_BUFFER, vbo);
-    /*
-    Particle* vboData = (Particle*)glMapBuffer(GL_ARRAY_BUFFER, GL_READ_ONLY);
-    if (vboData) {
-        printf("First particle position: (%f, %f, %f)\n", vboData[0].x, vboData[0].y, vboData[0].z);
-        glUnmapBuffer(GL_ARRAY_BUFFER);
-    } else {
-        printf("Failed to map VBO for reading.\n");
-    }
-    */
-   /*
-    for(int i = 0; i < 5; ++i) { // Check first 5 particles
-        printf("Particle %d: Position after upload: (%f, %f, %f)\n", i, h_particles[i].x, h_particles[i].y, h_particles[i].z);
-        //glUnmapBuffer(GL_ARRAY_BUFFER);
-    }
-   */
-    glBindBuffer(GL_ARRAY_BUFFER, 0);
-
-    // Register VBO with CUDA
-    registerVBOWithCUDA(vbo);
-
-    // Set point size for rendering particles (optional, can be controlled in shader)
-    //glPointSize(PARTICLE_SIZE);
-    glEnable(GL_PROGRAM_POINT_SIZE);
-
-    // Enable point sprites (for compatibility profile, optional in core profile)
-    glEnable(GL_POINT_SPRITE);
-
-    // Enable depth testing if needed (for 3D)
-    glEnable(GL_DEPTH_TEST);
-
-
-}
-
-// Cleanup resources
-void cleanup()
-{
-    unregisterVBOFromCUDA();
-    glDeleteBuffers(1, &vbo);
-}
-
-int main(int argc, char** argv)
-{
-    // Initialize GLUT
-    glutInit(&argc, argv);
-    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB /*| GLUT_DEPTH*/);
-    glutInitWindowSize(800, 600);
-    glutCreateWindow("CUDA-OpenGL Particle Simulation");
-
-    // Initialize OpenGL and CUDA
-    initialize();
-
-    // Register display callback
-    glutDisplayFunc(display);
-    glutIdleFunc(display);
-
-    // Cleanup on exit
-    atexit(cleanup);
-
-    // Start GLUT main loop
-    glutMainLoop();
-
-    return 0;
-}
\ No newline at end of file
diff --git a/src/main3.cpp b/src/main3.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b4ed3f8025dd745ec460e9e4a30ce2dd27f5835f
--- /dev/null
+++ b/src/main3.cpp
@@ -0,0 +1,416 @@
+// main.cpp: This example uses FreeGLUT for OpenGL and CUDA interoperation.
+// It creates a simple particle simulation and uses CUDA to update particle positions with gravity.
+// main.cpp
+#include <vector>
+#include <cstdio>
+#include <ctime>
+#include <cmath> // For std::pow
+#include <chrono>
+#include <cstdlib> // For exit()
+
+#ifdef _WIN32
+    #include <GL/glew.h>  // Include GLEW before other OpenGL headers
+#endif
+
+
+//#include <GL/freeglut.h>
+#include "MicroGlut.h"
+
+#include "GL_utilities.h"
+#include "kernels_cuda.h" // Header for CUDA function declarations
+
+//#include <glm/glm.hpp>
+//#include <glm/gtc/matrix_transform.hpp> // For transformation functions
+//#include <glm/gtc/type_ptr.hpp>
+
+#define MAIN
+#include "VectorUtils4.h"
+
+// Number of particles
+const int NUM_PARTICLES = 64000; // Should match a cubic grid for simplicity
+const float GRID_SPACING = 0.25f;  // Spacing between particles in the grid
+
+const float pointRadius = 6.0f;  // Point size in world space
+const float pointScale = 16.0f;
+
+auto previousTime = std::chrono::high_resolution_clock::now();
+float deltaTime = 0.016f;
+
+// 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;
+GLuint positionAttribLoc;
+// Vertex Buffer Object (VBO) for particle positions
+GLuint vbo;
+
+// Uniform locations
+GLuint modelMatrixLoc;
+GLuint viewMatrixLoc;
+GLuint projectionMatrixLoc;
+GLuint pointRadiusLoc;
+GLuint pointScaleLoc;
+
+// 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]);
+    }
+}
+ */
+
+void printFrameStats() {
+    // FPS and deltaTime accumulation
+    static float timeAccumulator = 0.0f;
+    static int frameCount = 0;
+
+    timeAccumulator += deltaTime;
+    frameCount++;
+
+    // Print FPS and average deltaTime every second
+    if (timeAccumulator >= 1.0f)
+    {
+        float averageDeltaTime = timeAccumulator / frameCount;
+        float fps = frameCount / timeAccumulator;
+
+        printf("FPS: %.2f, Average DeltaTime: %.1f ms\n", fps, averageDeltaTime*1000);
+
+        // Reset accumulators
+        timeAccumulator = 0.0f;
+        frameCount = 0;
+    }
+}
+
+// Display callback for OpenGL
+void display()
+{
+
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+    glEnable(GL_BLEND);
+    glDepthMask(GL_FALSE);
+   
+    auto currentTime = std::chrono::high_resolution_clock::now();
+    deltaTime = std::chrono::duration<float>(currentTime - previousTime).count();
+    previousTime = currentTime;
+
+     // Clamp deltaTime to prevent instability
+    if (deltaTime > 0.1f) // Max 100 ms
+        deltaTime = 0.1f;
+    if (deltaTime < 0.0001f) // Min 0.1 ms
+        deltaTime = 0.0001f;
+
+    // Update particle positions using CUDA
+    updateParticles(deltaTime);
+
+    // Use the shader program
+    glUseProgram(shaderProgram);
+
+
+    // Todo testing modifying stuff
+    mat4 projectionMatrix = perspective(fov, aspectRatio, zNear, zFar);
+    glUniformMatrix4fv(projectionMatrixLoc, 1, GL_TRUE, projectionMatrix.m);
+
+    mat4 viewMatrix = lookAt(campos, VectorAdd(campos, forward), upVector);
+    glUniformMatrix4fv(viewMatrixLoc, 1, GL_TRUE, viewMatrix.m);
+
+
+
+    // Bind VBO and set up vertex attribute for position
+    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 particles as points
+    glDrawArrays(GL_POINTS, 0, NUM_PARTICLES);
+
+    // Disable vertex attribute
+    glDisableVertexAttribArray(positionAttribLoc);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+    glutSwapBuffers();
+
+    glDepthMask(GL_TRUE);
+    glDisable(GL_BLEND);
+
+    printFrameStats();
+
+}
+
+void initParticles()
+{
+    // Seed the random number generator
+    srand(static_cast<unsigned int>(time(0)));
+
+    // Calculate particles per axis
+    int particlesPerAxis = static_cast<int>(ceil(pow(NUM_PARTICLES, 1.0f / 3.0f)));
+
+    int index = 0;
+    for (int i = 0; i < particlesPerAxis && index < NUM_PARTICLES; ++i) {
+        for (int j = 0; j < particlesPerAxis && index < NUM_PARTICLES; ++j) {
+            for (int k = 0; k < particlesPerAxis && index < NUM_PARTICLES; ++k) {
+                // Center particles around the origin
+                float x = (i - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING;
+                float y = (j - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + 2.0f;
+                float z = (k - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING;
+
+                // Add random perturbation within [-0.1*dx, 0.1*dx]
+                float perturbation = 0.1f * GRID_SPACING;
+                float rand_x = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
+                float rand_y = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
+                float rand_z = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
+
+                h_positions[index] = make_float3(
+                    x + rand_x,
+                    y + rand_y,
+                    z + rand_z
+                );
+
+                // Initialize velocities
+                h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f);
+
+                ++index;
+            }
+        }
+    }
+
+    // If NUM_PARTICLES is not a perfect cube, initialize remaining particles
+    for (; index < NUM_PARTICLES; ++index) {
+        h_positions[index] = make_float3(0.0f, 10.0f, 0.0f); // Default position
+        h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f); // Default velocity
+    }
+}
+
+// Initialize OpenGL and CUDA
+void initialize()
+{
+    // Initialize constants in CUDA
+    initializeConstants();
+
+    // Set glewExperimental to enable modern OpenGL features
+    /*
+    glewExperimental = GL_TRUE;
+
+    // Initialize GLEW
+    GLenum err = glewInit();
+    if (err != GLEW_OK) {
+        fprintf(stderr, "Error initializing GLEW: %s\n", glewGetErrorString(err));
+        exit(1);
+    }
+    */
+
+    // Load and compile shaders
+    //shaderProgram = loadShaders("../../shaders/particle.vert", "../../shaders/particle.frag");
+    shaderProgram = loadShaders("../shaders/particle.vert", "../shaders/particle.frag");
+    glUseProgram(shaderProgram);
+
+    // Query the attribute location for position
+    positionAttribLoc = glGetAttribLocation(shaderProgram, "aPosition");
+    if (positionAttribLoc == -1) {
+        fprintf(stderr, "Could not find attribute 'aPosition'\n");
+    } else {
+        printf("Attribute 'aPosition' location: %d\n", positionAttribLoc);
+    }
+
+    // Camera setup
+ 
+    mat4 viewMatrix = lookAt(campos,  VectorAdd(campos, forward), upVector);
+
+    // Create Projection Matrix using glm::perspective
+    mat4 projectionMatrix = perspective(fov, aspectRatio, zNear, zFar);
+
+    // Model Matrix (Identity matrix)
+    mat4 modelMatrix = IdentityMatrix();
+    //printMat4(modelMatrix);
+
+    // Query uniform locations
+    modelMatrixLoc = glGetUniformLocation(shaderProgram, "modelMatrix");
+    viewMatrixLoc = glGetUniformLocation(shaderProgram, "viewMatrix");
+    projectionMatrixLoc = glGetUniformLocation(shaderProgram, "projectionMatrix");
+    pointRadiusLoc = glGetUniformLocation(shaderProgram, "pointRadius");
+    pointScaleLoc = glGetUniformLocation(shaderProgram, "pointScale");
+
+
+    // Set uniform matrices
+    glUniformMatrix4fv(modelMatrixLoc, 1, GL_TRUE, modelMatrix.m);
+    glUniformMatrix4fv(viewMatrixLoc, 1, GL_TRUE, viewMatrix.m);
+    glUniformMatrix4fv(projectionMatrixLoc, 1, GL_TRUE, projectionMatrix.m);
+
+    // Set point size uniforms
+    if (pointRadiusLoc != -1) {
+        glUniform1f(pointRadiusLoc, pointRadius);
+    } else {
+        fprintf(stderr, "Failed to find 'pointRadius' uniform location!\n");
+    }
+
+    if (pointScaleLoc != -1) {
+        glUniform1f(pointScaleLoc, pointScale);
+    } else {
+        fprintf(stderr, "Failed to find 'pointScale' uniform location!\n");
+    }
+
+    // Initialize particle data
+    initParticles();
+
+    // 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);
+
+    // Register VBO with CUDA
+    registerVBOWithCUDA(vbo);
+
+    // Initialize device arrays
+    initializeDeviceArrays(
+        h_positions.data(),
+        h_velocities.data(),
+        NUM_PARTICLES
+    );
+
+    // Enable point size for rendering particles
+    glEnable(GL_PROGRAM_POINT_SIZE);
+    glEnable(GL_POINT_SPRITE);
+
+    // Enable depth testing
+    glEnable(GL_DEPTH_TEST);
+}
+
+// Cleanup resources
+void cleanup()
+{
+    unregisterVBOFromCUDA();
+    cleanupDeviceArrays();
+    glDeleteBuffers(1, &vbo);
+}
+
+
+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 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
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB /*| GLUT_DEPTH*/);
+    glutInitWindowSize(800, 600);
+    glutCreateWindow("CUDA-OpenGL Particle Simulation");
+    glClearColor(0.9f, 0.9f, 0.9f, 1.0f);
+
+    // Initialize OpenGL and CUDA
+    initialize();
+
+    // Register display callback
+    glutDisplayFunc(display);
+    glutIdleFunc(display);
+    glutKeyboardFunc(Key);
+    glutReshapeFunc(Reshape);
+
+    
+
+    // Cleanup on exit
+    atexit(cleanup);
+
+    // Start GLUT main loop
+    glutMainLoop();
+
+    return 0;
+}