
function V=vandermonde_subinterval(deg,omega,theta)

%--------------------------------------------------------------------------
% OBJECT.
%--------------------------------------------------------------------------
% Computation of the Vandermonde matrix on the points "theta",of the
% orthonormal basis of the space of trigonometric polynomials of degree
% "deg", relatively to the subinterval [-omega,omega] of [0,2*pi].
%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% deg  : degree of the space of trigonometric polynomials.
% omega: determined the interval [-omega,omega] in which the basis must be
%        evaluated.
% theta: column vector of points in which the orthonormal trig. basis of
%        degree "deg" is evaluated.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% V    : Vandermonde matrix w.r.t. orthonormal trigonometric basis in the 
%        subinterval [-omega,omega] of [0,2*pi]
%--------------------------------------------------------------------------
% REFERENCE PAPER.
%--------------------------------------------------------------------------
% Subperiodic Trigonometric Hyperinterpolation
% ⇥ "Festschrift for the 80th Birthday of Ian Sloan" (invited paper), 
% J. Dick, F.Y. Kuo, H. Wozniakowski Eds., Springer, 2018, pp. 283-304.
% » Joint paper with G. Da Fies and M. Vianello.
%--------------------------------------------------------------------------
% Copyright
%--------------------------------------------------------------------------
%% Copyright (C) 2020-
%% Gaspare Da Fies, 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:
%% Alvise Sommariva, Marco Vianello.
%%
%% Date: 2020
%% Date: DECEMBER 30, 2025
%--------------------------------------------------------------------------

xx=sin(theta/2)/sin(omega/2);

% ....................... Evaluation of Vodd  .............................

%---------------------------------
% CASE 1. beta=+1/2.
%---------------------------------

% Three terms recorsion coefficients
ab_p=r_trigauss_general(2*deg+1,omega,+1/2); % Recurr. coeff.

% Computing normalization terms
ab_p_sqrt2=sqrt(ab_p(:,2));
Ip=zeros(1,2*deg);
for ii=1:2*deg
    Ip(ii)=prod(ab_p_sqrt2(1:ii));
end
norm_parms_p=Ip(2:2:end);

% Evaluation of odd columns
V_p=orthpol_eval(ab_p(1:2*deg-1,:),xx);
VV_p=V_p(:,2:2:end);

Dp=diag(1./norm_parms_p(1:deg));

Vodd=VV_p*Dp;
cs=cos(theta/2);
Dpl=diag(cs);
Vodd=Dpl*Vodd;

% ....................... Evaluation of Veven  ............................

%---------------------------------
% CASE 2. beta=-1/2.
%---------------------------------

% Three terms recorsion coefficients
ab_m=r_trigauss_general(2*deg+1,omega,-1/2);

% Computing normalization terms
ab_m_sqrt2=sqrt(ab_m(:,2));

Im=zeros(1,2*deg+1);
for ii=1:2*deg+1
    Im(ii)=prod(ab_m_sqrt2(1:ii));
end
norm_parms_m=Im(1:2:end);

V_m=orthpol_eval(ab_m(1:2*deg,:),xx);
VV_m=V_m(:,1:2:end);
Dm=diag(1./norm_parms_m(1:deg+1));

Veven=VV_m*Dm;

% ....................... Assembling Vandermonde matrix ...................

V=Veven(:,1);
for ii=1:deg
    V=[V Vodd(:,ii) Veven(:,ii+1)];
end











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

function [ab,mom]=r_trigauss_general(n,omega,beta)


% PIESSENS RECURSION.
a_square=-1/(sin(omega/2))^2;
x_square=0;

[av,bv,cv,dv]=moments_recurrs(n,a_square,x_square,beta);

gamma=sin(omega/2);

switch beta
    case 1/2
        mod_factor=-2*1i*gamma^2; % beta=0.5;
        mom0=2*(omega/2+sin(omega)/2); % beta=1/2
        momfin=quadgk(@(x) cos(2*n*acos(x))*...
            (2*sin(omega/2)).*sqrt(1-(sin(omega/2))^2*x.^2),-1,1,...
            'MaxIntervalCount',50000);

    case -1/2
        mod_factor=2*1i;
        mom0=2*omega;
        momfin=quadgk(@(x) cos(2*n*acos(x))*...
            (2*sin(omega/2))./sqrt(1-(sin(omega/2))^2*x.^2),-1,1,...
            'MaxIntervalCount',50000);
    otherwise
        mod_factor=2*gamma*(-gamma)^(2*beta); % to be tested.
        % mom0e=quadgk(@(x) (2*sin(omega/2)).*...
        %     (1-(sin(omega/2))^2*x.^2).^beta,-1,1,'MaxIntervalCount',50000);
        momfin=quadgk(@(x) cos(2*n*acos(x))*...
            (2*sin(omega/2)).*(1-(sin(omega/2))^2*x.^2).^beta,-1,1,...
            'MaxIntervalCount',50000);
end

dv_mod=mod_factor*dv;
d=dv_mod;
d(1)=d(1)-cv(1)*mom0;
d(end)=d(end)-av(end)*momfin;

z(1,1)=mom0;
z(2:n,1)=tridisolve(cv(2:n-1),bv(1:n-1),av(1:n-2),d(1:n-1));
z(n+1)=momfin;

mom=zeros(2*n+2,1);
mom(1:2:2*n+1)=z(1:n+1);

% 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)=2*sin(omega/2)*(1-sin^2(omega/2)*x^2)^beta
% by the modified Chebyshev algorithm
ab=chebyshev(n+1,mom',abm);









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)

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









function [ab,normsq]=chebyshev(N,mom,abm)

% CHEBYSHEV Modified Chebyshev algorithm.
%
%    Given a weight function w encoded by its first 2n modified
%    moments, stored in the (row) vector mom, relative to monic
%    polynomials defined by the (2n-1)x2 array abm of their
%    recurrence coefficients, [ab,normsq]=CHEBYSHEV(n,mom,abm)
%    generates the array ab of the first n recurrence coefficients
%    of the orthogonal polynomials for the weight function w, and
%    the vector normsq of their squared norms. 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,normsq]=CHEBYSHEV(n,mom) does the same, but using the
%    classical Chebyshev algorithm. If n is larger than the sizes
%    of mom and abm warrant, then n is reduced accordingly.
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.

if N<=0, error('N out of range'), end
if N>size(mom,2)/2, N=size(mom,2)/2; end
if nargin<3, abm=zeros(2*N-1,2); end
if N>(size(abm,1)+1)/2, N=(size(abm,1)+1)/2; end
ab(1,1)=abm(1,1)+mom(2)/mom(1); ab(1,2)=mom(1);
if N==1, normsq(1)=mom(1); return, end
sig(1,1:2*N)=0; sig(2,:)=mom(1:2*N);
for n=3:N+1
    for m=n-1:2*N-n+2
        sig(n,m)=sig(n-1,m+1)-(ab(n-2,1)-abm(m,1))*sig(n-1,m) ...
            -ab(n-2,2)*sig(n-2,m)+abm(m,2)*sig(n-1,m-1);
    end
    ab(n-1,1)=abm(n-1,1)+sig(n,n)/sig(n,n-1)-sig(n-1,n-1)/ ...
        sig(n-1,n-2);
    ab(n-1,2)=sig(n,n-1)/sig(n-1,n-2);
end

normsq=zeros(N,1);

for n=1:N
    normsq(n)=sig(n+1,n);
end
normsq=normsq';









function [av,bv,cv,dv]=moments_recurrs(n,a_square,x_square,beta)

% See Da Fies thesis page 16.
%
% Taken from
% R. Piessens and M. Branders, 
% "Numerical solution of integral equations of mathematical physics using 
% Chebyshev polynomials", 
% Journal of computational physics 21 (1976),178-196.

k=2*(1:n-1)';
av=(1/4)+(beta+1)./(2*(k+1));
bv=(1/2)+a_square+x_square-(beta+1)./(k.^2-1);
cv=(1/4)-(beta+1)./(2*(k-1));
num_dv=2*(1+a_square+x_square).^(beta+1); 
den_dv=k.^2-1; 
dv=-num_dv./den_dv;









function V=orthpol_eval(ab,x)

% x column vector.

N=size(ab,1);

V=[zeros(size(x)) ones(size(x))];

aa=ab(:,1); bb=ab(:,2);
for ii=1:N
    Vn=(x-aa(ii)).*V(:,end)-bb(ii)*V(:,end-1);
    V=[V Vn];
end

V=V(:,2:end);
