conformal map from disc to funny domain I am working on numerical experiment about the Riemann mapping theorem.
In order to test my program, I would like to have some test function on highly non-symmetric and non convex shape, funnier than this one :

Ideally, something like the standard test domain, like the rabbit:

But I need an exact conformal map like a rational function:
$$\left(\sum_{0}^n a_k z^k\right)/\left(\sum_{0}^m b_k z^k\right), $$
or at least something I can give to python.
Of course for the rabbit like domain, I am not expected to easily get something  perfectly resembling with a rational function, but at least something that is roughly comparable.
Does anyone knows where I can find this can kind of conformal parametrization? I have google many key words without any good result.
Thank in advance!
Edit:
Thanks to Matt E., I get the desired shape, an I am also able to interpolate, but the code for the interpolation using Bezier curves, so the coefficients correspond to splines functions. How can get some trigonometric ones with python? Then I have to hope that the harmonic extension is bijective in the interior of the disc.

 A: If you want a specific shape, I would suggest using the Desmos graphing calculator to approximate points on your desired boundary. Follow this link here to see an example of what I mean: https://www.desmos.com/calculator/pkdtcwkthr
You can upload the picture you want, then eyeball enough points on the boundary to approximate the shape of a bunny. This will give you a finite collection of complex points, say $\{z_n\}_{n=1}^N$. From here, you want to find a parameterization $$p:[0,1] \to \mathbb{C}$$ such that for some collection $\{t_n\}_{n=1}^N$, we have $p(t_n) = z_n$.
Here's my attempt at finding an interpolation method, though it is unfinished. I would encourage you to verify that this method works for interpolating $\{z_1,z_2\}$ or $\{z_1,z_2,z_3\}$.
For $\{z_k\}_{k=1}^n \subset\mathbb{C}$, $\{t_k\}_{k=1}^n \subset (0,1)$ given, set $$p(t) = \sum_{k=1}^n z_k l_k(t), \quad l_k(t) 
= \frac{ \exp\left((2\pi i)\omega_n\displaystyle\prod_{j\neq k} (t - t_j)\right) -1}{\exp\left((2\pi i)\omega_n\displaystyle\prod_{j\neq k} (t_k - t_j)\right) -1},$$
where $\omega_n$ is chosen so that $p(0) = p(1)$.
Notice that $l_k(t_k) = 1$ and $l_k(t_j) = 0$ when $j\neq k$. Hence, $p(t_k) = \sum z_k l(t_k) = z_k$ for $1 \leq k \leq n$. It remains to be shown that one can choose some $\omega_n$ to satisfy $p(0) = p(1)$.
