Ellipsoid surface area

Question

Related: Ellipse circumference

Introduction

An ellipsoid (Wikipedia / MathWorld) is a 3D object analogous to an ellipse on 2D. Its shape is defined by three principal semi-axes \$a,b,c\$:

$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1 $$

Just like an ellipse, the volume of an ellipsoid is easy, but its surface area does not have an elementary formula. Even Ramanujan won't save you here.

The basic formula is given as the following:

$$ S = 2\pi c^2 + \frac{2\pi ab}{\sin\varphi} \left( E(\varphi,k) \sin^2\varphi + F(\varphi,k) \cos^2\varphi \right) \\ \text{where }\cos\varphi = \frac{c}{a},\quad k^2 = \frac{a^2(b^2-c^2)}{b^2(a^2-c^2)},\quad a \ge b \ge c $$

\$F\$ and \$E\$ are incomplete elliptic integral of the first kind and second kind respectively. Note that the formula does not work for a sphere.

A good approximation can be found on this archived page, where Knud Thomsen developed a symmetrical formula of

$$ S \approx 4\pi \left(\frac{a^p b^p + b^p c^p + c^p a^p} {3 - k\left(1-27abc/\left(a+b+c\right)^3\right)}\right)^{\frac{1}{p}} $$

with empirical values of \$p=\frac{\ln 2}{\ln (\pi/2)}\$ and \$k=3/32\$.

Challenge

Given the three principal semi-axes \$a,b,c\$ of an ellipsoid, compute its surface area.

All three input values are guaranteed to be positive, and you can use any reasonable representation of a real number for input. Also, you may assume the three values are given in a certain order (increasing or decreasing).

The result must be within 0.1% (=10^-3) relative error for the given test cases. You can go for the exact formula (if your language has the necessary built-ins) or Thomsen's approximation, or you can go for numerical integration (extra brownie points if you succeed in this way).

Test cases

The true answer was calculated by feeding the corresponding ellipsoid equation into WolframAlpha.

a  b  c  => answer
------------------
1  1  1  => 12.5664
1  1  2  => 21.4784
1  2  2  => 34.6875
1  1  10 => 99.151
1  2  3  => 48.8821
1  10 10 => 647.22
1  3  10 => 212.00

Answer 1

JavaScript (ES7),  86 83  82 bytes

Saved 3 bytes thanks to @xnor
Saved 1 byte thanks to @VarunVejalla

An approximation of Thomsen's approximation.

(a,b,c)=>(519*((a*b)**(p=1.535)+c**p*(b**p+a**p))/(31+27*a*b*c/(a+b+c)**3))**(1/p)

Try it online!

Answer 2

Wolfram Language (Mathematica), 27 bytes

Shamelessly copying from the 2D solution by @J42161217:

SurfaceArea[#~Ellipsoid~#]&

Try it online!

Answer 3

Ruby 2.7, ⁸⁶ 83 bytes

Saved three bytes, thanks to Dingus!

->a{(519*a.zip(a.rotate).sum{(_1*_2)**1.535}/(31+27*eval(a*?*)/a.sum**3))**0.65147}

Try it online!

Expects an array of \$a,b,c\$. TIO uses an older version of Ruby, whereas in Ruby 2.7, we've numbered parameters, which saves three bytes.

Answer 4

PowerShell, 325 bytes

Takes input as a descending list of semi-axis lengths; +5 bytes if this is not allowed, or +0 bytes if this is not allowed but we need only handle ascending-sorted lists as in the test cases. Spheroid handling takes up 26 bytes of this solution. This solution uses the formula for the surface area of an ellipsoid in the question, which does not hold for spheroids (a=b=c) and requires that a>=b>=c.

The integral calculated by this program is very rough to avoid taking too long to output; can easily improve the accuracy of results by reducing $d, and especially by calculating the area of the trapezoid between $x and $x+$d, rather than the rectangle of width $d at $x

$a,$b,$c=$args;$u,$p,$z={param($t,$g)for($d=1e-4;$x-lt$t){$i+=(&$g ($x+=$d))*$d}$i},($m=[Math])::PI,$c*$c;if($a-eq$c){4*$p*$z}else{$e=$c/$a;$h=$m::Acos($e);$s=$m::Sin($h);$j={$m::Sqrt(1-($a*$a*($b*$b-$z))/($b*$b*($a*$a-$z))*($v=$m::Sin($args[0]))*$v)};2*$p*$z+2*$p*$a*$b/$s*((&$u $h $j)*$s*$s+(&$u $h {1/(&$j @args)})*$e*$e)}

Try it online!

Ungolfed and Explanation

#Set a, b, and c to the elements in the input array
$a,$b,$c=$args;
#Declare a function for finding integrals of a provided function, for 0 to $to
$integralFunction={
    #take parameters $to (the value to calculate the integral to, from 0)
    #and $function (the function to calculate the integral of)
    param($to,$function)
    #set $dx to .0001 and loop until $x >= $to
    for($dx=1e-4;$x-lt$to){
        #increase x by dx
        $x+=$dx
        #calculate the area of the right-side rectangular sliver for this step (f(x)dx)
        $integral+=(&$function $x)*$dx
    }
    #return the integral of the function
    $integral
}
#pi!
$pi=[Math]::PI
#if this is a spheroid (due to the assumption a>=b>=c, a=c follows that a=b=c)
if($a-eq$c){
    #return the surface area of the sphere
    4*$pi*$c*$c
#if this is an ellipsoid
}else{
    #cosine of Phi equals the shortest side over the longest side
    $cosPhi=$c/$a;
    #we need Phi for the integral, so calculate the arccos of cosPhi
    $Phi=[Math]::Acos($cosPhi);
    #we need the sin of Phi for the final calculation
    $sinPhi=[Math]::Sin($Phi);
    #set k^2 to this value, because the formula says so
    #we don't need k, as the elliptical integrals use k^2
    $kSquared= ($a*$a*($b*$b-$c*$c))/($b*$b*($a*$a-$c*$c))
    #since both incomplete elliptical integrals have the term 
    #sqrt(1-k^2*sin^2(theta)), we assign this calculation it's own function
    $ellipticalIntegralFunctionBase={
        param($theta)
        [Math]::Sqrt(1-$kSquared*($sinTheta=[Math]::Sin($theta))*$sinTheta)
    };
    #the first elliptical integral (F) is the integral from 0 to Phi of 
    #(1/sqrt(1-k^2*sin^2(theta))) with respect to theta
    $FirstEllipticalIntegral = (&$integralFunction -to $Phi -function {param($theta) 1/(&$ellipticalIntegralFunctionBase -theta $theta)})
    #the second elliptical integral (E) is the integral from 0 to Phi of
    #sqrt(1-k^2*sin^2(theta) with respect to theta
    $SecondEllipticalIntegral = (&$integralFunction -to $Phi -function $ellipticalIntegralFunctionBase)
    #we have all of our intermediate calculations done, so now we plug those numbers into
    #(2*pi*c^2)+((2*pi*a*b)/sin(Phi))*(E*sin^2(Phi)+F*cos^2(Phi))
    #we return the result of this calculation implicitly
    (2*$pi*$c*$c)+((2*$pi*$a*$b)/$sinPhi)*($SecondEllipticalIntegral*$sinPhi*$sinPhi+$FirstEllipticalIntegral*$cosPhi*$cosPhi)
}

Answer 5

BQN, 50 bytes

{(519×(+´(⊢×1⊸⌽)𝕩⋆p)÷31+(27××´𝕩)÷3⋆˜+´𝕩)⋆÷p←1.535}

Try it!

An anonymous function which uses Arnauld's formula.

Input given as a single array.