
function rule=cub_interval_clenshawcurtis(ade,card_choice)

%--------------------------------------------------------------------------
% Object
%--------------------------------------------------------------------------
% Clenshaw-Curtis quadrature rule on (-1,1) with ADE equal to "ade".
% The rule is w.r.t. Legendre weight function.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% ade: algebraic degree of precision of the rule, i.e. number of points *1.
%      Tested with n <= 14 million nodes.
% card_choice: the Clenshaw-Curtis formula varies ade depending on its 
%       cardinality "n" ("ade" is always odd). Fixing "ade" we choose as
%       cardinality "n" following "card_choice"
%          1: n=ade+1, 
%          2: n=2*ceil((ade-1)/2)+1 (minimal).
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% rule: (n+1) x 2 matrix of nodes and weights.
%--------------------------------------------------------------------------
% Reference:
%--------------------------------------------------------------------------
% "Fast Construction of Fejer and Clenshaw-Curtis rules for general weight
% functions".
% A. Sommariva
%--------------------------------------------------------------------------
% Note:
%--------------------------------------------------------------------------
% If "n" is the cardinality of the rule, then ADE is 2*floor((n-1)/2)+1.
% Converse, a rule with fixed ADE must have n=ceil((ade-1)/2)+1 nodes.
% This has been checked numerically.
%
% The user may decide to use an interpolatory rule with ADE equal to "ade"
% and cardinality "n=ade+1" (see first lines of the actual code below).
%--------------------------------------------------------------------------
% License:
%--------------------------------------------------------------------------
% Copyright (C) 2025 Alvise Sommariva.
%
% 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 3 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:
%
% Alvise Sommariva <alvise@math.unipd.it>
%
% Date: December 21, 2025
%--------------------------------------------------------------------------

if nargin < 2
    card_choice=1; % 1: n=ade+1, 2: n=2*ceil((ade-1)/2)+1 (minimal).
end

if card_choice
    n=ade+1;
else
    n=2*ceil((ade-1)/2)+1; % example: ade=12 -> card=13, ade=13 -> card=13.
end

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];

rule=[x w];





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

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 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;

