bht_ellipse(t,der,varargin)

function y = bht_ellipse(t,der,varargin)
% BHT_ELLIPSE
% 
% Define an ellipse-shaped boundary.
% 
% Syntax
% 
%     y = BHT_ELLIPSE(t,der,settings)
% 
% Description
% 
%     This function defines an ellipse-shaped boundary as a parameterized
%     couterclockwise oriented closed line over [0,2*pi[. The input
%     variable t (the parameter) must be a vector with all of its elements
%     within [0,2*pi[. The function returns, depending on the flag der:
% 
%           o der = 0: the coordinates of the points corresponding to the
%             parameters t. 
%           o der = 1: the coordinates of the tangent vectors
%             (the first derivatives of the parameterization). 
%           o der = 2: the
%             second derivatives of the parameterization.
% 
%     The input variable settings is a row vector with five elements: 
%       settings(1:2) : ellipse center (x,y) coordinates
%       settings(3:4) : semimajor and semiminor axis' lengths
%       settings(5)   : a flag (-1 or 1) for the orientation of the boundary: 
%                       1 means counterclockwise orientation
%                      -1 means clockwise orientation. 
%
%     We recall that if the ellipse is a link of an articulated body, it must be
%     centered at the origin.
% 
% See also BHT_BOUNDARY_CHECK 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                      BIOHYDRODYNAMICS MATLAB TOOLBOX                    %
%                                                                         %
%                         A. MUNNIER and B. PINCON                        %
%                                                                         %
% [e-mail]        [e-mail]   %
% http://www.iecn.u-nancy.fr/~munnier  http://www.iecn.u-nancy.fr/~pincon %
%                                                                         %
%                                                                         %
%                      INSTITUT ELIE CARTAN, NANCY 1                      %
%                       http://www.iecn.u-nancy.fr/                       %
%                                                                         %
%                      INRIA Lorraine, Projet CORIDA                      %
%                    http://www.iecn.u-nancy.fr/~corida/                  %
%                                                                         %  
%                                                                         %  
%                                                                         %
%                                                                         %
%                            August 15th 2008                             %
%                                                                         %
%                                               GNU GPL v3 licence        %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%==========================================================================
if isempty(varargin)
    error('Boundary''s settings (row vector of five elements) are missing.')
else
settings = varargin{1};
end

if length(settings) ~= 5
    error('Boundary''s settings vector must have five elements');
end

% =========================================================================
% Ellipse's size
x1 = settings(1);  % x-coordinate of the center
x2 = settings(2);  % y-coordinate of the center
a = settings(3);   % semimajor axis
b = settings(4);   % semiminor axis
% orientation
orient = settings(5);
if orient^2 ~=1
    error('Last setting must be 1 (fluid outside) or -1 (fluid inside)')
end

b = b*orient;
%==========================================================================
n = length(t); % number of points
t = reshape(t,1,n); 

%==========================================================================
if der == 0 
  y = [a*cos(t)+x1; b*sin(t)+x2];
elseif der == 1
  y = [-a*sin(t); b*cos(t)];
elseif der == 2
  y = [-a*cos(t); -b*sin(t)];
else
  error('Flag der must be 0,1 or 2')
end