
function xyw = cub_multibubble(deg,cent,r)

%--------------------------------------------------------------------------
% Object:
% computes a quadrature formula of polynomial exactness degree deg on the
% union of a small number of planar disks (say up to 10-12 disks) by
% the inclusion-exclusion principle and subperiodic Gaussian quadrature on
% disk intersections.
%--------------------------------------------------------------------------
% 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.
%--------------------------------------------------------------------------
% Written by Alvise Sommariva and Marco Vianello, June 2015.
% Last update: January 04, 2021.
%--------------------------------------------------------------------------

m=length(r);xyw=[ ];good=(1:m);

for k=1:m
    if isempty(good)
        break
    end
    c=nchoosek(good,k);
    nw=[ ];good1=[ ];
    for j=1:length(c(:,1))
        nnww = gqintdisk(cent(c(j,:),:),r(c(j,:)),deg);
        if length(nnww(:,1))>1
            good1=[good1 c(j,:)];
            nw=[nw;nnww(:,1) nnww(:,2) (-1)^(k-1)*nnww(:,3)];
        end
    end
    
    xyw=[xyw;nw];
    good=unique(good1);
    
end



