function tw = trigquad(n,alpha,beta);

% by Gaspare Da Fies and Marco Vianello, University of Padova
% 8 Nov 2011

% computes the angles theta and weights w of a Fejer-like quadrature 
% formula exact for trigonometric polynomials of degree not greater 
% than n on the angular interval [alpha,beta] with beta-alpha<=2*pi

% input
% n: trigonometric degree of exactness 
% alpha,beta: endpoints of the angular interval of length <= 2*pi  

% output
% tw: 2-column array of the 2n+1 angular nodes and weights

% the formula integrates the canonical trigonometric basis with accuracy 
% from about 10^(-16) (small omega) to about 10^(-14) (omega-->pi) 
% up to n=300 


% omega: half-length of the angular interval, 0<omega<=pi
omega=(beta-alpha)/2;

% quadrature angles
k=[0:1:2*n];
x=cos((pi*k+pi/2)/(2*n+1));
tw(:,1)=2*asin(sin(omega/2)*x)+(alpha+beta)/2;

% even modified Chebyshev moments by recurrence
z(1)=2*omega;
z(n+1)=quadgk(@(t)cos(2*n*acos(sin(t/2)/sin(omega/2))),-omega,omega,'MaxIntervalCount',5000);
temp=(2:2:2*n-1);
dl=1/4 -1./(4*(temp-1));
dc=1/2 -1/sin(omega/2)^2 -1./(2*(temp.^2-1));
du=1/4 +1./(4*(temp+1));
d=4*cos(omega/2)/sin(omega/2)./(temp.^2-1)';
d(n-1)=d(n-1)-du(n-1)*z(n+1);
z(2:n)=tridisolve(dl(2:n-1),dc(1:n-1),du(1:n-2),d(1:n-1));
mom=z(1:n+1);

% quadrature weights
T=cos((2*(1:n))'*acos(x));
tw(:,2)=((mom(1)+2*(mom(2:n+1))*T)/(2*n+1))';

end


function x = tridisolve(a,b,c,d)    
%   TRIDISOLVE  Solve tridiagonal system of equations.
% From Cleve Moler's Matlab suite
% http://www.mathworks.it/moler/ncmfilelist.html 

%     x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
%     b(1)*x(1) + c(1)*x(2) = d(1),
%     a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
%     a(n-1)*x(n-1) + b(n)*x(n) = d(n).
%
%   The algorithm does not use pivoting, so the results might
%   be inaccurate if abs(b) is much smaller than abs(a)+abs(c).
%   More robust, but slower, alternatives with pivoting are:
%     x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%     x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)

x = d;
n = length(x);
for j = 1:n-1
   mu = a(j)/b(j);
   b(j+1) = b(j+1) - mu*c(j);
   x(j+1) = x(j+1) - mu*x(j);
end
x(n) = x(n)/b(n);
for j = n-1:-1:1
   x(j) = (x(j)-c(j)*x(j+1))/b(j);
end
end


