diff --git a/lenstronomy/LightModel/Profiles/lineprofile.py b/lenstronomy/LightModel/Profiles/lineprofile.py new file mode 100644 index 000000000..5b49d369a --- /dev/null +++ b/lenstronomy/LightModel/Profiles/lineprofile.py @@ -0,0 +1,69 @@ +import numpy as np + + +__all__ = ["LineProfile"] + + +class LineProfile(object): + """Horizontal line segment class. + + The line extends `length` arcseconds from + (`start_x`, `start_y`) at an angle `angle` degrees to the horizontal. Line `width` + is centered in the perpendicular direction, e.g. a profile with 1 arcsecond width + and `angle=0` will span -0.5 to 0.5 in the y-direction. Surface brightness is + constant and given by `amp`. + """ + + param_names = ["amp", "angle", "length", "width", "start_x", "start_y"] + lower_limit_default = { + "amp": 0, + "angle": -180, + "length": 0.01, + "width": 0.01, + "start_x": -100, + "start_y": -100, + } + upper_limit_default = { + "amp": 10, + "angle": 180, + "length": 10, + "width": 5, + "start_x": 100, + "start_y": 100, + } + + def __init__(self): + pass + + def function(self, x, y, amp, angle, length, width, start_x=0, start_y=0): + """Surface brightness per angular unit. + + :param x: x-coordinate on sky + :param y: y-coordinate on sky + :param amp: constant surface brightness of line + :param angle: angle of line to the horizontal (degrees) + :param length: length of line (arcseconds) + :param width: width of line (arcseconds), line width extends symmetrically + :param start_x: ra coordinate of start of line + :param start_y: dec-coordinate of start of line + :return: surface brightness, raise as definition is not defined + """ + ang = -np.deg2rad(angle) + x_ = np.cos(ang) * (x - start_x) + np.sin(ang) * (y - start_y) + y_ = np.cos(ang) * (y - start_y) - np.sin(ang) * (x - start_x) + flux = np.zeros_like(x_) + flux[(x_ >= 0) * (x_ <= length) * (abs(y_) <= width / 2)] = amp + return flux + + def total_flux(self, amp, angle, length, width, start_x=0, start_y=0): + """Integrated flux of the profile. + + :param amp: constant surface brightness of line + :param angle: angle of line to the horizontal (degrees) + :param length: length of line (arcseconds) + :param width: width of line (arcseconds), line width extends symmetrically + :param start_x: ra coordinate of start of line + :param start_y: dec-coordinate of start of line + :return: total flux + """ + return amp * length * width diff --git a/lenstronomy/LightModel/light_model_base.py b/lenstronomy/LightModel/light_model_base.py index 0582557b2..aaccaf213 100644 --- a/lenstronomy/LightModel/light_model_base.py +++ b/lenstronomy/LightModel/light_model_base.py @@ -41,6 +41,7 @@ "SLIT_STARLETS_GEN2", "LINEAR", "LINEAR_ELLIPSE", + "LINE_PROFILE", ] @@ -193,6 +194,10 @@ def __init__(self, light_model_list, smoothing=0.001, sersic_major_axis=None): from lenstronomy.LightModel.Profiles.linear import LinearEllipse self.func_list.append(LinearEllipse()) + elif profile_type == "LINE_PROFILE": + from lenstronomy.LightModel.Profiles.lineprofile import LineProfile + + self.func_list.append(LineProfile()) else: raise ValueError( "No light model of type %s found! Supported are the following models: %s" @@ -286,6 +291,7 @@ def total_flux(self, kwargs_list, norm=False, k=None): "GAUSSIAN_ELLIPSE", "MULTI_GAUSSIAN", "MULTI_GAUSSIAN_ELLIPSE", + "LINE_PROFILE", ]: kwargs_new = kwargs_list_standard[i].copy() if norm is True: diff --git a/lenstronomy/LightModel/linear_basis.py b/lenstronomy/LightModel/linear_basis.py index 2a3eb79fc..bee589ed4 100644 --- a/lenstronomy/LightModel/linear_basis.py +++ b/lenstronomy/LightModel/linear_basis.py @@ -77,6 +77,7 @@ def functions_split(self, x, y, kwargs_list, k=None): "ELLIPSOID", "LINEAR", "LINEAR_ELLIPSE", + "LINE_PROFILE", ]: kwargs_new = kwargs_list[i].copy() new = {"amp": 1} @@ -156,6 +157,7 @@ def num_param_linear_list(self, kwargs_list): "ELLIPSOID", "LINEAR", "LINEAR_ELLIPSE", + "LINE_PROFILE", ]: n_list += [1] elif model in ["MULTI_GAUSSIAN", "MULTI_GAUSSIAN_ELLIPSE"]: @@ -215,6 +217,7 @@ def update_linear(self, param, i, kwargs_list): "ELLIPSOID", "LINEAR", "LINEAR_ELLIPSE", + "LINE_PROFILE", ]: kwargs_list[k]["amp"] = param[i] i += 1 diff --git a/test/test_LightModel/test_Profiles/test_lineprofile.py b/test/test_LightModel/test_Profiles/test_lineprofile.py new file mode 100644 index 000000000..91a2ddd74 --- /dev/null +++ b/test/test_LightModel/test_Profiles/test_lineprofile.py @@ -0,0 +1,35 @@ +from lenstronomy.LightModel.Profiles.lineprofile import LineProfile +import numpy as np +import numpy.testing as npt + + +class TestLineProfile(object): + def setup_method(self): + self.lineprofile = LineProfile() + + def test_function(self): + amp = 1 + length = 1 + width = 0.01 + angle = 57 + x = np.array([0, 1, 0.5 * np.cos(np.deg2rad(-angle))]) + y = np.array([0, 0, 0.5 * np.sin(np.deg2rad(-angle))]) + x_single = 0.2 * np.cos(np.deg2rad(-angle)) + y_single = 0.2 * np.sin(np.deg2rad(-angle)) + flux_true = np.array([amp, 0, amp]) + flux = self.lineprofile.function(x, y, amp, angle, length, width) + single_flux_true = amp + single_flux = self.lineprofile.function( + x_single, y_single, amp, angle, length, width + ) + npt.assert_equal(flux_true, flux) + npt.assert_equal(single_flux_true, single_flux) + + def test_total_flux(self): + amp = 1 + length = 1 + width = 0.01 + angle = 57 + total_flux = self.lineprofile.total_flux(amp, angle, length, width) + total_flux_true = length * width * amp + npt.assert_equal(total_flux_true, total_flux) diff --git a/test/test_LightModel/test_light_model.py b/test/test_LightModel/test_light_model.py index 6eb4454ba..764bf5702 100644 --- a/test/test_LightModel/test_light_model.py +++ b/test/test_LightModel/test_light_model.py @@ -31,6 +31,7 @@ def setup_method(self): "INTERPOL", "SHAPELETS_POLAR_EXP", "ELLIPSOID", + "LINE_PROFILE", ] phi_G, q = 0.5, 0.8 e1, e2 = param_util.phi_q2_ellipticity(phi_G, q) @@ -129,6 +130,14 @@ def setup_method(self): "center_x": 0, "center_y": 0, }, # 'ELLIPSOID' + { + "amp": 1, + "length": 1.0, + "width": 0.01, + "angle": 57, + "start_x": 0, + "start_y": 0, + }, # 'LINE_PROFILE' ] self.LightModel = LightModel( @@ -176,7 +185,7 @@ def test_param_name_list_latex(self): def test_num_param_linear(self): num = self.LightModel.num_param_linear(self.kwargs, list_return=False) - assert num == 19 + assert num == 20 num_list = self.LightModel.num_param_linear(self.kwargs, list_return=True) assert num_list[0] == 1 @@ -199,6 +208,7 @@ def test_total_flux(self): "GAUSSIAN_ELLIPSE", "MULTI_GAUSSIAN", "MULTI_GAUSSIAN_ELLIPSE", + "LINE_PROFILE", ] kwargs_list = [ { @@ -247,6 +257,14 @@ def test_total_flux(self): "center_x": 0, "center_y": 0, }, # 'MULTI_GAUSSIAN_ELLIPSE' + { + "amp": 1, + "length": 1.0, + "width": 0.01, + "angle": 57, + "start_x": 0, + "start_y": 0, + }, # 'LINE_PROFILE' ] lightModel = LightModel(light_model_list=light_model_list) total_flux_list = lightModel.total_flux(kwargs_list)