  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);
 
% The slope and size of the ellipse is determined by four numbers,
% the coefficients of cVec and sVec in the initial x and y vectors.
% All other components of x and y are damped out by the power iteration.

% So let's get these four numbers using the facts we discovered
% about cVec and sVec...

% First, compute alfa and beta so

%          x = alfa*cVec + beta*sVec + (vector orthogonal to cVec and sVec)

  alfa = (cVec'*x)/(n/2);     % Recall n/2 = cVec'*cVec
  beta = (sVec'*x)/(n/2);     % Recall n/2 = sVec'*sVec
  
% Likewise         

%          y = gammma*cVec + delta*sVec + (vector orthogonal to cVec and sVec)

  gamma = (cVec'*y)/(n/2);    % Recall n/2 = cVec'*cVec
  delta = (sVec'*y)/(n/2);    % Recall n/2 = sVec'*sVec
  
% Let xf be a unit vector in the direction of  alfa*cVec +  beta*sVec 
% Let yf be a unit vector in the direction of gamma*cVec + delta*sVec
% These vectors define points that are on the limiting ellipse. We need
% to discover the semiaxes of this ellipse and its tilt. (We believe the
% tilt is 45 degrees.)

% The key is to parameterize xf and yf in a "trigonometric way.'

  [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);
  
% Since 
%            [alfa;beta]    = Gx'*[zx;0] = [cx;sx]*zx
%            [delta;-gamma] = Gy'*[zy;0] = [cy;sy]*zy 

% it follows that
 
%     xf = unit vector in direction of alfa*cVec + beta*sVec
%        = (cx*cVec + sx*sVec)/sqrt(n/2)
%     yf = unit vector in direction of gamma*cVec + delta*sVec
%        = (-sy*cVec + cy*sVec)/sqrt(n/2)

% Let theta = (0:n-1)'*(2pi/n), the vector of angles in cVec and sVec.
% Let theta_x = angle associated with cx and sx
% Let theta_y = angle associated with cy and sy
% It follows that
 
%      xf = (cos(theta_x)cos(theta) + sin(theta)sin(theta_x))/sqrt(n/2) 
%         = cos(theta-theta_x)/sqrt(n/2)

%      yf = (-sin(theta_y)cos(theta) + cos(theta_y)sin(theta))/sqrt(n/2) 
%         = sin(theta-theta_y)/sqrt(n/2)
%         = cos(theta-theta_y - pi/2)/sqrt(n/2)

% Since we conjecture that these points lie of a 45-degree tilted ellipse,
% let us rotate the point set by -45 degrees to get an axis-aligned
% ellipse. This means applying the rotation [1 1;-1 1]/sqrt to each point 
% giving

%   xfR = ( cos(theta-theta_x) + cos(theta - theta_y - pi/2) )/sqrt(n)
%   yfR = (-cos(theta-theta_x) + cos(theta - theta_y - pi/2) )/sqrt(n)

% We now use some trig identities which are based on familiar formulas for
% cos(a+b) and sin(a+b). In particular, if u and v are given, w = (u+v)/2,
% and del = (v-u)/2, then

%      cos(u) + cos(v) =  2cos(w)cos(del)
%     -cos(u) + cos(v) = -2sin(w)sin(del)

% Setting u = theta - theta_x, v = theta - theta_y - pi/2 we get

%    xfR = 2*cos((u+v)/2)*cos((v-u)/2)/sqrt(n) = a*cos(phi)
%    yfR = -2sin((u+v)/2)*sin((v-u)/2)/sqrt(n) = b*sin(phi)

% where 
%             a =  2*cos((theta_x - theta_y)/2 - pi/4)/sqrt(n)
%             b = -2*sin((theta_x - theta_y)/2 - pi/4)/sqrt(n)
%           phi = theta - (theta_x + theta_y)/2 - pi/4

% It follows that the points defined by xfR and yfR are on an axis-aligned
% ellipse with semiaxis in the x-direction given by abs(a) and semi-axis in
% the y-direction given by abs(b). "Unrotating", this means that the points
% defined by xf and yf are on an ellispe with semiaxis abs(a) along the 45
% degree line and semiaxis abs(b) along the 135-degree line

% Let's get a and b from cx = cos(theta_x), sx = sin(theta_x), 
% cy = cos(theta_y), and sy = sin(theta_y)...

   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);
   
