diff --git a/advection/src/RK4AdvectionKernel.cpp b/advection/src/RK4AdvectionKernel.cpp
index f103339..d9eb928 100644
--- a/advection/src/RK4AdvectionKernel.cpp
+++ b/advection/src/RK4AdvectionKernel.cpp
@@ -6,30 +6,30 @@ using namespace std;
 RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid) { }
 
 std::pair<double, double> RK4AdvectionKernel::advect(int time, double latitude, double longitude) const {
-    auto [v1, u1] = bilinearinterpolation(*grid, time, latitude, longitude);
+    auto [u1, v1] = bilinearinterpolation(*grid, time, latitude, longitude);
 //    lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt);
-    double lon1 = longitude + u1*.5*DT;
-    double lat1 = latitude + v1*.5*DT;
+    double lon1 = longitude + v1 * .5*DT;
+    double lat1 = latitude + u1*.5*DT;
 
 //    (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle]
-    auto [v2, u2] = bilinearinterpolation(*grid, time + 0.5*DT, lat1, lon1);
+    auto [u2, v2] = bilinearinterpolation(*grid, time + 0.5*DT, lat1, lon1);
 
 //    lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt)
-    double lon2 = longitude + u2 * 0.5 * DT;
-    double lat2 = latitude + v2 * 0.5 * DT;
+    double lon2 = longitude + v2 * 0.5 * DT;
+    double lat2 = latitude + u2 * 0.5 * DT;
 
 //    (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle]
-    auto [v3, u3] = bilinearinterpolation(*grid, time + 0.5 * DT, lat2, lon2);
+    auto [u3, v3] = bilinearinterpolation(*grid, time + 0.5 * DT, lat2, lon2);
 
 //    lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt)
-    double lon3 = longitude + u3 * DT;
-    double lat3 = latitude + v3 * DT;
+    double lon3 = longitude + v3 * DT;
+    double lat3 = latitude + u3 * DT;
 
 //    (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
-    auto [v4, u4] = bilinearinterpolation(*grid, time + DT, lat3, lon3);
+    auto [u4, v4] = bilinearinterpolation(*grid, time + DT, lat3, lon3);
 
-    double lonFinal = longitude + (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT;
-    double latFinal = latitude + (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT;
+    double lonFinal = longitude + (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT;
+    double latFinal = latitude + (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT;
 
     return {latFinal, lonFinal};
 }
diff --git a/advection/src/UVGrid.cpp b/advection/src/UVGrid.cpp
index 6673237..63f9079 100644
--- a/advection/src/UVGrid.cpp
+++ b/advection/src/UVGrid.cpp
@@ -34,6 +34,11 @@ UVGrid::UVGrid() {
 }
 
 const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const {
+    if(timeIndex < 0 or timeIndex >= timeSize
+        or latIndex < 0 or latIndex >= latSize
+        or lonIndex < 0 or lonIndex >= lonSize) {
+        throw std::out_of_range("Index out of bounds");
+    }
     size_t index = timeIndex * (latSize * lonSize) + latIndex * lonIndex + lonIndex;
     return uvData[index];
 }
diff --git a/advection/src/main.cpp b/advection/src/main.cpp
index 84ab7e8..f7e2e52 100644
--- a/advection/src/main.cpp
+++ b/advection/src/main.cpp
@@ -5,33 +5,37 @@
 #include <ranges>
 #include <chrono>
 
+#define NotAKernelError "Template parameter T must derive from AdvectionKernel"
+
 using namespace std;
 
+template <typename AdvectionKernelImpl>
+void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, double latstart, double lonstart) {
+    static_assert(std::is_base_of<AdvectionKernel, AdvectionKernelImpl>::value, NotAKernelError);
+    double lat1 = latstart, lon1 = lonstart;
+    for(int time = 100; time <= 10000; time += DT) {
+        cout << "lat = " << lat1 << " lon = " << lon1 << endl;
+        auto [templat, templon] = kernel.advect(time, lat1, lon1);
+        lat1 = templat;
+        lon1 = templon;
+    }
+}
+
+
 int main() {
     std::shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
 //    uvGrid->streamSlice(cout, 100);
 
     EulerAdvectionKernel kernelEuler = EulerAdvectionKernel(uvGrid);
+
     RK4AdvectionKernel kernelRK4 = RK4AdvectionKernel(uvGrid);
 
-    double latstart = 54.860801, lonstart = 4.075492;
+    double latstart = 52.881770, lonstart = 3.079979;
 
-    double lat1 = latstart, lon1 = lonstart;
     cout << "======= Euler Integration =======" << endl;
-    for(int time = 100; time <= 10000; time += DT) {
-        cout << "lat = " << lat1 << " lon = " << lon1 << endl;
-        auto [templat, templon] = kernelEuler.advect(time, lat1, lon1);
-        lat1 = templat;
-        lon1 = templon;
-    }
+    advectForSomeTime(*uvGrid, kernelEuler, latstart, lonstart);
     cout << "======= RK4 Integration =======" << endl;
-    lat1 = latstart, lon1 = lonstart;
-    for(int time = 100; time <= 10000; time += DT) {
-        cout << "lat = " << lat1 << " lon = " << lon1 << endl;
-        auto [templat, templon] = kernelRK4.advect(time, lat1, lon1);
-        lat1 = templat;
-        lon1 = templon;
-    }
+    advectForSomeTime(*uvGrid, kernelRK4, latstart, lonstart);
 
     return 0;
 }