## Copyright (C) 2001, James B. Rawlings and John W. Eaton
%%
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License as
## published by the Free Software Foundation; either version 3, or (at
## your option) any later version.
%%
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
%%
## You should have received a copy of the GNU General Public License
## along with this program; see the file COPYING.  If not, write to
## the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
## MA 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} ellipse (@var{a}, @var{level})
## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{shift})
## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{shift}, @var{n})
## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{shift}, @var{n}, @dots{})
## @deftypefnx{Function File} address@hidden =} ellipse (@dots{})
## @deftypefnx{Function File} address@hidden, @var{y}, @var{major}, @var{minor}, @
## @var{bbox}] =} ellipse (@dots{})
##
## Plot an ellipse with principal axes given the Eigen vectors of a 2 by 2 matrix.
##
## The Eigen vectors of @var{a} multiplied by the corresponding Eigen values and
## @var{level} are used as the principal axes of the plotted ellipse. By default
## @var{level} is 1. If it is a vector, several ellipses are plotted.
##
## The optional argument @var{shift} controls the center of the ellipse. By default
## this is a two-vector of zeros.
##
## The optional argument @var{n} controls the number of points used to plot the
## ellipse. By default 100 points are used.
##
## The rest of the arguments are passed to @code{plot}, and can be used to control
## the visual style of the plot.
##
## If one output argument is requested, a graphics handle of the plot is returned.
##
## If more than one output argument is requested, no plot is generated. Instead
## the data used to generate the plot is returned. @var{x} and @var{y} are
## @var{n}-vectors containing the points on the ellipse. @var{major} and @var{minor}
## contains the axes of the ellipse, and @var{bbox} contains the bounding box of
## the ellipse.
##
## The following example generates non-isotropic normally distributed data, and
## plots the ellipses with Mahalanobis' distance of 1 and 3 to the center.
##
## @example
## ## Generate data
## X = randn (100, 2);
## scale = 2;
## X (:, 1) *= scale;
## theta = 0.4;
## R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
## X = X * R;
##
## ## Estimate mean and covariance matrix
## mu = mean (X);
## C = cov (X);
##
## ## Plot data and ellipses
## figure
## plot (X (:, 1), X (:, 2), 'ro');
## hold on, ellipse (inv (C), [1, 3], mu, :, "k"); hold off
## axis equal
## @end example
##
## @noindent
## Note the use of @code{:} to denote the default value of @var{n}.
## @seealso{plot}
## @end deftypefn

function varargout = ellipse (amat, level = 1, shift = [0, 0], n = 100, varargin)
  ## Check input
  if (nargin < 1)
    error ("ellipse: not enough input arguments");
  endif
  
  if (!isnumeric (amat) || !isequal (size (amat), [2, 2]))
    error ("ellipse: first input argument must be a 2x2 matrix");
  endif
  
  if (!isvector (level))
    error ("ellipse: second input argument must be a vector");
  endif
  
  if (!isnumeric (shift) || numel (shift) != 2)
    error ("ellipse: third input argument must be a two-vector");
  endif
  shift = shift (:).';
  
  if (!isscalar (n) || n < 0 || n != round (n))
    error ("ellipse: fourth input argument must be a positive integer");
  endif
  
  ## We're ready to do the actual work
  [v, l] = eig (amat);
  
  dl = diag (l);
  if (any (imag (dl)) || any (dl <= 0))
    error ("ellipse: first input argument must be positive definite");
  endif

  ## Generate contour data.
  a = 1 / sqrt (dl (1));
  b = 1 / sqrt (dl (2));

  a = a * level (:);
  b = b * level (:);
    
  t = linspace (0, 2*pi, n);

  xt = a * cos (t);  xt = xt.';
  yt = b * sin (t);  yt = yt.';

  ## Rotate the contours.
  ra = atan2 (v (2, 1), v (1, 1));

  cos_ra = cos (ra);
  sin_ra = sin (ra);

  x = xt * cos_ra - yt * sin_ra + shift (1);
  y = xt * sin_ra + yt * cos_ra + shift (2);

  ## Endpoints of the major and minor axes.
  [mx, ii] = max (level);
    
  minor = (v * diag ([a(ii), b(ii)])).';
  major = minor;

  major (2, :) = -major (1,:);
  minor (1, :) = -minor (2,:);

  t = [1; 1] * shift;

  major = major + t;
  minor = minor + t;

  ## Bounding box for the ellipse using magic formula.
  ainv = inv (amat);
  xbox = level (ii) * sqrt (ainv (1,1));
  ybox = level(ii) * sqrt (ainv (2,2));

  bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];

  t = [1; 1; 1; 1; 1] * shift;
  bbox = bbox + t;
  
  ## Should we plot the ellipse, or return it?
  if (nargout <= 1)
    handle = plot (x, y, varargin {:}, major (:,1), major (:,2), varargin {:},
                   minor (:,1), minor (:,2), varargin {:});
    if (nargout == 1)
      varargout {1} = handle;
    endif
  else
    varargout {1} = x;
    varargout {2} = y;
    varargout {3} = major;
    varargout {4} = minor;
    varargout {5} = bbox;
  endif
endfunction

%!demo
%! ## Generate data
%! X = randn (100, 2);
%! scale = 2;
%! X (:, 1) *= scale;
%! theta = 0.4;
%! R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
%! X = X * R;
%!
%! ## Estimate mean and covariance matrix
%! mu = mean (X);
%! C = cov (X);
%!
%! ## Plot data and ellipses
%! figure
%! plot (X (:, 1), X (:, 2), 'ro');
%! hold on, ellipse (inv (C), [1, 3], mu, :, "k"); hold off
%! axis equal