diff --git a/source/MulensModel/binarylens.py b/source/MulensModel/binarylens.py index 4c77205bd..d708de5a4 100644 --- a/source/MulensModel/binarylens.py +++ b/source/MulensModel/binarylens.py @@ -4,7 +4,7 @@ from MulensModel.binarylensimports import ( _vbbl_wrapped, _adaptive_contouring_wrapped, - _vbbl_binary_mag_dark, _vbbl_binary_mag_0, + _vbbl_binary_mag_dark, _vbbl_binary_mag_finite, _vbbl_binary_mag_point, _vbbl_SG12_5, _adaptive_contouring_linear, _solver) import MulensModel as mm @@ -335,9 +335,7 @@ def _point_source_magnification_VBBL(self, source_x, source_y): """ args = [self.separation, self.mass_2/self.mass_1, source_x, source_y] - args += [0., 0., 0.] # THIS IS A HACK THAT SHOULD BE REMOVED BY PROPER IMPORT OF VBBL - - return _vbbl_binary_mag_0(*[float(arg) for arg in args]) + return _vbbl_binary_mag_point(*[float(arg) for arg in args]) def _point_source_magnification(self, source_x, source_y): """ @@ -626,24 +624,31 @@ def vbbl_magnification(self, source_x, source_y, rho, if not _vbbl_wrapped: raise ValueError('VBBL was not imported properly') + args = [float(self.separation), float(self.mass_2 / self.mass_1), + float(source_x), float(source_y), float(rho), float(accuracy)] + + u_limb_darkening = self._get_u(gamma, u_limb_darkening) + if u_limb_darkening is None: + _vbbl_function = _vbbl_binary_mag_finite + else: + args += [u_limb_darkening] + _vbbl_function = _vbbl_binary_mag_dark + + return _vbbl_function(*args) + + def _get_u(self, gamma, u_limb_darkening): + """ + Check if one one parameter is defined, extract the u value, + and make sure it's a float. + """ if gamma is not None and u_limb_darkening is not None: raise ValueError('Only one limb darkening parameters can be set' + ' in BinaryLens.vbbl_magnification()') elif gamma is not None: - u_limb_darkening = float(mm.Utils.gamma_to_u(gamma)) + out = float(mm.Utils.gamma_to_u(gamma)) elif u_limb_darkening is not None: - u_limb_darkening = float(u_limb_darkening) + out = float(u_limb_darkening) else: - u_limb_darkening = float(0.0) + out = None - s = float(self.separation) - q = float(self.mass_2 / self.mass_1) - x = float(source_x) - y = float(source_y) - rho = float(rho) - accuracy = float(accuracy) - - magnification = _vbbl_binary_mag_dark( - s, q, x, y, rho, u_limb_darkening, accuracy) - - return magnification + return out diff --git a/source/MulensModel/binarylensimports.py b/source/MulensModel/binarylensimports.py index 6da802e7c..585c9d583 100644 --- a/source/MulensModel/binarylensimports.py +++ b/source/MulensModel/binarylensimports.py @@ -49,13 +49,17 @@ def _import_compiled_VBBL(): _get_path_2('VBBL', "VBBinaryLensingLibrary_wrapper.so"), "VBBL") _vbbl_wrapped = (vbbl is not None) if not _vbbl_wrapped: - return (_vbbl_wrapped, None, None, None, None) + return (_vbbl_wrapped, None, None, None, None, None, None) - vbbl.VBBinaryLensing_BinaryMagDark.argtypes = 7 * [ctypes.c_double] - vbbl.VBBinaryLensing_BinaryMagDark.restype = ctypes.c_double + def _set_in_out(function, n_double): + """set input to n_double doubles and output to double""" + function.argtypes = n_double * [ctypes.c_double] + function.restype = ctypes.c_double - vbbl.VBBinaryLensing_BinaryMag0.argtypes = 7 * [ctypes.c_double] - vbbl.VBBinaryLensing_BinaryMag0.restype = ctypes.c_double + _set_in_out(vbbl.VBBinaryLensing_BinaryMagDark, 7) + _set_in_out(vbbl.VBBinaryLensing_BinaryMagFinite, 6) + _set_in_out(vbbl.VBBinaryLensing_BinaryMagPoint, 4) + _set_in_out(vbbl.VBBinaryLensing_BinaryMagPointShear, 7) vbbl.VBBL_SG12_5.argtypes = 12 * [ctypes.c_double] vbbl.VBBL_SG12_5.restype = np.ctypeslib.ndpointer( @@ -66,8 +70,11 @@ def _import_compiled_VBBL(): dtype=ctypes.c_double, shape=(18,)) return (_vbbl_wrapped, - vbbl.VBBinaryLensing_BinaryMagDark, vbbl.VBBL_SG12_5, - vbbl.VBBinaryLensing_BinaryMag0, vbbl.VBBL_SG12_9) + vbbl.VBBinaryLensing_BinaryMagDark, + vbbl.VBBinaryLensing_BinaryMagFinite, + vbbl.VBBinaryLensing_BinaryMagPoint, + vbbl.VBBinaryLensing_BinaryMagPointShear, + vbbl.VBBL_SG12_5, vbbl.VBBL_SG12_9) def _import_compiled_AdaptiveContouring(): @@ -87,16 +94,21 @@ def _import_compiled_AdaptiveContouring(): # Check import and try manually compiled versions. if _vbbl_wrapped: _vbbl_binary_mag_dark = mm_vbbl.VBBinaryLensing_BinaryMagDark - _vbbl_binary_mag_0 = mm_vbbl.VBBinaryLensing_BinaryMag0 + _vbbl_binary_mag_finite = mm_vbbl.VBBinaryLensing_BinaryMagFinite + _vbbl_binary_mag_point = mm_vbbl.VBBinaryLensing_BinaryMagPoint + _vbbl_binary_mag_point_shear = mm_vbbl.VBBinaryLensing_BinaryMagPointShear _vbbl_SG12_5 = mm_vbbl.VBBL_SG12_5 _vbbl_SG12_9 = mm_vbbl.VBBL_SG12_9 else: out = _import_compiled_VBBL() _vbbl_wrapped = out[0] _vbbl_binary_mag_dark = out[1] - _vbbl_SG12_5 = out[2] - _vbbl_binary_mag_0 = out[3] - _vbbl_SG12_9 = out[4] + _vbbl_binary_mag_finite = out[2] + _vbbl_binary_mag_point = out[3] + _vbbl_binary_mag_point_shear = out[4] + _vbbl_SG12_5 = out[5] + _vbbl_SG12_9 = out[6] + if not _vbbl_wrapped: _solver = 'numpy' diff --git a/source/MulensModel/binarylenswithshear.py b/source/MulensModel/binarylenswithshear.py index 7b7ef1d22..19fb79f56 100644 --- a/source/MulensModel/binarylenswithshear.py +++ b/source/MulensModel/binarylenswithshear.py @@ -3,8 +3,8 @@ from math import sqrt from MulensModel.binarylens import BinaryLens -from MulensModel.binarylensimports import (_vbbl_wrapped, - _vbbl_binary_mag_0, _vbbl_SG12_9) +from MulensModel.binarylensimports import ( + _vbbl_wrapped, _vbbl_binary_mag_point_shear, _vbbl_SG12_9) import MulensModel as mm @@ -568,7 +568,7 @@ def point_source_magnification(self, source_x, source_y, vbbl_on=True): x = float(source_x) y = float(source_y) - magnification = _vbbl_binary_mag_0( + magnification = _vbbl_binary_mag_point_shear( s, q, x, y, self.convergence_K, self.shear_G.real, self.shear_G.imag) else: diff --git a/source/MulensModel/tests/test_BinaryLensImports.py b/source/MulensModel/tests/test_BinaryLensImports.py new file mode 100644 index 000000000..f4833ca58 --- /dev/null +++ b/source/MulensModel/tests/test_BinaryLensImports.py @@ -0,0 +1,111 @@ +from numpy.testing import assert_almost_equal +import warnings + +import MulensModel as mm + + +# The test presented check if VBBL code was imported properly and +# they intentionally call private functions from mm.binarylensimports. + +# The values for tests were copied from other tests. + + +def test_VBBL_dark(): + """ + Directly (hence, calling private function) test imported VBBL function: + _vbbl_binary_mag_dark() + """ + if not mm.binarylensimports._vbbl_wrapped: + warnings.warn("VBBL not imported", UserWarning) + return + + out = mm.binarylensimports._vbbl_binary_mag_dark( + 1.35, 0.00578, 0.6598560303179819, -0.05280389642190758, 0.01, + 0.001, 0.0) + assert_almost_equal(out, 1.635980468652538) + + +def test_VBBL_finite(): + """ + Directly (hence, calling private function) test imported VBBL function: + _vbbl_binary_mag_finite() + """ + if not mm.binarylensimports._vbbl_wrapped: + warnings.warn("VBBL not imported", UserWarning) + return + + out = mm.binarylensimports._vbbl_binary_mag_finite( + 0.8, 0.1, 0.01, 0.01, 0.01, 0.001) + assert_almost_equal(out, 18.283392940574107) + + +def test_VBBL_point(): + """ + Directly (hence, calling private function) test imported VBBL function: + _vbbl_binary_mag_point() + """ + if not mm.binarylensimports._vbbl_wrapped: + warnings.warn("VBBL not imported", UserWarning) + return + + out = mm.binarylensimports._vbbl_binary_mag_point(0.8, 0.1, 0.01, 0.01) + assert_almost_equal(out, 18.185448359975954) + + +def test_VBBL_shear(): + """ + Directly (hence, calling private function) test imported VBBL function: + _vbbl_binary_mag_point_shear() + """ + if not mm.binarylensimports._vbbl_wrapped: + warnings.warn("VBBL not imported", UserWarning) + return + + out = mm.binarylensimports._vbbl_binary_mag_point_shear( + 1.1, 0.2, 0.19090909090909103, 0.0, 0.1, 0.1, -0.1) + assert_almost_equal(out, 7.797177952275903) + + +def test_VBBL_SG12_5(): + """ + Directly (hence, calling private function) test imported VBBL function: + _vbbl_SG12_5() + """ + if not mm.binarylensimports._vbbl_wrapped: + warnings.warn("VBBL not imported", UserWarning) + return + + out = mm.binarylensimports._vbbl_SG12_5( + -0.17073496203535826, -0.2960054288981844, 0.2635469112719409, + 0.8257865631242759, 0.8806819361918513, 0.31179330108200043, + 0.029491785441276515, 0.17299410121242445, 0.2364363674161051, + 0.13879110194896355, 0.056792767235166075, -0.003964567527886535) + expected = [ + -0.58000512, -1.6727006, -0.57653394, 0.56004431, -0.5526021, + 0.90541606, -0.16774112, -0.76401705, -0.14860692, -0.04307995] + assert_almost_equal(out, expected) + + +def test_VBBL_SG12_9(): + """ + Directly (hence, calling private function) test imported VBBL function: + _vbbl_SG12_9() + """ + if not mm.binarylensimports._vbbl_wrapped: + warnings.warn("VBBL not imported", UserWarning) + return + + out = mm.binarylensimports._vbbl_SG12_9( + -0.0006162037037037004, 0.018455861111111117, 0.06027712777777768, + -0.08050824074074037, -0.30724011207070717, -0.5752337759963271, + -0.524794951286975, -0.2288742865013774, -0.0350081818181818, 0.0, + 0.0006162037037037004, 0.005693722222222227, 0.006298813888888853, + -0.0035301525925927266, 0.063350366010101, 0.1906118155555555, + 0.09094615694687937, 0.029791669421487667, 0.03019545454545458, + 0.0158) + expected = [ + -1.50992662, 2.47560471, -1.61711916, -0.4453668, -0.75676062, + -0.14298556, 0.36097464, -0.29933664, 0.02381135, -0.74906129, + -3.63807783, 0.90118702, 1.29719497, 0.60196247, -0.63649524, + 0.0565556, -0.0138864, -0.03508702] + assert_almost_equal(out, expected) diff --git a/source/VBBL/VBBinaryLensingLibrary_wrapper.cpp b/source/VBBL/VBBinaryLensingLibrary_wrapper.cpp index d00107852..8e9cd7acd 100644 --- a/source/VBBL/VBBinaryLensingLibrary_wrapper.cpp +++ b/source/VBBL/VBBinaryLensingLibrary_wrapper.cpp @@ -3,15 +3,33 @@ #include "VBBinaryLensingLibrary.h" extern "C" { - double VBBinaryLensing_BinaryMagDark(double a, double q, double y1, double y2, double RSv, double a1, double Tol) { + double VBBinaryLensing_BinaryMagDark(double a, double q, double y1, double y2, double RSv, double tolerance, double a1) { static VBBinaryLensing VBBL; - - return VBBL.BinaryMagDark(a, q, y1, y2, RSv, a1, Tol); + + return VBBL.BinaryMagDark(a, q, y1, y2, RSv, a1, tolerance); + } +} + +extern "C" { + double VBBinaryLensing_BinaryMagFinite(double a, double q, double y1, double y2, double RSv, double tolerance) { + static VBBinaryLensing VBBL; + + VBBL.Tol = tolerance; + + return VBBL.BinaryMag2(a, q, y1, y2, RSv); + } +} + +extern "C" { + double VBBinaryLensing_BinaryMagPoint(double a, double q, double y1, double y2) { + static VBBinaryLensing VBBL; + + return VBBL.BinaryMag0(a, q, y1, y2); } } extern "C" { - double VBBinaryLensing_BinaryMag0(double a,double q,double y1,double y2, double K, double G, double Gi) { + double VBBinaryLensing_BinaryMagPointShear(double a,double q,double y1,double y2, double K, double G, double Gi) { static VBBinaryLensing_shear VBBL; return VBBL.BinaryMag0_shear(a, q, y1, y2, K, G, Gi); @@ -43,15 +61,6 @@ extern "C" double* VBBL_SG12_5(double p0, double p1, return roots; } -extern "C" { - double VBBL_BinaryMag(double a, double q, double y1, double y2) { - static VBBinaryLensing VBBL; - - return VBBL.BinaryMag0(a, q, y1, y2); - } -} - - extern "C" double* VBBL_SG12_9(double p0, double p1, double p2, double p3, double p4, double p5, double p6, double p7, double p8, double p9, double p10, double p11, diff --git a/source/VBBL/vbbl_wrapper.cpp b/source/VBBL/vbbl_wrapper.cpp index 93a8eb1c3..d8594dfd0 100644 --- a/source/VBBL/vbbl_wrapper.cpp +++ b/source/VBBL/vbbl_wrapper.cpp @@ -9,18 +9,44 @@ Based on code written by Przemek Mroz static PyObject * VBBinaryLensing_BinaryMagDark_wrapper(PyObject *self, PyObject *args) { - double a, q, y1, y2, RSv, a1, Tol, mag; + double a, q, y1, y2, RSv, a1, tolerance, mag; static VBBinaryLensing VBBL; - if (!PyArg_ParseTuple(args, "ddddddd", &a, &q, &y1, &y2, &RSv, &a1, &Tol)) return NULL; + if (!PyArg_ParseTuple(args, "ddddddd", &a, &q, &y1, &y2, &RSv, &tolerance, &a1)) return NULL; - mag = VBBL.BinaryMagDark(a, q, y1, y2, RSv, a1, Tol); + mag = VBBL.BinaryMagDark(a, q, y1, y2, RSv, a1, tolerance); return Py_BuildValue("d", mag); } -static PyObject * -VBBinaryLensing_BinaryMag0_wrapper(PyObject *self, PyObject *args) { +static PyObject * +VBBinaryLensing_BinaryMagFinite_wrapper(PyObject *self, PyObject *args) { + double a, q, y1, y2, RSv, tolerance, mag; + static VBBinaryLensing VBBL; + + if (!PyArg_ParseTuple(args, "dddddd", &a, &q, &y1, &y2, &RSv, &tolerance)) return NULL; + + VBBL.Tol = tolerance; + + mag = VBBL.BinaryMag2(a, q, y1, y2, RSv); + + return Py_BuildValue("d", mag); +} + +static PyObject * +VBBinaryLensing_BinaryMagPoint_wrapper(PyObject *self, PyObject *args) { + double a, q, y1, y2, mag; + static VBBinaryLensing VBBL; + + if (!PyArg_ParseTuple(args, "dddd", &a, &q, &y1, &y2)) return NULL; + + mag = VBBL.BinaryMag0(a, q, y1, y2); + + return Py_BuildValue("d", mag); +} + +static PyObject * +VBBinaryLensing_BinaryMagPointShear_wrapper(PyObject *self, PyObject *args) { double a, q, y1, y2, K, G, Gi, mag; static VBBinaryLensing_shear VBBL; @@ -111,7 +137,9 @@ VBBL_SG12_9_wrapper(PyObject *self, PyObject *args) { static PyMethodDef VBBLMethods[] = { {"VBBinaryLensing_BinaryMagDark", VBBinaryLensing_BinaryMagDark_wrapper, METH_VARARGS, "some notes here"}, - {"VBBinaryLensing_BinaryMag0", VBBinaryLensing_BinaryMag0_wrapper, METH_VARARGS, "some notes here"}, + {"VBBinaryLensing_BinaryMagFinite", VBBinaryLensing_BinaryMagFinite_wrapper, METH_VARARGS, "some notes here"}, + {"VBBinaryLensing_BinaryMagPoint", VBBinaryLensing_BinaryMagPoint_wrapper, METH_VARARGS, "some notes here"}, + {"VBBinaryLensing_BinaryMagPointShear", VBBinaryLensing_BinaryMagPointShear_wrapper, METH_VARARGS, "some notes here"}, {"VBBL_SG12_5", VBBL_SG12_5_wrapper, METH_VARARGS, "some notes here"}, {"VBBL_BinaryMag", VBBinaryLensing_BinaryMag_wrapper, METH_VARARGS, "some notes here"}, {"VBBL_SG12_9", VBBL_SG12_9_wrapper, METH_VARARGS, "some notes here"},