  function PolyAve()
% Displays the averaging process N times for the each of the problem sizes
% specifed by probSizes
clc
N=3; 
probSizes = [20 30 40];
for eg=1:N
  for n = probSizes
    % Initializations...
    close all
    % Generate and display a random  polygon
    [x,y] = GenEG(n);
    [u,v] = D2Ellipse(x,y);
    d = 1.6*max(max(u),max(v));
    % Show the initial polygon...
    figure
    set(gcf,'position',[250 200 900 600])
    hold on
    plot([x;x(1)],[y;y(1)],'r',x,y,'.r','Markersize',20)
    axis([-d d -d d])
    axis  square off
    title(sprintf('A Random Polygon  (n = %2d )',n),'FontWeight','bold','FontSize',24)
    pause(2)
    % Show the limiting ellipse...
    plot(u,v,'k','linewidth',2')
    hold off
    % Repeatedly smooth and display the polygon...
    
    % Approximate damping factor...
    rho_n = 1-(3/2)*(pi/n)^2;
    % rho_n^m <=.01  m*log(rho) <= -3 
    m = ceil(2./(-log10(rho_n)));
    for k=1:m
        [x,y] = Smooth(x,y);
        plot(u,v,'k','linewidth',2')
        hold on
        plot([x;x(1)],[y;y(1)],'r',x,y,'.r','Markersize',20)
        hold off
        axis([-d d -d d])
        axis  square off
        title(sprintf('n = %2d   Averagings = %1d',n,k),'FontWeight','bold','FontSize',24)
        shg
    end
    hold on manual
    d = .9*d;
    plot([-d d],[-d d],'--k','linewidth',2)
    plot([-d d],[d -d],'--k','linewidth',2)
    shg
    pause(2)
  end
end


  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.

% Initializations...
n = length(x); xNew = zeros(n,1); yNew = zeros(n,1);
for i=1:n-1
    % i-th new vertex is a midpoint...
    xNew(i) = (x(i) + x(i+1))/2;
    yNew(i) = (y(i) + y(i+1))/2;
end
% Special case...
xNew(n) = (x(n)+x(1))/2; 
yNew(n) = (y(n)+y(1))/2; ;
xNew = xNew - mean(xNew); 
xNew = xNew/norm(xNew);
yNew = yNew - mean(yNew); 
yNew = yNew/norm(yNew);

  function [x,y] = GenEG(n)
% n>=3 and x and y are column n-vectors with unit 2-norm that define the 
% vertices of a random polygon with centroid (0,0). 
% 
x = rand(n,1); 
x = x-mean(x); 
x = x/norm(x);
y = rand(n,1); 
y = y-mean(y); 
y = y/norm(y);

function [u,v] = D2Ellipse(x0,y0)
% plot(u,v) will disply the limiting ellipse
n = length(x0);
tau = (2*pi/n)*(0:n-1)'; c = sqrt(2/n)*cos(tau); s = sqrt(2/n)*sin(tau);
cos_theta_u = (c'*x0)/sqrt((c'*x0)^2+(s'*x0)^2);
sin_theta_u = (s'*x0)/sqrt((c'*x0)^2+(s'*x0)^2);
cos_theta_v = (c'*y0)/sqrt((c'*y0)^2+(s'*y0)^2);
sin_theta_v = (s'*y0)/sqrt((c'*y0)^2+(s'*y0)^2);

t = linspace(0,2*pi,300);
M = sqrt(2/n)*[cos_theta_u sin_theta_u ; cos_theta_v sin_theta_v]*...
    [cos(t);sin(t)];
u = M(1,:)';
v = M(2,:)';
