
function rule = cub_interval_fejer(n,ftype)

%--------------------------------------------------------------------------
% Object
%--------------------------------------------------------------------------
% Clenshaw-Curtis quadrature rule on (-1,1) with ADE equal to "ade".
% The rule is w.r.t. Legendre weight function.
%
% Here "n" is the cardinality of the interpolatory rule described by
% the variable "ftype".
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% n: cardinality of the cubature rule
% ftype: 1: Fejer I, 2: Fejer II, 3: Clenshaw-Curtis.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% rule: n x 2 matrix of nodes "rule(:,1)" and weights "rule(:,2)".
%--------------------------------------------------------------------------
% Reference:
%--------------------------------------------------------------------------
% J. Waldvogel
% "Fast Construction of Fejer and Clenshaw-Curtis Quadrature rules".
% BIT, 2003, Vol. 43, No. 1, pp. 001–018
%--------------------------------------------------------------------------
% 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 18, 2025
%--------------------------------------------------------------------------
% Note
%--------------------------------------------------------------------------
% Based on a routine by J. Waldvogel.
%--------------------------------------------------------------------------

if nargin < 2, ftype=2; end

% .............................. Nodes/weights  ...........................


switch ftype

    case 1  % Fejer1

        N=(1:2:n-1)'; l=length(N); m=n-l; K=(0:m-1)';
        k=1:n;
        x = cos((2*k-1)'*pi/(2*n));

        v0=[2*exp(1i*pi*K/n)./(1-4*K.^2); zeros(l+1,1)];
        v1=v0(1:end-1)+conj(v0(end:-1:2)); w=ifft(v1);

    case 2 % Fejer2

        n=n+1;
        k=(1:n-1)';
        x=cos(k*pi/n);

        N=(1:2:n-1)'; l=length(N); m=n-l; 
        v0=[2./N./(N-2); 1/N(end); zeros(m,1)];
        v2=-v0(1:end-1)-v0(end:-1:2); wf2=ifft(v2);

        w=wf2(2:end);


    case 3 % Clenshaw-Curtis

        n=n-1;

        N=(1:2:n-1)'; l=length(N); m=n-l; K=(0:m-1)';

        v0=[2./N./(N-2); 1/N(end); zeros(m,1)];
        v2=-v0(1:end-1)-v0(end:-1:2);

        g0=-ones(n,1); g0(1+l)=g0(1+l)+n; g0(1+m)=g0(1+m)+n;
        g=g0/(n^2-1+mod(n,2)); wcc=ifft(v2+g);

        w=[wcc; wcc(1)];

        NN=length(w)-1; theta=(0:NN)*pi;
        x=cos(theta'/NN);

end


rule=[x w];
