
function xyw = cub_circzone(n,alpha,beta,r)

%--------------------------------------------------------------------------
% Object:
% The routine computes the nodes and weights of a product gaussian formula
% on a circular zone of a disk centered at the origin namely, the
% portion of the disk between the lines x=r*cos(alpha) and x=r*cos(beta).
%
% Note: in the special case alpha=0 the zone is a circular segment.
%--------------------------------------------------------------------------
% Input:
% n: algebraic degree of exactness
% [alpha,beta]: angular interval, 0<=alpha<beta<=pi
% r: radius of the disk
%--------------------------------------------------------------------------
% Output:
% xyw: (ceil((n+1)/2) x (n+3)) 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: November 8, 2011.
% Last update: January 4. 2020.
%--------------------------------------------------------------------------

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

% algebraic gaussian formula on [-1,1]
ab=r_jacobi(ceil((n+1)/2),0,0);
xw=gauss(ceil((n+1)/2),ab);

% creating the polar grid
[t,theta]=meshgrid(xw(:,1),tw(:,1));
[w1,w2]=meshgrid(xw(:,2),tw(:,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(:);




