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 2015

% 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); ic=0;

% inclusion test (checks whether the intersection is one of the disks)
for i=1:m
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;
if sum(incl) == m
ic=i;
break;
end;
end;

if ic>0
% one disk is icluded in all the others
vert=[[cent(ic,:)+r(ic)*[1 0] ic ic];[cent(ic,:)+r(ic)*[-1 0] ic ic]];

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));
pts=[pts;xout' yout' [i j;i j]];
end;
end;

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

% tiny imaginary parts may seldom appear
vert(:,1:2)=real(vert(:,1:2));

% eliminating multiple vertices
if ~isempty(vert)
[v,is]=singlepts(vert(:,1:2),ceil(abs(log10(tol))));
vert=vert(is,:);
end;

if isempty(vert) || length(vert(:,1))==1
% zero measure intersection
xyw(1,1:2)=cent(1,:);xyw(1,3)=0;
    
else
% labelling each vertex with the corresponding pair of circles
vertord(1,:)=vert(1,:);
vert=setdiff(vert,vert(1,:),'rows');
k=1;
while ~(isempty(vert))
ind=find((vertord(k,3)==vert(:,3)));
if isempty(ind)
ind=find((vertord(k,4)==vert(:,4)));
end;
if isempty(ind)
ind=find((vertord(k,3)==vert(:,4)));
end;
if isempty(ind)
ind=find((vertord(k,4)==vert(:,3)));
end;
k=k+1;
vertord(k,:)=vert(ind,:);
vert=setdiff(vert,vert(ind,:),'rows');

end;

% ordering the vertices
if length(vertord(:,1))==2
ord=[1;2;1];  
else
ord=convhull(vertord(:,1),vertord(:,2));
end;
vertord=vertord(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


