function V=vandermonde_koornwinder_classII(deg,x,y,alpha)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Vandermonde matrix of degree "deg" of the Koornwinder Class II orthogonal
% polynomials for the unit disk.
%
% The evaluation points (x,y) must belong to the unit-disk B((0,0),1).
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% deg  : degree of the Vandermonde matrix;
% x,y  : column vectors of points (x(k),y(k)) of the unit-disk defining the
%      Vandermonde matrix (Koornwinder class II).
%--------------------------------------------------------------------------
% Outputs:
%--------------------------------------------------------------------------
% V : Vandermonde matrix (w.r.t. Koornwinder class II basis).
%--------------------------------------------------------------------------
% Routines
%--------------------------------------------------------------------------
% 1. vandermonde_jacobi (attached below)
%--------------------------------------------------------------------------
% Reference.
%--------------------------------------------------------------------------
% [1] 
% T.H. Koornwinder
% Two variable analogues of the classical orthogonal polynomials
% R. Askey (Ed.), Theory and Application of Special Functions, Academic 
% Press (1975), pp. 435-495.
%
% [2]
% Lidia Fernández, Teresa E. Pérez, Miguel A. Piñar
% On Koornwinder classical orthogonal polynomials in two variables
% Journal of Computational and Applied Mathematics
% Volume 236, Issue 15, September 2012, Pages 3817-3826
%
% [3]
% L.Bos, A. Sommariva, M. Vianello,
% Least-squares polynomial approximation on weakly admissible meshes
% Journal of Computational and Applied Mathematics
% Volume 235, Issue 3, 1 December 2010, Pages 660-668
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%% Copyright (C) 2009- 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
%%
%% Author:  
%%          Alvise Sommariva <alvise@math.unipd.it>
%%          Marco Vianello   <marcov@math.unipd.it>
%%
%% Date    : January 3, 2026
%% Modified: January 3, 2026
%--------------------------------------------------------------------------

if (size(x,1) == 1)
    x=x'; y=y';
end

% We refer to definition of Koornwinder polynomials given in [2], section
% "5.2. Class II: orthogonal polynomials on the unit disk".
%
% The Jacobi-polynomials use the standardized condition 
% P_n^{(alpha,beta)}(1)=nchoosek(n+alpha,alpha).
%
% In the attached subroutine "vandermonde_jacobi" they correspond to those
% computed setting "basis_type='trad'".

t1=( (1-x.^2) == 0 ); % useful for avoid singularity at "x=-1", "x=1".
t=( ((1-x.^2)+t1).^(-1/2) ).*y;
Vdx=vandermonde_jacobi(deg,t,alpha,alpha,'trad'); 

V=[];
degs=[]; % If degs(m,:)=(i,j) then the m-th column corresponds to 
         % "P_(i,j) (x,y)".
for n=0:deg
    for k=0:n 
        a=alpha+k+1/2;
        Vsx=vandermonde_jacobi(n-k,x,a,a,'trad');
        V(:,end+1)=Vsx(:,end).*( (1-x.^2).^(k/2) ).*Vdx(:,k+1);
        degsL=[n-k k];
        degs=[degs; degsL];
    end
end









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



function V=vandermonde_jacobi(N,t,a,b,basis_type)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Vandermonde matrix of Jacobi polynomial with exponents "a", "b", i.e.
% w.r.t. the function "w(t)=(1-t)^a * (1+t)^b" in "(-1,1)" at the points
% "t" (column vector of points). The degree of the polynomial is "N".
%
% The variable "basis_type" is a string and can be
%
% * 'monic': Jacobi monic polynomials (default);
% * 'orthn': Jacobi orthonormal polynomials;
% * 'trad' : traditional polynomials (see https://dlmf.nist.gov/18.2#iii,
%            section "Standardizations");
% * 'gegcl': classical Gegenbauer (ultraspherical) polynomials (see
%            Abramowitz-Stegun); for representing the polynomials  w.r.t.
%            the weight function (1-x^2)^(lambda-1/2), set the input
%            variable "a" equal to "lambda" and "b=[]"; it is normalized so
%            that C^lambda_n(1)=nchoosek(n+2*lambda-1,n);
% * 'ch1cl': classical Chebyshev polynomials of the first kind, that is
%            T_n(x)=cos(n*acos(x)) if "x" belongs to "[-1,1]". If this
%            option is choosen then set "a=-1/2" and "b=-1/2";
% * 'ch2cl': classical Chebyshev polynomials of the second kind, that is
%            U_n(x)=sin((n+1)*acos(x))/sin(acos(x)) if "x" belongs to 
%            "[-1,1]". If this option is choosen then set "a=1/2" and 
%             "b=1/2".
%--------------------------------------------------------------------------
% REFERENCE
%--------------------------------------------------------------------------
% [1] "Orthogonal Polynomials, Computation and Approximation"
%      Walter Gautschi
% [2] "Digital Library of Mathematical Functions",
%      https://dlmf.nist.gov/18.3#T1
% [3] "Handbook of Mathematical Functions with Formulas, Graphs, and
%      Mathematical Tables", Tenth Printing, December 1972.
%      Milton Abramowitz e Irene Stegun
% 
% This routine has been introduced as part of
%
% [a] "ORTHOCUB: integral and  differential cubature rules by orthogonal
%      moments", 2025.
%      by L. Rinaldi, A. Sommariva, M. Vianello.
%--------------------------------------------------------------------------
% License:
%--------------------------------------------------------------------------
% Copyright (C) 2025 Laura Rinaldi, 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 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:
%
% Laura Rinaldi    <laura.rinaldi@unipd.it>
% Alvise Sommariva <alvise@math.unipd.it>
% Marco Vianello   <marcov@math.unipd.it>
%
% Date: November 21, 2025
%--------------------------------------------------------------------------

if nargin < 2
    error('vandermonde_jacobi: too few variables in input.');
end
if nargin < 3, a=0; end
if nargin < 4, b=0; end
if nargin < 5, basis_type='monic'; end




% ......................... Monic polynomials  ............................

if strcmp(basis_type,'monic')
    % See [1] at page 9.
    V=[zeros(length(t),1) ones(length(t),1)];

    % Computing Vandermonde matrix at degree "n".
    if N >= 1
        ab=r_jacobi(N,a,b);
        alpha=ab(:,1); beta=ab(:,2);
        % p_{k+1}(t)=(t-alpha_k) p_k(t) - beta_k p_{k-1}(t)
        for k=1:N % The variable "k" is the degree of the polynomial.
            V(:,k+2)=(t-alpha(k)).*V(:,k+1)-beta(k)*V(:,k);
        end
    end

    % The first column is regarding degree "-1" and must not be considered.
    V=V(:,2:end);
end





% ........................ Orthonormal polynomials  .......................

%--------------------------------------------------------------------------
% Referring to
% See [1], i.e. Gautschi monography, at page 12.
% 
% Orthonormal Vandermonde Jacobi polynomials are computed. Notice that
% increasing the degree the orthonormality tends to deteriorate.
%--------------------------------------------------------------------------
% Example:
%--------------------------------------------------------------------------
%
% deg: 10  a: -4.00000e-01 b: -4.00000e-01 inf.err.: 2.618e-14
% deg: 30  a: -4.00000e-01 b: -4.00000e-01 inf.err.: 3.906e-14
% deg: 50  a: -4.00000e-01 b: -4.00000e-01 inf.err.: 1.577e-13
% deg: 100 a: -4.00000e-01 b: -4.00000e-01 inf.err.: 1.046e-12
%
% Here "inf.err." is the error between the numerical evaluation of the 
% matrix ((phi_i,phi_j)_2,w)_{i,j} and the identity matrix.
%--------------------------------------------------------------------------

if strcmp(basis_type,'orthn')

    ab=r_jacobi(N+1,a,b);
    alpha=ab(:,1); beta=ab(:,2);
    V=[zeros(length(t),1) ones(length(t),1)/sqrt(beta(1))];

    % Computing orthonormal Vandermonde matrix at degree "n".
    for k=1:N % The variable "k" is the degree of the polynomial.
        V(:,k+2)=((t-alpha(k)).*V(:,k+1)-sqrt(beta(k))*V(:,k))/...
            sqrt(beta(k+1));
    end

    % The first column is regarding degree "-1" and must not be considered.
    V=V(:,2:end);

end




% ................. Matlab and Mathematica standard polynomials  ..........

%--------------------------------------------------------------------------
% Referring to
% "https://www.mathworks.com/help/symbolic/sym.jacobip.html".
%
% The formula below gives the family implemented in Matlab by "jacobiP"
% and by Mathematica routine "jacobiP".
% See also "https://mathworld.wolfram.com/JacobiPolynomial.html".
%
% They verify the standardization proposed by Abramowicz-Stegun [3] at
% page 774. The polynomial of degree "n" of this family, say "P(n,a,b)" is
% such that "P(n,a,b)(1)=nchoosek(n+a,n)". We have verified this property.
%--------------------------------------------------------------------------
% Example:
%--------------------------------------------------------------------------
% >> vandermonde_jacobi(3,1,0.3,0.4,'trad')
% ans =
%  1.000000000000000 1.300000000000000 1.495000000000001 1.644500000000001
% >> jacobiP(3,0.3,0.4,1)
% ans =
%  1.644500000000000
% >> % nchoosek via gamma functions (we use that "n!=gamma(n+1)").
% >> gamma((3+0.3)+1)/( gamma(3+1)*gamma(0.3+1) )
% ans =
%  1.644500000000000
%--------------------------------------------------------------------------


if strcmp(basis_type,'trad')
    if N == 0
        V=ones(length(t),1);
    else
        V=[ones(length(t),1) zeros(length(t),N)];
        V(:,2)=(a-b)/2+(1+(a+b)/2)*t; % degree 1
        nV=0:2*N;
        cV=nV+a+b; % the component "c_n" is "cV(n+1)"
        for n=2:N
            term1=2*n*cV(n+1).*cV(2*n-1);
            term2=cV(2*n)*(cV(2*n-1).*cV(2*n+1)*t+a^2-b^2);
            term3=-2*(n-1+a)*(n-1+b).*cV(2*n+1);
            V(:,n+1)=(term2./term1).*V(:,n)+(term3./term1).*V(:,n-1);
        end
    end

end




% ................. Matlab and Mathematica standard Gegenbauer ............

%--------------------------------------------------------------------------
% See "https://www.mathworks.com/help/symbolic/sym.gegenbauerc.html".
% This routine is an open source version of Matlab/Mathematica named
% "gegenbauerC".
% As normalization, it is used that the polynomial C^a_n orthogonal with
% respect to the weight function (1-x^2)^(a-1/2) in (-1,1), a > -1/2 is
% such that C^a_n(1)=nchoosek(n+2a-1,n) if a is not 0, otherwise
% C^a_0(1) = 1 and C^0_n(1)=2/n.
% See Abramovicz-Stegun p.774. This property has been verified.
%
% Note: the parameter "a" correspond to the Gegenbauer parameter "lambda"
%       and not to the Jacobi exponent.
%--------------------------------------------------------------------------
% Example:
%--------------------------------------------------------------------------
% >> V=vandermonde_jacobi(5,1,1/2,1/2,'ch2cl')
% V =
%      1     2     3     4     5     6
% >> chebyshevU(5,1)
% ans =
%      6
% >> V=vandermonde_jacobi(3,0.3,1/2,1/2,'ch2cl')
% >> V =
% 1.000000000000000 0.600000000000000 -0.640000000000000 -0.984000000000000
% >> sin((3+1)*acos(x))/sin(acos(x))
% ans =
% -0.984000000000000
%--------------------------------------------------------------------------

if strcmp(basis_type,'gegcl')
    if N == 0
        V=ones(length(t),1);
    else
        V=[ones(length(t),1) zeros(length(t),N)];
        V(:,2)=2*a*t; % degree 1
        for n=2:N
            term1=2*t*(n+a-1)/n;
            term2=-(n+2*a-2)/n;
            V(:,n+1)=term1*V(:,n)+term2*V(:,n-1);
        end
    end
end





% ................. Matlab and Mathematica standard chebyshevT ............

%--------------------------------------------------------------------------
% See "https://it.mathworks.com/help/symbolic/sym.chebyshevu.html".
% This routine is an open source version of Matlab/Mathematica named
% "gegenbauerC".
% They are the Chebyshev polynomials that in the interval "[-1,1]" are of
% the form "T_n(x)=cos(n*acos(x)".
%
% As normalization, it is used that T_n(1)=1, see Abramovicz-Stegun p.774.
%--------------------------------------------------------------------------
% Example:
%--------------------------------------------------------------------------
% >> V=vandermonde_jacobi(3,1,-1/2,-1/2,'ch1cl')
% V =
%      1     1     1     1
% >> V=vandermonde_jacobi(3,0.5,-1/2,-1/2,'ch1cl')
% V =
% 1.000000000000000 0.500000000000000 -0.500000000000000 -1.000000000000000
% >> n=0:3; cos(n*acos(0.5))
% ans =
% 1.000000000000000 0.500000000000000 -0.500000000000000 -1.000000000000000
%--------------------------------------------------------------------------

if strcmp(basis_type,'ch1cl')
    if N == 0
        V=ones(length(t),1);
    else
        V=[ones(length(t),1) zeros(length(t),N)];
        V(:,2)=t; % degree 1
        for n=2:N
            term1=2*t;
            term2=-1;
            V(:,n+1)=term1*V(:,n)+term2*V(:,n-1);
        end
    end
end




% ................. Matlab and Mathematica standard chebyshevU ............

%--------------------------------------------------------------------------
% See "https://it.mathworks.com/help/symbolic/sym.chebyshevu.html".
% This routine is an open source version of Matlab/Mathematica named
% "chebyshevU".
% They are the Chebyshev polynomials that in the interval "[-1,1]" are of
% the form "U_n(x)=sin((n+1)*acos(x))/sin(acos(x))".
%
% As normalization, it is used that _n(1)=n+1, see Abramovicz-Stegun
% p.774. We have verified this property.
%--------------------------------------------------------------------------
% Example:
%--------------------------------------------------------------------------
% >> V=vandermonde_jacobi(3,1,1/2,1/2,'ch2cl')
% V =
%      1     2     3     4
% >> chebyshevU(3,1)
% ans =
%      4
% >> V=vandermonde_jacobi(3,0.2,1/2,1/2,'ch2cl')
% V =
% 1.000000000000000 0.400000000000000 -0.840000000000000 -0.736000000000000
% >> chebyshevU(3,0.2)
% ans =
% -0.736000000000000
% >> x=0.2; n=3; sin((n+1)*acos(x))/sin(acos(x))
% ans =
% -0.736000000000000
%--------------------------------------------------------------------------

if strcmp(basis_type,'ch2cl')

    if N == 0
        V=ones(length(t),1);
    else
        V=[ones(length(t),1) zeros(length(t),N)];
        V(:,2)=2*t; % degree 1
        for n=2:N
            term1=2*t;
            term2=-1;
            V(:,n+1)=term1*V(:,n)+term2*V(:,n-1);
        end
    end

end





% 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.


function ab=r_jacobi(N,a,b)
if nargin<2, a=0; end;  if nargin<3, b=a; end
if((N<=0)||(a<=-1)||(b<=-1))
    error('parameter(s) out of range')
end
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;
A=[nu (b^2-a^2)*ones(1,N)./(nab.*(nab+2))];
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));
ab=[A' [mu; B1; B']];
