diff --git a/source/MulensModel/binarylens.py b/source/MulensModel/binarylens.py
index b506870e..66d5eb1d 100644
--- a/source/MulensModel/binarylens.py
+++ b/source/MulensModel/binarylens.py
@@ -95,7 +95,10 @@ def source_x(self):
         if self._source_x is None:
             print('traj.x', self.trajectory.x)
             print('type', type(self.trajectory.x))
-            self._source_x = float(self.trajectory.x)
+            if isinstance(self.trajectory.x, np.ndarray):
+                self._source_x = float(self.trajectory.x[0])
+            else:
+                self._source_x = float(self.trajectory.x)
 
         return self._source_x
 
@@ -112,7 +115,10 @@ def source_y(self):
         magnification calculation.
         """
         if self._source_y is None:
-            self._source_y = float(self.trajectory.y)
+            if isinstance(self.trajectory.x, np.ndarray):
+                self._source_y = float(self.trajectory.y[0])
+            else:
+                self._source_y = float(self.trajectory.y)
 
         return self._source_y
 
@@ -354,7 +360,10 @@ def source_x(self):
             print('traj.x', self.trajectory.x)
             print('type(traj.x)', type(self.trajectory.x))
             print('x_shift', x_shift)
-            self._source_x = float(self.trajectory.x + x_shift)
+            if isinstance(self.trajectory.x, np.ndarray):
+                self._source_x = float(self.trajectory.x[0] + x_shift)
+            else:
+                self._source_x = float(self.trajectory.x + x_shift)
 
         return self._source_x
 
@@ -815,7 +824,10 @@ def source_x(self):
         # AdaptiveContouring uses different coordinates conventions,
         # so we have to transform the coordinates below.
         if self._source_x is None:
-            self._source_x = float(-self.trajectory.x)
+            if isinstance(self.trajectory.x, np.ndarray):
+                self._source_x = float(-self.trajectory.x[0])
+            else:
+                self._source_x = float(-self.trajectory.x)
 
         return self._source_x
 
@@ -829,6 +841,9 @@ def source_y(self):
         # AdaptiveContouring uses different coordinates conventions,
         # so we have to transform the coordinates below.
         if self._source_y is None:
-            self._source_y = float(-self.trajectory.y)
+            if isinstance(self.trajectory.y, np.ndarray):
+                self._source_y = float(-self.trajectory.y[0])
+            else:
+                self._source_y = float(-self.trajectory.y)
 
         return self._source_y
diff --git a/source/MulensModel/binarylenswithshear.py b/source/MulensModel/binarylenswithshear.py
index c70ba97f..6d6384fd 100644
--- a/source/MulensModel/binarylenswithshear.py
+++ b/source/MulensModel/binarylenswithshear.py
@@ -606,6 +606,6 @@ def get_magnification(self):
     @property
     def source_x(self):
         if self._source_x is None:
-            self._source_x = float(self.trajectory.x)
+            self._source_x = float(self.trajectory.x[0])
 
         return self._source_x
diff --git a/source/MulensModel/elliputils.py b/source/MulensModel/elliputils.py
index 6c833190..1bbe23b5 100644
--- a/source/MulensModel/elliputils.py
+++ b/source/MulensModel/elliputils.py
@@ -1,6 +1,7 @@
 import os.path
 import numpy as np
 from scipy.interpolate import interp1d, interp2d
+from scipy.interpolate import RegularGridInterpolator as RGI
 
 import MulensModel as mm
 
@@ -42,7 +43,11 @@ def _read_elliptic_files(self):
                 if line[:3] == "# Y":
                     yy = np.array([float(t) for t in line.split()[2:]])
         pp = np.loadtxt(self.file_3)
-        EllipUtils._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic')
+        try:
+            EllipUtils._interpolate_3 = RGI((xx, yy), pp, method='cubic',
+                                            bounds_error=False)
+        except ValueError:
+            EllipUtils._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic')
         EllipUtils._interpolate_3_min_x = np.min(xx)
         EllipUtils._interpolate_3_max_x = np.max(xx)
         EllipUtils._interpolate_3_min_y = np.min(yy)