function xyw = gqellblend(n,A,B,C,alpha,beta)
% by Gaspare Da Fies, Alvise Sommariva and Marco Vianello
% University of Padova
% 2 June 2011

% computes the nodes and weights of a product gaussian formula
% exact on total-degree bivariate polynomials of degree <=n 
% on the planar region R obtained by linear blending (convex combination) 
% of two trigonometric arcs with parametric equations 
% P(theta)=A1*cos(theta)+B1*sin(theta)+C1
% Q(theta)=A2*cos(theta)+B2*sin(theta)+C2
% namely
% R = {(x,y)=t*P(theta)+(1-t)*Q(theta), t in [0,1], theta in [alpha,beta], 
% 0<beta-alpha<=2*pi} 

% uses the routines:
%
% r_jacobi.m, gauss.m from
% www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
%
% trigauss.m 
% http://www.math.unipd.it/~marcov/mysoft/trigauss.m
% this will be soon substituted by an optimized version 

% input: 
% n: algebraic degree of exactness
% A,B,C: 2x2 matrices of the parametric arc coefficients:
% A1=A(1,:), B1=B(1,:), C1=C(1,:)
% A2=A(2,:), B2=B(2,:), C2=C(2,:)
% [alpha,beta]: angular interval, 0<beta-alpha<=2*pi

% output:
% xyw: 3 columns array of (xnodes,ynodes,weights) 


% computing the algebraic and trigonometric degree increment 
S1=abs((A(1,1)-A(2,1))*(B(1,2)-B(2,2))+(A(1,2)-A(2,2))*(B(2,1)-B(1,1)))>10*eps;
S2=abs((C(1,1)-C(2,1))*(B(1,2)-B(2,2))+(C(1,2)-C(2,2))*(B(2,1)-B(1,1)))>10*eps;
S3=abs((A(1,1)-A(2,1))*(C(1,2)-C(2,2))+(A(1,2)-A(2,2))*(C(2,1)-C(1,1)))>10*eps;

if (S1 || S2 || S3)
    h=1;
else
    h=0;
end

S4=abs(A(1,2)*A(2,1)-A(1,1)*A(2,2)-B(1,2)*B(2,1)+B(1,1)*B(2,2))>10*eps;
S5=abs(A(1,2)*B(2,1)-A(1,1)*B(2,2)+B(1,2)*A(2,1)-B(1,1)*A(2,2))>10*eps;
S6=abs(B(2,1)*(C(1,2)-C(2,2))-B(2,2)*(C(1,1)-C(2,1)))>10*eps;
S7=abs(A(2,1)*(C(1,2)-C(2,2))-A(2,2)*(C(1,1)-C(2,1)))>10*eps;
S8=abs((C(1,1)-C(2,1))*(B(1,2)-B(2,2))+(C(1,2)-C(2,2))*(B(2,1)-B(1,1)))>10*eps;
S9=abs((A(1,1)-A(2,1))*(C(1,2)-C(2,2))+(A(1,2)-A(2,2))*(C(2,1)-C(1,1)))>10*eps;

if (S4 || S5)
    k=2;
elseif (S6 || S7 || S8 || S9)
    k=1;
else
    k=0;
end
    
% trigonometric gaussian formula on the arc
tw=trigauss(n+k,alpha,beta);
           
% algebraic gaussian formula on [0,1] 
ab=r_jacobi(ceil((n+h+1)/2),0,0);
xw=gauss(ceil((n+h+1)/2),ab);
xw(:,1)=xw(:,1)/2+1/2;
xw(:,2)=xw(:,2)/2;

% creating the grid
[t,theta]=meshgrid(xw(:,1),tw(:,1));
[w1,w2]=meshgrid(xw(:,2),tw(:,2));

% nodal cartesian coordinates and weights  
s=sin(theta(:));
c=cos(theta(:));
p1=A(1,1)*c+B(1,1)*s+C(1,1);
p2=A(1,2)*c+B(1,2)*s+C(1,2);
q1=A(2,1)*c+B(2,1)*s+C(2,1);
q2=A(2,2)*c+B(2,2)*s+C(2,2);
dp1=-A(1,1)*s+B(1,1)*c;% plot(xyw(:,1)+trasl,xyw(:,2),'.','MarkerSize',6);
dp2=-A(1,2)*s+B(1,2)*c;
dq1=-A(2,1)*s+B(2,1)*c;
dq2=-A(2,2)*s+B(2,2)*c;

xyw(:,1)=p1.*t(:)+q1.*(1-t(:));
xyw(:,2)=p2.*t(:)+q2.*(1-t(:));
xyw(:,3)=abs((p1-q1).*(dp2.*t(:)+dq2.*(1-t(:))) - ... 
(p2-q2).*(dp1.*t(:)+dq1.*(1-t(:)))).* w1(:).*w2(:);
         

% this part plots the region R and the cubature nodes        
% plot(xyw(:,1),xyw(:,2),'.','MarkerSize',6);
% axis equal;
% hold on;
% theta=linspace(alpha,beta,500);
% s=sin(theta);
% c=cos(theta);
% p1=A(1,1)*c+B(1,1)*s+C(1,1);
% p2=A(1,2)*c+B(1,2)*s+C(1,2);
% plot(p1,p2);
% q1=A(2,1)*c+B(2,1)*s+C(2,1);
% q2=A(2,2)*c+B(2,2)*s+C(2,2);
% plot(q1,q2);
% t=linspace(0,1,100);
% seg11=(A(1,1)*cos(alpha)+B(1,1)*sin(alpha)+C(1,1))*t + ...
% (A(2,1)*cos(alpha)+B(2,1)*sin(alpha)+C(2,1))*(1-t);
% seg12=(A(1,2)*cos(alpha)+B(1,2)*sin(alpha)+C(1,2))*t + ...
% (A(2,2)*cos(alpha)+B(2,2)*sin(alpha)+C(2,2))*(1-t);
% plot(seg11,seg12);
% seg21=(A(1,1)*cos(beta)+B(1,1)*sin(beta)+C(1,1))*t + ...
% (A(2,1)*cos(beta)+B(2,1)*sin(beta)+C(2,1))*(1-t);
% seg22=(A(1,2)*cos(beta)+B(1,2)*sin(beta)+C(1,2))*t + ...
% (A(2,2)*cos(beta)+B(2,2)*sin(beta)+C(2,2))*(1-t);
% plot(seg21,seg22);

         
end