In paper Shadows of rotating black holes approximated by Durer-Pascal limacons, the author tries to find an approximation to the equations
$$ x=\frac{r ∆+r Q^2 −M(r^2 −a^2)}{a(r −M)\sin θ} $$ and $$ y^2=\frac{4r^2∆}{(r −M)^2}−(x+a\sinθ)^2 $$ where $∆ = r ^2 − 2Mr + a^ 2 + Q^ 2$
Note: $r$ does not represent the radius here, i.e $r\neq\sqrt{x^{2}+y^{2}}$.
The author shows that the above equations can be approximated using limacon equation $$ (x ^2 +y^ 2 −Ax)^ 2 =B ^2 (x ^2 +y^ 2 ) $$ for $Q=0$ and $\theta=\frac{\pi}{2}$.
Then it was stated (without proof) that for different values of $\theta$ the "paramters A and B are given by the following simple Fourier series"
$$ A_a(θ) = A_a \sin θ +0.2 \,a \sin^ 3 θ \cos ^2 θ $$ $$ B_a(θ) = B_a +0.23M\left( 1− \sqrt{ 1-\frac{ a^4}{ M^4}}\right) \cos^4 θ $$
My question is how did the author use the fourier series to arrive at the above equations for the paramters A and B?
The statement of the original paper, which is admittedly rather vague, is that for $Q=0$, the orbit of $x(\theta)$ and $y(\theta)$ as given in your question, can be approximated by the limacon $(x^2+y^2-A x)^2 = B^2(x^2+y^2)$. This approximation is not analytical, as far as I can see. Rather, to my mind, the author uses Table 11 in the original paper as data to determine $A$ and $B$. If the orbit was a real limacon, of course $A$ and $B$ would be constant. Here, the author fits a Fourier series for both $A$ and $B$ to the available data, numerically. Additional evidence of this procedure is the appearance of numerical (i.e. inexact) quantities '$0.2$' and '$0.23$' in the expansions for $A$ and $B$: if these were derived analytically, from the original equations, no such numerical approximations would have been necessary.