function curve =  primitive(sType, varargin)
% curve = PRIMITIVE(shape, ...) creates a curve object centered at the origin.
% shape is either 'circle', 'ellipse', or 'cookie'.
% When shape is 'circle', the next argument is the radius
% When shape is 'ellipse', the arguments are the X-semiaxis and the ratio of 
%   the Y-semiaxis to the X-semiaxis.
% When shape is 'cookie', the shape is defined in polar coordinates 
%   by the formula: r(phi) = R*(1+m*cos(n*phi))
%   The order of arguments is R, m, and n.
switch lower(sType)
    case 'circle'
        nPoints = 30;
        phi = linspace(0, 2*pi - 2*pi/nPoints, nPoints);
        phi = phi+pi/nPoints;
        R = varargin{1};
        x = R*cos(phi);
        y = R*sin(phi);
        x(end+1) = x(1);
        y(end+1) = y(1);
    case 'ellipse'
        nPoints = 50;
        phi = linspace(0, 2*pi - 2*pi/nPoints, nPoints);
        phi = phi+pi/nPoints;
        R = varargin{1};
        ell = varargin{2};
        x = R*cos(phi);
        y = ell*R*sin(phi);
        x(end+1) = x(1);
        y(end+1) = y(1);
    case 'cookie'
        nPoints = 50;
        phi = linspace(0, 2*pi - 2*pi/nPoints, nPoints);
        phi = phi+pi/nPoints;
        R = varargin{1};
        mod_amp = varargin{2};
        n = varargin{3};
        r = R*(1+mod_amp*cos(n*phi));
        x = r.*cos(phi);
        y = r.*sin(phi);
        x(end+1) = x(1);
        y(end+1) = y(1);
    otherwise
        error('primitive: Unknown primitive type')
        return
end
points = [x; y];
cs = cscvn(points);
N = 500;
vT = linspace(cs.breaks(1), cs.breaks(end), N);
mF = ppval(cs, vT);
xline = mF(1,:);
yline = mF(2,:);
dt = (vT(2)-vT(1));
dydt = diff(yline)/dt;
dxdt = diff(xline)/dt;
curve.length = sum(sqrt(dydt.^2 + dxdt.^2))*dt;
curve.iCurve = 1;
curve.cs = cs;
curve.max_x = max(xline);
curve.max_y = max(yline);
curve.min_x = min(xline);
curve.min_y = min(yline);
curve.xc = 0;
curve.yc = 0;