function [qx,qy]=cPolygon(px,py,D,N)
% cPolygon <Curved Polygon>
% Usage:: [qx,qy]=cPolygon(px,py,D[inf],N[1])
% px,py uncurved polygon verticies
% D diameter of side curvature
% N number of points on curved side
% sign(D) determines convex/concave
% sign(N) determines long/short
%

% revision history:
% 02/28/09 Mark D. Shattuck <mds> cPolygon.m
% 06/25/21 mds added comments

%% Parse Input
Ns=length(px);
if(~exist('D','var') || isempty(D))
  D=inf*ones(1,Ns);
end

if(~exist('N','var') || isempty(N))
  N=ones(1,Ns);
end

%% Main
% extract sign of D
sD=sign(D);
D=abs(D);

% extract sign of N
sN=N<0;
N=abs(N);

% convert to complex for simplicity
p=px+1i*py;

l=p([2:end 1])-p;                 % vector sides
L=abs(l);                         % length of sides
n=(imag(l)-1i*real(l))./L;        % unit normal to sides
th=sN.*2*pi+sD.*2.*asin(L./D);    % angle of curve
S=sqrt(D.^2-L.^2);                % 2*distance to center of curvature 
nn=[0 cumsum(N)];                 % indices for new polygon
c=(l-sD.*S.*n)/2;                 % vector to center

q(1,nn(end))=1i;                  % initialize output
for ns=1:Ns
  s=(0:N(ns)-1)/N(ns);            % equal spacing for side ns
  if(isinf(D(ns)))                % special case for D->inf    
    q(nn(ns)+1:nn(ns+1))=p(ns)+s*l(ns);
  else
    q(nn(ns)+1:nn(ns+1))=p(ns)+(1-exp(1i*s*th(ns)))*c(ns); % rotate radial
  end                                                      % vector  
end

% convert back to x,y
qx=real(q);
qy=imag(q);

