function tw=trigauss_clenshawcurtis_2013(tde,alpha,beta,rule_type)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Computation of trigonometric Clenshaw-Curtis/Fejer rules on the unit arc
% "[alpha,beta]" with trigonometric degree of exactness "tde".
%
% In this routine this is done by weighted Clenshaw-Curtis/Fejer rule.
%--------------------------------------------------------------------------
% Note:
%--------------------------------------------------------------------------
% This routine requires Signal Processing Toolbox.
%--------------------------------------------------------------------------
% Inputs:
%--------------------------------------------------------------------------
% tde  : the rule computes computes the n+1 angles and weights of a
%       trigonometric Clenshaw-Curtis/Fejer quadrature formula on
%       [alpha,beta], where
%                          0 < beta-alpha <= 2*pi;
% alpha, beta: arc angles for trigonometric Clenshaw-Curtis rules on the
%        unit arc from "alpha" to "beta".
% rule_type: 1: Fejer I, 2: Fejer II, 3: Clenshaw-Curtis.
%--------------------------------------------------------------------------
% Outputs:
%--------------------------------------------------------------------------
% tw   : (tde+1) x 2 matrix, where
%        a) "tw(:,1)" determines the nodes,
%        b) "tw(:,2)" determines the weights.
%--------------------------------------------------------------------------
% Reference:
%--------------------------------------------------------------------------
% Fast Construction of Fejer and Clenshaw-Curtis/Fejer rules for general
% weight functions,
% Computers Mathematics with Applications,
% Volume 65, Issue 4, February 2013, Pages 682-693.(Open Access)
% A. Sommariva
%--------------------------------------------------------------------------
% 1. First version: 2013
%    By A. Sommariva
%
% 2. Modified on: December 30, 25.
%--------------------------------------------------------------------------
% Examples
%--------------------------------------------------------------------------
% >> % FEJER I
% >> tde=5; alpha=pi/4; beta=pi/2; rule_type=1;
% >> tw=trigauss_clenshawcurtis_2013(tde,alpha,beta,rule_type);
% >> tw
% tw =
%     1.5667    0.0141
%     1.5349    0.0483
%     1.4741    0.0729
%     1.3894    0.0950
%     1.2881    0.1065
%     1.1781    0.1120
%     1.0681    0.1065
%     0.9668    0.0950
%     0.8821    0.0729
%     0.8213    0.0483
%     0.7894    0.0141
% >> length(tw)
% ans =
%     11
% >> % FEJER II
% >> tde=5; alpha=pi/4; beta=pi/2; rule_type=2;
% >> tw=trigauss_clenshawcurtis_2013(tde,alpha,beta,rule_type);
% >> tw
% tw =
% 1.5513    0.0376
% 1.4951    0.0734
% 1.4079    0.0996
% 1.2987    0.1170
% 1.1781    0.1224
% 1.0575    0.1170
% 0.9482    0.0996
% 0.8611    0.0734
% 0.8049    0.0376
% 0.7854    0.0040
% >> length(tw)
% ans =
%     11
% >> % CLENSHAW-CURTIS
% >> tde=5; alpha=pi/4; beta=pi/2; rule_type=3;
% >> tw=trigauss_clenshawcurtis_2013(tde,alpha,beta,rule_type);
% >> tw
% tw =
%     1.5708    0.0040
%     1.5513    0.0376
%     1.4951    0.0734
%     1.4079    0.0996
%     1.2987    0.1170
%     1.1781    0.1224
%     1.0575    0.1170
%     0.9482    0.0996
%     0.8611    0.0734
%     0.8049    0.0376
%     0.7854    0.0040
% >> length(tw)
% ans =
%     11
% >> 
%--------------------------------------------------------------------------
%% Copyright (C) 2013-
%% Gaspare Da Fies, Gerard Meurant, Alvise Sommariva, Marco Vianello.
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Authors:
%% Gaspare Da Fies, Gerard Meurant, Alvise Sommariva, Marco Vianello.
%%
%% Date: 2013
%% Date: DECEMBER 20, 2025
%--------------------------------------------------------------------------

% .......................... troubleshooting  .............................

if nargin < 1, tde=10; end
if nargin < 2, beta=pi; alpha=-beta; end
if nargin < 3, beta=alpha; alpha=-beta; end
if nargin < 4, rule_type=3; end

n=2*tde+1;

omega=(beta-alpha)/2;

switch rule_type
    case 1
        moms=moms_cheb1type_trig(n-1,omega);
        [x,w]=fejerI(n,moms);

    case 2
        moms=moms_cheb1type_trig(n-1,omega);
        [x,w]=fejerII(n,moms);

    otherwise
        moms=moms_cheb1type_trig(n-1,omega);
        [x,w]=clenshaw_curtis(n,moms);

end

t=(alpha+beta)/2+2*asin(sin(omega/2)*x);

tw=[t w];







%--------------------------------------------------------------------------
% Subroutines
%--------------------------------------------------------------------------

function [x,w]=fejerI(n,moms)

% OBJECT: FEJER RULE. TYPE I.

% INPUTS:
% n: NUMBER OF POINTS. TESTED WITH n <= 14 MILLION NODES.
% moms: WEIGHT FUNCTION MODIFIED MOMENTS (w.r.t. CHEB. POLY. TYPE I).

% OUTPUTS:
% x: QUADRATURE NODES IN THE INTERVAL (-1,1). COLUMN VECTOR.
% w: QUADRATURE WEIGHTS. COLUMN VECTOR.

% REFERENCE:
% ALVISE SOMMARIVA
% "Fast Construction of Fejer and Clenshaw-Curtis rules for general weight
% functions".

% AUTHOR:
% ALVISE SOMMARIVA, JULY 23, 2012.

if nargin < 2
    moms=moms_cheb1type_legendre(n-1);
end

moms(1)=sqrt(1/2)*moms(1);
x=cos((2*(1:n)'-1)*pi/(2*n));
w=sqrt(2/n)*idct(moms);









function [x,w]=fejerII(n,moms)

% OBJECT: FEJER RULE. TYPE II.

% INPUTS:
% n: NUMBER OF POINTS. TESTED WITH n <= 14 MILLION NODES.
% moms: WEIGHT FUNCTION MODIFIED MOMENTS (w.r.t. CHEB. POLY. TYPE I).

% OUTPUTS:
% x: QUADRATURE NODES IN THE INTERVAL (-1,1). COLUMN VECTOR.
% w: QUADRATURE WEIGHTS. COLUMN VECTOR.

% REFERENCE:
% ALVISE SOMMARIVA
% "Fast Construction of Fej?er and Clenshaw-Curtis rules for general weight
% functions".

% AUTHOR:
% ALVISE SOMMARIVA, JULY 23, 2012.

if nargin < 2
    moms=moms_cheb1type_legendre(n-1);
end

momsII=compute_moments_IIw(n-1,moms);
theta=(1:n)*pi/(n+1); theta=theta'; x=cos(theta);
w=(2*sin(theta)/(n+1)).*dst(momsII);









function [x,w]=clenshaw_curtis(n,moms)

% OBJECT: CLENSHAW-CURTIS RULE.

% INPUTS:
% n: NUMBER OF POINTS. TESTED WITH n <= 14 MILLION NODES.
% moms: WEIGHT FUNCTION MODIFIED MOMENTS (w.r.t. CHEB. POLY. TYPE I).

% OUTPUTS:
% x: QUADRATURE NODES IN THE INTERVAL [-1,1]. COLUMN VECTOR.
% w: QUADRATURE WEIGHTS. COLUMN VECTOR.

% REFERENCE:
% ALVISE SOMMARIVA
% "Fast Construction of Fej?er and Clenshaw-Curtis rules for general weight
% functions".

% AUTHOR:
% ALVISE SOMMARIVA, JULY 23, 2012.

if nargin < 2
    moms=moms_cheb1type_legendre(n-1);
end

momsIIcc=compute_moments_IIwcc(n-3,moms);
theta=(1:n-2)'*pi/(n-1); xx=cos(theta);
w=((2*sin(theta)/(n-1)).*dst(momsIIcc))./(1-xx.^2);
w1=(2*sum(moms)-moms(1)-moms(end))/(2*(n-1));
wn=moms(1)-w1-sum(w);
x=[1;xx;-1]; w=[w1;w;wn];









function [momsII,momsI]=compute_moments_IIwcc(n,momsI)

% COMPUTATION OF MOMENTS W.R.T. CHEBYSHEV SECOND TYPE FROM CHEBYSHEV FIRST
% TYPE.

momsI=momsI(1:n+3);
momsII=0.5*(momsI(1:end-2)-momsI(3:end));









function momsII=compute_moments_IIw(n,momsI)

% COMPUTATION OF MOMENTS W.R.T. CHEBYSHEV SECOND TYPE FROM CHEBYSHEV FIRST
% TYPE.

momsI=momsI(1:n+1);

momsIe=momsI(1:2:end); momsIo=momsI(2:2:end);
momsII=zeros(size(momsI));
momsII(1:2:end)=2*cumsum(momsIe)-momsI(1);
momsII(2:2:end)=2*cumsum(momsIo);









function moms=moms_cheb1type_legendre(n)

moms=zeros(n+1,0);
m=0:2:n; m=m';
moms_even_degree=2./(1-m.^2);
moms(m+1,1)=moms_even_degree;









function moms=moms_cheb1type_trig(n,omega)

% fprintf('\n \t SUBRANGE');

moms=zeros(1,n+1);
moms(1)=2*omega; % FIRST MOMENT.

if(n>=2)

    if(omega<=1/4*pi)
        l=10;
    elseif(omega<=1/2*pi)
        l=20;
    elseif(omega<=3/4*pi)
        l=40;
    else
        if omega == pi
            l=2*ceil(10*pi);
        else
            l=2*ceil(10*pi/(pi-omega));
        end
    end


    temp=(2:2:n+2*l-2); % AUXILIAR VECTORS.
    temp2=temp.^2-1;

    dl=1/4 -1./(4*(temp-1)); % DIAGONALS.
    dc=1/2 -1/sin(omega/2)^2 -1./(2*temp2);
    du=1/4 +1./(4*(temp+1));

    d=4*cos(omega/2)/sin(omega/2)./temp2'; % COMPUTING KNOWN TERM.
    d(end)=d(end);                         % PUT LAST MOMENT NULL.

    z=tridisolve(dl(2:end),dc,du(1:end-1),d); % SOLVE SYSTEM.
    moms(3:2:n+1)=z(1:floor(n/2)); % SET ODD MOMENTS.

end

moms=moms';

normalized = 1;

if normalized == 0
    M=length(moms);
    kk=2.^(-((1:2:M)-2))'; kk(1)=1;
    v=ones(M,1);
    v(1:2:M)=kk;
    moms=v.*moms;
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)

% optimized version
x = d;
n = length(x);
bi = zeros(n,1);

for j = 1:n-1
    bi(j) = 1 / b(j);
    mu = a(j) * bi(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)) * bi(j);
end