h=powel(varargin)

function h=powel(varargin)
%   POWEr ELlipse: a mirror-symmetrical x-y line defined as 
%       (|x/radx|)^powx+(|y/rady|)^powy = 1
% Combines ellipse, oval, rectangle, caro and boat shapes
%Call:
%       h=powel([radpows,][colnump])
%Input: 
%		radpows = [radx rady powx powy]
%       radx = x-radial size, default 1
%       rady = y-radial size, default radx
%       powx = x-power, default 2
%       powy = y-power, default powx
%
%		colnump = (default {'k',10})
%        LINCOL (char) or  
%        {LINCOL} or 
%        {NUMP} ,   class(NUMP)= 'double', or 
%        {LINCOL,NUMP} or {NUMP,LINCOL} or {NUMP;LINCOL} or {LINCOL;NUMP}
%
%       NUMP > 1 , default 10. 2*NUMP-1 points is in 1-st quadrant
%       Total NUMber of Points is 8*NUMP-7
%
%       LINCOL(default 'k') = LINe/COLor specification for CLINE
%
%       If only NUMP present, then the curve is not plotted, 
%       and the output is 8*NUMP-7 of x-y vertices. 
%Output:
%		h = handle to the object plotted, or vertices [x y]
%
%   Examples
%       powel    - unit black circle (73 points) at the axis origin
%       h=powel([2 3],'r2') ellipse at the origin with halfaxes 2 and 3, 
%       drawn by CLINE(...,'r2'); h = handle
%       h=powel([2 3 11 20],{'b',21}) = blue rectoid
%       h=powel([2 3 1])   black romb
%       powel([2 3 .5],'g') green sharp diamond
%       powel([2 3 .9 1.3],'r') red sharp-obtuse diamond
%       powel([2 3 .5 3])(dataaspectratio 1 1 1)  black y-boat (smooth-sharp diamond)
%       powel([5 1 2 .5])  x-boat
%       xy=powel([5 1 2 .5],{50}) vertices of x-boat, no plot
%       

%	Vassili Pastushenko	 Nov 2007
%==========================================================
%Defaults:
radpows=1; %unit circle
colnump={10,'k'};
%Parse input
NAR=numel(varargin);    
        for i=1:NAR
            locvar=varargin{i};
            locID=class(locvar);
            if strmatch(locID,'double','exact')
             radpows=locvar;
            end
            if strmatch(locID,'cell','exact')
             colnump=locvar;
            elseif strmatch(locID,'char','exact')
             colnump={locvar};
            end
        end
     
 %Default settings
 nrad=numel(radpows);
 switch nrad
     case 1
         radpows=[radpows*[1 1] 2 2]; %circle
     case 2
         radpows=[radpows 2 2]; %ellipse
     case 3
         radpows=radpows([1:3,3]);%rectoid
 end
 %Check powx & powy
 if any(~isreal(radpows(3:4)))|any(radpows(3:4))<=0
     error('powx and powy should be real and positive')
 end
 
 %parse/deal colnump
 ncol=numel(colnump);
  for i=1:ncol
     locvar=colnump{i};
            locID=class(locvar);
            if strmatch(locID,'double','exact')
                NUMP=locvar;
            end
            if strmatch(locID,'char','exact')
                col = locvar;
            end
  end
 
 NUM=exist('NUMP');
 COL=exist('col');
 NOPLOT=NUM&~COL;
  
 if ~NUM
     NUMP=10;
 end
 
 if ~COL
     col='5.-k'; %default
 end
   
   %deal radpows
   rx=radpows(1);
   ry=radpows(2);
   pows=radpows(3:4);
   powx=pows(1);
   powy=pows(2);
 

NUMP=2*round(abs(NUMP))-1;
theta=linspace(0,pi/2,NUMP)';   

xbas=cos(theta);%.^(1/sqrt(powx*powy));
proxX=xbas(1:2:NUMP);
proyY=xbas(2:2:NUMP);
proyX=(1-proxX.^powx).^(1/powy);
proxY=(1-proyY.^powy).^(1/powx);
z=sortrows([-abs(rx)*[proxX;proxY],abs(ry)*[proyX;proyY]],1);
x=-z(:,1);
y=z(:,2);

ind=1:NUMP-1;
ind2=2:NUMP;
X=[x;-x(NUMP-1:-1:2);-x;x(NUMP-1:-1:1)];
Y=[y;y(NUMP-1:-1:2);-y;-y(NUMP-1:-1:1)];
z=[X Y];

if  NOPLOT
    h=z;
    return
else
h=cline(z,col);
end