Inverting a function in the parametrisation of Riemann's minimal surface While writing code to render Riemann's minimal surface – which is now available here as part of my Malibu STL surface generator – the following problem came up. Consider the following function:
$$y(\theta)=\operatorname{Re}\left(\sqrt{\frac{1-e^{2i\theta}}{(\sigma-1)e^{2i\theta}+(\sigma+1)e^{i\theta}}}\right)\qquad\theta\in[0,\pi],\sigma>0$$
This is (after constant factors have been taken out) the $y$-coordinate of a point on the surface corresponding to some complex number $w$, which itself is the image under a conformal map of a point on the unit circle in the upper half-plane, hence the appearance of $e^{i\theta}$. $\sigma$ is a given constant; for all values of $\sigma$, $y$ is strictly increasing up to a point and strictly decreasing after that.
In order for neighbouring patches to match up exactly, right down to the vertices in the output STL file, I needed to invert $y$, i.e. given $y(\theta)$ find the two possibilities for $\theta$. Failing that, I needed to find for a given $\theta_1$ another angle $\theta_2$ such that $y(\theta_1)=y(\theta_2)$. I could not do either analytically and had to resort to numerical methods.
If there was no taking of the real component I would be able to solve this rather easily, but $y$ is the result of using the Weierstrass–Enneper parametrisation and there's no way around it here. The $\operatorname{Re}$ causes SymPy to output a very complicated expression for $y$ that does not use $\operatorname{Re}$ but uses $\sin\theta$ and $\cos\theta$, and I made no headway this way.

Is it possible to solve the "preimage" or "second-preimage" problems for $y$ described above analytically?

Edit: My full parametrisation of a fundamental patch of the Riemann surface is as follows. Let
$$m=\frac1{1+\sigma},m'=1-m=\sigma m,c=\frac{1+\sigma}{\sqrt\sigma}$$
Then
$$x(w)=c\operatorname{Im}(2(m'F(\varphi=\sin^{-1}\sqrt{1-w},m)-E(\varphi,m))+p)$$
$$y(w)=c\operatorname{Re}\left(p=\sqrt{\frac{(\sigma+w)(1-w)}{w(1+\sigma)}}\right)$$
$$z(w)=2\operatorname{Im}(F(\varphi,m))$$
where $w$ lies in the semicircle with diameter $1$ and $\sigma$ in the upper half-plane punctured at $(0,0)$. The conformal map I am using is the Möbius transformation
$$w=\frac{2\sigma v}{(\sigma-1)v+(\sigma+1)}$$
where $v$ lies in the upper half-unit disc. Substituting this map in the $y(w)$ equation and then restricting $v$ to lie on the unit circle ($v=e^{i\theta}$) gives the equation at the top of this question.
The translation vector between unit cells is given by $(x(-\sigma), 0, z(-\sigma))=\left(\frac2{\sqrt\sigma}((1+\sigma)E(m')-K(m')),0,2K(m')\right)$.
 A: I tried a different tack. Let $x=\cos\theta$ so that $x\in[-1,1]$ and $e^{i\theta}=x+i\sqrt{1-x^2}$, then after explicitly determining the real and imaginary parts of the expression under the square root in the formula for $y$ I got the following equivalent representation:
$$2y^2(x)=(1-\sigma)L+\sqrt{2L}\qquad L=\frac{1-x^2}{(x+1)\sigma^2-x+1}$$
Since I only care about finding matching angles that give the same result under $y$, I don't need the square and doubling on the LHS, so from now on define $y(x)=(1-\sigma)L+\sqrt{2L}$. This $y$ is easily inverted, being the composition of two quadratic expressions. Suppose $k=y(x)$ is known, then
$$L=\frac{k(1-\sigma)+1\pm\sqrt{2k(1-\sigma)+1}}{(1-\sigma)^2}$$
$$x=\frac12\left(L(1-\sigma^2)\pm\sqrt{(L(1+\sigma^2)-2)^2-(2L\sigma)^2}\right)$$
Finally $\cos^{-1}x=\theta$.

The previous version of my code in Malibu had assumed that "for all values of $\sigma$, $y$ is strictly increasing up to a point and strictly decreasing after that". This is actually only true for $\sigma\le3$. What is true for all $\sigma$ is that $y(\cos\theta)$ has a local extremum at
$$\theta_0=\cos^{-1}\frac{1-\sigma}{1+\sigma}\qquad\left[y(\cos\theta_0)=\left(\frac2{1+\sigma}\right)^2\right]$$
and this is a maximum for $\sigma\le3$ but a minimum for $\sigma>3$. In the latter case there are two local maxima, one on each side of $\theta_0$:
$$\theta_\pm=\cos^{-1}\frac{1+\sigma\pm\sqrt{3(\sigma-3)(3\sigma-1)}}{4(1-\sigma)}\qquad\left[y(\cos\theta_\pm)=\frac1{2(\sigma-1)}\right]$$
Thus there can be up to four inverses for a given $k$, which corresponds to the patch boundary caving in at $\sigma>3$ ($\sigma=10$ is pictured):

The following table shows how the sign choices produce the different inverses:
$$\begin{array}{c|cccc}
\text{sign choice for }(L,x)&(-,+)&(+,+)&(+,-)&(-,-)\\
\theta\text{ range if }\sigma\le3&[0,\theta_0]& & &[\theta_0,\pi]\\
\theta\text{ range if }\sigma>3&[0,\theta_-]&[\theta_-,\theta_0]&[\theta_0,\theta_+]&[\theta_+,\pi]\\
\end{array}$$
