function xyw = gqintdisk(deg,cent,r)

% computes a cubature formula of polynomial exactness degree deg
% on the intersection of an arbitrary number of planar disks

% by Alvise Sommariva and Marco Vianello, May 29 2015

% WARNING: the algorithm may fail when two "vertices" of the intersection
% can be connected clockwise or counterclockwise by more than one arc

% INPUT
% deg: degree of polynomial exactness
% cent: 2-column array of the centers coordinates
% r: 1-column array of the radii

% OUTPUT
% xyw: 3 columns array of the nodes coordinates and weights


% FUNCTION BODY

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

% 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));

if m==1
% one disk is included in all the others
xyw = gqcircsegm(deg,pi,r(1));
xyw(:,1)=xyw(:,1)+cent(1,1);
xyw(:,2)=xyw(:,2)+cent(1,2);

else
if m==2
% there are only two distinct disks
xyw = gqlens(deg,cent(1,1),cent(1,2),r(1),cent(2,1),cent(2,2),r(2));

else
% computing pairwise disk intersection points
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) i j];
else
pts=[pts;xout' yout' [i j;i j]];
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;

end;

if isempty(vert) || length(vert(:,1))==1
% zero measure intersection
xyw(1,1:2)=cent(1,:);xyw(1,3)=0;

else

% ordering the vertices
% ord=convhull(vertord(:,1),vertord(:,2));
% vertord=vertord(ord,:);
ord=convhull(vert(:,1),vert(:,2));
vertord=vert(ord,:);

% barycenter of the vertices
L=length(vertord(1:end-1,1));
bary=sum(vertord(1:end-1,1:2))/L;

xyw=[ ];
for i=1:length(vertord(:,1))-1
% circle arc connecting two consecutive vertices
if (vertord(i,3)==vertord(i+1,3) || vertord(i,3)==vertord(i+1,4))
    j=vertord(i,3);
end;
if (vertord(i,4)==vertord(i+1,4) || vertord(i,4)==vertord(i+1,3))
    j=vertord(i,4);
end;

% angular interval corresponding to the consecutive vertices
theta1=angle(sqrt(-1)*(vertord(i,2)-cent(j,2))+(vertord(i,1)-cent(j,1)));
theta2=angle(sqrt(-1)*(vertord(i+1,2)-cent(j,2))+(vertord(i+1,1)-cent(j,1)));
alpha=theta1;
if theta1 < theta2
    beta=theta2;
else
    beta=2*pi+theta2;
end;

% computation of the nodes and weights via union of asymmetric sectors
A=[r(j) 0;0 0];
B=[0 r(j);0 0];
C=[cent(j,:);bary];
xyw = [xyw;gqellblend(deg,A,B,C,alpha,beta)];

end;
end;
end;
end;

end


