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

% computes the n+1 angles and weights of a trigonometric gaussian
% quadrature formula on [alpha,beta], 0<beta-alpha<=pi with respect to
% the weight function w(theta)=abs(sin(theta))

% uses the routines chebyshev.m, gauss.m from
% www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html 
% we suggest to put the following statements
% ab = zeros(N,2); sig = zeros(N+1,2*N);
% at the beginning of the body of chebyshev.m to speed-up execution   

% input:
% n: trigonometric degree of exactness
% [alpha,beta]: angular interval, 0<beta-alpha<=pi

% output:
% tw: (n+1) x 2 array of (angles,weights)

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



% modified Chebyshev moments
mom=(quadvgk(@(t)chebpolys_mod(t,2*n+2,omega),[-1;1],2*n+2))';

% normalization of the moments (monic polynomials)
k=(3:length(mom));
mom(3:end)=exp((2-k)*log(2)).*mom(3:end);

% recurrence coeffs of the monic Chebyshev polynomials
abm(:,1)=zeros(2*n+1,1);
abm(:,2)=0.25*ones(2*n+1,1); abm(1,2)=pi; abm(2,2)=0.5; 


% recurrence coeffs for the monic OPS w.r.t. the weight function
% w(x)=4*(sin(omega/2))^2*abs(x) (it will
% by the modified Chebyshev algorithm 
[ab,normsq]=chebyshev(n+1,mom,abm);

% Gaussian formula for the weight function above
xw=gauss(n+1,ab);

% angles and weights for the trigonometric gaussian formula w.r.t. the
% weight function w(theta)=abs(sin(theta))
tw(:,1)=2*asin(sin(omega/2)*xw(:,1));
tw(:,2)=xw(:,2);

end

function Y = chebpolys_mod(x,m,omega)
k=(1:m);
Y=cos((k'-1)*acos(x)).*(ones(m,1)*abs(x).*4.*(sin(omega/2))^2);
end