Consider $$f(x,y)=\sum_{i=1}^n\dfrac{\sin(\omega_ix)}{\sin(\omega_iy)}r_i$$ where $n,\omega_i,r_i>0$ are known parameters. Restrict to domains where $ \sin(\omega_iy)\neq 0$, and by symmetry $x,y>0$.
There is in general no closed-form solution. But assume I know one (an approximation obtained with Newton-Raphson for example) denoted by $(x_0,y_0)$.
Using the implicit function theorem, I can calculate the Taylor expansion of $x\mapsto f(x,\varphi(x))$:
to the first order: $$f(x,\varphi(x))=f(x_0,\varphi(x_0)+(x-x_0)(\partial_1 f(x_0,\varphi(x_0))+\partial_2 f(x_0,\varphi(x_0)) \varphi'(x_0)) +o((x-x_0)^2)$$ Since by definition of $\varphi$ this is equal to $0$ and $f(x_0,y_0)=0$, the only unknown is$\varphi'(x_0)$ (recall that $\varphi(x_0)=y_0$ which is known). We can solve this first order linear equation and get $\varphi'(x_0)$.
then to the second order: $$f(x,\varphi(x))=(x-x_0)(...)+(x-x_0)^2(...)+o(x^3)$$
Here, using the value of $\varphi'(x_0)$ for the previous step, we get the value of the only unknown $\varphi''(x_0)$.
- etc. up to the desired order $N$.
It is then possible to reconstruct the Taylor expansion of $\varphi$ in the neighbourhood of $x_0$ using the identified $\varphi^{(i)}(x_0)$.
The obtained result, when choosing values for $n,w_i,r_i$ is represented below: blue curves are numerical approximation of roots of $f$ in the $(x,y)$ plane, red dot is $(x_0,y_0)$ and red curve is the graph of $(x,\varphi(x))$.
Note that the blue lines in fact correspond to the roots of $f(x,y)\Pi_{i=1}^n \sin(\omega_iy)$ to get rid of the singularities which are not of interest here.
As you can see, the approximation is OK but very limited: the radius of convergence is much smaller than the blue curve delimited by vertical tangents. Here the order $N$ is 15 and getting above that results in too much computation time.
The question is, would there be a way to get an analytical (with explicit dependence on $n,\omega_i,r_i$) expression corresponding to a better approximation? Possibilities could include Chebyshev series, Padé approximations, or why not a nonlinear change of coordinates...). Ideally, I would get an approximation of the whole closed curve---but I'm not sure if it is possible. I'm open to any comments.
Possible hints
Reducing $f(x,y)$ to a common denominator yields $$f(x,y)=\dfrac{\sum_{i=1}^n \big(\prod_{j\neq i} \sin(\omega_j y) \big)\sin(\omega_i x)}{\prod_{i=1}^n \sin(\omega_j y)}$$
So for a fixed $y$, $f$ is a linear combination of $\sin(\omega_1x),\dots,\sin(\omega_n x)$.
Also, with $X_i=\sin(\omega_i x)$ and $Y_i=\sin(\omega_i y)$ the numerator writes $$r_1 Y_2\dots Y_n X_1 + r_2 Y_1Y_3\dots Y_n X_2 + \dots + r_n Y_1\dots Y_{n-1} X_n=0$$ so problem can also be seen as solving a polynomial equation.
But I am not sure this helps...