  function ShowAveraging()
% Illustrates functions use to examine repeated polygon smoothing.
close all
n = 20;                         % The number of polygon edges
speedy = 1;                     % Setting this to zero slows the simulation
its = 0;                        % Keeps track of the number of averagings
[x,y,u,v,d,itMax] = GenEG(n);   % Generate an example
Show(x,y,its,d,u,v)             % Display the polygon. 
pause(1)
close
[x,y] = AverageOnce(x,y);       % Average once
its = its+1;
Show(x,y,its,d,u,v)
pause(1)

close
[x,y] = AverageOnce(x,y);       % Average again
its = its+1;
Show(x,y,its,d,u,v)
pause(1)

[x,y] = Simulate(x,y,d,its,itMax,speedy,u,v);  % Repeatedly average

% To suppress the display of the limiting ellipse, change the each
% Show(x,y,its,d,u,v)  to Show(x,y,its,d) and change  the last line
% to [x,y] = Simulate(x,y,d,its,itMax,speedy);
 
% For your convenience, the following functions are also available in
% UntangleGUI.m.
% It is recommended that you study how AverageOnce, Simulate, and Show
% work. GenEg involves more advanced math--do not spend time on how
% it works.

  function [x,y] = Simulate(x0,y0,d,k,itMax,speedy,u,v)
% x0 and y0 are column n-vectors that define the vertices of a normalized
%    polygon P.
% k is a positive integer that represents the number of averagings that
%    have been performed on some original polygon to obtain P.
% itMax a positive integer >=k.
% x and y are column n-vectors that define the vertices of a normalized
%    polygon that is the result of itMax averagings of the same original
%    polygon.
% u and v are vectors with the property that plot(u,v) displays the 
%    limiting ellipse.
% speedy has the value of 1 if a fast simulation is desired, 0 otherwise.
% d is a real number with the property that the simulation takes place in
%    the square with vertices (d,d), (-d,d), (-d,-d), (d,-d).

% Set up the time delay based on the value of speedy...
if speedy
    tDelay = 0.01;
else
    tDelay = .15;
end
% The simulation...
x = x0;
y = y0;
for it=k+1:itMax
   % Average and display..
   [x,y] = AverageOnce(x,y);
   if nargin==6
      % Do not include the limiting ellipse...
      Show(x,y,it,d)
   else
      % Include the limiting ellipse...
      Show(x,y,it,d,u,v)
   end
   pause(tDelay)
end
  


  function Show(x,y,k,d,u,v)
% Displays the polygon whose vertices are defined by the column n-vectors
%    x and y. 
% k is an integer, typically the number of completed averaging steps.
% d is a positive real. The x and y ranges will be -d to +d.
% If nargin==6, then the u and v are vector of the limiting ellipse and it
%   will be displayed. 
% A call of the form Show(x,y,k,d) is also valid but does not display the
%   limiting ellipse.
          
plot([x;x(1)],[y;y(1)],'r',x,y,'.r','Markersize',20)
axis([-d d -d d])
axis  square off
title(sprintf('n = %2d       Number of Averagings = %1d',length(x),k),'FontSize',18)
hold on
plot([d d -d -d d],[-d d d -d -d],'linewidth',10)
if nargin==6
   plot(u,v,'k','linewidth',2')
   plot([-d d],[-d d],'--k','linewidth',2)
   plot([-d d],[d -d],'--k','linewidth',2)
end
hold off
    

  function [xNew,yNew] = AverageOnce(x,y)
% x and y are column n-vectors that define the vertices of a normalized 
%    polygon P.
% xNew and yNew are column n-vectors that define the vertices of a normalized
%    polygon Q obtained by connecting the midpoints of P's edges.
n = length(x); 
xNew = (x + [x(2:n);x(1)])/2; xNew = xNew - mean(xNew); xNew = xNew/norm(xNew);
yNew = (y + [y(2:n);y(1)])/2; yNew = yNew - mean(yNew); yNew = yNew/norm(yNew);

  function [x,y,u,v,d,itMax] = GenEG(n)
% n is an integer >= 4. 
% x and y are column n-vectors that define the vertices of a normalized polygon.
% u and v are vectors with the property that plot(u,v) displays the
%    limiting ellipse.
% d is a positive real with the property that the simulation takes
%    place in the square with vertices (d,d), (-d,d), (-d,-d), (d,-d).
% itMax is a positive integer with the property that with that many
%    iterations, the polygon has visiably converged to the limiting ellipse.

% x and y vectors define the normalized random polygon P...
x = rand(n,1); x = x-mean(x); x = x/norm(x);
y = rand(n,1); y = y-mean(y); y = y/norm(y);
% u and v vectors...
tau = (2*pi/n)*(0:n-1)'; 
c = sqrt(2/n)*cos(tau); 
s = sqrt(2/n)*sin(tau);
cos_theta_u = (c'*x)/sqrt((c'*x)^2+(s'*x)^2);
sin_theta_u = (s'*x)/sqrt((c'*x)^2+(s'*x)^2);
cos_theta_v = (c'*y)/sqrt((c'*y)^2+(s'*y)^2);
sin_theta_v = (s'*y)/sqrt((c'*y)^2+(s'*y)^2);
R =  sqrt(2/n)*[cos_theta_u sin_theta_u ; cos_theta_v sin_theta_v];
t = linspace(0,2*pi);
M = R*[cos(t);sin(t)];
u = M(1,:)';
v = M(2,:)';
% Simulation parameters...
d = 1.6*max(max(u),max(v));
itMax = ceil(3./(-log10(  1-(3/2)*(pi/n)^2 )));