A PHP Function to Compute the Circumference or Perimeter of an Ellipse or Circle
by Evaluating the Complete Elliptic Integral of the Second Kind Using an Infinite
Power Series Summation Based on Legendre Polynomials.

MathJax equations display requires Java Script to be enabled.
PHP Function by Jay Tanner - 2017

The Complete Elliptic Integral Of The Second Kind

The complete elliptic integral of the second kind may be symbolized as:

$\begin{align*} (Equation~1):~~~~~~~~ \int_0^{\large\frac{\pi}{2}}{\sqrt{1-\epsilon^2\cdot\sin^2(\phi)}\;\;d\phi} \end{align*}$ =
The integral defines the length of the elliptic arc from $0$ to $2\pi$ (or $0$ to $90$°) or only $\frac{1}{4}$ of the full $360$° span of the ellipse.  It is a line integral in the sense that it defines the length of a line, in this case, one quarter of the ellipse circumference.  Numerically solving the integral is difficult enough as it is, so solving only for the first quarter and then multiplying the result by $4$ seems like the simplest way to go about it.

Notice that the angular arc-span is measured positively counter-clockwise.

The Basic Ellipse

To make use of the integral, we need to adapt it to our purpose, which in this case, is to 'simply' compute the circumference of an ellipse. 

Generally, the idealized equations equate the major radius $(R)$ to unity or ($R = 1$) and this is often left invisible and not clearly indicated.  In the real world, we need to be able to use any units of measure we wish, so we will use $R$ rather than an implied invisible $1$.  Thus, $R$ can be interpreted as any implied units, meters, kilometers, miles, etc.

Given an ellipse with parameters $(R \ge r)$   =  

The elliptic eccentricity $(\epsilon)$, may be found from:

$ (Eq.~2):~~~~~~~~~~~~ \begin{align*} \epsilon~^2 ~=~ 1 - \frac{\Large r^2}{R^2} \end{align*} $

from which
$ (Eq.~3):~~~~~~~~ \begin{align*} \epsilon ~=~ \sqrt{1 - \frac{\Large r^2}{R^2}} \end{align*} $

The eccentricity tells us how much the form of the ellipse deviates from a perfect circle.  Zero means the form is a perfect circle and the greater the value of $\large\epsilon$, within the contstraints $(0 ~\le~ \large\epsilon ~\lt~ 1)$, the 'flatter' the ellipse. 

The symbol $\large(\epsilon)$ (Greek lower-case letter epsilon) is used here for the eccentricity instead of $\large(e)$, so as not to confuse it with Euler's number $\large(e = 2.71828...)$, used as the base of natural logarithms.

To be able to use any desired units of length or distance measure, simply specify the $(R)$ value in any units $($m, km, mi, etc.$)$, so that the resulting circumference $(C)$ will be in the same units as $(R)$.

$\begin{align*} (Eq.~3):~~~~~ \frac{C}{4} ~=~ R \times \int_0^{\large\frac{\pi}{2}}{\sqrt{1-\epsilon^2\cdot\sin^2(\phi)}\;\;d\phi} \end{align*}$  = =

Since the integral equates to $\frac{1}{4}$ the length defined by the complete elliptic integral, we must multiply it by $4$ to obtain the full circumference.

$\begin{align*} (Eq.~4):~~~~~ C ~=~ 4R~\times \int_0^{\large\frac{\pi}{2}}{\sqrt{1-\epsilon^2\cdot sin^2(\phi)}\;\;d\phi} \end{align*}$  =  $4$ $\times$ =

The numerical computation of the elliptic integral will involve evaluating the infinite power series summation below until a given precision limit is reached.  This series is based on Legendre polynomials.  Other series are possible, but this series was chosen for the simplicity of programming it as well as its fast convergence for many practical problems.

The integral on the left formally defines the problem and the summation equation on the right side provides one of several possible numerical solutions adopted to solve the problem, defined purely in terms of the basic ellipse parameters $(R,~r)$.

$\begin{align*} (Eq.~5):~~~~~ \frac{C}{4} ~=~ R \int_0^{\large\frac{\pi}{2}}{\sqrt{1-\epsilon^2\cdot sin^2(\phi)}\;\;d\phi} ~~=~~ \frac{\pi R}{2}\left[ 1 - \sum_{n=1}^\infty \left(\frac{(2n-1)!!}{(2n)!!}\right)^2 \cdot \frac{{\left(\sqrt{1 - \Large\frac{r^2}{R^2}}\right) ^ {2n}}}{2n-1}\right] \end{align*}$

As stated above, since equation $(4)$ only solves for the length of $\frac{1}{4}$ of the ellipse, we have to multply it by $4$ to obtain the full circumference.  Doing this for the working summation gives:

$\begin{align*} (Eq.~6):~~~~~~~ C ~~=~~ 2\pi{R}\left[ 1 - \sum_{n=1}^\infty \left(\frac{(2n-1)!!}{(2n)!!}\right)^2 \cdot \frac{{\left(\sqrt{1 - \Large\frac{r^2}{R^2}}\right) ^ {2n}}}{2n-1}\right] \end{align*}$

Now, let's rewrite a more simplified equation by combining all the integer values in the numerator and denominator into an integer fraction of the form $\left(\frac{Y}{X}\right)$, where

$ \begin{align*} (Eq.~7):~~~~~~~ C ~~=~~ 2\pi{R}\left[ 1 - \sum_{n=1}^\infty \left(\frac{(2n-1)!!~^2}{(2n)!!~^2\cdot (2n-1)}\right) \cdot {\left(1 - \large\frac{r^2}{R^2}\right) ^ n}~\right] \end{align*} $

and setting

$ \begin{align*} (Eq.~8):~~~~~~ Y ~=~ (2n-1)!!~^2 ~=~ Integer~numerator \end{align*} $

$ \begin{align*} (Eq.~9):~~~~~~ X ~=~ (2n)!!~^2\cdot(2n-1) ~=~ Integer~denominator \end{align*} $

puts the equation into the form:

$\begin{align*} (Eq.~10):~~~~~~ C ~~=~~ 2\pi{R}\left[ 1 - \sum_{n=1}^\infty \left(\frac{Y}{X}\right) \cdot {\left(1 - \large\frac{r^2}{R^2}\right) ^ n}~\right] \end{align*}$

which can be simplified even more by recalling $(Eq.~2):~~~~ \large\epsilon^2 = 1 - \frac{\Large r^2}{R^2}$

and substituting $\large(\epsilon^2)$ for $\left(1 - \large\frac{\Large r^2}{R^2}\right)$ in $(Eq.~10)$, as a pre-computed constant, to obtain:

$ \begin{align*} (Eq.~11):~~~~~~~~ C ~=~ 2\pi{R}\left[ 1 - \sum_{n=1}^\infty \left(\frac{Y}{X}\right) \cdot (\epsilon^2)^n\right] ~=~ 2\pi{R}~\left[ 1 - Summation\right] \end{align*} $

If the ellipse is the cross-section of an ellipsoid, then the $(R, ~r)$ parameters respectively represent the equatorial and polar radii, sometimes symbolized as $(a, ~b)$, where $(a \ge b)$.

In the special case of a circle, where ($R = r$), $(\epsilon = 0)$, then $[Summation ~=~ 0]$, and the equation for the circumference reduces to the classical circular circumference equation:

$\begin{align*} (Eq.~12):~~~~~~~~ C ~=~ 2\pi{R}~\left[ 1 - Summation\right] ~=~ 2\pi{R}~\left[ 1 - 0\right] ~=~ 2\pi R \end{align*}$

When $(R \ne r)$, then $(0 ~<~ \epsilon ~<~ 1)$ and the form is an ellipse, and computing the circumference will require solving the elliptic integral via a power series summation like the one used here. 

Below is a listing of a PHP function designed to compute the high-precision general circumference of an ellipse.  It is based on $(Eq.~11)$ above and makes use of arbitrary-precision arithmetic to compute the EXACT values of the extremely large integer fractions involved, so it could take several seconds in some extreme cases.

For most practical applications, such as geodesic and meridian computations, the function should converge sufficiently fast and accurately to $16$ decimals  Also, the function could be easily tweaked for higher accuracy, if actually needed. 

Since the ellipse circumference function calls another function to compute the double factorial, the code for that function is also included.

For any integer $N$, the general equation for the double factorial of $N$ can be expressed as:

$\begin{align*} (Eq.~13):~~~~~~~~ N~!! ~=~ \prod_{n ~=~ 0}^{\lceil\frac{N}{2}\rceil-1} (N - 2n) ~=~ N\cdot(N-2)\cdot(N-4)\cdot(N-6)\cdot ~...~ \end{align*}$

The code for both PHP functions is included in the demonstration program listed below.

Double-Click to select ALL code

Double-Click to select ALL code

What the functions above do is delineated by the first few terms of the infinite series worked out
numerically below.  Due to the extremely large numbers, arbitrary-precision is used to maintain
accuracy within the program.

Notice how quickly the numbers become very large in the series.  Although very large at first, the integer
fractions can subsequently be reduced to lowest terms, which reduces the sizes of the numbers used in
the terms, as demonstrated below.

As a visual and mental aid, the first few terms of the summation series are worked out numerically below.

This is the summation that is evaluated to 16 decimals by the PHP functions defined above.

Summation Series Form 1:

The initial symbolic series, before full numerical expansion, may be written as:

$ \begin{align*} C ~=~ 2\pi R\left[1 - \left[\left(\frac{1!!^2}{2!!^2} \cdot \frac{\epsilon^{2}}{1}\right) + \left(\frac{3!!^2}{4!!^2} \cdot \frac{\epsilon^{4}}{3}\right) + \left(\frac{5!!^2}{6!!^2} \cdot \frac{\epsilon^{6}}{5}\right) + \left(\frac{7!!^2}{8!!^2} \cdot \frac{\epsilon^{8}}{7}\right) + \left(\frac{9!!^2}{10!!^2} \cdot \frac{\epsilon^{10}}{9}\right) + \left(\frac{11!!^2}{12!!^2} \cdot \frac{\epsilon^{12}}{11}\right) + \left(\frac{13!!^2}{14!!^2} \cdot \frac{\epsilon^{14}}{13}\right) + \left(\frac{15!!^2}{16!!^2} \cdot \frac{\epsilon^{16}}{15}\right) + \left(\frac{17!!^2}{18!!^2} \cdot \frac{\epsilon^{18}}{17}\right) + \left(\frac{19!!^2}{20!!^2} \cdot \frac{\epsilon^{20}}{19}\right) + \left(\frac{(2n-1)!!^2}{(2n)!!^2} \cdot \frac{\epsilon^{2n}}{2n-1}\right) + ~ ... ~\right]\right] \end{align*} $

Summation Series Form 2:

Computing the values of the double factorials gives:

$ \begin{align*} C ~=~ 2\pi R\left[1 - \left[\left(\frac{1^2}{2^2} \cdot \frac{\epsilon^{2}}{1}\right) + \left(\frac{3^2}{8^2} \cdot \frac{\epsilon^{4}}{3}\right) + \left(\frac{15^2}{48^2} \cdot \frac{\epsilon^{6}}{5}\right) + \left(\frac{105^2}{384^2} \cdot \frac{\epsilon^{8}}{7}\right) + \left(\frac{945^2}{3840^2} \cdot \frac{\epsilon^{10}}{9}\right) + \left(\frac{10395^2}{46080^2} \cdot \frac{\epsilon^{12}}{11}\right) + \left(\frac{135135^2}{645120^2} \cdot \frac{\epsilon^{14}}{13}\right) + \left(\frac{2027025^2}{10321920^2} \cdot \frac{\epsilon^{16}}{15}\right) + \left(\frac{34459425^2}{185794560^2} \cdot \frac{\epsilon^{18}}{17}\right) + \left(\frac{654729075^2}{3715891200^2} \cdot \frac{\epsilon^{20}}{19}\right) + ~...~ \right] \right] \end{align*} $

Summation Series Form 3:

Computing the integer values in full gives the series in the form:

$ \begin{align*} C ~=~ 2\pi R\left[1 - \left[\left(\frac{1}{4} \cdot \frac{\epsilon^{2}}{1}\right) + \left(\frac{9}{64} \cdot \frac{\epsilon^{4}}{3}\right) + \left(\frac{225}{2304} \cdot \frac{\epsilon^{6}}{5}\right) + \left(\frac{11025}{147456} \cdot \frac{\epsilon^{8}}{7}\right) + \left(\frac{893025}{14745600} \cdot \frac{\epsilon^{10}}{9}\right) + \left(\frac{108056025}{2123366400} \cdot \frac{\epsilon^{12}}{11}\right) + \left(\frac{18261468225}{416179814400} \cdot \frac{\epsilon^{14}}{13}\right) + \left(\frac{4108830350625}{106542032486400} \cdot \frac{\epsilon^{16}}{15}\right) + \left(\frac{1187451971330625}{34519618525593600} \cdot \frac{\epsilon^{18}}{17}\right) + \left(\frac{428670161650355625}{13807847410237440000} \cdot \frac{\epsilon^{20}}{19}\right) + ~ ... ~ \right]\right] \end{align*} $

Summation Series Form 4:

Combining the three integer values in the numerators and denominators of each term in the

series into integer fractions of the form $\left(\frac{Y}{X}\right)$, separate from the powers of $(\epsilon)$, gives:

$ \begin{align*} C ~=~ 2\pi R\left[1 - \left[\left(\frac{1}{4}\right)\cdot{\epsilon^{2}} + \left(\frac{9}{192}\right)\cdot{\epsilon^{4}} + \left(\frac{225}{11520}\right)\cdot{\epsilon^{6}} + \left(\frac{11025}{1032192}\right)\cdot{\epsilon^{8}} + \left(\frac{893025}{132710400}\right)\cdot{\epsilon^{10}} + \left(\frac{108056025}{23357030400}\right)\cdot{\epsilon^{12}} + \left(\frac{18261468225}{5410337587200}\right)\cdot{\epsilon^{14}} + \left(\frac{4108830350625}{1598130487296000}\right)\cdot{\epsilon^{16}} + \left(\frac{1187451971330625}{586833514935091200}\right)\cdot{\epsilon^{18}} + \left(\frac{428670161650355625}{262349100794511360000}\right)\cdot{\epsilon^{20}} + ~ ... ~ \right]\right] \end{align*} $

Summation Series Form 5:

Reducing all of the integer $\left(\frac{Y}{X}\right)$ fractions to their lowest terms to obtain the series in the smallest possible integer values gives:

$ \begin{align*} C ~=~ 2\pi R\left[1 - \left[\left(\frac{1}{4}\right)\cdot{\epsilon^{2}} + \left(\frac{3}{64}\right)\cdot{\epsilon^{4}} + \left(\frac{5}{256}\right)\cdot{\epsilon^{6}} + \left(\frac{175}{16384}\right)\cdot{\epsilon^{8}} + \left(\frac{441}{65536}\right)\cdot{\epsilon^{10}} + \left(\frac{4851}{1048576}\right)\cdot{\epsilon^{12}} + \left(\frac{14157}{4194304}\right)\cdot{\epsilon^{14}} + \left(\frac{2760615}{1073741824}\right)\cdot{\epsilon^{16}} + \left(\frac{8690825}{4294967296}\right)\cdot{\epsilon^{18}} + \left(\frac{112285459}{68719476736}\right)\cdot{\epsilon^{20}} + ~ ... ~ \right]\right] \end{align*} $

It might be noted that in this final reduced fractions series, the fractional denominators all equate to
even powers of $2$, so the summation series could also be expressed in the somewhat abbreviated
form shown below:

Abbreviated Powers-of-2 Series:

Abbreviating the powers of $2$ puts the series into the semi-symbolic form:

$ \begin{align*} C = 2\pi R\left[1 - \left[\left(\frac{1}{2^{2}}\right)\cdot{\epsilon^{2}} + \left(\frac{3}{2^{6}}\right)\cdot{\epsilon^{4}} + \left(\frac{5}{2^{8}}\right)\cdot{\epsilon^{6}} + \left(\frac{175}{2^{14}}\right)\cdot{\epsilon^{8}} + \left(\frac{441}{2^{16}}\right)\cdot{\epsilon^{10}} + \left(\frac{4851}{2^{20}}\right)\cdot{\epsilon^{12}} + \left(\frac{14157}{2^{22}}\right)\cdot{\epsilon^{14}} + \left(\frac{2760615}{2^{30}}\right)\cdot{\epsilon^{16}} + \left(\frac{8690825}{2^{32}}\right)\cdot{\epsilon^{18}} + \left(\frac{112285459}{2^{36}}\right)\cdot{\epsilon^{20}} + \left(\frac{370263621}{2^{38}}\right)\cdot{\epsilon^{22}} + \left(\frac{19870814327}{2^{44}}\right)\cdot{\epsilon^{24}} + \left(\frac{67607800225}{2^{46}}\right)\cdot{\epsilon^{26}} + \left(\frac{931331941875}{2^{50}}\right)\cdot{\epsilon^{28}} + \left(\frac{3241035157725}{2^{52}}\right)\cdot{\epsilon^{30}} + \left(\frac{2913690606794775}{2^{62}}\right)\cdot{\epsilon^{32}} + \left(\frac{10313859829588425}{2^{64}}\right)\cdot{\epsilon^{34}} + \left(\frac{147068001273760875}{2^{68}}\right)\cdot{\epsilon^{36}} + \left(\frac{527570807893408125}{2^{70}}\right)\cdot{\epsilon^{38}} + \left(\frac{30451387031607516975}{2^{76}}\right)\cdot{\epsilon^{40}} + \left(\frac{110412172026168752025}{2^{78}}\right)\cdot{\epsilon^{42}} + \left(\frac{1608732721339962891075}{2^{82}}\right)\cdot{\epsilon^{44}} + \left(\frac{5884494925884363316125}{2^{84}}\right)\cdot{\epsilon^{46}} + \left(\frac{1382856307582825379289375}{2^{92}}\right)\cdot{\epsilon^{48}} + \left(\frac{5095548922181194957605489}{2^{94}}\right)\cdot{\epsilon^{50}} + \left(\frac{75347791458762166858320219}{2^{98}}\right)\cdot{\epsilon^{52}} + \left(\frac{279375967507591408803895133}{2^{100}}\right)\cdot{\epsilon^{54}} + \left(\frac{16620019291523039931905190055}{2^{106}}\right)\cdot{\epsilon^{56}} + \left(\frac{61954530890516920554723865425}{2^{108}}\right)\cdot{\epsilon^{58}} + \left(\frac{926013721710259572557939375219}{2^{112}}\right)\cdot{\epsilon^{60}} + \left(\frac{3467974385468495527196694913021}{2^{114}}\right)\cdot{\epsilon^{62}} + \left(\frac{13327425563355428311016898550739703}{2^{126}}\right)\cdot{\epsilon^{64}} + \left(\frac{50115525878733222161261891244517065}{2^{128}}\right)\cdot{\epsilon^{66}} + \left(\frac{755201090663955648831472444186407675}{2^{132}}\right)\cdot{\epsilon^{68}} + \left(\frac{2850036442562830175141140497529602189}{2^{134}}\right)\cdot{\epsilon^{70}} + \left(\frac{172374426322411173185388238239475569431}{2^{140}}\right)\cdot{\epsilon^{72}} + \left(\frac{652605297026338283871342029799270910417}{2^{142}}\right)\cdot{\epsilon^{74}} + \left(\frac{9897545709748482283090298097371213946075}{2^{146}}\right)\cdot{\epsilon^{76}} + \left(\frac{37579438838788616163607147608362104233125}{2^{148}}\right)\cdot{\epsilon^{78}} + \left(\frac{9143829058254046084928891156066667202003975}{2^{156}}\right)\cdot{\epsilon^{80}} + \left(\frac{34807473018303177214431870617293636779074025}{2^{158}}\right)\cdot{\epsilon^{82}} + \left(\frac{530636374381070885289400150022823401509557075}{2^{162}}\right)\cdot{\epsilon^{84}} + \left(\frac{2024683407927774524454687971017316980881517125}{2^{164}}\right)\cdot{\epsilon^{86}} + \left(\frac{123739948773767707506962128476636851848089414375}{2^{170}}\right)\cdot{\epsilon^{88}} + \left(\frac{473144900422362152704398894219555132770250042225}{2^{172}}\right)\cdot{\epsilon^{90}} + \left(\frac{7243857369604368761347687418306572817214092801475}{2^{176}}\right)\cdot{\epsilon^{92}} + \left(\frac{27752270221349829256353770312869409575410985685325}{2^{178}}\right)\cdot{\epsilon^{94}} + \left(\frac{27243478600625082386653951190466803733195117614427375}{2^{188}}\right)\cdot{\epsilon^{96}} + \left(\frac{104560039693777648560189987596897791087627242322760625}{2^{190}}\right)\cdot{\epsilon^{98}} + \left(\frac{1606544097886954814597607121428815180503175052840752451}{2^{194}}\right)\cdot{\epsilon^{100}} + ~...~ \right]\right] \end{align*} $

In all of the final reduced fractions, the numerator is always an odd integer and the denominator is always
an even power of $2$.

When $(\epsilon ~=~ 0)$, then the form is a perfect circle and $(Summation = 0)$.

When $(0 ~<~ \epsilon ~<~ 1)$, then the form is an ellipse and $(0 < Summation < 1)$.

The closer that the eccentricity $\large(\epsilon)$ gets to $1$, the longer the series will take to converge. However, for the
intended geophysical and astronomical applications, the elliptic eccentricities don't tend to be so extreme.

One can see that the integers become very large, while the numerator to denominator ratio diminishes
rather slowly, which translates to the series converging more slowly as the eccentricity, or the elliptical
deviation from a perfect circle, increases. 

For small values of eccentricity $\large(\epsilon)$, sufficient approximation formulas could possibly be derived from
the numbers in first few terms of the series.