Merge pull request #137 from rapoliveira/allow-q-gt-1

Allowed q > 1 and added unit tests to check magnification

source/MulensModel/__init__.py

@@ -26,7 +26,15 @@
 
 from .version import __version__
 
-__all__ = ['mulensobjects', 'MODULE_PATH', 'DATA_PATH', 'BinaryLens']
+__all__ = ['mulensobjects', 'MODULE_PATH', 'DATA_PATH', 'BinaryLens',
+           'BinaryLensWithShear', 'Caustics', 'CausticsPointWithShear',
+           'CausticsWithShear', 'Coordinates', 'Event', 'FitData', 'Horizons',
+           'LimbDarkeningCoeffs', 'MagnificationCurve', 'Model',
+           'ModelParameters', 'which_parameters', 'MulensData', 'Lens',
+           'Source', 'MulensSystem', 'orbits', 'PointLens',
+           'get_pspl_magnification', 'PointLensWithShear',
+           'PointLensFiniteSource', 'SatelliteSkyCoord', 'Trajectory',
+           'UniformCausticSampling', 'MAG_ZEROPOINT',  'Utils', '__version__']
 
 MODULE_PATH = path.abspath(__file__)
 for i in range(3):

source/MulensModel/modelparameters.py

@@ -878,8 +878,8 @@ def _check_valid_parameter_values(self, parameters):
 
         for name in ['q']:
             if name in parameters.keys():
-                if parameters[name] <= 0. or parameters[name] >= 1.:
-                    msg = "Parameter {:} has to be in (0, 1) range, not {:}"
+                if parameters[name] <= 0.:
+                    msg = "Parameter {:} has to be larger than 0, not {:}"
                     raise ValueError(msg.format(name, parameters[name]))
 
         for name in ['xi_eccentricity']:
@@ -1234,8 +1234,8 @@ def q(self):
 
     @q.setter
     def q(self, new_q):
-        if new_q < 0. or new_q > 1.:
-            raise ValueError('mass ratio q has to be between 0 and 1')
+        if new_q <= 0.:
+            raise ValueError('mass ratio q has to be larger than 0')
         self.parameters['q'] = new_q
         self._update_sources('q', new_q)
 

source/MulensModel/tests/test_ModelParameters.py

@@ -30,7 +30,7 @@ def test_wrong_type_of_parameters(self):
                                 'shear_G': 0.1, 'alpha': 123.})
         with self.assertRaises(ValueError):
             mm.ModelParameters({'t_0': 123., 'u_0': 1, 't_E': 10., 's': 1.2,
-                               'alpha': 34.56, 'q': 1.5})
+                               'alpha': 34.56, 'q': -0.5})
 
     def test_init_for_2_sources(self):
         """
@@ -150,6 +150,70 @@ def test_positive_t_E():
     assert params.t_E == params.t_eff / abs(params.u_0)
 
 
+def test_q_gt_1_is_good():
+    """
+    Check if the magnification is reproduced by transforming q -> 1/q and
+    alpha -> alpha +/- 180, including for q > 1. See issue #84.
+    """
+    t_0 = 3583.
+    u_0 = 0.3
+    t_E = 12.
+    s = 1.65
+    q = 0.25
+    alpha = 339.0
+    rho = 0.001
+
+    planet = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': s, 'q': q,
+                       'alpha': alpha, 'rho': rho})
+    planet_2 = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': s, 'q': 1/q,
+                         'alpha': alpha+180., 'rho': rho})
+    planet_3 = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': s, 'q': 1/q,
+                         'alpha': alpha-180., 'rho': rho})
+    list_of_methods = [3588., 'VBBL', 3594., 'hexadecapole', 3598.0]
+    planet.set_magnification_methods(list_of_methods)
+    planet_2.set_magnification_methods(list_of_methods)
+    planet_3.set_magnification_methods(list_of_methods)
+    t_checks = [3580, 3589, 3590, 3592, 3593, 3595]
+    magnifications_1 = planet.get_magnification(time=t_checks)
+    magnifications_2 = planet_2.get_magnification(time=t_checks)
+    magnifications_3 = planet_3.get_magnification(time=t_checks)
+
+    assert max(magnifications_1 - magnifications_2) < 1e-10
+    assert max(magnifications_1 - magnifications_3) < 1e-10
+
+
+def test_q_gt_1_is_smooth():
+    """
+    Check that there is a smooth transition between q = 0.97, 0.99 and 1.01.
+    In this case, a trajectory with caustic approaching is adopted.
+    """
+    t_0 = 3583.
+    u_0 = 0.3
+    t_E = 12.
+    s = 2.18
+    q = 0.99
+    alpha = 310.
+    rho = 0.001
+    planet = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': s, 'q': q,
+                       'alpha': alpha, 'rho': rho})
+    q_min = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': s, 'q': q-0.02,
+                      'alpha': alpha, 'rho': rho})
+    q_max = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E, 's': s, 'q': q+0.02,
+                      'alpha': alpha, 'rho': rho})
+
+    planet.set_magnification_methods([3580., 'VBBL', 3595.])
+    q_min.set_magnification_methods([3580., 'VBBL', 3595.])
+    q_max.set_magnification_methods([3580., 'VBBL', 3595.])
+    t_checks = [3571, 3583, 3585.5, 3586, 3586.5, 3592.5]
+    magnification = planet.get_magnification(time=t_checks)
+    diff_min = magnification - q_min.get_magnification(time=t_checks)
+    diff_max = magnification - q_max.get_magnification(time=t_checks)
+    limits = np.array([0.01, 0.01, 0.01, 0.56, 0.01, 0.03])
+
+    assert np.all(abs(diff_min) < limits)
+    assert np.all(abs(diff_max) < limits)
+
+
 def test_rho_t_e_t_star():
     """check if conversions between rho, t_E, and t_star work ok"""
     t_0 = 2450000.

source/MulensModel/version.py

@@ -1 +1 @@
-__version__ = "2.23.0"
+__version__ = "2.24.0"