11import numpy as np
22from typing import Tuple
3+ from autogalaxy .convert import multipole_comps_from , multipole_k_m_and_phi_m_from
34
45
56from autogalaxy .ellipse .ellipse .ellipse import Ellipse
@@ -37,6 +38,132 @@ class representing the multipole of an ellispe with, which is used to perform el
3738 self .m = m
3839 self .multipole_comps = multipole_comps
3940
41+ def get_shape_angle (
42+ self ,
43+ ellipse : Ellipse ,
44+ ) -> float :
45+ """
46+ The shape angle is the offset between the angle of the ellipse and the angle of the multipole,
47+ this defines the shape that the multipole takes.
48+
49+ In the case of the m=4 multipole, angles of 0 indicate pure diskiness, angles +- 45
50+ indicate pure boxiness.
51+
52+ Parameters
53+ ----------
54+ ellipse
55+ The base ellipse profile that is perturbed by the multipole.
56+
57+ Returns
58+ -------
59+ The angle between the ellipse and the multipole, in degrees between +- 180/m.
60+ """
61+
62+ angle = (
63+ ellipse .angle ()
64+ - multipole_k_m_and_phi_m_from (self .multipole_comps , self .m )[1 ]
65+ )
66+ if angle < - 180 / self .m :
67+ angle += 360 / self .m
68+ elif angle > 180 / self .m :
69+ angle -= 360 / self .m
70+
71+ return angle
72+
73+ def points_perturbed_from (
74+ self , pixel_scale , points , ellipse : Ellipse , n_i : int = 0
75+ ) -> np .ndarray :
76+ """
77+ Returns the (y,x) coordinates of the input points, which are perturbed by the multipole of the ellipse.
78+
79+ Parameters
80+ ----------
81+ pixel_scale
82+ The pixel scale of the data that the ellipse is fitted to and interpolated over.
83+ points
84+ The (y,x) coordinates of the ellipse that are perturbed by the multipole.
85+ ellipse
86+ The ellipse that is perturbed by the multipole, which is used to compute the angles of the ellipse.
87+
88+ Returns
89+ -------
90+ The (y,x) coordinates of the input points, which are perturbed by the multipole.
91+ """
92+ symmetry = 360 / self .m
93+ k_orig , phi_orig = multipole_k_m_and_phi_m_from (self .multipole_comps , self .m )
94+ comps_adjusted = multipole_comps_from (
95+ k_orig ,
96+ symmetry
97+ - 2 * phi_orig
98+ + (symmetry - (ellipse .angle () - phi_orig )), # Re-align light to match mass
99+ self .m ,
100+ )
101+
102+ # 1) compute cartesian (polar) angle
103+ theta = np .arctan2 (points [:, 0 ], points [:, 1 ]) # <- true polar angle
104+
105+ # 2) multipole in that same frame
106+ delta_theta = self .m * (theta - ellipse .angle_radians ())
107+ radial = comps_adjusted [1 ] * np .cos (delta_theta ) + comps_adjusted [0 ] * np .sin (
108+ delta_theta
109+ )
110+
111+ # 3) perturb along the true radial direction
112+ x = points [:, 1 ] + radial * np .cos (theta )
113+ y = points [:, 0 ] + radial * np .sin (theta )
114+
115+ return np .stack ((y , x ), axis = - 1 )
116+
117+
118+ class EllipseMultipoleScaled (EllipseMultipole ):
119+ def __init__ (
120+ self ,
121+ m = 4 ,
122+ scaled_multipole_comps : Tuple [float , float ] = (0.0 , 0.0 ),
123+ major_axis = 1.0 ,
124+ ):
125+ """
126+ class representing the multipole of an ellipse, which is used to perform ellipse fitting to
127+ 2D data (e.g. an image). This multipole is fit with its strength held relative to an ellipse with a
128+ major_axis of 1, allowing for a set of ellipse multipoles to be fit at different major axes but with
129+ the same scaled strength k/a.
130+
131+ The scaled_multipole_comps (for all ellipses) are converted to a k value, which is then reset to
132+ its `true' value for a multipole at the given major axis value, which is then used to perturb an ellipse
133+ as per the normal `EllipseMultipole' class and below.
134+
135+ The multipole is added to the (y,x) coordinates of an ellipse that are already computed via the `Ellipse` class.
136+
137+ The addition of the multipole is performed as follows:
138+
139+ :math: r_m = \sum_{i=1}^{m} \left( a_i \cos(i(\t heta - \phi)) + b_i \sin(i(\t heta - \phi)) \r ight)
140+ :math: y_m = r_m \sin(\t heta)
141+ :math: x_m = r_m \cos(\t heta)
142+
143+ Where:
144+
145+ m = The order of the multipole.
146+ r = The radial coordinate of the ellipse perturbed by the multipole.
147+ \phi = The angle of the ellipse.
148+ a = The amplitude of the cosine term of the multipole.
149+ b = The amplitude of the sine term of the multipole.
150+ y = The y-coordinate of the ellipse perturbed by the multipole.
151+ x = The x-coordinate of the ellipse perturbed by the multipole.
152+ """
153+
154+ self .scaled_multipole_comps = scaled_multipole_comps
155+ k , phi = multipole_k_m_and_phi_m_from (
156+ multipole_comps = scaled_multipole_comps , m = m
157+ )
158+ k_adjusted = k * major_axis
159+
160+ specific_multipole_comps = multipole_comps_from (k_adjusted , phi , m )
161+
162+ super ().__init__ (m , specific_multipole_comps )
163+
164+ self .specific_multipole_comps = specific_multipole_comps
165+ self .m = m
166+
40167 def points_perturbed_from (
41168 self , pixel_scale , points , ellipse : Ellipse , n_i : int = 0
42169 ) -> np .ndarray :
@@ -56,17 +183,36 @@ def points_perturbed_from(
56183 -------
57184 The (y,x) coordinates of the input points, which are perturbed by the multipole.
58185 """
186+ symmetry = 360 / self .m
187+ k_orig , phi_orig = multipole_k_m_and_phi_m_from (
188+ self .specific_multipole_comps , self .m
189+ )
190+ comps_adjusted = multipole_comps_from (
191+ k_orig ,
192+ symmetry - 2 * phi_orig + (symmetry - (ellipse .angle () - phi_orig )),
193+ self .m ,
194+ )
59195
60- angles = ellipse .angles_from_x0_from (pixel_scale = pixel_scale , n_i = n_i )
196+ # 1) compute cartesian (polar) angle
197+ theta = np .arctan2 (points [:, 0 ], points [:, 1 ]) # <- true polar angle
61198
62- radial = np .add (
63- self .multipole_comps [1 ]
64- * np .cos (self .m * (angles - ellipse .angle_radians ())),
65- self .multipole_comps [0 ]
66- * np .sin (self .m * (angles - ellipse .angle_radians ())),
199+ # 2) multipole in that same frame
200+ delta_theta = self .m * (theta - ellipse .angle_radians ())
201+ radial = comps_adjusted [1 ] * np .cos (delta_theta ) + comps_adjusted [0 ] * np .sin (
202+ delta_theta
67203 )
68204
69- x = points [:, 1 ] + (radial * np .cos (angles ))
70- y = points [:, 0 ] + (radial * np .sin (angles ))
205+ # Old code, delete in fuure but keep for debugging for now:
206+
207+ # radial = np.add(
208+ # self.multipole_comps[1]
209+ # * np.cos(self.m * (angles - ellipse.angle_radians())),
210+ # self.multipole_comps[0]
211+ # * np.sin(self.m * (angles - ellipse.angle_radians())),
212+ # )
213+
214+ # 3) perturb along the true radial direction
215+ x = points [:, 1 ] + radial * np .cos (theta )
216+ y = points [:, 0 ] + radial * np .sin (theta )
71217
72218 return np .stack (arrays = (y , x ), axis = - 1 )
0 commit comments