diff --git a/README.md b/README.md
index 45045d6a0150404b2ce52c677649ef1017fb8cc5..6b4104a6959bc4bc54298059a91cab6e5b58aa9e 100644
--- a/README.md
+++ b/README.md
@@ -38,7 +38,6 @@ These values are adjustable within the code and can be fine-tuned based on desir
    sudo apt install -y build-essential cmake libglew-dev freeglut3-dev libglm-dev
    ```
    Install the **CUDA Toolkit** from NVIDIA’s official site [here](https://developer.nvidia.com/cuda-toolkit), which includes both CUDA and Thrust libraries.
-    ```
 
 2. **Build the Project with CMake**
 
diff --git a/src/constants.h b/src/constants.h
index 4665a6f462825ab6bc2ad39ecb9cbe5437b0c774..5a5d500f5927e1d018798d4f320c2ecd2ecee56f 100644
--- a/src/constants.h
+++ b/src/constants.h
@@ -2,17 +2,17 @@
 #ifndef CONSTANTS_H
 #define CONSTANTS_H
 
-const float h_BOUND_X_MIN = -5.0f;
-const float h_BOUND_X_MAX = 5.0f;
-const float h_BOUND_Y_MIN = -12.0f;
+const float h_BOUND_X_MIN = -6.0f;
+const float h_BOUND_X_MAX = 12.0f;
+const float h_BOUND_Y_MIN = 0.0f;
 const float h_BOUND_Y_MAX = 12.0f;
-const float h_BOUND_Z_MIN = -5.0f;
-const float h_BOUND_Z_MAX = 5.0f;
+const float h_BOUND_Z_MIN = -4.0f;
+const float h_BOUND_Z_MAX = 4.0f;
 
-const float h_GAS_CONSTANT = 0.02f;
-const float h_VISCOSITY = 3.2f;
+const float h_GAS_CONSTANT = 300.0f;
+const float h_VISCOSITY = 0.8f;
 
 const float h_PI = 3.14159265358979323846f;
-const float h_H = 0.15f; // Smoothing radius
+const float h_H = 0.4f; // Smoothing radius
 
 #endif // CONSTANTS_H
diff --git a/src/kernels_cuda.cu b/src/kernels_cuda.cu
index 2dd90f96df3979012d60d05e95a59d2059bc15ed..3a1717753fdd45399331e7b297be8aac7db31a94 100644
--- a/src/kernels_cuda.cu
+++ b/src/kernels_cuda.cu
@@ -23,8 +23,8 @@ __constant__ float POLY6_CONST;
 __constant__ float SPIKY_GRAD_CONST;
 __constant__ float VISC_LAP_CONST;
 
-__constant__ float PARTICLE_MASS = 0.0156f;
-__constant__ float REST_DENSITY = 20.0f;
+__constant__ float PARTICLE_MASS = 8.0f;
+__constant__ float REST_DENSITY = 1000.0f;
 __constant__ float GAS_CONSTANT;
 __constant__ float VISCOSITY;
 
@@ -135,6 +135,7 @@ __global__ void reorderData(
 }
 
 // Find cell start and end indices
+// can bli reeplaceed med atomicAdd för varje ceell index??
 __global__ void findCellStartEnd(
     unsigned int* cellIndices,
     unsigned int* cellStart,
@@ -219,9 +220,9 @@ __global__ void computeDensityPressure(
         }
 
         densities[idx] = density;
-        //pressures[idx] = GAS_CONSTANT * fmaxf((density - REST_DENSITY), 0.0);
-        //pressures[idx] = GAS_CONSTANT * fmaxf((density - REST_DENSITY), 0.0);
-        pressures[idx] = GAS_CONSTANT * (density - REST_DENSITY);
+        pressures[idx] = GAS_CONSTANT * fmaxf((density - REST_DENSITY), 0.0);
+        //pressures[idx] = GAS_CONSTANT * min((density - REST_DENSITY), 0.0);
+        //pressures[idx] = GAS_CONSTANT * (density - REST_DENSITY);
     }
 }
 
@@ -249,7 +250,7 @@ __global__ void computeForces(
 
         float3 fPressure = make_float3(0.0f, 0.0f, 0.0f);
         float3 fViscosity = make_float3(0.0f, 0.0f, 0.0f);
-        float3 fGravity = make_float3(0.0f, -1.0f * density_i, 0.0f);
+        float3 fGravity = make_float3(0.0f, -16.0f * density_i, 0.0f);
 
         unsigned int cellIndex = cellIndices[idx];
 
@@ -280,7 +281,6 @@ __global__ void computeForces(
 
                             float3 pos_j = positions[j];
                             float3 r = (pos_i - pos_j);
-                            //float3 r = (pos_j - pos_i); not sure if it should be i -j or j-i
                             
                             float distSq = r.x * r.x + r.y * r.y + r.z * r.z;
 
@@ -300,7 +300,8 @@ __global__ void computeForces(
                                 fPressure += pressureTerm*(r/dist); // r/dist gives direction
                                 */
 
-                                // Symmetric pressure
+                                // Symmetric pressure -- can cahnge to inverse dennsity in density/pressure calc to avoid devide by density here
+                                // Add particle_mass to pressure term according to equations !!! Added down below in aceleration calcualtion
                                 float pressureTerm = - (pressure_i / (density_i * density_i) + pressure_j / (density_j * density_j)) * spikyGrad;
                                 float3 pressureForce = pressureTerm * (r / dist); // Directional force
                                 fPressure += pressureForce;
@@ -317,7 +318,11 @@ __global__ void computeForces(
         }
 
         // Total acceleration
-        accelerations[idx] = (fPressure + fViscosity + fGravity) / density_i;
+        // According to paper v_i(t+dt) = v_(i+1) + (dt/m_i)*Fp_i
+        // While external and viscocisty force is calculated before But mass and nnot deeensity is used to deevidee force
+        // Deviding by mass, looks a lot better
+        accelerations[idx] = (fPressure * PARTICLE_MASS + fViscosity + fGravity) / PARTICLE_MASS;
+        //accelerations[idx] = (fPressure * PARTICLE_MASS + fViscosity + fGravity) / density_i;
     }
 }
 
@@ -359,7 +364,7 @@ __global__ void enforceBoundary(
             velocities[idx].x *= -damping;
         }
         if (positions[idx].x > BOUND_X_MAX) {
-            positions[idx].x = BOUND_X_MAX;
+            positions[idx].x = BOUND_X_MAX; // Should account for particle radius!!
             velocities[idx].x *= -damping;
         }
 
@@ -481,7 +486,7 @@ void updateParticles(float deltaTime)
     cudaDeviceSynchronize();
 
     // Integrate
-    float timeStep = 0.017f;
+    float timeStep = 0.008f;
     integrate<<<blocksPerGrid, threadsPerBlock>>>(
         d_positions,
         d_velocities,
@@ -492,7 +497,7 @@ void updateParticles(float deltaTime)
     cudaDeviceSynchronize();
 
     // Enforce boundary conditions
-    float damping = 0.95f;
+    float damping = 0.7f;
     enforceBoundary<<<blocksPerGrid, threadsPerBlock>>>(
         d_positions,
         d_velocities,
diff --git a/src/main.cpp b/src/main.cpp
index d72bdf539ddbd49a5f7383d7564de0d070a7170f..601bbe11b0d682608496d11c89a37b3cc2f131fd 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -23,11 +23,11 @@
 #include "GL_utilities.h"
 
 // Number of particles
-const int NUM_PARTICLES = 64000; // Should match a cubic grid for simplicity
-const float GRID_SPACING = 0.125f;  // Spacing between particles in the grid
+const int NUM_PARTICLES = 8000; // Should match a cubic grid for simplicity
+const float GRID_SPACING = 0.5f;  // Spacing between particles in the grid
 
-const float pointRadius = 2.5f;  // Point size in world space
-const float pointScale = 10.0f;
+const float pointRadius = 8.0f;  // Point size in world space
+const float pointScale = 16.0f;
 
 auto previousTime = std::chrono::high_resolution_clock::now();
 float deltaTime = 0.016f;
@@ -149,7 +149,7 @@ void initParticles()
             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 + 10.0f;
+                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]
@@ -208,8 +208,8 @@ void initialize()
     }
 
     // Camera setup
-    glm::vec3 eyePosition(2.0f, 25.0f, 30.0f);
-    glm::vec3 targetPosition(0.0f, 0.0f, 0.0f);
+    glm::vec3 eyePosition(5.0f, 15.0f, 20.0f);
+    glm::vec3 targetPosition(0.0f, 4.0f, 0.0f);
     glm::vec3 upVector(0.0f, 1.0f, 0.0f);
 
     // Create View Matrix using glm::lookAt