Skip to content

Commit

Permalink
bugs in finite_source_uniform_WittMao94 and finite_source_LD_WittMao9…
Browse files Browse the repository at this point in the history
…4 methods
  • Loading branch information
rpoleski committed Jan 18, 2025
1 parent 453acb6 commit 0be8dce
Show file tree
Hide file tree
Showing 4 changed files with 24 additions and 24 deletions.
4 changes: 2 additions & 2 deletions source/MulensModel/elliputils.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,9 @@ def _read_elliptic_files(self):

values = np.loadtxt(self.file_3)
try:
EllipUtils._interpolate_3 = RGI((xx, yy), values.T, method='cubic', bounds_error=False)
EllipUtils._interpolate_3 = RGI((xx, yy), values, method='cubic', bounds_error=False)
except ValueError:
EllipUtils._interpolate_3 = interp2d(xx, yy, values.T, kind='cubic')
EllipUtils._interpolate_3 = interp2d(xx, yy, values, kind='cubic')

EllipUtils._interpolate_3_min_x = np.min(xx)
EllipUtils._interpolate_3_max_x = np.max(xx)
Expand Down
11 changes: 8 additions & 3 deletions source/MulensModel/pointlens.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
from scipy import integrate
from scipy.special import ellipk, ellipe
# These are complete elliptic integrals of the first and the second kind.
from scipy.interpolate import RegularGridInterpolator as RGI
from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3

import MulensModel as mm
Expand Down Expand Up @@ -675,11 +676,12 @@ def get_magnification(self):

return self._magnification

def _get_magnification_WM94(self, u):
def _get_magnification_WM94(self, u, rho=None):
"""
Get point-lens finite-source magnification without LD.
"""
rho = self.trajectory.parameters.rho
if rho is None:
rho = self.trajectory.parameters.rho

if u == rho:
u2 = u**2
Expand Down Expand Up @@ -740,7 +742,10 @@ def _get_ellip3(self, n, k):
cond_4 = (k <= self._ellip_data._interpolate_3_max_y)

if cond_1 and cond_2 and cond_3 and cond_4:
return self._ellip_data._interpolate_3(n, k)[0]
if isinstance(self._ellip_data._interpolate_3, RGI):
return float(self._ellip_data._interpolate_3((n, k)).T)
else:
return self._ellip_data._interpolate_3(n, k)[0]

return ellip3(n, k)

Expand Down
31 changes: 13 additions & 18 deletions source/MulensModel/tests/test_MagnificationCurve.py
Original file line number Diff line number Diff line change
Expand Up @@ -101,46 +101,42 @@ def test_Lee09_and_WittMao94():
t_vec = np.array([3.5, 2., 1., 0.5, 0.])

# The values below were calculated using code developed by P. Mroz.
expected_0 = np.array([1.01084060513, 1.06962639343, 1.42451408166,
2.02334097551, 2.13919086656])
expected_1 = np.array([1.01110609638, 1.07461016241, 1.57232954942,
2.21990790526, 2.39458814753])
# expected_2 = np.array([1.0110829794, 1.07404148634, 1.55620547462,
# 2.24809136704, 2.44503143812])
expected_0 = np.array([1.01084060513, 1.06962639343, 1.42451408166, 2.02334097551, 2.13919086656])
expected_1 = np.array([1.01110609638, 1.07461016241, 1.57232954942, 2.21990790526, 2.39458814753])
# expected_2 = np.array([1.0110829794, 1.07404148634, 1.55620547462, 2.24809136704, 2.44503143812])
# The last values are for 2-parameter LD with same settings and lambda=0.3.
# Correction is:
# -lambda*(1-1.25*sqrt(costh))
# and for 1-parameter LD we used:
# 1-gamma*(1-1.5*costh)

# Test uniform source first.
params_0 = mm.ModelParameters(
{'t_0': 0., 'u_0': 0.5, 't_E': 1., 'rho': 1.})
params_0 = mm.ModelParameters({'t_0': 0., 'u_0': 0.5, 't_E': 1., 'rho': 1.})
mag_curve_0 = mm.MagnificationCurve(times=t_vec, parameters=params_0)
methods_0 = [-5., 'finite_source_uniform_Lee09', 5.]
mag_curve_0.set_magnification_methods(methods_0, 'point_source')
results_0 = mag_curve_0.get_point_lens_magnification()
np.testing.assert_almost_equal(expected_0, results_0, decimal=4)

# Then test 1-parameter limb-darkening.
params_1 = mm.ModelParameters(
{'t_0': 0., 'u_0': 0.1, 't_E': 1., 'rho': 1.})
mag_curve_1 = mm.MagnificationCurve(times=t_vec, parameters=params_1,
gamma=0.5)
params_1 = mm.ModelParameters({'t_0': 0., 'u_0': 0.1, 't_E': 1., 'rho': 1.})
mag_curve_1 = mm.MagnificationCurve(times=t_vec, parameters=params_1, gamma=0.5)
methods_1 = [-5., 'finite_source_LD_Lee09', 5.]
mag_curve_1.set_magnification_methods(methods_1, 'point_source')
results_1 = mag_curve_1.get_point_lens_magnification()
np.testing.assert_almost_equal(expected_1, results_1, decimal=3)

# Tests for Witt & Mao 1994 start here
methods_2 = [-5., 'finite_source_uniform_WittMao94', 5.]
mag_curve_0.set_magnification_methods(methods_2, 'point_source')
results_2 = mag_curve_0.get_point_lens_magnification()
mag_curve_2 = mm.MagnificationCurve(times=t_vec, parameters=params_0)
mag_curve_2.set_magnification_methods(methods_2, 'point_source')
results_2 = mag_curve_2.get_point_lens_magnification()
np.testing.assert_almost_equal(expected_0, results_2, decimal=4)

methods_3 = [-5., 'finite_source_LD_WittMao94', 5.]
mag_curve_1.set_magnification_methods(methods_3, 'point_source')
results_3 = mag_curve_1.get_point_lens_magnification()
mag_curve_3 = mm.MagnificationCurve(times=t_vec, parameters=params_1, gamma=0.5)
mag_curve_3.set_magnification_methods(methods_3, 'point_source')
results_3 = mag_curve_3.get_point_lens_magnification()
np.testing.assert_almost_equal(expected_1, results_3, decimal=3)


Expand All @@ -152,8 +148,7 @@ def test_PSPL_for_binary():
t_E = 20.
u_0 = 1.
t_vec = np.array([10., 100.]) * t_E + t_0
params = mm.ModelParameters({
't_0': t_0, 'u_0': u_0, 't_E': t_E, 's': 1.2, 'q': 0.1, 'alpha': 0.})
params = mm.ModelParameters({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': 1.2, 'q': 0.1, 'alpha': 0.})
mag_curve = mm.MagnificationCurve(times=t_vec, parameters=params)
mag_curve.set_magnification_methods(None, 'point_source_point_lens')
u2 = u_0**2 + ((t_vec - t_0) / t_E)**2
Expand Down
2 changes: 1 addition & 1 deletion source/MulensModel/version.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = "3.0.2"
__version__ = "3.0.3"

0 comments on commit 0be8dce

Please sign in to comment.