Perimeter of Ellipse

[Joe Horn]

Although simple formulae for the perimeter of an ellipse exist, they are only approximations. Exact formulae are complicated. The following program uses a converging iteration technique instead. Ported from a QBASIC program by Gérard P. Michon.

Syntax: EllipsePerimeter(a,b), where a & b are the lengths of the axes (order doesn't matter)
Output: perimeter of ellipse

Examples:
EllipsePerimeter(0.5, 0.5) --> 3.14159265359
EllipsePerimeter(3, 4) --> 22.1034921607

Code:

gk(h);
cayley(x);

EXPORT EllipsePerimeter(a,b)
BEGIN
LOCAL x,h,P;
a:=ABS(a);
b:=ABS(b);
IF a < b THEN x:=a; a:=b; b:=x; END;
IF b < 0.28*a THEN
 P := 4*a*cayley((b/a)^2);
ELSE
 h := ((a-b)/(a+b))^2;
 P := pI*(a+b)*gk(h);
END;
RETURN P;
END
;

gk(h)
BEGIN
LOCAL z,x,n;
z := 0; x := 1; n := 0;
WHILE z + x <> z
 DO
 n := n + 1;
 x := h * x * ((n-1.5)/n)^2;
 z := z + x;
END;
RETURN 1 + z;
END;

cayley(x)
BEGIN
LOCAL y,t,n,z,u,v,w;
y := LOG(16 / x) - 1;
t := x / 4;
n := 1;
z := 0;
u := t * y;
v := (n - .5) / n;
w := .5 / ((n - .5) * n);
WHILE z <> z + u
 DO
 z := z + u;
 n := n + 1;
 t := x * t * v;
 v := (n - .5) / n;
 t := t * v;
 y := y - w;
 w := .5 / ((n - .5) * n);
 y := y - w;
 u := t * y;
END;
RETURN 1 + z;
END;