Ellipsoid surface area

[Albert Chan]

This version use Carlson symmetric elliptic integrals, to avoid integration.

\(F\!\left(\phi, m\right) = \sin(\phi) R_F\!\left(\cos^{2}\!\left(\phi\right), 1 - m \sin^{2}\!\left(\phi\right), 1\right)\)

\(F\!\left(\phi, m\right) - E\!\left(\phi, m\right) = \frac{1}{3} m \sin^{3}\!\left(\phi\right) R_D\!\left(\cos^{2}\!\left(\phi\right), 1 - m \sin^{2}\!\left(\phi\right), 1\right)\)

We force R arguments to get close, with Duplication theorems.
Assume x,y,z > 0 (x or y can be zero, but not both)
Let λ = √(xy) + √(yz) + √(xz) > 0:

RF(x,y,z) = RF(x+λ, y+λ, z+λ)*2
RD(x,y,z) = RD(x+λ, y+λ, z+λ)*2 + 3/(√(z)*(z+λ))

If arguments close enough, we can stop.

RF(x,x,x) = x^(-1/2)
RD(x,x,x) = x^(-3/2)

From AM-GM inequality, we have (x+y)/2 ≥ √(xy)

λ = √(xy) + √(yz) + √(xz) ≤ (x+y)/2 + (y+z)/2 + (x+z)/2 = x+y+z

Code assumed λ*(1+1e-12) > x+y+z, we consider x,y,z close enough.

Code:

#cas
RFD(x,y,z)  // return RF(x,y,z), RD(x,y,z)
BEGIN
    LOCAL rx,ry,rz,k;
    rx := sqrt(approx(x));
    ry := sqrt(approx(y));
    rz := sqrt(approx(z));
    k := rx*(ry+rz) + ry*rz;
    rx := k*(1+1e-12) > x+y+z; // x,y,z close ?
    IF rx THEN rx:=sqrt(3/k); return rx, 3*rx/k END 
    rx, ry := RFD(x+k, y+k, z+k);
    return rx*2, ry*2 + 3/(rz*(z+k));
END;

ellipsoid_area(a,b,c)
BEGIN
    LOCAL rf, rd, k;
    rf := approx(b/a)^2;
    rd := approx(b/c)^2;
    k := (1-rf)*(1-rd)/3;
    rf, rd := RFD(rf,rd,1);
    return 2.*pi*(b*b+(rf-k*rd)*a*c);
END;    
#end

>>> mp.pretty = True
>>> elliprf(1,2,3), elliprd(1,2,3) # see https://mpmath.org/doc/current/functions/elliptic.html
(0.726945935468908, 0.290460281028991)

CAS> RFD(1,2,3)                     → [0.726945935469, 0.290460281029]
CAS> ellipsoid_area(1,2,3)       → 48.8821463026

Numerical Recipes, Elliptic Integrals: https://core.ac.uk/download/pdf/211378009.pdf
HP-41 Module: Elliptic Functions and Orthogonal Polynomials, ellipsoid area example (page 9), using RG