
function tw=trigauss_gausslegendre_2017(n,alpha,beta)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Computation of trigonometric gaussian rules on the unit arc
% "[alpha,beta]" with trigonometric degree of precision "n".
%--------------------------------------------------------------------------
% Inputs:
%--------------------------------------------------------------------------
% n    : the rule computes computes the n+1 angles and weights of a
%        trigonometric gaussian quadrature formula on [alpha,beta], with
%                          0 < beta-alpha <= 2*pi;
% alpha, beta: arc angles for trigonometric gaussian rules on the unit arc
%       from "alpha" to "beta".
%--------------------------------------------------------------------------
% Outputs:
%--------------------------------------------------------------------------
% tw   : (n+1) x 2 matrix, where
%        a) "tw(:,1)" determines the nodes,
%        b) "tw(:,2)" determines the weights.
%--------------------------------------------------------------------------
% Reference:
%--------------------------------------------------------------------------
% Subperiodic trigonometric subsampling: a numerical approach
% A. Sommariva and Marco Vianello
% Appl. Anal. Discrete Math. 11 (2017), 470--483
%--------------------------------------------------------------------------
% 1. First version: 2017
%    By A. Sommariva, M. Vianello.
%
% 2. Modified on: December 30, 25.
%--------------------------------------------------------------------------
%% Copyright (C) 2021-
%% 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: MAY 18, 2013
%% Date: DECEMBER 20, 2025
%--------------------------------------------------------------------------

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

if nargin < 1, n=10; end
if nargin < 2, beta=pi; alpha=-beta; end
if nargin < 3, beta=alpha; alpha=-beta; end

if beta-alpha==2*pi
    tw=circle_trigquad(n,alpha,beta);
else
    NN=degree_trigGL(n,alpha,beta);

    ab=r_jacobi(NN,0,0);
    xw=gauss(NN,ab); xGL=xw(:,1); wGL=xw(:,2);
    t=(beta+alpha)/2+(beta-alpha)*xGL/2;
    w=(beta-alpha)*wGL/2;
    tw=[t w];

end









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

function NN=degree_trigGL(n,alpha,beta)

% this routine computes the number of nodes that Gauss-Legendre needs so to
% match the trigonometric moments up to 10^(-14).

theta=(beta-alpha)/2;

if n < 500

    u=[1 (25:25:500)];

    L=[8 29    45    60    75    89   103   119   132   146   158   173  ...
        185   199 212   225   240   253   264   279   291];

    s=spline(u,L);

    NN=ceil(ppval(s,n*theta));

else

    u=n*theta;
    s=0.54*u+21;
    NN=ceil(s);

end









function xw=gauss(N,ab)

% GAUSS Gauss quadrature rule.
%
%    Given a weight function w encoded by the nx2 array ab of the
%    first n recurrence coefficients for the associated orthogonal
%    polynomials, the first column of ab containing the n alpha-
%    coefficients and the second column the n beta-coefficients,
%    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
%    the n-point Gauss quadrature rule for the weight function w.
%    The nodes, in increasing order, are stored in the first
%    column, the n corresponding weights in the second column, of
%    the nx2 array xw.
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.

N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
    J(n,n-1)=sqrt(ab(n,2));
    J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];









function ab=r_jacobi(N,a,b)

% R_JACOBI Recurrence coefficients for monic Jacobi polynomials.
%
%    ab=R_JACOBI(n,a,b) generates the first n recurrence
%    coefficients for monic Jacobi polynomials with parameters
%    a and b. These are orthogonal on [-1,1] relative to the
%    weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
%    are stored in the first column, the n beta-coefficients in
%    the second column, of the nx2 array ab. The call ab=
%    R_JACOBI(n,a) is the same as ab=R_JACOBI(n,a,a) and
%    ab=R_JACOBI(n) the same as ab=R_JACOBI(n,0,0).
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.

nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1
    ab=[nu mu]; return
end

N=N-1;
n=1:N;
nab=2*n+a+b;
nuadd=(b^2-a^2)*ones(1,N)./(nab.*(nab+2));
A=[nu nuadd];
n=2:N;
nab=nab(n);
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
abadd=[mu; B1; B'];
ab=[A' abadd];