Uniform distribution on a toroidal surface

[Barnak]

I'm not sure this is the right thread to post my problem :

I'm trying to define a uniform distribution on the toroidal surface associated to a dipolar magnetic field (or electric). More specifically, the surface (in 3D euclidian space) is parametrised as this, using the usual polar coordinates :

$x(\theta, \phi) = \sin^3{\theta} \; \cos{\phi},$
$y(\theta, \phi) = \sin^3{\theta} \; \sin{\phi},$
$z(\theta) = \sin^2{\theta} \; \cos{\theta}.$

The surface element is this :

$dS(\theta, \phi) = \sin^7{\theta} \; d\theta \; d\phi.$

For a simple sphere, we get

$dS_{sphere}(\theta, \phi) = \sin{\theta} \; d\theta \; d\phi = du \; d\phi,$

where $u = \cos \theta$ is the natural variable to define the uniform distribution on the sphere.

In the case of my toroidal surface defined above, the "natural" variable (if I'm not doing a mistake) is really complicated :

$u = \cos{\theta} - \cos^3{\theta} + \tfrac{3}{5} \cos^5{\theta} - \tfrac{1}{7}\cos^7{\theta},$

so $dS(u, \phi) = du \; d\phi$. This is the variable I should use to define an uniform distribution of points on the surface.
However, how should I define the three parametric coordinates $x(u, \phi)$, $y(u, \phi)$, $z(u)$ ? I'm unable to invert the function above to give $\cos \theta = f(u) = ?$

Help please !

I'm using Mathematica to do my calculations.

 

[fresh_42]

$0=-u + \cos \theta \cdot \left( 1 - (\cos^2 \theta) + \frac{3}{5} (\cos^2 \theta)^2 -\frac{1}{7}(\cos^2 \theta)^3\right)$
$= -u + \cos \theta \cdot (1-v+\frac{3}{5}v^2-\frac{1}{7}v^3)$ which can basically be solved for $v=\cos^2 \theta$ and then for $\cos \theta$.