I would like to calculate and plot the points of a 2D Cassinian curve in a Cartesian plane.
The Cassinian curves are those curves for which the product of the distances of the points belonging to the curve from 2 other points, called foci, is constant:
$$r_1 = |p-a|$$ $$r_2 = |p-b|$$ $$r_1 r_2 = k$$
The curves are shown in the figure. The outer curves have a bigger value of k.
My idea, given one curve hence one $k$ is to iterate using from $0$ to $2\pi$ and calculate at each step the x-value(s) and y-value(s) but I do not know which calculations I should perform (particularly when the result is 2 closed non-intersecting curves).
I could find something similar but for the ellipsis here

Let $O$ be the center of the cassinian, $S$ and $T$ its foci, $P$ a point on the curve, and set:
$$a=OS=OT$$ $$r=OP$$ $$\theta=\angle SOP$$
From the cosine rule applied to triangle $SOP$ we get: $$PS^2=a^2+r^2-2ar\cos\theta$$
and from the same rule applied to triangle $TOP$ we get $$PT^2=a^2+r^2+2ar\cos\theta$$
Multiplying together those two equalities, and remembering that $PS\cdot PT=k$, we obtain the equation: $$ k^2=(a^2+r^2)^2-4a^2r^2\cos^2\theta, $$ which can be solved for $r^2$ to yield the polar equation of the curve: $$ r^2=a^2\left(\cos2\theta\pm\sqrt{{k^2\over a^4}-\sin^2 2\theta} \right). $$
If $k/a^2\ge1$ then one value of $r^2$ is negative and must be discarded, so you get a single value of $r$ for each $\theta$.
If $k/a^2<1$ both values are positive but they exist only for $|\sin2\theta|\le k/a^2$: in this case you have two separate branches.
Once you have computed $r$ you can the find the coordinates of point $P$ as usual: $$ x=r\cos\theta,\quad y=r\sin\theta. $$