diff --git a/README.md b/README.md
index ee39093..005af4b 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,83 @@
 **University of Pennsylvania, CIS 5650: GPU Programming and Architecture,
 Project 1 - Flocking**
 
-* (TODO) YOUR NAME HERE
-  * (TODO) [LinkedIn](), [personal website](), [twitter](), etc.
-* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab)
+* CARLOS LOPEZ GARCES
+  * [LinkedIn](https://www.linkedin.com/in/clopezgarces/)
+  * [Personal website](https://carlos-lopez-garces.github.io/)
+* Tested on: Windows 11, 13th Gen Intel(R) Core(TM) i9-13900HX @ 2.20 GHz, RAM 32GB, NVIDIA GeForce RTX 4060, personal laptop.
 
-### (TODO: Your README)
+# Boid flocking simulation
 
-Include screenshots, analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+Simulation of flocking behavior of boids (bird-oid particles), accelerated using CUDA.
+
+At each time step of the simulation, the position and velocity of each boid is updated according to 3 rules: (1) *cohesion*, according to which a boid will head towards the center of mass of nearby boids; (2) *separation*, according to which a boid will keep some distance from nearby boids; and (3) *alignment*, according to which a boid will head in the same direction and with the same velocity as nearby boids. 
+
+## 3 implementations
+
+The "nearby boids" are determined based on a fixed radius, which may be different for each rule.
+
+The search for nearby boids has 3 implementations, a naive one, one that partitions the space using a uniform grid, and one that optimizes memory access using that same uniform grid. Details follow.
+
+### Naive
+
+To search for nearby boids, each boid compares distance against every other boid against each of the 3 rules' thresholds.
+
+GIF shows the default 5000 boids, using default block size 128.
+
+![](images/naive_5000_boids.gif)
+
+### Uniform grid
+
+To avoid comparing distance with every other boid, space is partitioned into a uniform grid and each boid only compares distance against boids in up to 2 cells along each dimension (x, y, z): since each cell is twice the maximum neighborhood radius, each boid only needs to compare distance against boids in its own cell and up to 1 extra adjacent cell in each dimension (the neighborhood can only every overlap uo to 2 adjacent cells in each dimension).
+
+GIF shows the default 5000 boids, using default block size 128.
+
+![](images/uniform_grid_5000_boids.gif)
+
+### Coherent uniform grid
+
+To avoid using an intermediary index buffer to find position and velocity of boids, these arrays are rearranged so that they can be found at the same index as the boid index / cell index array. This allows for coalesced memory accesses.
+
+GIF shows the default 5000 boids, using default block size 128.
+
+![](images/coherent_grid_5000_boids.gif)
+
+
+## Performance graphs
+
+The uniform grid optimization outperforms the naive one because each boid is compared against fewer other boids. The coherent grid on top of the uniform grid optimization also boosts framerate because coalesced memory access allows for efficient use of the cache.
+
+![](images/FPS_as_Boid_Count_Increases_with_Visualization.png)
+
+As a result of turning visualization off, FPS naturally goes up.
+
+![](images/FPS_as_Boid_Count_Increases_without_Visualization.png)
+
+Block sizes of 32, 64, 128, 256, 512, and 1024 are compared for 5000 boids.
+
+![](images/FPS_as_Block_Size_Increases_without_Visualization.png)
+
+
+## Questions
+
+***For each implementation, how does changing the number of boids affect performance? Why do you think this is?***
+
+In the naive implementation, fps sharply decreases because the number of neighborhood checks ***per boid*** increases linearly with the number of total boids.
+
+In the two uniform grid implementations, fps decreases less sharply than in the naive implementation; this is because, while the number of boids to check in adjacent cells varies, on average there are fewer boids there than the total.
+
+With the coherent grid optimization, fps decreases even less sharply. The memory access pattern improves and, as a result, the cache is used more efficiently, decreasing time per frame.
+
+***For each implementation, how does changing the block count and block size affect performance? Why do you think this is?***
+
+I compared block sizes of 32, 64, 128, 256, 512, and 1024 for 5000 boids. In the naive implementation, increasing block size led to decreased framerate. In the other two implementations, average framerate increases as block size increases from 32 to 64, and then appears to decrease as it increases to 128, 256, and 512, to finally increase when block size is 1024. A fluctuation is observed, but it's not very drastic.
+
+There may be many reasons behind these observations. One explanation for the increased fps from 32 to 64, is that fewer blocks (with more threads) are used and these threads are able to leverage cache loads from other threads in the same block. A reason for seing decreased fps afterwards is increased register pressure that prevents parallel scheduling.
+
+***For the coherent uniform grid: did you experience any performance improvements with the more coherent uniform grid? Was this the outcome you expected? Why or why not?***
+
+Yes. Framerate decreases less sharply with this optimization. This is likely the result of using the cache more efficiently, since we are now accessing contiguous array elements more frequently (better memory access pattern that leads to reduced memory latency), as opposed to scattered elements.
+
+***Did changing cell width and checking 27 vs 8 neighboring cells affect performance? Why or why not? Be careful: it is insufficient (and possibly incorrect) to say that 27-cell is slower simply because there are more cells to check!***
+
+Decreasing cell width appears to improve performance because the set of nearby adjacent cells to check surrounds nearby boids more tightly; this results in fewer boids being checked (fewer computations per boid and fewer memory accesses). When cells are larger, there may be boids in the outermost checked cell that may not be within the neighborhood radius (and still they are checked).
\ No newline at end of file
diff --git a/images/FPS_as_Block_Size_Increases_without_Visualization.png b/images/FPS_as_Block_Size_Increases_without_Visualization.png
new file mode 100644
index 0000000..17260c1
Binary files /dev/null and b/images/FPS_as_Block_Size_Increases_without_Visualization.png differ
diff --git a/images/FPS_as_Boid_Count_Increases_with_Visualization.png b/images/FPS_as_Boid_Count_Increases_with_Visualization.png
new file mode 100644
index 0000000..46656c6
Binary files /dev/null and b/images/FPS_as_Boid_Count_Increases_with_Visualization.png differ
diff --git a/images/FPS_as_Boid_Count_Increases_without_Visualization.png b/images/FPS_as_Boid_Count_Increases_without_Visualization.png
new file mode 100644
index 0000000..4d3bb34
Binary files /dev/null and b/images/FPS_as_Boid_Count_Increases_without_Visualization.png differ
diff --git a/images/coherent_grid_5000_boids.gif b/images/coherent_grid_5000_boids.gif
new file mode 100644
index 0000000..7c26ee6
Binary files /dev/null and b/images/coherent_grid_5000_boids.gif differ
diff --git a/images/naive_5000_boids.gif b/images/naive_5000_boids.gif
new file mode 100644
index 0000000..19e1ecf
Binary files /dev/null and b/images/naive_5000_boids.gif differ
diff --git a/images/uniform_grid_5000_boids.gif b/images/uniform_grid_5000_boids.gif
new file mode 100644
index 0000000..d342732
Binary files /dev/null and b/images/uniform_grid_5000_boids.gif differ
diff --git a/src/kernel.cu b/src/kernel.cu
index 74dffcb..85f2ca7 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <cuda.h>
 #include <cmath>
+#include <utility>
 #include <glm/glm.hpp>
 #include "utilityCore.hpp"
 #include "kernel.h"
@@ -49,7 +50,7 @@ void checkCUDAError(const char *msg, int line = -1) {
 #define rule2Scale 0.1f
 #define rule3Scale 0.1f
 
-#define maxSpeed 1.0f
+#define maxSpeed 2.5f
 
 /*! Size of the starting area in simulation space. */
 #define scene_scale 100.0f
@@ -66,14 +67,22 @@ dim3 threadsPerBlock(blockSize);
 // Consider why you would need two velocity buffers in a simulation where each
 // boid cares about its neighbors' velocities.
 // These are called ping-pong buffers.
+// Since each boid thread will be updating its velocity and reading the velocities
+// of other boids concurrently with the other boid threads, we need all threads to
+// avoid writing to one of the buffers (where neighbor velocities are to be read)
+// in a given simulation step and only write to one of them and read from the other;
+// in the next simulation step, the roles of these buffers will swap.
 glm::vec3 *dev_pos;
+glm::vec3 *dev_posReshuffled;
 glm::vec3 *dev_vel1;
+glm::vec3 *dev_vel1Reshuffled;
 glm::vec3 *dev_vel2;
 
 // LOOK-2.1 - these are NOT allocated for you. You'll have to set up the thrust
 // pointers on your own too.
 
 // For efficient sorting and the uniform grid. These should always be parallel.
+// These 2 will be sorted according to uniform grid cell index.
 int *dev_particleArrayIndices; // What index in dev_pos and dev_velX represents this particle?
 int *dev_particleGridIndices; // What grid cell is this particle in?
 // needed for use with thrust
@@ -169,6 +178,29 @@ void Boids::initSimulation(int N) {
   gridMinimum.z -= halfGridWidth;
 
   // TODO-2.1 TODO-2.3 - Allocate additional buffers here.
+
+  cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int));
+  checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!");
+  // Wrap the device pointer to the array with a thrust device pointer.
+  dev_thrust_particleArrayIndices = thrust::device_pointer_cast(dev_particleArrayIndices);
+
+  cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int));
+  checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!");
+  // Wrap the device pointer to the array with a thrust device pointer.
+  dev_thrust_particleGridIndices = thrust::device_pointer_cast(dev_particleGridIndices);
+
+  cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int));
+  checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!");
+
+  cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int));
+  checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!");
+
+  cudaMalloc((void**)&dev_posReshuffled, N * sizeof(glm::vec3));
+  checkCUDAErrorWithLine("cudaMalloc dev_posReshuffled failed!");
+
+  cudaMalloc((void**)&dev_vel1Reshuffled, N * sizeof(glm::vec3));
+  checkCUDAErrorWithLine("cudaMalloc dev_vel1Reshuffled failed!");
+
   cudaDeviceSynchronize();
 }
 
@@ -223,17 +255,181 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities)
 * stepSimulation *
 ******************/
 
+__device__ glm::vec3 computeVelocityChangeDueToRule1Cohesion(int N, int iSelf, const glm::vec3 *pos) {
+  glm::vec3 perceivedCenter(0.0f, 0.0f, 0.0f);
+  glm::vec3 selfPos = pos[iSelf];
+  int nearbyBoidsCount = 0;
+
+  for (int b = 0; b < N; ++b) {
+    glm::vec3 bPos = pos[b];  
+    if (b != iSelf && distance(bPos, selfPos) < rule1Distance) {
+        perceivedCenter += bPos;
+        nearbyBoidsCount++;
+    }
+  }
+
+  if (nearbyBoidsCount > 0) {
+    // Center of mass of nearby boids.
+    perceivedCenter /= nearbyBoidsCount;
+
+    return (perceivedCenter - selfPos) * rule1Scale;
+  }
+
+  // No nearby boids, no change in velocity.
+  return glm::vec3(0.0f, 0.0f, 0.0f);
+}
+
+__device__ glm::vec3 computeVelocityChangeDueToRule2Separation(int N, int iSelf, const glm::vec3 *pos) {
+  glm::vec3 offsetToKeepDistance(0.0f, 0.0f, 0.0f);
+  glm::vec3 selfPos = pos[iSelf];
+
+  for (int b = 0; b < N; ++b) {
+    glm::vec3 bPos = pos[b];  
+    if (b != iSelf && distance(bPos, selfPos) < rule2Distance) {
+      offsetToKeepDistance -= bPos - selfPos;
+    }
+  }
+
+  return offsetToKeepDistance * rule2Scale;
+}
+
+__device__ glm::vec3 computeVelocityChangeDueToRule3Alignment(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) {
+  glm::vec3 perceivedVelocity(0.0f, 0.0f, 0.0f);
+  glm::vec3 selfPos = pos[iSelf];
+  int nearbyBoidsCount = 0;
+
+  for (int b = 0; b < N; ++b) {
+    if (b != iSelf && distance(pos[b], selfPos) < rule3Distance) {
+        perceivedVelocity += vel[b];
+        nearbyBoidsCount++;
+    }
+  }
+
+  if (nearbyBoidsCount > 0) {
+    // Average velocity.
+    perceivedVelocity /= nearbyBoidsCount;
+  }
+
+  // If there were no nearby boids, nearbyBoidsCount = 0, then perceivedVelocity = 0 too.
+  return perceivedVelocity * rule3Scale;
+}
+
+__device__ glm::vec3 computeVelocityChangeDueToRule1Cohesion(int iSelf, int startIndex, int endIndex, int *particleArrayIndices, const glm::vec3 *pos) {
+  glm::vec3 perceivedCenter(0.0f, 0.0f, 0.0f);
+  glm::vec3 selfPos = pos[iSelf];
+  int nearbyBoidsCount = 0;
+
+  for (int i = startIndex; i <= endIndex; ++i) {
+    int b = particleArrayIndices[i];
+    glm::vec3 bPos = pos[b];  
+    if (b != iSelf && distance(bPos, selfPos) < rule1Distance) {
+        perceivedCenter += bPos;
+        nearbyBoidsCount++;
+    }
+  }
+
+  if (nearbyBoidsCount > 0) {
+    // Center of mass of nearby boids.
+    perceivedCenter /= nearbyBoidsCount;
+
+    return (perceivedCenter - selfPos) * rule1Scale;
+  }
+
+  // No nearby boids, no change in velocity.
+  return glm::vec3(0.0f, 0.0f, 0.0f);
+}
+
+__device__ glm::vec3 computeVelocityChangeDueToRule2Separation(int iSelf, int startIndex, int endIndex, int *particleArrayIndices, const glm::vec3 *pos) {
+  glm::vec3 offsetToKeepDistance(0.0f, 0.0f, 0.0f);
+  glm::vec3 selfPos = pos[iSelf];
+
+  for (int i = startIndex; i <= endIndex; ++i) {
+    int b = particleArrayIndices[i];
+    glm::vec3 bPos = pos[b];  
+    if (b != iSelf && distance(bPos, selfPos) < rule2Distance) {
+      offsetToKeepDistance -= bPos - selfPos;
+    }
+  }
+
+  return offsetToKeepDistance * rule2Scale;
+}
+
+__device__ glm::vec3 computeVelocityChangeDueToRule3Alignment(int iSelf, int startIndex, int endIndex, int *particleArrayIndices, const glm::vec3 *pos, const glm::vec3 *vel) {
+  glm::vec3 perceivedVelocity(0.0f, 0.0f, 0.0f);
+  glm::vec3 selfPos = pos[iSelf];
+  int nearbyBoidsCount = 0;
+
+  for (int i = startIndex; i <= endIndex; ++i) {
+    int b = particleArrayIndices[i];
+    if (b != iSelf && distance(pos[b], selfPos) < rule3Distance) {
+        perceivedVelocity += vel[b];
+        nearbyBoidsCount++;
+    }
+  }
+
+  if (nearbyBoidsCount > 0) {
+    // Average velocity.
+    perceivedVelocity /= nearbyBoidsCount;
+  }
+
+  // If there were no nearby boids, nearbyBoidsCount = 0, then perceivedVelocity = 0 too.
+  return perceivedVelocity * rule3Scale;
+}
+
 /**
 * LOOK-1.2 You can use this as a helper for kernUpdateVelocityBruteForce.
 * __device__ code can be called from a __global__ context
 * Compute the new velocity on the body with index `iSelf` due to the `N` boids
 * in the `pos` and `vel` arrays.
 */
-__device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) {
-  // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves
-  // Rule 2: boids try to stay a distance d away from each other
-  // Rule 3: boids try to match the speed of surrounding boids
-  return glm::vec3(0.0f, 0.0f, 0.0f);
+__device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) { 
+  glm::vec3 selfPos = pos[iSelf];
+
+  glm::vec3 rule1PerceivedCenter(0.0f, 0.0f, 0.0f);
+  glm::vec3 rule2OffsetToKeepDistance(0.0f, 0.0f, 0.0f);
+  glm::vec3 rule3PerceivedVelocity(0.0f, 0.0f, 0.0f);
+  
+  int rule1NearbyBoidsCount = 0;
+  int rule3NearbyBoidsCount = 0;
+
+  for (int b = 0; b < N; ++b) {
+    if (b == iSelf) continue;
+
+    glm::vec3 bPos = pos[b];
+
+    float selfToBDist = distance(bPos, selfPos);
+
+    if (selfToBDist < rule1Distance) {
+      rule1PerceivedCenter += bPos;
+      rule1NearbyBoidsCount++;
+    }
+
+    if (selfToBDist < rule2Distance) {
+      rule2OffsetToKeepDistance -= bPos - selfPos;
+    }
+
+    if (selfToBDist < rule3Distance) {
+      rule3PerceivedVelocity += vel[b];
+      rule3NearbyBoidsCount++;
+    }
+  }
+
+  // Velocity change due to rule 2: separation.
+  glm::vec3 velChange = rule2OffsetToKeepDistance * rule2Scale;
+
+  // Add velocity change due to rule 1: cohesion.
+  if (rule1NearbyBoidsCount > 0) {
+    // From self position to center of mass of nearby boids.
+    velChange += ((rule1PerceivedCenter/ (float)rule1NearbyBoidsCount) - selfPos) * rule1Scale;
+  }
+
+  // Add velocity change due to rule 3: alignment.
+  if (rule3NearbyBoidsCount > 0) {
+    // Average velocity of nearby boids, scaled.
+    velChange += (rule3PerceivedVelocity / (float)rule3NearbyBoidsCount) * rule3Scale;
+  }
+
+  return velChange;
 }
 
 /**
@@ -245,6 +441,23 @@ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos,
   // Compute a new velocity based on pos and vel1
   // Clamp the speed
   // Record the new velocity into vel2. Question: why NOT vel1?
+
+  int iSelf = blockIdx.x*blockDim.x + threadIdx.x;
+  if (iSelf >= N) {
+    // This is one of the extra threads in the last block.
+    return;
+  }
+
+  glm::vec3 newVel = vel1[iSelf] + computeVelocityChange(N, iSelf, pos, vel1);
+
+  float speed = glm::length(newVel);
+  if (speed > maxSpeed) {
+    // Normalizes the newVel vector so that newVel has same direction but with
+    // maxSpeed length.
+    newVel = (newVel / speed) * maxSpeed;
+  }
+
+  vel2[iSelf] = newVel;
 }
 
 /**
@@ -278,17 +491,43 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) {
 //          for(x)
 //            for(y)
 //             for(z)? Or some other order?
+// Most memory efficient method is for(z) for(y) for(x).
 __device__ int gridIndex3Dto1D(int x, int y, int z, int gridResolution) {
   return x + y * gridResolution + z * gridResolution * gridResolution;
 }
 
-__global__ void kernComputeIndices(int N, int gridResolution,
-  glm::vec3 gridMin, float inverseCellWidth,
-  glm::vec3 *pos, int *indices, int *gridIndices) {
-    // TODO-2.1
-    // - Label each boid with the index of its grid cell.
-    // - Set up a parallel array of integer indices as pointers to the actual
-    //   boid data in pos and vel1/vel2
+__global__ void kernComputeIndices(
+  int N,
+  int gridResolution,
+  glm::vec3 gridMin,
+  float inverseCellWidth,
+  glm::vec3 *pos,
+  int *indices,
+  int *gridIndices
+) {
+  // TODO-2.1
+  // - Label each boid with the index of its grid cell.
+  // - Set up a parallel array of integer indices as pointers to the actual
+  //   boid data in pos and vel1/vel2
+
+  int boidIndex = blockIdx.x*blockDim.x + threadIdx.x;
+
+  glm::ivec3 gridIndex3D = (pos[boidIndex] - gridMin) * inverseCellWidth;
+
+  int gridIndex1D = gridIndex3Dto1D(
+    gridIndex3D.x,
+    gridIndex3D.y,
+    gridIndex3D.z,
+    gridResolution
+  );
+
+  // This thread records its corresponding boid index in the indices array
+  // and its corresponding grid cell index. The array index is irrelevant
+  // because these arrays will be sorted later; what matters is that the
+  // recorded boidIndex and gridIndex1D are stored at the same array index
+  // in their corresponding arrays.
+  indices[boidIndex] = boidIndex;
+  gridIndices[boidIndex] = gridIndex1D;
 }
 
 // LOOK-2.1 Consider how this could be useful for indicating that a cell
@@ -300,12 +539,47 @@ __global__ void kernResetIntBuffer(int N, int *intBuffer, int value) {
   }
 }
 
+__global__ void kernReshuffleBuffers(
+  int N,
+  int *particleArrayIndices,
+  glm::vec3 *particleData1,
+  glm::vec3 *particleData2,
+  glm::vec3 *particleDataReshuffled1,
+  glm::vec3 *particleDataReshuffled2
+) {
+  int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+  if (index < N) {
+    particleDataReshuffled1[index] = particleData1[particleArrayIndices[index]];
+    particleDataReshuffled2[index] = particleData2[particleArrayIndices[index]];
+  }
+}
+
+// Assumes that particleGridIndices has been sorted in ascending order.
 __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices,
   int *gridCellStartIndices, int *gridCellEndIndices) {
-  // TODO-2.1
   // Identify the start point of each cell in the gridIndices array.
   // This is basically a parallel unrolling of a loop that goes
-  // "this index doesn't match the one before it, must be a new cell!"
+  // "this index doesn't match the one before it, must be a new cell!".
+
+  int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+  if (index >= N) {
+    // This is one of the extra threads in the last block.
+    return;
+  }
+
+  // Index 0 in particleGridIndices is always the start of a cell.
+  if (index == 0 || particleGridIndices[index-1] < particleGridIndices[index]) {
+    // The boids of the grid cell with index particleGridIndices[index] start at
+    // index 'index' in particleGridIndices.
+    gridCellStartIndices[particleGridIndices[index]] = index;
+  }
+
+  // Index N-1 in particleGridIndices is always the end of a cell.
+  if (index == N-1 || particleGridIndices[index] < particleGridIndices[index+1]) {
+    // The boids of the grid cell with index particleGridIndices[index] end at
+    // index 'index' in particleGridIndices (inclusive).
+    gridCellEndIndices[particleGridIndices[index]] = index;
+  }
 }
 
 __global__ void kernUpdateVelNeighborSearchScattered(
@@ -313,22 +587,137 @@ __global__ void kernUpdateVelNeighborSearchScattered(
   float inverseCellWidth, float cellWidth,
   int *gridCellStartIndices, int *gridCellEndIndices,
   int *particleArrayIndices,
-  glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
-  // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce
-  // the number of boids that need to be checked.
-  // - Identify the grid cell that this particle is in
-  // - Identify which cells may contain neighbors. This isn't always 8.
-  // - For each cell, read the start/end indices in the boid pointer array.
-  // - Access each boid in the cell and compute velocity change from
-  //   the boids rules, if this boid is within the neighborhood distance.
-  // - Clamp the speed change before putting the new speed in vel2
+  glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2
+) {
+  // Update a boid's velocity using the uniform grid to reduce the number of 
+  // boids that need to be checked.
+
+  int iSelf = (blockIdx.x * blockDim.x) + threadIdx.x;
+  if (iSelf >= N) {
+    // This is one of the extra threads in the last block.
+    return;
+  }
+
+  glm::vec3 selfPos = pos[iSelf];
+
+  // Position of this boid relative to the origin of the grid, scaled according to
+  // the uniform grid coordinate system.
+  glm::vec3 selfPosInGridSpace = (pos[iSelf] - gridMin) * inverseCellWidth;
+
+  // Grid cell index that contains this boid. (Cast vec3 to ivec3 to obtain integer
+  // grid cell index from position relative to grid.)
+  glm::ivec3 gridIndex3D = selfPosInGridSpace;
+
+  // Data to be computed by iterating over nearby cells and the boids they contain.
+  // Center of mass of nearby boids, to be computed.
+  glm::vec3 rule1PerceivedCenter(0.0f, 0.0f, 0.0f);
+  glm::vec3 rule2OffsetToKeepDistance(0.0f, 0.0f, 0.0f);
+  // Average velocity of nearby boids, to be computed.
+  glm::vec3 rule3PerceivedVelocity(0.0f, 0.0f, 0.0f);
+  int rule1NearbyBoidsCount = 0;
+  int rule3NearbyBoidsCount = 0;
+
+  // Iterate over the grid cells that overlap with the sphere centered around this
+  // boid and whose radius is max(rule1Distance, rule2Distance, rule3Distance). This
+  // way, we narrow down the search to only the grid cells that contain boids that
+  // could satisfy at least one of the neighborhood distance rules.
+  //
+  // A cell's width/height/depth is twice the maximum neighborhood distance in simulation
+  // space. Since a cell's width is 1 in grid space, the maximum neighborhood distance in 
+  // grid space is 0.5.
+  //
+  // Since a cell's width/height/depth is twice the maximum neighborhood size, the neighborhood
+  // can only occupy 1 or 2 adjacent grid cells in each dimension (x, y, z). Determine which.
+  float maxNeighborhoodDistInGridSpace = 0.5f;
+  int minX = imax(int(selfPosInGridSpace.x - maxNeighborhoodDistInGridSpace), 0);
+  int maxX = imin(int(selfPosInGridSpace.x + maxNeighborhoodDistInGridSpace), gridResolution-1);
+  int minY = imax(int(selfPosInGridSpace.y - maxNeighborhoodDistInGridSpace), 0);
+  int maxY = imin(int(selfPosInGridSpace.y + maxNeighborhoodDistInGridSpace), gridResolution-1);
+  int minZ = imax(int(selfPosInGridSpace.z - maxNeighborhoodDistInGridSpace), 0);
+  int maxZ = imin(int(selfPosInGridSpace.z + maxNeighborhoodDistInGridSpace), gridResolution-1);
+  // This is the most memory efficient order for the iteration. Because of the gridIndex3Dto1D 
+  // mapping of gridCellStartIndices/gridCellEndIndices, if we fix z and y, all (x, y, z) grid
+  // coordinates map to contiguous elements in these arrays.
+  for (int z = minZ; z <= maxZ; ++z) {
+    for (int y = minY; y <= maxY; ++y) {
+      for (int x = minX; x <= maxX; ++x) {
+        glm::ivec3 neighGridIndex3D = glm::ivec3(x, y, z);
+
+        int neighGridIndex1D = gridIndex3Dto1D(
+          neighGridIndex3D.x,
+          neighGridIndex3D.y,
+          neighGridIndex3D.z,
+          gridResolution
+        );
+
+        int neighGridStartIndex = gridCellStartIndices[neighGridIndex1D];
+        if (neighGridStartIndex == -1) {
+          // -1 is a special value for grid cells that don't contain boids.
+          continue;
+        }
+        // Inclusive.
+        int neighGridEndIndex = gridCellEndIndices[neighGridIndex1D];
+
+        for (int i = neighGridStartIndex; i <= neighGridEndIndex; ++i) {
+          int b = particleArrayIndices[i];
+
+          if (b == iSelf) continue;
+
+          glm::vec3 bPos = pos[b];
+
+          float selfToBDist = distance(bPos, selfPos);
+
+          if (selfToBDist < rule1Distance) {
+            rule1PerceivedCenter += bPos;
+            rule1NearbyBoidsCount++;
+          }
+
+          if (selfToBDist < rule2Distance) {
+            rule2OffsetToKeepDistance -= bPos - selfPos;
+          }
+
+          if (selfToBDist < rule3Distance) {
+            rule3PerceivedVelocity += vel1[b];
+            rule3NearbyBoidsCount++;
+          }
+        }
+      }
+    }
+  }
+
+  // Velocity change due to rule 2: separation.
+  glm::vec3 velChange = rule2OffsetToKeepDistance * rule2Scale;
+
+  // Add velocity change due to rule 1: cohesion.
+  if (rule1NearbyBoidsCount > 0) {
+    // From self position to center of mass of nearby boids.
+    velChange += ((rule1PerceivedCenter/ (float)rule1NearbyBoidsCount) - selfPos) * rule1Scale;
+  }
+
+  // Add velocity change due to rule 3: alignment.
+  if (rule3NearbyBoidsCount > 0) {
+    // Average velocity of nearby boids, scaled.
+    velChange += (rule3PerceivedVelocity / (float)rule3NearbyBoidsCount) * rule3Scale;
+  }
+
+  glm::vec3 newVel = vel1[iSelf] + velChange;
+
+  float speed = glm::length(newVel);
+  if (speed > maxSpeed) {
+    // Normalizes the newVel vector so that newVel has same direction but with
+    // maxSpeed length.
+    newVel = (newVel / speed) * maxSpeed;
+  }
+
+  vel2[iSelf] = newVel;
 }
 
 __global__ void kernUpdateVelNeighborSearchCoherent(
   int N, int gridResolution, glm::vec3 gridMin,
   float inverseCellWidth, float cellWidth,
   int *gridCellStartIndices, int *gridCellEndIndices,
-  glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
+  glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2
+) {
   // TODO-2.3 - This should be very similar to kernUpdateVelNeighborSearchScattered,
   // except with one less level of indirection.
   // This should expect gridCellStartIndices and gridCellEndIndices to refer
@@ -341,6 +730,126 @@ __global__ void kernUpdateVelNeighborSearchCoherent(
   // - Access each boid in the cell and compute velocity change from
   //   the boids rules, if this boid is within the neighborhood distance.
   // - Clamp the speed change before putting the new speed in vel2
+
+  // Update a boid's velocity using the uniform grid to reduce the number of 
+  // boids that need to be checked.
+
+  int iSelf = (blockIdx.x * blockDim.x) + threadIdx.x;
+  if (iSelf >= N) {
+    // This is one of the extra threads in the last block.
+    return;
+  }
+
+  glm::vec3 selfPos = pos[iSelf];
+
+  // Position of this boid relative to the origin of the grid, scaled according to
+  // the uniform grid coordinate system.
+  glm::vec3 selfPosInGridSpace = (pos[iSelf] - gridMin) * inverseCellWidth;
+
+  // Grid cell index that contains this boid. (Cast vec3 to ivec3 to obtain integer
+  // grid cell index from position relative to grid.)
+  glm::ivec3 gridIndex3D = selfPosInGridSpace;
+
+  // Data to be computed by iterating over nearby cells and the boids they contain.
+  // Center of mass of nearby boids, to be computed.
+  glm::vec3 rule1PerceivedCenter(0.0f, 0.0f, 0.0f);
+  glm::vec3 rule2OffsetToKeepDistance(0.0f, 0.0f, 0.0f);
+  // Average velocity of nearby boids, to be computed.
+  glm::vec3 rule3PerceivedVelocity(0.0f, 0.0f, 0.0f);
+  int rule1NearbyBoidsCount = 0;
+  int rule3NearbyBoidsCount = 0;
+
+  // Iterate over the grid cells that overlap with the sphere centered around this
+  // boid and whose radius is max(rule1Distance, rule2Distance, rule3Distance). This
+  // way, we narrow down the search to only the grid cells that contain boids that
+  // could satisfy at least one of the neighborhood distance rules.
+  //
+  // A cell's width/height/depth is twice the maximum neighborhood distance in simulation
+  // space. Since a cell's width is 1 in grid space, the maximum neighborhood distance in 
+  // grid space is 0.5.
+  //
+  // Since a cell's width/height/depth is twice the maximum neighborhood size, the neighborhood
+  // can only occupy 1 or 2 adjacent grid cells in each dimension (x, y, z). Determine which.
+  float maxNeighborhoodDistInGridSpace = 0.5f;
+  int minX = imax(int(selfPosInGridSpace.x - maxNeighborhoodDistInGridSpace), 0);
+  int maxX = imin(int(selfPosInGridSpace.x + maxNeighborhoodDistInGridSpace), gridResolution-1);
+  int minY = imax(int(selfPosInGridSpace.y - maxNeighborhoodDistInGridSpace), 0);
+  int maxY = imin(int(selfPosInGridSpace.y + maxNeighborhoodDistInGridSpace), gridResolution-1);
+  int minZ = imax(int(selfPosInGridSpace.z - maxNeighborhoodDistInGridSpace), 0);
+  int maxZ = imin(int(selfPosInGridSpace.z + maxNeighborhoodDistInGridSpace), gridResolution-1);
+  // This is the most memory efficient order for the iteration. Because of the gridIndex3Dto1D 
+  // mapping of gridCellStartIndices/gridCellEndIndices, if we fix z and y, all (x, y, z) grid
+  // coordinates map to contiguous elements in these arrays.
+  for (int z = minZ; z <= maxZ; ++z) {
+    for (int y = minY; y <= maxY; ++y) {
+      for (int x = minX; x <= maxX; ++x) {
+        glm::ivec3 neighGridIndex3D = glm::ivec3(x, y, z);
+
+        int neighGridIndex1D = gridIndex3Dto1D(
+          neighGridIndex3D.x,
+          neighGridIndex3D.y,
+          neighGridIndex3D.z,
+          gridResolution
+        );
+
+        int neighGridStartIndex = gridCellStartIndices[neighGridIndex1D];
+        if (neighGridStartIndex == -1) {
+          // -1 is a special value for grid cells that don't contain boids.
+          continue;
+        }
+        // Inclusive.
+        int neighGridEndIndex = gridCellEndIndices[neighGridIndex1D];
+
+        for (int b = neighGridStartIndex; b <= neighGridEndIndex; ++b) {
+          if (b == iSelf) continue;
+
+          glm::vec3 bPos = pos[b];
+
+          float selfToBDist = distance(bPos, selfPos);
+
+          if (selfToBDist < rule1Distance) {
+            rule1PerceivedCenter += bPos;
+            rule1NearbyBoidsCount++;
+          }
+
+          if (selfToBDist < rule2Distance) {
+            rule2OffsetToKeepDistance -= bPos - selfPos;
+          }
+
+          if (selfToBDist < rule3Distance) {
+            rule3PerceivedVelocity += vel1[b];
+            rule3NearbyBoidsCount++;
+          }
+        }
+      }
+    }
+  }
+
+  // Velocity change due to rule 2: separation.
+  glm::vec3 velChange = rule2OffsetToKeepDistance * rule2Scale;
+
+  // Add velocity change due to rule 1: cohesion.
+  if (rule1NearbyBoidsCount > 0) {
+    // From self position to center of mass of nearby boids.
+    velChange += ((rule1PerceivedCenter/ (float)rule1NearbyBoidsCount) - selfPos) * rule1Scale;
+  }
+
+  // Add velocity change due to rule 3: alignment.
+  if (rule3NearbyBoidsCount > 0) {
+    // Average velocity of nearby boids, scaled.
+    velChange += (rule3PerceivedVelocity / (float)rule3NearbyBoidsCount) * rule3Scale;
+  }
+
+  glm::vec3 newVel = vel1[iSelf] + velChange;
+
+  float speed = glm::length(newVel);
+  if (speed > maxSpeed) {
+    // Normalizes the newVel vector so that newVel has same direction but with
+    // maxSpeed length.
+    newVel = (newVel / speed) * maxSpeed;
+  }
+
+  vel2[iSelf] = newVel;
 }
 
 /**
@@ -349,21 +858,96 @@ __global__ void kernUpdateVelNeighborSearchCoherent(
 void Boids::stepSimulationNaive(float dt) {
   // TODO-1.2 - use the kernels you wrote to step the simulation forward in time.
   // TODO-1.2 ping-pong the velocity buffers
+
+  // numObjects was initialized in Boids::initSimulation with the total number
+  // of boids in the simulation. y and z components of dim3 default to 1.
+  dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+  
+  kernUpdateVelocityBruteForce<<<fullBlocksPerGrid, blockSize>>>(
+    numObjects,
+    // Initialized in Boids::initSimulation with random positions.
+    dev_pos,
+    dev_vel1,
+    dev_vel2
+  );
+
+  // dev_vel2 contains the latest velocities.
+  kernUpdatePos<<<fullBlocksPerGrid, blockSize>>>(numObjects, dt, dev_pos, dev_vel2);
+
+  // Boids::copyBoidsToVBO uses dev_vel1 to pass velocities to the visualization
+  // code, which will find the latest velocities in dev_vel1 after the swap.
+  std::swap(dev_vel1, dev_vel2);
 }
 
 void Boids::stepSimulationScatteredGrid(float dt) {
-  // TODO-2.1
+  int N = numObjects;
   // Uniform Grid Neighbor search using Thrust sort.
-  // In Parallel:
-  // - label each particle with its array index as well as its grid index.
-  //   Use 2x width grids.
-  // - Unstable key sort using Thrust. A stable sort isn't necessary, but you
-  //   are welcome to do a performance comparison.
-  // - Naively unroll the loop for finding the start and end indices of each
-  //   cell's data pointers in the array of boid indices
-  // - Perform velocity updates using neighbor search
-  // - Update positions
-  // - Ping-pong buffers as needed
+
+  // Label each particle with its array index as well as its grid index.
+  // TODO: use 2x width grids???
+  dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize);
+  kernComputeIndices<<<fullBlocksPerGrid, blockSize>>>(
+    N,
+    gridSideCount,
+    gridMinimum,
+    gridInverseCellWidth,
+    dev_pos,
+    dev_particleArrayIndices,
+    dev_particleGridIndices
+  );
+
+  // Sort dev_particleGridIndices and dev_particleArrayIndices simultaneously
+  // in ascending order according to grid cell index (dev_particleGridIndices).
+  thrust::sort_by_key(
+    dev_thrust_particleGridIndices,
+    dev_thrust_particleGridIndices + N,
+    dev_thrust_particleArrayIndices
+  );
+
+  // Naively unroll the loop for finding the start and end indices of each cell's 
+  // data pointers in the array of boid indices.
+  
+  // A grid with at least gridCellCount threads.
+  dim3 fullBlocksPerGridForGrid((gridCellCount + blockSize - 1) / blockSize);
+  // Set cell start index and cell end index arrays to -1. That way, cells that
+  // don't contain boids will continue to be -1 after kernIdentifyCellStartEnd and
+  // we can avoid trying to process them.
+  //
+  // This kernel is launched with at least gridCellCount because that's the number
+  // of elements in dev_gridCellStartIndices and dev_gridCellEndIndices. Each thread
+  // resets the value for a single grid cell.
+  kernResetIntBuffer<<<fullBlocksPerGridForGrid, blockSize>>>(gridCellCount, dev_gridCellStartIndices, -1);
+  kernResetIntBuffer<<<fullBlocksPerGridForGrid, blockSize>>>(gridCellCount, dev_gridCellEndIndices, -1);
+  // This kernel is launched with at least N = <number of boids> threads because that's
+  // the number of elements in dev_particleGridIndices. Each thread determines if its
+  // boid is the start/end of that boid's grid cell.
+  kernIdentifyCellStartEnd<<<fullBlocksPerGrid, blockSize>>>(
+    N,
+    dev_particleGridIndices,
+    dev_gridCellStartIndices,
+    dev_gridCellEndIndices
+  );
+
+  // Update current velocities (dev_vel1) according to 3 rules and store them in dev_vel2.
+  kernUpdateVelNeighborSearchScattered<<<fullBlocksPerGrid, blockSize>>>(
+    N,
+    gridSideCount,
+    gridMinimum,
+    gridInverseCellWidth,
+    gridCellWidth,
+    dev_gridCellStartIndices,
+    dev_gridCellEndIndices,
+    dev_particleArrayIndices,
+    dev_pos,
+    dev_vel1,
+    dev_vel2
+  );
+
+  // Update positions according to updated velocities found in dev_vel2.
+  kernUpdatePos<<<fullBlocksPerGrid, blockSize>>>(N, dt, dev_pos, dev_vel2);
+
+  // Ping-pong velocity buffers.
+  std::swap(dev_vel1, dev_vel2);
 }
 
 void Boids::stepSimulationCoherentGrid(float dt) {
@@ -382,14 +966,107 @@ void Boids::stepSimulationCoherentGrid(float dt) {
   // - Perform velocity updates using neighbor search
   // - Update positions
   // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+
+  int N = numObjects;
+  // Uniform Grid Neighbor search using Thrust sort.
+
+  // Label each particle with its array index as well as its grid index.
+  // TODO: use 2x width grids???
+  dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize);
+  kernComputeIndices<<<fullBlocksPerGrid, blockSize>>>(
+    N,
+    gridSideCount,
+    gridMinimum,
+    gridInverseCellWidth,
+    dev_pos,
+    dev_particleArrayIndices,
+    dev_particleGridIndices
+  );
+
+  // Sort dev_particleGridIndices and dev_particleArrayIndices simultaneously
+  // in ascending order according to grid cell index (dev_particleGridIndices).
+  thrust::sort_by_key(
+    dev_thrust_particleGridIndices,
+    dev_thrust_particleGridIndices + N,
+    dev_thrust_particleArrayIndices
+  );
+
+  // Naively unroll the loop for finding the start and end indices of each cell's 
+  // data pointers in the array of boid indices.
+  
+  // A grid with at least gridCellCount threads.
+  dim3 fullBlocksPerGridForGrid((gridCellCount + blockSize - 1) / blockSize);
+  // Set cell start index and cell end index arrays to -1. That way, cells that
+  // don't contain boids will continue to be -1 after kernIdentifyCellStartEnd and
+  // we can avoid trying to process them.
+  //
+  // This kernel is launched with at least gridCellCount because that's the number
+  // of elements in dev_gridCellStartIndices and dev_gridCellEndIndices. Each thread
+  // resets the value for a single grid cell.
+  kernResetIntBuffer<<<fullBlocksPerGridForGrid, blockSize>>>(gridCellCount, dev_gridCellStartIndices, -1);
+  kernResetIntBuffer<<<fullBlocksPerGridForGrid, blockSize>>>(gridCellCount, dev_gridCellEndIndices, -1);
+  // This kernel is launched with at least N = <number of boids> threads because that's
+  // the number of elements in dev_particleGridIndices. Each thread determines if its
+  // boid is the start/end of that boid's grid cell.
+  kernIdentifyCellStartEnd<<<fullBlocksPerGrid, blockSize>>>(
+    N,
+    dev_particleGridIndices,
+    dev_gridCellStartIndices,
+    dev_gridCellEndIndices
+  );
+
+  kernReshuffleBuffers<<<fullBlocksPerGrid, blockSize>>>(
+    N, 
+    dev_particleArrayIndices, 
+    dev_pos,
+    dev_vel1,
+    dev_posReshuffled,
+    dev_vel1Reshuffled
+  );
+
+  // Update current velocities (dev_vel1) according to 3 rules and store them in dev_vel2.
+  kernUpdateVelNeighborSearchCoherent<<<fullBlocksPerGrid, blockSize>>>(
+    N,
+    gridSideCount,
+    gridMinimum,
+    gridInverseCellWidth,
+    gridCellWidth,
+    dev_gridCellStartIndices,
+    dev_gridCellEndIndices,
+    dev_posReshuffled,
+    dev_vel1Reshuffled,
+    dev_vel2
+  );
+
+  // Update positions according to updated velocities found in dev_vel2.
+  kernUpdatePos<<<fullBlocksPerGrid, blockSize>>>(N, dt, dev_posReshuffled, dev_vel2);
+
+  // Ping-pong velocity buffers.
+  std::swap(dev_pos, dev_posReshuffled);
+  std::swap(dev_vel1, dev_vel2);
 }
 
 void Boids::endSimulation() {
   cudaFree(dev_vel1);
+  checkCUDAErrorWithLine("cudaFree dev_vel1 failed!");
   cudaFree(dev_vel2);
+  checkCUDAErrorWithLine("cudaFree dev_vel2 failed!");
   cudaFree(dev_pos);
+  checkCUDAErrorWithLine("cudaFree dev_pos failed!");
 
   // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+  cudaFree(dev_particleArrayIndices);
+  checkCUDAErrorWithLine("cudaFree dev_particleArrayIndices failed!");
+  cudaFree(dev_particleGridIndices);
+  checkCUDAErrorWithLine("cudaFree dev_particleGridIndices failed!");
+  cudaFree(dev_gridCellStartIndices);
+  checkCUDAErrorWithLine("cudaFree dev_gridCellStartIndices failed!");
+  cudaFree(dev_gridCellEndIndices);
+  checkCUDAErrorWithLine("cudaFree dev_gridCellEndIndices failed!");
+  cudaFree(dev_posReshuffled);
+  checkCUDAErrorWithLine("cudaFree dev_posReshuffled failed!");
+  cudaFree(dev_vel1Reshuffled);
+  checkCUDAErrorWithLine("cudaFree dev_vel1Reshuffled failed!");
 }
 
 void Boids::unitTest() {
diff --git a/src/main.cpp b/src/main.cpp
index fe657ed..bc00324 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -17,8 +17,8 @@
 
 // LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID
 #define VISUALIZE 1
-#define UNIFORM_GRID 0
-#define COHERENT_GRID 0
+#define UNIFORM_GRID 1
+#define COHERENT_GRID 1
 
 // LOOK-1.2 - change this to adjust particle count in the simulation
 const int N_FOR_VIS = 5000;