Related Links (Outside this Site)

(2001-10-23)
What is the surface area of an ellipsoid? [spheroid]

There are simple formulas for the surface area of an ellipsoid of revolution,  but when the 3 semiaxes (a, b, c) are distinct, the formula isn't elementary:

The surface area of an ellipsoid of equation   (x/a)²+(y/b)²+(z/c)²=1   is:

The above was originally posted here to provide a correct version of a flawed formula given in the Mathematica 4 documentation  [where "EllipticE" and "EllipticF" are interchanged, as David W. Cantrell first pointed out].  As of 2005, the typo has only been corrected in the Mathematica 5 documentation.

Knud Thomsen  (Denmark, 2004-04-26; e-mail)
The following symmetrical formula seems to give the surface area of a general ellipsoid with a relative error < 1.42%.  Could it be new?   [...]
S   »   2p  ( a^pb^p + a^pc^p + b^pc^p ) ^1/p   where   p = lg(3) = ln(3)/ln(2)

The above (correct) graphic rendition won't match the previously quoted formula if your web browser fails to provide proper legacy support for the Symbol font.  The "bad boys", including some versions of Safari, Opera and Firefox,  display "p" instead pf "π"  and render other kinds of unreliable gibberish in all scientific web pages whose authors make explicit use of the "Symbol font", following the relevant W3C recommendation  (issued in 1997 and universally supported until 2005 or so).  To be safe, you are urged to switch to a "nice guy", like  Internet Explorer  or  Google Chrome  (otherwise, you have to use one of the temporary solutions discussed elsewhere).

Somebody must have thought of this formula before but it seems unpublished.  If you know better, please tell me and I'll post updates here.  Meanwhile, I shall refer to the above as Thomsen's formula  (2004-04-26).

For a very long prolate ellipsoid of revolution, Thomsen's expression is  2^1+1/p / p  times the true surface area  (a relative error of about -1.41544 %).

This formula is also discussed in the next article.

On 2004-05-13, Knud Thomsen (Denmark) wrote:   [edited summary] The expression   S   »   4p  [ ( apbp + apcp + bpcp ) / 3 ] 1/p approximates the true surface area of the ellipsoid with the least relative error (± 1.061 %  worst case) when   p » 1.6075 Best regards, Knud Thomsen

This second expression is exact for a sphere regardless of the value of p,  whereas the previous one is always exact for a degenerate ellipsoid (c = 0).  Both coincide for the value  p = lg 3 = 1.5849625...  first given by Thomsen.

Knud Thomsen (2004-05-14) and David Cantrell (>2004-05-16) independently mentioned the obvious generalization to n-dimensional ellipsoids (with semiaxes a₁, ... a_n )  which may be expressed in terms of the Hölder mean (H) of the n products of  n-1  semiaxes ( H = 0  when two or more semiaxes are zero):

H   =   (a₁a₂ ... a_n ) [ (a₁^-p + a₂^-p + ... + a_n^-p ) / n ] ^1/p   [nonzero semiaxes]
H   =   (a₁a₂ ... a_n-1 ) / n ^1/p     [when one semiaxis, say a_n, is zero]

The hypervolume (V) and hyperarea (S) of an n-dimensional hypersphere are:

V   =   Rⁿ p^n/2 / G(1+n/2)   [R is the hypersphere radius]
S   =   2 R^n-1 p^n/2 / G(n/2)   =   2 H p^n/2 / G(n/2)

Scaling considerations translate into the exact formula for the hypervolume (V) of an ellipsoid (replace Rⁿ  by a₁a₂...a_n ).  The n-1 dimensional case gives half the hyperarea of the [two sided] degenerate n-dimensional ellipsoid, which is thus:

S'   =   2 ( n^1/p H ) p^[n-1]/2 / G(1+[n-1]/2)

The expressions for S (hypersphere) and S' (degenerate hyperellipsoid) are identical functions of H  (generalizing the YNOT formula for the ellipse)  if:

n^1/p   =   Öp  G(n/2+1/2) / G(n/2)

p   =	Log (n)
	Log (A)

A  =  Öp  G(n/2+1/2) / G(n/2)   has different elementary expressions in terms of  double-factorials  (A006882)  depending on the parity of n :

• If n = 2k, then   A   =   p (2k-1)!! / 2^k (k-1)!   =   (p/2) (n-1)!! / (n-2)!!
• If n = 2k+1, then   A   =   2^k k! / (2k-1)!!   =   (n-1)!! / (n-2)!!

The limit of this value of p is 2 for [very] large values of n...

Thomsen and Cantrell both expressed some doubts about the future popularity of such approximate formulas for the hyperareas of hyper-ellipsoids...

On 2004-05-14, David W. Cantrell wrote:   [edited summary] I just now noticed [the recent update(s) to this page].  Congratulations, Knud, on your approximation of 2004-04-26 and your more recent one with   p = 1.6075...   I would prefer   p = 8/5   which makes the latter algebraic and optimal for nearly spherical ellipsoids.  The relative error is still never worse than -1.178 % Best regards, David W. Cantrell

David posted other approximations shortly after that  fast  feedback.

(2004-05-02)   Approximate Surface Area of a Scalene Ellipsoid
Ellipsoids of revolution and flat ones may be considered trivial cases...

An ellipsoid of semiaxes a, b and c is called scalene if these three quantities are distinct  (or when no two of them are known to coincide,  per the inclusive viewpoint which is preferred in a mathematical context).

Let a be the largest semiaxis and c the smallest one.  Let  e = Ö(1-c²/a²).

The surface area (S) of the ellipsoid has a simple expression in 3 special cases:  for an oblate or prolate ellipsoid of revolution, and for a degenerate ellipsoid  (namely, a flat spheroid whose surface consists of the two sides of an ellipse) :

• If a = b, then   S = 2p [ a² + c² atanh(e)/e ]    Oblate ellipsoid (M&M's).
• If b = c, then   S = 2p [ c² + ac arcsin(e)/e ]   Prolate ellipsoid (cigar).
• If c = 0, then   S = 2p ab 

The above symmetrical formula, proposed in 2004 by the Danish geologist Knud Thomsen, is exact for either a sphere  ( a = b = c )  or a flat spheroid.  It's to the ellipsoid area what the YNOT formula is to the ellipse perimeter.

Dropping the simplicity and symmetry of Thomsen's formula, it's easy to devise expressions that are correct in all three of the above trivial cases...  For example,  consider the following expression: (where atanh(x)/x and arcsin(x)/x are defined by continuity to be equal to 1 when x is zero):

S   »   2p  [ a (b-c) + ac arcsin(u)/u + c² atanh(v)/v ]
where   u = Ö(1-b²/a²)   and   v = Ö(1-c²/b²)

This turns out to be a poor kind of  "improvement"  over Thomsen's formula except, of course, when the ellipsoid is a solid of revolution, or nearly so.

The worst discrepancy occurs when  a  is much larger than the other two semiaxes and those lesser axes are in a particular ratio:

c   =   lb       with       l   =   (p/2-1)^1/(p-1)

In that case,  Thomsen's value exceed this expression by a factor of

(1+l^p )^1-1/p   »   1.07578

As Thomsen's expression is a -1½ % underestimate, the above is -9 % off.  There's a unacceptable price to pay for breaking up the natural symmetry of the approximation,  as asymmetrical error terms no longer cancel out.

This remark also applies to another approximation (< 2.1 %) proposed in 2003 by Dr. Andreas Dieckmann (of Bonn University) in terms of   r = arcsin(e)/e :

This gives the correct area  S = 2pc (c+ar)  for prolate spheroids.

On 2004-05-17, we received the first attempt at optimizing a symmetrical formula by Thomsen, who investigated the following expression, featuring a second parameter (k) generalizing his earlier formula (the case k = 0):

S   »   4p  [ ( a^pb^p + a^pc^p + b^pc^p ) / ( 3 - k { 1 - 27abc / (a+b+c)³ } ) ] ^1/p

The formula is designed to be correct in the case of a sphere for all values of p and k.  Knud Thomsen selects

p   =   ln(2) / ln(p/2)

He then claims optimal results when k is around 0.0942 and suggests the convergent 3/32 (0.9375) which is good enough to yield a relative error between  -0.204 %  and  +0.187 %  [the next convergents would be 5/53 and 8/85].  This represents an improvement of one order of magnitude over his original formula (k=0).

Earlier work of Achim Flammenkamp (of Universität Bielefeld, Germany) along similar lines was brought to our attention on 2004-06-14 by his colleague Torsten Sillke:  Building on an article of Klamkin, Flammenkamp investigated several approximations to the surface area of an ellipsoid around 1990, including the following two expressions.  The first one has only a 10% accuracy, whereas a worst relative error of 2.09% is claimed for the second formula.

S   »   2p  [ ( ab + ac + bc ) - 3 abc / (a+b+c) ]
 
S   »   p  [ (1-1/Ö3) ( ab + ac + bc )  +  (1+1/Ö3) ( a²b² + a²c² + b²c² )^½ ]

From 1994 notes (2006-10-20)   The  nautical mile  (1852 m)
A lesson in averages...  over the surface of an oblate spheroid.

Originally, one nautical mile was meant to correspond to the distance between two points on the Ocean separated by a  geocentric  angle of one minute  (1°/60).  Taking into account the oblateness of the Earth, the nautical mile should be defined, more precisely, as an "average" minute of  latitude.

Following a 1929 resolution of the International Hydrographic Conference, the  international nautical mile (NM) is now officially defined as exactly equal to 1852 m.  (This conventional definition had been legal in France since 1906.)  The so-called "British Admiralty" nautical mile also has a conventional value:  6080 ft, or exactly 1853.184 m.

If the meter was still defined in reference to the Meridian, there would be exactly 10 000 000 meters in a quarter meridian, corresponding to a change in latitude of 5400 minutes.  A "typical" single minute of latitude would then be equal to 100000/54 or about 1851.852 m.  To the nearest whole number of meters, this gives the above definition of the international nautical mile.

The nautical mile formerly used in the United States had a different definition:  It was  1/21600  of the meridian of the so-called  Clark spheroid  (which is a sphere having the same surface as the Earth).  Modern values for the Reference Ellipsoid (see below) would put the radius of Clark's spheroid at about 6371007.181 m and the "conventional" value of the US nautical mile at 1853.250866 m (6080.219 ft, or 6080.207 "US Survey" ft).

The true shape of the Earth (the so-called "geoid") is defined as an equipotential surface  (i.e., a surface which is everywhere perpendicular to the local vertical indicated by a plumb line).  Such a shape turns out to be quite complicated, but its irregularities are best charted in comparison with a close "regular" approximation, the Reference Ellipsoid, whose exact shape and size have been defined in 1979, in Canberra (IUGG 1980).  The Reference Ellipsoid is  defined to have the following characteristics:

• Equatorial radius  (a)  exactly  equal to 6378137 m.
• Oblateness  f  =  1 / 298.257222101   =   1 - b/a

The polar radius  b  =  a (1- f )  is the distance of either pole to the center.  The Meridian is an ellipse of eccentricity  e = Ö(2f-f ² ) = 0.0818191910428.  The length of the reference Meridian is thus:

40007862.91692186(12) m

The "uncertainty" stated (between parentheses) is in units of the last significant digit shown and corresponds to an "uncertainty" of half a unit in the least significant digit of the value for the oblateness  f  specified above.  That's a ludicrous precision of  120 nm.  Divide this by 21600 and you obtain:

1852.215875783419(5)m

Forsaking most of the ludicrous precision, we may state this result in terms of the aforementioned international nautical mile (NM) henceforth equal to 1852 m.

1852.216 m   =   1.000117 NM

Now, none of the above computations yield the proper "average minute of latitude" on that spheroid.  Actually, proper averaging yields subtler and  simpler  results  with expressions involving only elementary functions and shunning all the  arcana  related to the perimeter of an ellipse.

Like many questions about averages, this one can only be properly addressed if probabilistic assumptions are made explicit:  Let's do our "averaging" by giving equal weight to any square inch on the surface of the Globe.  This is the simplest natural approach which refers only to the spheroid's geometry, shunning  ad hoc  alternatives for the real Earth (like considering only open oceans, or even weighing maritime routes according to their estimated traffic).

Thus, a "random" point on the Earth is more likely to be at a low latitude than a high one.  Loosely speaking, there are "more points" in the vicinity of the equator, at latitude 0°, than in the North Pole's vicinity, at latitude +90°.

On a sphere, for any infinitesimal positive quantity d, a random
latitude would be between q and q+d with probability  ½ cos(q) d.

The correct "weighing" for an oblate spheroid is somewhat more complicated and yields the following (exact) formula, in terms of the radius of the equator (a) and the eccentricity of the meridian (e):

Average Minute of Geocentric Latitude
For the Earth, this is :   1853.256 m = 1.000678 NM

2p a   2 [ 1 - e2 ( 1 - 2 e2 / 15) ] 21600 Ö(1-e2 ) (1 + (1-e2 ) atanh(e) / e )

The correction factor is close to unity for small values of e, but it can become large when e approaches one.  Technically, the value of an average minute of latitude is indeed  infinite  for a flat disk  (e = 1).

The difference with the "British Admiralty nautical mile" of 6080 ft is less than 3 inches, and the former definition (6080.2 ft) of the US nautical mile is only a centimeter off the mark.  These two are indeed good approximations to the "average minute of latitude", while the more naive "new" nautical mile of 1852 m is over a meter short!

This is not the end of the story, though.  So far, the "latitude" we spoke about has been the "geocentric" one:  the angle between the plane of the equator and a line going through the center of the Earth.  This is not the "geodetic" latitude which is really seen by navigators who measure the angle between the horizon (or a vertical line, which does  not  go through the center of the Earth) and a fixed direction, like the celestial pole.  If the Earth was a sphere, there would not be any difference whatsoever.  Also, since there are still 21600 "geodetic" minutes of latitude in a full meridian; this yields the same "naive" average of 1852.216 m as for the "geocentric" minute.  Things change when we average "correctly" over the surface of the globe: First, we observe that distance changes along a meridian with geodetic latitude at a rate (per radian) equal to the meridian's radius of curvature —this is one way the radius of curvature can be defined, for any planar curve.  Thus, at a point where the radius of curvature of the meridian is R, the geodetic minute of latitude is equal to Rp/10800.  The average (over the ellipsoid's surface) of R is 6356828.89 m, which gives a low value (1849.12657 m) for the average "geodetic" mile, almost 3 meters short of the international nautical mile.  The corresponding formula is very similar to the "geocentric" one:

Average Minute of Geodetic Latitude
For the Earth, this is :   1849.127 m = 0.998448 NM

2p a   2 [ 1 - e2 ( 4/3 - 8 e2 / 15) ] 21600 Ö(1-e2 ) (1 + (1-e2 ) atanh(e) / e )

This is smaller than the geocentric expression for roundish spheroids but becomes higher when the polar radius is below  a/Ö6, or about 40.825% of the equatorial radius.

Compared to the geodetic latitude, the geocentric value is between 33.333% lower  (for e near 1)  and 14.564% higher  (for e = Ö(9-Ö21)/Ö8) ).

Kaimbridge M. GoldChild  (2011-09-21)   Great Ellipses
What is the average length of a  great ellipse ?

A  great ellipse  on an ellipsoid is defined as the intersection of the surface of the ellipsoid with a plane that goes through its center.  Great ellipses  are not necessarily geodesics  (which is to say that the shortest path between two points on the surface of an ellipsoid isn't always an arc of a great ellipse).

If the center of the ellipsoid and two points on its surface are aligned, those two points are said to be  antipodal.  The  great ellipses  through two antipodal points may have different lengths.

Sidey P. Timmins, U.S. Census Bureau  (2012-01-09)
Curved Surface Area with a Given Boundary on the Surface of the Earth
What area is delimited by a closed curve on an oblate spheroid?

On an oblate ellipsoid of revolution with equatorial radius  a  and polar radius  b,  let's consider the curved surface area bounded by a "rectangle" consisting of the equator, two arcs of meridians and an arc of the parallel line of constant geodetic latitude  j.

We tally that surface area  positively  in the Northern Hemisphere and  negatively  in the Southern Hemisphere when going from West to East  (increasing longitudes)  and  vice-versa  westward.

That area is directly proportional to the difference between the longitudes of the two meridians.  The coefficient of proportionality is some  odd  function  f  (which shall be computed below)  of the aforementioned geodetic latitude.  In particular, for an infinitesimal longitude difference  dq, the above curved rectangle encloses the following infinitesimal surface area:

dA_o   =   - f (j)  dq

The negative sign makes the above consistent with the sign convention used for planar areas.  The equator oriented eastward plays the rôle of the x-axis and the northward meridian is analoguous to the y-axis.

Indeed, in the Northern hemisphere where  f  is positive, a positive dq corresponds to a clockwise boundary  (eastward at the given longitude and westward back on the equator)  which must encircle a negative infinitesimal area, by convention.

Along a  closed loop  (possibly specified by giving  [j(t),q(t)]  as a periodic function of a parameter  t )  the contour integration of the above would give the area  enclosed  by the curve  only  if it doesn't go around the polar axis.  To remove this inconvenient restriction  (which has no counterpart in planar geometry)  we shall consider instead the infinitesimal area swept by the arc of a meridian originating at the  North pole, counted positively eastward :

dA   =   [  f (p/2) - f (j) ]  dq

The integral of that is the correct area encircled by the curve  (defined only modulo the total surface area of the ellipsoid)  even if the curve goes around the polar axis many times.  (The area bounded by a self-crossing loop is tallied like in the planar case, as depicted at right.)

Signed Area Encircled by a Closed Curve Drawn on a Spheroid :

A   =   ò  [  f (p/2) - f (j) ]  dq   =   ò  [  f (p/2) - f ( j(t) )  ]  q' (t)  dt

The general expression of  f  can be obtained by means of cartesian coordinates  (well, cylindrical ones, actually).  The correspondence with geodetic coordinates is given elsewhere on this site).  The basic infinitesimal rectangular surface area is equal to the displacement  r dq  along the parallel multiplied into the (perpendicular) displacement  ds  along the meridian:

df dq   =   r dq ds   =   r dq dz [ 1 + (dr/dz)² ]^½

Using the equation of the meridian   (r/a)² + (z/b)² = 1   we obtain:

df   =   r dz [ 1 + (dr/dz)² ]^½       with   dr/dz = -(a/b)² z/r

Therefore,

( df / dz ) ²   =   r² + (a/b)⁴ z²   =   a² - (a/b)² z² + (a/b)⁴ z²
=   a² [ 1 + ( a e z / b² ) ² ]
 
df   =   a [ 1 + ( a e z / b² ) ² ] ^½  dz

Introducing  u = a e z / b²   we have:   df   =   ( b² / e )  [ 1 + u ² ] ^½  du   which is readily integrated  (positing  f = 0  when  z = 0) :

f   =   ( b ²/ 2e )  [ asinh (u)  +  u ( 1 + u² )^½ ]

Using the expression of z, we may express that result in geodetic terms:

Defining Geodetic Latitude (j)

e   =   ( 1 - b 2/ a 2 )½ u   =   e a  sin j   /   ( a2 cos2 j  +  b2 sin2 j ) ½ f (j)   =   ( b 2/ 2e )  [ asinh (u)  +  u ( 1 + u2 )½ ] where   asinh (u)   =   Log [ u  +  ( 1 + u2 )½ ]

Let's work that out, piece by piece:

a² cos² j  +  b² sin² j   =   a²  [ 1 - e² sin² j ]
1 + u²   =   1 + e² sin² j / [ 1 - e² sin² j ]   =   1 / [ 1 - e² sin² j ]
u  ( 1 + u² )^½   =   e sin j / [ 1 - e² sin² j ]
 
[ u  +  ( 1 + u² )^½ ] ²   =   [ 1 + e sin j ] ² / [ 1 - e² sin² j ]
  =   (1+e sin j) / (1-e sin j)
Log [ u + ( 1 + u² )^½ ]   =   ½ Log [(1+e sin j) / (1-e sin j)]   =   atanh (e sin j)

Plugging both of those results into our previous expression, we obtain:

Weighing function   f (j)   for an oblate spheroid of eccentricity  e

f (j)   =   ( b 2/ 2 )  [ atanh (e sin j) / e  +  sin j / (1 - e2 sin2 j) ]

On 2012-01-23, Sidey Timmins wrote:   [edited summary] It works!  A very simple formula - wow!   [...]   Thank you.

Total Surface Area of an Oblate Spheroid :

The surface area of the Northern hemisphere being  2 p f (p/2)  the surface area of the entire oblate spheroid is thus shown to be   4 p f (p/2)   which boils down to the very nice formula advertised elsewhere:

Surface Area of an Oblate Spheroid

A   =   2 p [ a 2  +  b 2 atanh(e)/e ]

For a sphere  ( a = b  and e = 0 )  that's Archimedes' formula:  4 p a ²

Area Bounded by a Closed Polygonal Line :

Sidey Timmins' original question was actually about the case where the boundary curve is a "polygonal" line.  The value of the above integral can be obtained through the following  trapezoidal approximation  (which does turn into an  exact formula  when the special type of "straight" lines  investigated below is used to join one vertex of the polygon to the next).

Area Bounded by a Noncircumnavigating Polygonal Line

- A   =   ò  f (j)  dq   =   Si   f ( ji )  ( qi+1 - qi-1 ) / 2

In this, there are n vertices in the polygon and the indices involved are understood to be modulo n  (the points of indices 0 and n coincide).

(2012-01-22)   "Pseudo-Straight Lines"  on an Oblate Spheroid
Curves that sweep out areas varying  quadratically  with longitude.

When used as the sides of polygonal boundaries, such curves turn the above trapezoidal approximation into an  exact formula.

With the notations of the previous section, those curves feature an increase in  f  proportional to the increase in longitude:

d f   =   k dq
f (j)   =   f_o  +  k q

A parametrical family of curves is thus explicitly defined which does not appear to have any particularly enlightening characteristics...

By the  Gauss-Bonnet theorem, those curves would actually be  geodesics  on a surface of revolution of constant Gaussian curvature  (sphere or pseudo-sphere).  This is not the case for an oblate ellipsoid, but it's close enough for practical purposes when the vertices are not widely separated  (so the variations in curvature from one vertex to the next are negligible).