12
\$\begingroup\$

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
\$\endgroup\$
7
  • 2
    \$\begingroup\$ @CommandMaster It's fine if it terminates with 100% chance and all possible outputs are within the error margin. \$\endgroup\$ – Bubbler Jan 4 at 5:47
  • 1
    \$\begingroup\$ Your exact formula has an issue; the exact formula requires the semi-axes to follow a>=b>=c, which is not true of the test cases provided. If we are allowed to order the semi-axes to work for the formula provided, that should be specified in the question, or the test cases should be changed to reflect the formula provided. Also, F and E should be the incomplete elliptical integral of the first and second kind, respectively; you have them backwards in the question. \$\endgroup\$ – Zaelin Goodman Jan 4 at 15:36
  • 1
    \$\begingroup\$ Also worth noting that the exact formula provided does not hold for a sphere, as in the first test case. As it stands, the provided exact formula cannot solve the test cases on it's own; it will have to be supplemented with special case handling for spheres, or replaced by a different formula. \$\endgroup\$ – Zaelin Goodman Jan 4 at 16:57
  • 1
    \$\begingroup\$ @ZaelinGoodman Sorry, I missed some details (and misread the E and F) when I was copying the formula from Wikipedia. You're also right about the sphere because \$k\$ would be undefined. I decided to copy the Wikipedia one because the ones on MathWorld were way too complex. (And if I'm understanding it correctly, equations 18 and 23 have the same problem of being undefined for a sphere.) \$\endgroup\$ – Bubbler Jan 4 at 23:13
  • 1
    \$\begingroup\$ @Bubbler no worries, that's a lot to transcribe, and Wikipedia did you dirty by not mentioning the spherical edge case! I just noticed those as I was developing my powershell implementation of that function, figured it'd be helpful for anybody else who decides to go that route over the approximation :) Both 18 and 23 fail for spheres; however, the equations laid out in section 1.2 of This paper from the university of Auckland, New Zealand work for spheres, but the expansion in 1.2.1 converges slowly for highly eccentric ellipsoids \$\endgroup\$ – Zaelin Goodman Jan 5 at 0:38
7
\$\begingroup\$

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!

\$\endgroup\$
4
  • 1
    \$\begingroup\$ Looks like you can do (b*c)**p, etc \$\endgroup\$ – xnor Jan 4 at 5:53
  • 1
    \$\begingroup\$ You can replace (b*c)**p+(c*a)**p with c**p*(b**p+a**p) which saves a byte. \$\endgroup\$ – Varun Vejalla Jan 4 at 17:45
  • \$\begingroup\$ I think your try it online link points to a different answer. \$\endgroup\$ – Razetime Jan 6 at 6:11
  • 1
    \$\begingroup\$ @Razetime Thank you for reporting this. Now fixed. \$\endgroup\$ – Arnauld Jan 6 at 10:35
5
\$\begingroup\$

Wolfram Language (Mathematica), 27 bytes

Shamelessly copying from the 2D solution by @J42161217:

SurfaceArea[#~Ellipsoid~#]&

Try it online!

\$\endgroup\$
4
\$\begingroup\$

Ruby 2.7, 86 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.

\$\endgroup\$
1
  • \$\begingroup\$ It's maybe a bit less interesting but it turns out that porting Arnauld's answer saves 6 bytes. \$\endgroup\$ – Dingus Jan 4 at 21:55
2
\$\begingroup\$

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)
}
\$\endgroup\$
2
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.