From fe58801d1d729b0d60cfb0b4fa21e6bd5c73a6ce Mon Sep 17 00:00:00 2001
From: "Raphael A. P. Oliveira" <roliveira.astro@gmail.com>
Date: Thu, 7 Dec 2023 16:50:48 +0100
Subject: [PATCH] Fixed small bug and added scipy version, issue #74

---
 data/interpolate_elliptic_integral_3.py | 29 ++++++++++++++-----------
 source/MulensModel/pointlens.py         | 16 +++++++++-----
 2 files changed, 27 insertions(+), 18 deletions(-)

diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py
index 26d0dd79..1eb8d105 100644
--- a/data/interpolate_elliptic_integral_3.py
+++ b/data/interpolate_elliptic_integral_3.py
@@ -1,11 +1,12 @@
 """
 Calculates interpolation tables for elliptical integral of the third kind.
 """
-import math
+# import math
 import numpy as np
-from math import sin, cos, sqrt
-from scipy import integrate
-from scipy.interpolate import interp1d, interp2d
+# from math import sin, cos, sqrt
+import scipy
+# from scipy import integrate
+from scipy.interpolate import interp2d  # , interp1d
 from scipy.interpolate import RegularGridInterpolator as RGI
 from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3
 
@@ -46,8 +47,10 @@ def get_ellip(x, y):
     add_y = []
     p = get_ellip(x, y)
 
-    # interp_p = interp2d(x, y, p.T, kind='cubic')
-    interp_rgi = RGI((x, y), p, method='cubic', bounds_error=False)
+    if scipy.__version__ >= "1.9.0":
+        interp_p = RGI((x, y), p, method='cubic', bounds_error=False)
+    else:
+        interp_p = interp2d(x, y, p.T, kind='cubic')
 
     check_x = []
     for i in range(len(x)-1):
@@ -60,13 +63,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]
-                check_p[ix, iy] = float(interp_rgi((cx, cy)).T)
+        cond = np.array(check_y) > cx
+        check_p[ix, cond] = 1.
+        check_true_p[ix, cond] = 1.
+        if scipy.__version__ >= "1.9.0":
+            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 ec7e36ce..b084e336 100644
--- a/source/MulensModel/pointlens.py
+++ b/source/MulensModel/pointlens.py
@@ -2,6 +2,7 @@
 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
@@ -93,9 +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')
-        PointLens._interpolate_3 = RGI((xx, yy), pp.T, method='cubic',
-                                       bounds_error=False)
+
+        if scipy.__version__ >= "1.9.0":
+            PointLens._interpolate_3 = RGI((xx, yy), pp, method='cubic',
+                                           bounds_error=False)
+        else:
+            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)
@@ -598,8 +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]
-            return float(PointLens._interpolate_3((n,k)).T)
+            if scipy.__version__ >= "1.9.0":
+                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):