
function xyw = cub_dbubble(n,x1,y1,r1,x2,y2,r2)

%--------------------------------------------------------------------------
% Object:
% The routine computes the nodes and weights of a product gaussian formula
% on a planar double bubble, namely the union of two disks with possibly
% different radii.
%--------------------------------------------------------------------------
% Input:
% n: algebraic degree of exactness;
% x1,y1: coordinates of the center of the first disk;
% r1: radius of the first disk;
% x2,y2: coordinates of the center of the second disk;
% r2: radius of the second disk.
%--------------------------------------------------------------------------
% Output:
% xyw: (2 x ceil((n+1)/2) x ceil((n+2)/2)) x 3 array of
%                        (xnodes,ynodes,weights)
%--------------------------------------------------------------------------
% Required routines:
% 1. r_jacobi.m (www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html)
% 2. gauss.m (www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html)
% 3. trigauss.m
%--------------------------------------------------------------------------
% Written by Gaspare Da Fies and Marco Vianello, University of Padova
% Date: December 2, 2011.
% Last update: January 4. 2020.
%--------------------------------------------------------------------------


% distance between the centers
d=sqrt((x1-x2)^2+(y1-y2)^2);
if d>=(r1+r2)
    
    % the disks are disjoint
    % trigonometric gaussian formula on the arc
    tw=trigauss(n+2,-pi,pi);
    
    % algebraic gaussian formula on [-1,1]
    ab=r_jacobi(ceil((n+1)/2),0,0);
    xw=gauss(ceil((n+1)/2),ab);
    
    % creating the grid
    [t,theta]=meshgrid(xw(:,1),tw((1:ceil((n+2)/2)),1));
    [w1,w2]=meshgrid(xw(:,2),tw((1:ceil((n+2)/2)),2));
    
    % nodal cartesian coordinates and weights
    c=cos(theta(:)); s=sin(theta(:));
    xyw1(:,1)=r1*c+x1;
    xyw1(:,2)=r1*t(:).*s+y1;
    xyw1(:,3)=r1^2*s.^2.*w1(:).*w2(:);
    xyw2(:,1)=r2*c+x2;
    xyw2(:,2)=r2*t(:).*s+y2;
    xyw2(:,3)=r2^2*s.^2.*w1(:).*w2(:);
    xyw=[xyw1; xyw2];
    
elseif d<=(max(r1,r2)-min(r1,r2))
    
    % one circle is included in the other
    % quadrature formula on the larger disk
    if(r1>=r2)
        r=r1;
        x=x1;
        y=y1;
    else
        r=r2;
        x=x2;
        y=y2;
    end
    % trigonometric gaussian formula on the arc
    tw=trigauss(n+2,-pi,pi);
    
    % algebraic gaussian formula on [-1,1]
    ab=r_jacobi(ceil((n+1)/2),0,0);
    xw=gauss(ceil((n+1)/2),ab);
    
    % creating the grid
    [t,theta]=meshgrid(xw(:,1),tw((1:ceil((n+2)/2)),1));
    [w1,w2]=meshgrid(xw(:,2),tw((1:ceil((n+2)/2)),2));
    
    % nodal cartesian coordinates and weights
    xyw(:,1)=r*cos(theta(:))+x;
    xyw(:,2)=r*t(:).*sin(theta(:))+y;
    xyw(:,3)=r^2*(sin(theta(:))).^2.*w1(:).*w2(:);
    
else
    
    % the union is a double bubble
    % angles of the circular segments
    omega1=pi-acos((d^2+r1^2-r2^2)/(2*d*r1));
    omega2=pi-acos((d^2+r2^2-r1^2)/(2*d*r2));
    % angle with the x-axis of the line connecting the centers
    if (y2-y1)==0
        phi1=acos((x1-x2)/d);
        phi2=acos((x2-x1)/d);
    else
        phi1=sign(y1-y2)*acos((x1-x2)/d);
        phi2=sign(y2-y1)*acos((x2-x1)/d);
    end
    
    % trigonometric gaussian formula on the arc
    tw1=trigauss(n+2,-omega1,omega1);
    % equal radii: symmetric bubble
    if (r1==r2)
        tw2=tw1;
    else
        tw2=trigauss(n+2,-omega2,omega2);
    end
    
    % algebraic gaussian formula on [-1,1]
    ab=r_jacobi(ceil((n+1)/2),0,0);
    xw=gauss(ceil((n+1)/2),ab);
    
    % creating the grid
    [t1,theta1]=meshgrid(xw(:,1),tw1((1:ceil((n+2)/2)),1));
    [w11,w12]=meshgrid(xw(:,2),tw1((1:ceil((n+2)/2)),2));
    [t2,theta2]=meshgrid(xw(:,1),tw2((1:ceil((n+2)/2)),1));
    [w21,w22]=meshgrid(xw(:,2),tw2((1:ceil((n+2)/2)),2));
    
    % nodal cartesian coordinates and weights
    c1=cos(theta1(:)); s1=sin(theta1(:));
    c2=cos(theta2(:)); s2=sin(theta2(:));
    xyw1(:,1)=r1*c1*cos(phi1)-r1*t1(:).*s1*sin(phi1)+x1;
    xyw1(:,2)=r1*c1*sin(phi1)+r1*t1(:).*s1*cos(phi1)+y1;
    xyw1(:,3)=r1^2*s1.^2.*w11(:).*w12(:);
    xyw2(:,1)=r2*c2*cos(phi2)-r2*t2(:).*s2*sin(phi2)+x2;
    xyw2(:,2)=r2*c2*sin(phi2)+r2*t2(:).*s2*cos(phi2)+y2;
    xyw2(:,3)=r2^2*s2.^2.*w21(:).*w22(:);
    xyw=[xyw1; xyw2];
    
end



