
function [polar_interval_output,disk_sequence_output,centers,rs] = ...
    determine_intersectiondisks_boundary(cent,r)

%----------------------------------------------------------------------------------------
% OBJECT
% Determines the boundary
% on the intersection of an arbitrary number of planar disks,
% based on subperiodic Gaussian quadrature
%----------------------------------------------------------------------------------------
% AUTHORS
% by Alvise Sommariva and Marco Vianello, June 3 2015
% modified on September 1, 2022.
%----------------------------------------------------------------------------------------
% INPUT
% cent: 2-column array of the centers coordinates
% r: 1-column array of the radii
% deg: degree of polynomial exactness
%----------------------------------------------------------------------------------------
% OUTPUT
% xyw: 3-column array of the nodes coordinates and weights
%----------------------------------------------------------------------------------------
% IMPORTANT: 
% cent must be an N x 2 matrix and r must be a N x 1 vector. 
%----------------------------------------------------------------------------------------

% FUNCTION BODY

r=r(:);

m=length(r); vert=[ ]; tol=10^(-10);

polar_interval_output=[];
disk_sequence_output=[];
centers=[];
rs=[];

% ... inclusion test ...
still=(1:m);i=1;
while length(still)>2 & (i<length(still) || i==length(still))
    theta=angle(sqrt(-1)*(cent(:,2)-cent(i,2))+(cent(:,1)-cent(i,1)));
    incl=(cent(:,1)-cent(i,1)+r(i)*cos(theta)).^2 + ...
        (cent(:,2)-cent(i,2)+r(i)*sin(theta)).^2 < r.^2+tol;
    ind=find(incl==1);
    ind=setdiff(ind,i);
    
    if length(ind)>0
        still=setdiff(still,ind);
        cent=cent(still,:);
        r=r(still);
        still=(1:length(still));
        i=1;
    else
        i=i+1;
    end
    
end

% ... number of distinct disks ...
m=length(cent(:,1));

centers=cent;
rs=r;

if (m == 1)
    polar_interval_output=[0 2*pi];
    disk_sequence_output=1;
    return;
end

% computing pairwise disk intersection points when m > 1.
pts=[ ];
for i=1:m-1
    for j=i+1:m
        [xout,yout]=...
            circcirc(cent(i,1),cent(i,2),r(i),cent(j,1),cent(j,2),r(j));
        if isnan(xout(1))
            pts=[ ];
            break
        end
        if xout(1)==xout(2) & yout(1)==yout(2)
            % the two circles are tangent
            pts=[pts;xout(1) yout(1)];
        else
            pts=[pts;xout' yout'];
        end
    end
end

if isempty(pts)
    vert=[ ];
    
else
    
    % tiny imaginary parts may seldom appear
    pts(:,1:2)=real(pts(:,1:2));
    
    % selecting the vertices of the circle intersection
    k=1;
    for i=1:length(pts(:,1))
        s=sum((pts(i,1)-cent(:,1)).^2+(pts(i,2)-cent(:,2)).^2 < r.^2+tol);
        if s==m
            vert(k,:)=pts(i,:);k=k+1;
        end
    end
    
    % eliminating multiple points
    vert = singlepts(vert,12);
    
end

if isempty(vert) || length(vert(:,1))==1
    % zero measure intersection
    % xyw(1,1:2)=cent(1,:);xyw(1,3)=0;
    
else
    
    % ordering counterclockwise the vertices
    if length(vert(:,1))==2
        vertord=[vert;vert(1,:)];
    else
        ord=convhull(vert(:,1),vert(:,2));
        vertord=vert(ord,:);
    end
    
    % circles passing through a vertex
    circles=zeros(length(vertord(:,1)),length(cent(:,1)));
    for i=1:length(vertord(:,1))
        circles(i,:)=abs(((vertord(i,1)-cent(:,1)).^2 + ...
            (vertord(i,2)-cent(:,2)).^2  - r.^2))<tol;
    end
    
    % barycenter of the vertices
    L=length(vertord(1:end-1,1));
    bary=sum(vertord(1:end-1,1:2))/L;
    
    % computation of the nodes and weights
    for i=1:length(vertord(:,1))-1
        % minimal length arc connecting (counterclockwise) two consecutive 
        % vertices
        ind=find(circles(i,:)==1 & circles(i+1,:)==1);
        theta1=angle(sqrt(-1)*(vertord(i,2)-cent(ind,2))+...
            (vertord(i,1)-cent(ind,1)));
        theta2=angle(sqrt(-1)*(vertord(i+1,2)-cent(ind,2))+ ...
            (vertord(i+1,1) - cent(ind,1)));
        ind2=find(theta2<theta1);
        alpha=theta1;
        beta=theta2;
        beta(ind2)=2*pi+beta(ind2);
        [y,k]=min((beta-alpha).*r(ind));
        j=ind(k);
        
        
        polar_interval_output=[polar_interval_output; alpha(k) beta(k)];
        disk_sequence_output=[disk_sequence_output; j];
        
    end
end





function spts = singlepts(pts,dig)

% eliminates from an array of points pts those having dig digits 
% in common with a subset of other points, keeping only one per subset  

m=10.^(floor(log10(abs(pts))-dig+1));
pts1=round(pts./m).*m;
pts1(find(abs(pts)<eps)) = 0;
spts=unique(pts1,'rows');

