Replaced interp2d in data folder for issue rpoleski#74, still checkin…

…g difference

data/interpolate_elliptic_integral_3.py

@@ -6,6 +6,7 @@
 from math import sin, cos, sqrt
 from scipy import integrate
 from scipy.interpolate import interp1d, interp2d
+from scipy.interpolate import RegularGridInterpolator as RGI
 from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3
 
 
@@ -45,7 +46,29 @@ def get_ellip(x, y):
     add_y = []
     p = get_ellip(x, y)
 
-    interp_p = interp2d(x, y, p.T, kind='cubic')
+    # 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):
@@ -63,7 +86,8 @@ def get_ellip(x, y):
                 check_p[ix, iy] = 1.
                 check_true_p[ix, iy] = 1.
             else:
-                check_p[ix, iy] = interp_p(cx, cy)[0]
+                # check_p[ix, iy] = interp_p(cx, cy)[0]
+                check_p[ix, iy] = float(interp_rgi((cx, cy)).T)
     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)