diff --git a/data/interpolate_elliptic_integral_1_2.py b/data/interpolate_elliptic_integral_1_2.py
index f74ee9f6..61263163 100644
--- a/data/interpolate_elliptic_integral_1_2.py
+++ b/data/interpolate_elliptic_integral_1_2.py
@@ -2,10 +2,7 @@
 Calculates interpolation tables for elliptical integral of
 the first and second kind.
 """
-import math
 import numpy as np
-from math import sin, cos, sqrt
-from scipy import integrate
 from scipy.interpolate import interp1d
 from scipy.special import ellipk, ellipe
 # These are complete elliptic integrals of the first and the second kind.
diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py
index 58a02653..f38276ca 100644
--- a/data/interpolate_elliptic_integral_3.py
+++ b/data/interpolate_elliptic_integral_3.py
@@ -1,11 +1,9 @@
 """
 Calculates interpolation tables for elliptical integral of the third kind.
 """
-import math
 import numpy as np
-from math import sin, cos, sqrt
-from scipy import integrate
-from scipy.interpolate import interp1d, interp2d
+from scipy.interpolate import interp2d  # , interp1d
+from scipy.interpolate import RegularGridInterpolator as RGI
 from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3
 
 
@@ -22,7 +20,6 @@
 
 
 def get_ellip(x, y):
-    p = []
     z = np.zeros((len(x), len(y)))
     for (i, x_) in enumerate(x):
         for (j, y_) in enumerate(y):
@@ -45,7 +42,10 @@ def get_ellip(x, y):
     add_y = []
     p = get_ellip(x, y)
 
-    interp_p = interp2d(x, y, p.T, kind='cubic')
+    try:
+        interp_p = RGI((x, y), p, method="cubic", bounds_error=False)
+    except ValueError:
+        interp_p = interp2d(x, y, p.T, kind="cubic")
 
     check_x = []
     for i in range(len(x)-1):
@@ -58,12 +58,13 @@ def get_ellip(x, y):
     check_true_p = get_ellip(check_x, check_y)
     check_p = np.zeros((len(check_x), len(check_y)))
     for (ix, cx) in enumerate(check_x):
-        for (iy, cy) in enumerate(check_y):
-            if cy > cx:
-                check_p[ix, iy] = 1.
-                check_true_p[ix, iy] = 1.
-            else:
-                check_p[ix, iy] = interp_p(cx, cy)[0]
+        cond = np.array(check_y) > cx
+        check_p[ix, cond] = 1.
+        check_true_p[ix, cond] = 1.
+        if isinstance(interp_p, RGI):
+            check_p[ix, ~cond] = interp_p((cx, np.array(check_y)[~cond]))
+        else:
+            check_p[ix, ~cond] = interp_p(cx, np.array(check_y)[~cond]).T[0]
     relative_diff_p = np.abs(check_p - check_true_p) / check_true_p
     index = np.unravel_index(relative_diff_p.argmax(), relative_diff_p.shape)
 
diff --git a/source/MulensModel/pointlens.py b/source/MulensModel/pointlens.py
index a187193b..eb541465 100644
--- a/source/MulensModel/pointlens.py
+++ b/source/MulensModel/pointlens.py
@@ -2,8 +2,10 @@
 import warnings
 import numpy as np
 from math import sin, cos, sqrt, log10
+import scipy
 from scipy import integrate
 from scipy.interpolate import interp1d, interp2d
+from scipy.interpolate import RegularGridInterpolator as RGI
 from scipy.special import ellipk, ellipe
 # These are complete elliptic integrals of the first and the second kind.
 from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3
@@ -92,7 +94,12 @@ def _read_elliptic_files(self):
                 if line[:3] == "# Y":
                     yy = np.array([float(t) for t in line.split()[2:]])
         pp = np.loadtxt(file_3)
-        PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic')
+
+        try:
+            PointLens._interpolate_3 = RGI((xx, yy), pp, method='cubic',
+                                           bounds_error=False)
+        except ValueError:
+            PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic')
         PointLens._interpolate_3_min_x = np.min(xx)
         PointLens._interpolate_3_max_x = np.max(xx)
         PointLens._interpolate_3_min_y = np.min(yy)
@@ -595,7 +602,10 @@ def _get_ellip3(self, n, k):
         cond_4 = (k <= PointLens._interpolate_3_max_y)
 
         if cond_1 and cond_2 and cond_3 and cond_4:
-            return PointLens._interpolate_3(n, k)[0]
+            if isinstance(PointLens._interpolate_3, RGI):
+                return float(PointLens._interpolate_3((n, k)).T)
+            else:
+                return PointLens._interpolate_3(n, k)[0]
         return ellip3(n, k)
 
     def get_point_lens_large_LD_integrated_magnification(self, u, gamma):
diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py
index 326b86ff..71fae1bc 100644
--- a/source/MulensModel/version.py
+++ b/source/MulensModel/version.py
@@ -1 +1 @@
-__version__ = "2.19.0"
+__version__ = "2.19.3"