function [area axes offset rot] = computeellipse(x,y,fvalue)
% COMPUTEELLIPSE Calculate the confidence ellipse area.
%   COMPUTEELLIPSE(X,Y) calculates the area enclosed by an ellipse that encloses 
%   the samples in X and Y with confidence corresponding to 95%, 
%   FVALUE (alpha = 0.05, 2, n-2) = 3.00. 
%
%   COMPUTEELLIPSE(X,Y,FVALUE) Supply an alternate FVALUE for a confidence level 
%   different than 95%. This value is determined by looking up the FVALUE for a 
%   bivariate distribution with a significantly large number of data points 
%   (n > 120). 
%   
%   [AREA AXES OFFSET ROT] = COMPUTEELLIPSE(...) additionally returns the major 
%   and minor axes, AXES(1) and AXES(2), respectively, an OFFSET from the 
%   origin, and a rotation matrix, ROT.
%
%   Input
%       X      | m x 1 | double | m is the number of samples
%       Y      | m x 1 | double |
%       FVALUE | 1 x 1 | double | optional
%
%   Output
%       AREA   | 1 x 1 | double | 
%       AXES   | 1 x 2 | double | [major minor] axis lengths
%       OFFSET | 1 x 2 | double | [x y] offsets from origin
%       ROT    | 2 x 2 | double | rotation matrix with axis scaling
%   
%   Example 
%       x = randn(500,1)-2;
%       y = 0.5*randn(500,1)+0.25*x;
%       [area axes offset rot] = computeellipse(x,y);
%       theta  = linspace(0,2*pi,100)';
%       ellipse = [cos(theta) sin(theta)]*rot;
%       ellipse = ellipse+repmat(offset,size(ellipse,1),1);
%       plot(x,y,'.',ellipse(:,1),ellipse(:,2));
%       axis equal;
%
%   Reference
%     [1] Prieto TE, Myklebust JB, Hoffmann RG, Lovett EG, Myklebust BM. 
%         Measures of postural steadiness: differences between healthy young and 
%         elderly adults. IEEE Trans Biomed Eng 1996;43(9): 956–66.
%
%   See also 

% Created by: Jesse Norris
%             Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://jesse.musculoskeletaldynamics.com
% Last Revision: 13 June 2006

% Check sizes of inputs
  [m n]   = size(x);
  [my ny] = size(y);    
  
% Verify inputs are suitable
  if any([n ny] ~= 1); error('Inputs should be column vectors'); end;
  if (m ~= my); error('Number of samples of <x> and <y> disagree'); end;    
  
% If FVALUE is not passed in, assume default of 3.0
  if ~exist('fvalue','var'); fvalue = 3.0; end;
  
% Mean zero data & store the offset
  offset = mean([x y],1);
  x = x-offset(1); 
  y = y-offset(2);     
  
% Compute covariance matrix   
  covar = [x y]'*[x y]/(m-1);
  
% Calculate eigenvectors & eigenvalues of covariance matrix
  [vecs,eigs] = eig(covar);
    
% Determine the axes of the ellipse
  axes = sqrt(2*fvalue*eigs);
  rot  = axes*vecs;
    
% Sort such that axes(1) = major, axes(2) = minor
  axes = sort(diag(axes),1,'descend')';               
  
  %{ 
    % ALTERNATE METHOD TO CALCULATE AXES
    
    % Assign covariance to meaningful variable names
      covxx = covar(1,1);
      covyy = covar(2,2);
      covxy = covar(1,2);              
  
    % Compute major and minor ellipse axes
      D = sqrt((covxx-covyy)^2+4*covxy^2);
      axes(1) = sqrt(fvalue*(covxx+covyy+D));
      axes(2) = sqrt(fvalue*(covxx+covyy-D));              
  %}
    
% Calculate area
  area = pi*axes(1)*axes(2);          

end