diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py
new file mode 100644
index 00000000..d98d58d6
--- /dev/null
+++ b/SimPEG/tests/test_interpolation.py
@@ -0,0 +1,81 @@
+import numpy as np
+import unittest
+from TestUtils import OrderTest
+
+MESHTYPES = ['uniformTensorMesh']
+call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
+call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
+cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
+cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
+cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
+cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
+cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
+cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
+
+
+LOCS = np.random.rand(50,3)*0.6+0.2
+
+class TestInterpolation(OrderTest):
+    name = "Interpolation"
+    meshTypes = MESHTYPES
+    meshDimension = 3
+
+    def getError(self):
+        funX = lambda x, y, z: np.cos(2*np.pi*y)
+        funY = lambda x, y, z: np.cos(2*np.pi*z)
+        funZ = lambda x, y, z: np.cos(2*np.pi*x)
+
+        if 'x' in self.type:
+            anal = call3(funX, LOCS)
+        elif 'y' in self.type:
+            anal = call3(funY, LOCS)
+        elif 'z' in self.type:
+            anal = call3(funZ, LOCS)
+
+        if 'F' in self.type:
+            Fc = cartF3(self.M, funX, funY, funZ)
+            grid = self.M.projectFaceVector(Fc)
+        elif 'E' in self.type:
+            Ec = cartE3(self.M, funX, funY, funZ)
+            grid = self.M.projectEdgeVector(Ec)
+
+        comp = self.M.getInterpolationMat(LOCS, self.type)*grid
+
+        err = np.linalg.norm((comp - anal), np.inf)
+        return err
+
+    def test_orderFx(self):
+        self.type = 'Fx'
+        self.name = 'Interpolation Fx'
+        self.orderTest()
+
+    def test_orderFy(self):
+        self.type = 'Fy'
+        self.name = 'Interpolation Fy'
+        self.orderTest()
+
+    def test_orderFz(self):
+        self.type = 'Fz'
+        self.name = 'Interpolation Fz'
+        self.orderTest()
+
+    def test_orderEx(self):
+        self.type = 'Ex'
+        self.name = 'Interpolation Ex'
+        self.orderTest()
+
+    def test_orderEy(self):
+        self.type = 'Ey'
+        self.name = 'Interpolation Ey'
+        self.orderTest()
+
+    def test_orderEz(self):
+        self.type = 'Ez'
+        self.name = 'Interpolation Ez'
+        self.orderTest()
+
+
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py
index 0d528e93..12eb8c47 100644
--- a/SimPEG/utils/interputils.py
+++ b/SimPEG/utils/interputils.py
@@ -7,22 +7,22 @@ def interpmat(x,y,z,xr,yr,zr):
 
     """ Local interpolation computed for each receiver point in turn """
 
-    nx = max(x.shape)
-    ny = max(y.shape)
-    nz = max(z.shape)
-    npts = max(xr.shape)
+    nx = x.size
+    ny = y.size
+    nz = z.size
+    npts = xr.shape[0]
 
     Q = sp.lil_matrix((npts, nx*ny*nz))
 
     for i in range(npts):
-        # in x-direction   
+        # in x-direction
         im = np.argmin(abs(x-xr[i]))
         if  xr[i] - x[im] >= 0:  # Point on the left
             ind_x1 = im
             ind_x2 = im+1
         elif  xr[i] - x[im] < 0:  # Point on the right
             ind_x1 = im-1
-            ind_x2 = im        
+            ind_x2 = im
         dx1 = xr[i] - x[ind_x1]
         dx2 = x[ind_x2] - xr[i]
         # in y-direction
@@ -32,9 +32,9 @@ def interpmat(x,y,z,xr,yr,zr):
             ind_y2 = im+1
         elif  yr[i] - y[im] < 0:  # Point on the right
             ind_y1 = im-1
-            ind_y2 = im        
+            ind_y2 = im
         dy1 = yr[i] - y[ind_y1]
-        dy2 = y[ind_y2] - yr[i]    
+        dy2 = y[ind_y2] - yr[i]
         # in z-direction
         im = np.argmin(abs(z-zr[i]))
         if  zr[i] - z[im] >= 0:  # Point on the left
@@ -42,9 +42,9 @@ def interpmat(x,y,z,xr,yr,zr):
             ind_z2 = im+1
         elif  zr[i] - z[im] < 0:  # Point on the right
             ind_z1 = im-1
-            ind_z2 = im        
+            ind_z2 = im
         dz1 = zr[i] - z[ind_z1]
-        dz2 = z[ind_z2] - zr[i]    
+        dz2 = z[ind_z2] - zr[i]
         dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
 
         Dx =  x[ind_x2] - x[ind_x1]
@@ -53,7 +53,7 @@ def interpmat(x,y,z,xr,yr,zr):
 
         # Get the row in the matrix
 
-        inds = sub2ind((nx,ny,nz),[ 
+        inds = sub2ind((nx,ny,nz),[
             ( ind_x1,  ind_y2,  ind_z1),
             ( ind_x1,  ind_y1,  ind_z1),
             ( ind_x2,  ind_y1,  ind_z1),
@@ -62,7 +62,7 @@ def interpmat(x,y,z,xr,yr,zr):
             ( ind_x1,  ind_y2,  ind_z2),
             ( ind_x2,  ind_y1,  ind_z2),
             ( ind_x2,  ind_y2,  ind_z2)])
-        
+
         vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz),
                 (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz),
                 (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz),
@@ -74,4 +74,4 @@ def interpmat(x,y,z,xr,yr,zr):
 
         Q[i, mkvc(inds)] = vals
     Q = Q.tocsr()
-    return Q
\ No newline at end of file
+    return Q