Simplifying the Willmore energy of an ellipsoid Willmore energy measures how "non-spherical" a smooth surface $S$ is. It is defined by
$$W(S)=\int_SH^2\,dA$$
where $H$ is the mean curvature.
For a torus of revolution with major and minor radii $a$ and $b$ respectively where $a>b$, if we let $p=b/a$ then its Willmore energy is easily shown to be $\frac{\pi^2}{p\sqrt{1-p^2}}$, which attains its minimum at $p=1/\sqrt2$. The (proved) Willmore conjecture states that the torus thus obtained has the minimum energy among all genus-$1$ surfaces of $2\pi^2$. I took Blender out for a ride and produced a render of this "perfect doughnut":

Now I want to calculate the Willmore energy of an ellipsoid. In this genus-$0$ case the extremal results are easy to get: $W(S)\ge4\pi$ and equality is attained iff $S$ is a sphere. But I still want numerical results for the fun of it. Using fundamental forms (and cross-checking with this), I found that for an ellipsoid $E$ with semi-axes $1,a,b$:
$$W(E)=\frac{a^2b^2}4\int_0^{2\pi}\int_0^\pi\frac{(a^2+b^2+(1-(a\cos u)^2-(b\sin u)^2)\sin^2v)^2\sin v}{((ab\cos v)^2+((a\sin u)^2+(b\cos u)^2)\sin^2v)^{5/2}}\,dv\,du$$
which I cannot seem to simplify further.

Does the double integral above have a simpler or even closed form?

Or am I doing it wrong, whereby I would have an easier time using implicit equations as mooted here?
Edit: Using Zhou's ellipsoidal coordinates as suggested by Jean Marie in the comments I have got an expression using only single integrals. For an ellipsoid $E$ with semi-axes $a>b>c>0$ let
$$R_\eta(k)=\int_{c^2}^{b^2}\frac{\eta^k}{\sqrt{(a^2-\eta)(b^2-\eta)(\eta-c^2)\eta}}\,d\eta$$
and
$$R_\zeta(k)=\int_{b^2}^{a^2}\frac{\zeta^k}{\sqrt{(a^2-\zeta)(\zeta-b^2)(\zeta-c^2)\zeta}}\,d\zeta$$
Then
$$W(E)=\frac{(abc)^2}2(R_\eta(-2)R_\zeta(1)+R_\eta(-1)R_\zeta(0)-R_\eta(0)R_\zeta(-1)-R_\eta(1)R_\zeta(-2))$$
Edit 2: The above expression can be simplified to
$$W(E)=\frac{(abc)^2}2(R_\eta(-2)R_\zeta(1)-R_\eta(1)R_\zeta(-2))+\pi$$
 A: For an ellipsoid $E$ with semi-axes $a\ge b\ge c>0$, define
$$A=a^2,B=b^2,C=c^2,\varphi=\cos^{-1}\frac ca$$
$$g=\sqrt{(A-C)B},m=\frac{(B-C)A}{(A-C)B}$$
Then
$$\color{red}{\begin{align}W(E)=\frac\pi{3ABC}&\Big(C(5AB+2AC+2BC)\\
&+2(AB+AC+BC)E(\varphi,m)g\\
&+BC(AB+AC+2BC-A^2)F(\varphi,m)/g\Big)\end{align}}$$
If $b=c$ (the spheroidal case) the above formula simplifies to
$$W(E)=\frac{\pi A\varphi}g+\frac\pi3\left(7+\frac{2B}A\right)$$
while if $a=b=c$ then $E$ is a sphere and $W(E)=4\pi$. For the spheroidal case, whether $a>b$ does not matter – the complex parts cancel out.
All elliptic integrals in this answer use the same argument interpretation as their implementations in Mathematica and mpmath.

This is the result of a long and tedious simplification of the eight integrals in the expression resulting from the use of Zhou's ellipsoidal coordinates. In the equations below, $g=\frac2{\sqrt{(A-C)B}}$, and where $m$ appears in a formula for $R_\zeta(\cdot)$ its value is $1-m=\frac{(A-B)C}{(A-C)B}$ and not $\frac{(B-C)A}{(A-C)B}$.
$$R_\eta(1)=Cg\Pi\left(\frac{B-C}B,m\right)$$
$$R_\zeta(1)=g\left((B-C)\Pi\left(\frac{A-B}{A-C},m\right)+CK(m)\right)$$
$$R_\eta(0)=gK(m)$$
$$R_\zeta(0)=gK(m)$$
$$R_\eta(-1)=\frac gA\left(K(m)+\frac{A-C}CE(m)\right)$$
$$R_\zeta(-1)=\frac gC\left(K(m)+\frac{C-A}AE(m)\right)$$
$$R_\eta(-2)=\frac g{3A^2BC^2}\begin{pmatrix}C(2BC+AB+AC-A^2)K(m)\\+2(A-C)(AB+AC+BC)E(m)\end{pmatrix}$$
$$R_\zeta(-2)=\frac g{3A^2BC^2}\begin{pmatrix}A(2AB+AC+BC-C^2)K(m)\\+2(C-A)(AB+AC+BC)E(m)\end{pmatrix}$$
The terms $R_\eta(-1)R_\zeta(0)-R_\eta(0)R_\zeta(-1)$ in the main expression simplify to just $\frac{2\pi}{ABC}$. The other two terms simplify to a more complicated expression, but still one not involving any integral signs.

Here is Python code to calculate the Willmore energy:
from mpmath import *

def W_ellipsoid(a, b, c):
    c, b, a = sorted([fabs(a), fabs(b), fabs(c)])
    if c == 0: return inf
    if a == c: return 4*pi
    A, B, C = a*a, b*b, c*c
    g = sqrt((A-C)*B)
    phi = acos(c/a)
    m = ((B-C)*A) / ((A-C)*B)
    return pi/(3*A*B*C) * ( C*(5*A*B + 2*A*C + 2*B*C)
                          + 2*(A*B + A*C + B*C)*ellipe(phi, m)*g
                          + B*C*(A*B + A*C + 2*B*C - A*A)*ellipf(phi, m)/g )

def W_spheroid(a, b):
    a, b = fabs(a), fabs(b)
    if a == 0 or b == 0: return inf
    if a == b: return 4*pi
    A, B = a*a, b*b
    g = sqrt((A-B)*B)
    phi = acos(b/a)
    return re(pi*(A*phi/g + (7+2*B/A)/3))


Edit: I got around to symmetrising the elliptic integrals and obtained the following formula for the Willmore energy, where $A,B,C$ keep their meanings from above:
$$W(E)=\frac\pi3\left(3+4\left(\frac1A+\frac1B+\frac1C\right)R_G(AB,CA,BC)-(A+B+C)R_F(AB,CA,BC)\right)$$
