From 8b8e4675c04e3b2805dafdc25a3e591e4531fe94 Mon Sep 17 00:00:00 2001 From: "Raphael A. P. Oliveira" Date: Mon, 4 Dec 2023 17:35:56 +0100 Subject: [PATCH] Replaced interp2d in pointlens.py, issue #74 --- data/interpolate_elliptic_integral_3.py | 21 --------------------- source/MulensModel/pointlens.py | 8 ++++++-- 2 files changed, 6 insertions(+), 23 deletions(-) diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py index f6a96fe9..26d0dd79 100644 --- a/data/interpolate_elliptic_integral_3.py +++ b/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) diff --git a/source/MulensModel/pointlens.py b/source/MulensModel/pointlens.py index a187193b..ec7e36ce 100644 --- a/source/MulensModel/pointlens.py +++ b/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):