diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
index 912fc154b3a442008fd8e2dc006e6b4288895f41..1b4a818913a4ebd869d600474e2455d82b8e400a 100644
--- a/src/ipglasma.cpp
+++ b/src/ipglasma.cpp
@@ -143,22 +143,20 @@ WilsonLine& IPGlasma::GetWilsonLine(double x, double y)
     ApplyPeriodicBoundaryConditions(q);
     x=q[0];
     y=q[1];
-    int xind = FindIndex(x, xcoords);
-    int yind = FindIndex(y, ycoords);
-    
+    std::vector<int> coords = LatticeCoordinates(x,y);
+
     // Handle edges
-    if (xind < 0)
-        xind = 0;
-    if (xind >= xcoords.size())
-        xind = xcoords.size()-1;
-
-    if (yind < 0)
-        yind = 0;
-    if (yind >= ycoords.size())
-        yind = ycoords.size()-1;
-    //cout << "Coordinates " << x << ", "  << y << " indeces " << xind << ", " << yind << endl;
+    if (coords[0] < 0)
+        coords[0] = 0;
+    if (coords[0] >= xcoords.size())
+        coords[0] = xcoords.size()-1;
+
+    if (coords[1] < 0)
+        coords[1] = 0;
+    if (coords[1] >= ycoords.size())
+        coords[1] = ycoords.size()-1;
     
-    return wilsonlines[ xind*xcoords.size() + yind];
+    return GetWilsonLine( WilsonLineCoordinate(coords[0],coords[1]));
     
 }
 
@@ -452,38 +450,23 @@ void IPGlasma::SetSchwinger(bool s, double rc)
     }
 }
 
-/* Returns index i for which
- * vec[i]<=val
- * Assumes that vec[i]<vec[i+1]
- * If such index can't be found, returns -1
- */
 
-int FindIndex(double val, std::vector<double> &vec)
+std::vector<int> IPGlasma::LatticeCoordinates(double x, double y)
 {
-    if (val < vec[0]) return -1;
-   
-    int ind=-1;
-   
-    uint start=0; uint end=vec.size()-1;
-    while(end-start>5)
-    {
-        int tmp = static_cast<int>((start+end)/2.0);
-    
-        if (vec[tmp]>=val)
-            end=tmp;
-        else
-            start=tmp;
-    }
-   
-   
-    for (uint i=start; i<=end; i++)
-    {
-        if (vec[i]<=val and vec[i+1]>val)
-        {
-            ind=i;
-            break;
-        }
-    }
-    if (ind == -1) return vec.size()-1;
-    return ind;
+    std::vector<int> ret;
+
+    // Note: My lattice is from -L/2 to L/2, so I need to shift the coordinates
+    x = x + xcoords[xcoords.size()-1];
+    y = y + ycoords[ycoords.size()-1];
+    double lattice_spacing = xcoords[1]-xcoords[0];
+
+    int ix = x/lattice_spacing; 
+    int iy = y/lattice_spacing; 
+
+    return std::vector<int> {ix, iy}; 
+}
+
+int IPGlasma::WilsonLineCoordinate(int  xind, int yind)
+{
+    return xcoords.size()*xind + yind; 
 }
diff --git a/src/ipglasma.hpp b/src/ipglasma.hpp
index 1d3c1e330e158b44536e6f88a961645f940c8d61..2f8c52288dccec5930d2d271d3465d7fcf97de07 100644
--- a/src/ipglasma.hpp
+++ b/src/ipglasma.hpp
@@ -32,7 +32,7 @@ public:
     double Amplitude(double xpom, double q1[2], double q2[2]);
     double AmplitudeImaginaryPart(double xpom, double q1[2], double q2[2] );
 
-    WilsonLine& GetWilsonLine( double x, double y); // Find Wilson line that corresponds to the coordinate
+    WilsonLine& GetWilsonLine( double x, double y); // Find Wilson line at the given point (x,y) [GeV^-1]
     
     std::string InfoStr();
     
@@ -48,6 +48,13 @@ public:
     void SetSchwinger(bool s, double rc=0);
     void ApplyPeriodicBoundaryConditions(double q[2]); 
 
+    std::vector<int> LatticeCoordinates(double x, double y);
+    int WilsonLineCoordinate(int  xind, int yind);
+    WilsonLine& GetWilsonLine(int i) { return wilsonlines[i]; }
+
+    double X(int ix) { return xcoords[ix]; }
+    double Y(int iy) { return ycoords[iy]; }
+
 private: