function xyw = gqsymmlens(n,omega,r)
% by Gaspare Da Fies and Marco Vianello, University of Padova
% 8 Nov 2011

% computes the nodes and weights of a product gaussian formula
% on a circular lens symmetric w.r.t. the y-axis, namely   
% the intersection of the two disks with radius r centered 
% at (-r*cos(omega),0) and (r*cos(omega),0), respectively 

% 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

% input: 
% n: algebraic degree of exactness
% omega: half-length of the angular interval, 0<omega<=pi/2 
% r: radius of the disk 

% output:
% xyw: (ceil((n+1)/2) x (n+3)) x 3 array of (xnodes,ynodes,weights) 


% trigonometric gaussian formula on the arc
tw=trigauss(n+2,-omega,omega);

% 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));
[w1,w2]=meshgrid(xw(:,2),tw(:,2));

% nodal cartesian coordinates and weights  
diffc=cos(theta(:))-cos(omega);
xyw(:,1)=r*t(:).*diffc;
xyw(:,2)=r*sin(theta(:));
xyw(:,3)=r^2*cos(theta(:)).*diffc.*w1(:).*w2(:);

end


