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$.