  function Demo5
% Shows the convergence of vertices to the limiting ellipse
% for a range of n-values
  for n=10:10:40
     for eg=1:3
        ShowAll(n)
     end
  end

  function ShowAll(n)
% Shows polygon smoothing for a random n-point problem
 
  nSmoothings = 10*n;
% Generate random example... 
  x = randn(n,1); x = x-sum(x)/n; x = x/norm(x);
  y = randn(n,1); y = y-sum(y)/n; y = y/norm(y);
  
% Get points on the limiting ellipse...
  [semiaxis45,semiaxis135] = Ellipse(x,y);
  xfR = semiaxis45 *cos(linspace(0,2*pi));
  yfR = semiaxis135*sin(linspace(0,2*pi));
  u = (xfR - yfR)/sqrt(2);
  v = (xfR + yfR)/sqrt(2);

% Display the original polygon...
  close all
  figure
  set(gcf,'position',[100 100 1000 500])
  subplot(1,2,1)
  a = .6;
  plot([x;x(1)],[y;y(1)],'k',x,y,'or',[-a a],[a -a],'m',[-a a],[-a a],'m')
  axis equal off
  title('Original Polygon','FontWeight','bold','FontSize',14)
  
  % Repeatedly smooth and display the polygon... 
  for k=1:nSmoothings
     subplot(1,2,2)
     [x,y] = Smooth(x,y);
     plot(u,v,'g','Linewidth',3)
     hold on
     plot([x ; x(1)],[y ; y(1)],'k',x,y,'or',[-a a],[a -a],'m',[-a a],[-a a],'m')
     hold off
     axis equal off
     title(sprintf('After %1d Smoothings',k),'FontWeight','bold','FontSize',14)
     if k<=20
         pause(.15)
     end
     shg
  end

  function [semiaxis45,semiaxis135] = Ellipse(x,y)
% Determines the semiaxes of the limiting ellipse when the
% polygon smoothing process is applied to the vectors x and y
% which are assumed to have unit length and centroid zero.
% The ellipse is tilted 45 degrees and the semiaxes along the
% 45 and 135 lines are returned in semiaxis45 and semiaxis135 resp.

  n = length(x);
  [cVec,sVec] = csVecs(n);
  alfa = (cVec'*x)/(n/2);     
  beta = (sVec'*x)/(n/2);     
  gamma = (cVec'*y)/(n/2);    
  delta = (sVec'*y)/(n/2);    
  [Gx,zx] = planerot([alfa;beta]);    cx = Gx(1,1); sx = Gx(1,2);
  [Gy,zx] = planerot([delta;-gamma]); cy = Gy(1,1); sy = Gy(1,2);
  theta_x = atan2(sx,cx);
  theta_y = atan2(sy,cy);
  a = 2*cos((theta_x - theta_y)/2 - pi/4)/sqrt(n);
  b = -2*sin((theta_x - theta_y)/2 - pi/4)/sqrt(n);
  semiaxis45  = abs(a);
  semiaxis135 = abs(b);
   


  function [xNew,yNew] = Smooth(x,y)
% x and y are column n-vectors that define the vertices of a polygon P. 
% xNew and yNew are column n-vectors that define a polygon Q
% whose edges are the midpoints of P's edges.

  n = length(x); xNew = zeros(n,1); yNew = zeros(n,1);
  for i=1:n-1
     xNew(i) = (x(i) + x(i+1))/2;
     yNew(i) = (y(i) + y(i+1))/2;
  end
  xNew(n) = (x(n)+x(1))/2;  xNew = xNew-sum(xNew)/n; xNew = xNew/norm(xNew);
  yNew(n) = (y(n)+y(1))/2;  yNew = yNew-sum(yNew)/n; yNew = yNew/norm(yNew);