How spherical is my ellipsoid?

Define the (unnormalised) Willmore energy of a surface as the integral of squared mean curvature over it: $$W=\int_SH^2\,dA$$ For surfaces topologically equivalent to a sphere $$\W\ge4\pi\$$, and $$\W=4\pi\$$ iff it is actually a sphere. Thus the energy quantifies how spherical a surface is.

In September 2020 I asked on MathsSE how to simplify the Willmore energy of an ellipsoid with semi-major axes $$\a\ge b\ge c\ge0\$$. Two days later I had done it myself:

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=\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}}

Here $$\F\$$ and $$\E\$$ are the elliptic integrals of the first and second kinds respectively. (Note that as with all my elliptic integral/function answers on MathsSE the argument convention is as in Mathematica and mpmath, where $$\m\$$ is the parameter.) Later on I derived a much cleaner and order-agnostic formula in terms of Carlson's symmetric integrals$$\A,B,C\$$ keep their meanings: $$W=\frac\pi3\left(3+4\left(\frac1A+\frac1B+\frac1C\right)R_G(AB,CA,BC)-(A+B+C)R_F(AB,CA,BC)\right)$$ I was inspired to write this question after writing answers to the ellipse perimeter and ellipsoid surface area questions that use the symmetric integrals. In the former case it actually saves one byte over an answer using the classic formula.

Given an ellipsoid's semi-axes $$\a,b,c\$$, which you may assume are sorted and positive, output the ellipsoid's Willmore energy with a relative error of at most $$\10^{-3}\$$. You may use either formula above or something else entirely like explicitly integrating the squared mean curvature.

This is ; fewest bytes wins.

Test cases

(1, 1, 1) 12.5663706143592
(2, 1, 1) 15.4516066443266
(2, 2, 1) 16.9023119660317
(3, 2, 1) 21.2244261324396
(2, 2, 2) 12.5663706143592
(6, 1, 1) 34.2162831541838
(6, 5, 1) 70.4793621781325
(13, 9, 6) 15.9643343585267


JavaScript (ES6), 195 bytes

Quite obviously not the right tool for the job, at least without any specialized library. The only math builtins used here are $$\\sin\$$, $$\\cos\$$ and $$\\operatorname{hypot}\$$.

The precision can be improved by using a smaller increment (the $$\i\$$ in the code). The version below uses $$\i=10^{-2}\$$ so that it doesn't time-out on TIO. But as a result, the relative error is $$\\approx10^{-1}\$$.

with(Math)f=(a,b,c)=>{i=1e-2;for(F=G=p=0;S=sin(p+=i),P=S*i*i/12,p<PI;)for(q=0;q<2*PI;G+=P*Q)F+=P/(Q=hypot(a*b*S*cos(q+=i),c*a*S*sin(q),b*c*cos(p)));return(4/a/a+4/b/b+4/c/c)*G-(a*a+b*b+c*c)*F+PI}


Try it online!

Formulas

This is derived from the second formula given in the challenge and the following formulas from this page for Carlson's symmetric integrals:

$$R_F(x,y,z)=\frac{1}{4\pi}\int_{0}^{2\pi}\int_{0}^{\pi}\frac{\sin\theta\:d\theta\:d\phi}{(x\sin^2\theta\cos^2\phi+y\sin^2\theta\sin^2\phi+z\cos^2\theta)^{1/2}}$$ $$R_G(x,y,z)=\frac{1}{4\pi}\int_{0}^{2\pi}\int_{0}^{\pi}(x\sin^2\theta\cos^2\phi+y\sin^2\theta\sin^2\phi+z\cos^2\theta)^{1/2}\sin\theta\:d\theta\:d\phi$$

Wolfram Language (Mathematica), 63 bytes

-1 byte and fixed a bug thanks to @att.

Pi/3(3.-Tr[s=#^2]CarlsonRF@@(t=##&@@s/s)+4Tr[1/s]CarlsonRG@@t)&


CarlsonRG and CarlsonRF are new functions in Mathematica 12.3, so this does not work on TIO.

• 1##&@@#/#->##&@@#/#
– att
Commented Feb 15, 2023 at 7:31
• I think you need to square the argument (since $A=a^2$ etc.)
– att
Commented Feb 15, 2023 at 7:33