Replaced interp2d in pointlens.py, issue rpoleski#74

data/interpolate_elliptic_integral_3.py

@@ -49,27 +49,6 @@ def get_ellip(x, y):
     # interp_p = interp2d(x, y, p.T, kind='cubic')
     interp_rgi = RGI((x, y), p, method='cubic', bounds_error=False)
 
-    # Testing the change between scipy's interp2d and RGI
-    # xnew = np.arange(min(x), max(x), 1e-2)
-    # ynew = np.arange(min(y), max(y), 1e-2)
-    # xxnew, yynew = np.meshgrid(xnew, ynew, indexing='ij', sparse=True)
-    # znew_interp2d = interp_p(xnew, ynew)
-    # znew_rgi = interp_rgi((xxnew, yynew)).T
-    # diff = abs(znew_interp2d - znew_rgi).flatten()
-    # diff_max, diff_mean = diff.max(), diff.mean()  # maximum: 1e-10, 1e-12
-
-    # Plotting to check difference
-    # import matplotlib.pyplot as plt
-    # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
-    # ax1.plot(x, p.T[0, :], 'ro-', xnew, znew_interp2d[0, :], 'b-')
-    # im = ax2.imshow(znew)
-    # plt.colorbar(im, ax=ax2)
-    # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
-    # ax1.plot(x, p.T[0, :], 'ro-', xnew, znew_rgi[0, :], 'b-')
-    # im = ax2.imshow(znew_rgi)
-    # plt.colorbar(im, ax=ax2)
-    # plt.show()
-
     check_x = []
     for i in range(len(x)-1):
         check_ = np.logspace(np.log10(x[i]), np.log10(x[i+1]), n_divide)

source/MulensModel/pointlens.py

@@ -4,6 +4,7 @@
 from math import sin, cos, sqrt, log10
 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 +93,9 @@ 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')
+        # PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic')
+        PointLens._interpolate_3 = RGI((xx, yy), pp.T, method='cubic',
+                                       bounds_error=False)
         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 +598,8 @@ 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]
+            # return PointLens._interpolate_3(n, k)[0]
+            return float(PointLens._interpolate_3((n,k)).T)
         return ellip3(n, k)
 
     def get_point_lens_large_LD_integrated_magnification(self, u, gamma):