  function [P,n] = EllipsePerimeter(a,b,tol)
% P is an approximation to the perimeter of the ellipse (x/a)^2 + (y/b)^2 = 1.
% Its relative error is less than or equal to tol, a positive number
%     which should be no smaller than 10^-10.
% n is the number of edges in the final approximating polygons.

  tol = max(tol,10^-10);
  n = 2; Inner = 4*a; Outer = inf;
  while ( abs(Inner-Outer)>=tol*Inner)
      n = 2*n;
%     The necessary sines and cosines...
      theta = linspace(0,2*pi,n+1); c = cos(theta); s = sin(theta);
%     The x and y coordinates of the vertices...
      gamma = (1 - c(2))/s(2);
      Px = a*c; Qx = a*(c - gamma*s);
      Py = b*s; Qy = b*(s + gamma*c);
%     The two perimeters..
      Inner = PolyLine(Px,Py); Outer = PolyLine(Qx,Qy);
  end     
  P = (Inner+Outer)/2;
  
  function L = PolyLine(u,v)
% u and v are column (n+1)-vectors with u(1) = u(n+1) and v(1) = v(n+1).
% L is the perimeter of the polygon with vertices 
% (u(1),v(1)),...,(u(n),v(n))
  n = length(u)-1;
  L = 0;
  for k=1:n
%    Add in the distance from (u(k),v(k)) to (u(k+1),v(k+1))...
     L = L + sqrt((u(k) - u(k+1))^2 + (v(k) - v(k+1))^2 );
  end