function wam = wam_lune(deg,x1,y1,r1,x2,y2,r2)

%--------------------------------------------------------------------------
% Object:
% This routine computes the nodes of a WAM (Weakly Admissible Mesh)
% for total-degree bivariate polynomials of degree "deg" on the circular 
% lune obtained as difference of the disk with center "(x1,y1)" and radius
% "r1", and the disk with center "(x2,y2)" and radius "r2".
%
% The formula is chosen to have the lowest possible cardinality
%--------------------------------------------------------------------------
% Input:
% deg: polynomial degree
% x1,y1: center coordinates of the first disk
% r1: radius of the first disk
% x2,y2: center coordinates of the second disk
% r2: radius of the second disk
%--------------------------------------------------------------------------
% Output:
% wam: 2-column array of the WAM point coordinates
%--------------------------------------------------------------------------
% Dates:
% Written by Alvise Sommariva and Marco Vianello, University of Padova
% October 2014
%
% Modified:
% November 5, 2020, by A. Sommariva 
%--------------------------------------------------------------------------

% uses the function wam_blend in the case of the second disk
% contained into the first


% distance between the centers
d=sqrt((x1-x2)^2+(y1-y2)^2);

% by translation and rotation we can always reduce the lune to the case
% disk(0,0,r1) setminus disk(-d,0,r2)

if d>=(r1+r2)
    
    disp('the disks do not overlap')
    
else
    
    if d<=abs(r1-r2)
        
        % the disks are contained one into the other
        
        if(r1<=r2)
            
            disp('empty difference')
            
        else
            
            % the difference is an asymmetric annulus
            
            A=[r1 0;r2 0]; B=[0 r1;0 r2]; C=[0 0;-d 0];
            wam=wam_blend(deg,A,B,C,-pi,pi);
            
        end
        
    else
        
        % the difference is a proper lune
        
        % angles of the two arcs
        omega2=pi-acos((d^2+r1^2-r2^2)/(2*d*r1));
        omega1=acos((d^2+r2^2-r1^2)/(2*d*r2));
        
        if omega1<=atan(2*(1-cos(omega2))/sin(omega2))
            
            % subperiodic Chebyshev angles
            t1=subpcheb(deg,-omega1,omega1);
            t2=subpcheb(deg,-omega2,omega2);
            t2=[t2(1:deg)];
            
            % creating the grid
            [phi,theta]=meshgrid(t1,t2);
            
            % nodal cartesian coordinates
            c1=cos(phi(:));
            s1=sin(phi(:));
            c2=cos(theta(:));
            s2=sin(theta(:));
            
            wam(:,1)=r1*(c2 + (1-c2)/(1-cos(omega2))*sin(omega2)/sin(omega1).* ...
                (c1-cos(omega1)));
            wam(:,2)=r1/sin(omega1)*s1.*s2;
            wam=[wam;[r1 0]];
            
        elseif (cos(omega2)+cos(omega2-omega1))^2<=4*cos(omega1)
            
            % subperiodic Chebyshev angles
            t1=subpcheb(deg,-omega1,omega1);
            t1=[t1(1:deg)];
            t2=subpcheb(deg,-omega2,omega2);
            
            % creating the grid
            [phi,theta]=meshgrid(t1,t2);
            
            % nodal cartesian coordinates
            c1=cos(phi(:));
            s1=sin(phi(:));
            c2=cos(theta(:));
            s2=sin(theta(:));
            
            wam(:,1)=r1*((c1-cos(omega1))*(cos(omega2)/(1-cos(omega1)) ...
                +sin(omega2)/sin(omega1)) + c2.*(1-c1)/(1-cos(omega1)));
            wam(:,2)=r1/sin(omega1)*s1.*s2;
            wam=[wam;[-d+r2 0]];
            
        else
            
            % subperiodic Chebyshev angles
            t1=subpcheb(deg,-omega1,omega1);
            t2=subpcheb(deg,omega1,omega2);
            
            % creating the grid
            [phi,theta]=meshgrid(t1,t2);
            
            % nodal cartesian coordinates
            c1=cos(phi(:));
            s1=sin(phi(:));
            c2=cos(theta(:));
            s2=sin(theta(:));
            
            wam(:,1)=r1/sin(omega1)*(sin(omega1)*c2-cos(omega1)*s2+s2.*c1);
            wam(:,2)=r1/sin(omega1)*(s1.*s2);
            
        end
        
    end
    
    % rotating and translating the nodes to go back to the original lune
    % disk(x1,y1,r1) setminus disk(x2,y2,r2)
    
    % rotation
    if(d~=0)
        
        R=[(x1-x2)/d (y1-y2)/d; (y1-y2)/d -(x1-x2)/d];
        wam=wam*R';
        
        % translation
        wam(:,1)=x1+wam(:,1);
        wam(:,2)=y1+wam(:,2);
        
    end
    
end

end









function angles=subpcheb(deg,alpha,beta)

% computes Chebyshev-like angles in [alpha,beta]

omega=(beta-alpha)/2;
k=[0:1:2*deg];
x=cos((pi*k+pi/2)/(2*deg+1));
angles=2*asin(sin(omega/2)*x)+(alpha+beta)/2;

end




