function [xx,yy]=calcExtHull(px,py,R,Lx,Ly,ift,jft,xx,yy,Ns,Nc)
% calcExtHull <Calculate exterior hull of Minkowski sum of circle with polygon (px,py).>
% Usage:: [xx,yy]=calcExtHull(px,py,R,[Lx],[Ly])
%

% revision history:
% 02/13/19 Mark D. Shattuck <mds> calcExtHull.m

%% Parse Input
if(~exist('Lx','var') || isempty(Lx))
  isper=false;
else
  isper=true;
end

if(~exist('ift','var') || isempty(ift))
  ift=[2:Ns 1];    % used to compute circular diff
end

if(~exist('jft','var') || isempty(ift))
  jft=[Ns 1:Ns-1]; % used to compute circular diff
end

if(~exist('xx','var') || isempty(xx))
  xx=zeros(2*Ns,Nc);
  yy=zeros(2*Ns,Nc);
end

if(~exist('Ns','var') || isempty(Ns))
  [Ns Nc]=size(px);
end

%% Main
lx=px(ift,:)-px;  % circular diff
ly=py(ift,:)-py;

if(isper)
  lx=lx-Lx*round(lx/Lx);
  ly=ly-Ly*round(ly/Ly);
end

lk=sqrt(lx.^2+ly.^2);

nx=ly./lk;
ny=-lx./lk;
kc=ny.*nx(jft,:)-nx.*ny(jft,:);

for m=1:Nc;
  xx(2:2:end,m)=nx(:,m);
  xx(1:2:end,m)=nx(jft,m);
  yy(2:2:end,m)=ny(:,m);
  yy(1:2:end,m)=ny(jft,m);
  
  ii=find(kc(:,m)<-1e-6);
  if(~isempty(ii))
    lxn1=lx(ii,m)./lk(ii,m);
    lyn1=ly(ii,m)./lk(ii,m);
    lxnN=lx(jft(ii),m)./lk(jft(ii),m);
    lynN=ly(jft(ii),m)./lk(jft(ii),m);
    C=lxn1.*lxnN+lyn1.*lynN;
    S=sqrt(1+max(C-.5));
    ll=1./sqrt((1-C))./S;
    xx(2*ii,m)=(lxn1-lxnN).*ll;
    xx(2*ii-1,m)=xx(2*ii,m);
    
    yy(2*ii,m)=(lyn1-lynN).*ll;
    yy(2*ii-1,m)=yy(2*ii,m);
  end
  xx(:,m)=px([1 1 fix(2:.5:Ns) Ns],m)+R(m).*xx(:,m);
  yy(:,m)=py([1 1 fix(2:.5:Ns) Ns],m)+R(m).*yy(:,m);
end

