
function xyw = cub_circsegm(n,omega,r)

%--------------------------------------------------------------------------
% Object:
% The routine computed the nodes and weights of a product gaussian formula
% on a circular segment of a disk centered at the origin
% with angles in [-omega,omega]
%--------------------------------------------------------------------------
% Input:
% n: algebraic degree of exactness
% omega: half-length of the angular interval, 0<omega<=pi
% r: radius of the disk
%--------------------------------------------------------------------------
% Output:
% xyw: (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.
%--------------------------------------------------------------------------

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

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


