
function [V,dVx,dVy]=proriol_dubiner(pts,deg)

%--------------------------------------------------------------------------
% PRORIOL-DUBINER basis on the unit triangle, i.e. the triangle with
% vertices [0,0], [1,0], [0,1] (as well as its derivatives dVx, dVy).
%
% The matrices V,dVx,dVy have size "#poly.space x #pts" (i.e. Vandermonde
% matrix suitable for cubature).
%
% This version can be speeded up by computing together V, dVx, dVy and not
% separately as in this code. Unfortunately the resulting code would become
% difficult to be understood, but some orthogonal polynomial evaluated less
% times.
%
% The polynomial basis is tested to be orthogonal but not orthonormal. 
%
% See:
% M.G. BLYTH and C.POZRIKIDIS
% "A Lobatto interpolation grid over the triangle".
%--------------------------------------------------------------------------
% Last update: February 20, 2018.
%--------------------------------------------------------------------------

dVx=[];
dVy=[];

V=determine_vand(pts,deg);

if nargout >= 2
    dVx=determine_dVx(pts,deg);
end

if nargout == 3
    dVy=determine_dVy(pts,deg);
end




function V=determine_vand(pts,deg)

% Proriol-Dubiner Vandermonde matrix of degree "deg", at "pts".
% The matrix has size "#poly.space x #pts" (i.e. Vandermonde matrix
% suitable for cubature).

x=pts(:,1); y=pts(:,2);

% Working with row vectors.
if size(x,1) > size(x,2)
    x=x';
end

if size(y,1) > size(y,2)
    y=y';
end

% Following Blyth-Pozrikidis notation, paying attention to singularities.

eps1=((2*x)./(1-y+ (1-abs(sign(1-y))).*eps))-1;
eta1=2*y-1;

term1=OP_eval(eps1,deg,0,0);
term2=MB_eval((1-eta1)/2,deg);
term_left=term1.*term2;

V=[];
for ii=0:deg
    term_rightL=OP_eval(eta1,deg-ii,2*ii+1,0);
    % note: term_left(ii+1,:) is the component at degree "ii".
    term_leftL=repmat(term_left(ii+1,:),size(term_rightL,1),1);
    VL=term_leftL.*term_rightL;
    V=[V; VL];
end








function dVx=determine_dVx(pts,deg)

% Derivative of Proriol-Dubiner Vandermonde basis w.r.t "x" variable.
%
% NB. Proriol-Dubiner Vandermonde matrix of degree "deg", at "pts".
% The matrix has size "#poly.space x #pts" (i.e. Vandermonde matrix
% suitable for cubature).

x=pts(:,1); y=pts(:,2);

% Working with row vectors.
if size(x,1) > size(x,2)
    x=x';
end

if size(y,1) > size(y,2)
    y=y';
end

% Following Blyth-Pozrikidis notation, paying attention to singularities.
% NB: if size(x)= 1 x N then size(eps1)=size(eta1)=size(deps1_dx)= 1 x N.
eps1=((2*x)./(1-y+ (1-abs(sign(1-y))).*eps))-1;
eta1=2*y-1;
deps1_dx=2./(1-y+ (1-abs(sign(1-y))).*eps);

% Here we have to pay attention to partial derivatives.
% About the sizes (important to match up dimensions):
% size(term1a)= deg x N
% size(term1b_pre)= deg x 1
% size(deps1_dx) = 1 x N
% size(term1b) = deg x N
term1a=diffOP_eval(eps1,deg,0,0);
term1b=repmat(deps1_dx,size(term1a,1),1);
term1=term1a.*term1b;
term2=MB_eval((1-eta1)/2,deg);
term_left=term1.*term2;

dVx=[];
for ii=0:deg
    term_rightL=OP_eval(eta1,deg-ii,2*ii+1,0);
    % note: term_left(ii+1,:) is the component at degree "ii".
    term_leftL=repmat(term_left(ii+1,:),size(term_rightL,1),1);
    dVxL=term_leftL.*term_rightL;
    dVx=[dVx; dVxL];
end

term_2_1=OP_eval(eps1,deg,0,0);







function dVy=determine_dVy(pts,deg)

% Derivative of Proriol-Dubiner Vandermonde basis w.r.t "y" variable.
%
% NB. Proriol-Dubiner Vandermonde matrix of degree "deg", at "pts".
% The matrix has size "#poly.space x #pts" (i.e. Vandermonde matrix
% suitable for cubature).
%
% This derivative is the more complicated, by effects of products and chain
% rules.

x=pts(:,1); y=pts(:,2);

% Working with row vectors.
if size(x,1) > size(x,2)
    x=x';
end

if size(y,1) > size(y,2)
    y=y';
end

% Following Blyth-Pozrikidis notation, paying attention to singularities.
eps1=((2*x)./(1-y+ (1-abs(sign(1-y))).*eps))-1;
eta1=2*y-1;

deps1_dy=+2*x./(1-y+ (1-abs(sign(1-y))).*eps).^2;

% Here we have to pay attention to partial derivatives. In particular the
% first one of Legendre polynomials (i.e. degree 0). It is not the
% derivative of a constant but just a constant.
% About the sizes (important to match up dimensions):
% size(term1a)= deg x N
% size(term1b_pre)= deg x 1
% size(deps1_dy) = 1 x N
% size(term1b) = deg x N



% First term of the sums of dVy.
Aterm1a=diffOP_eval(eps1,deg,0,0);
Aterm1b=repmat(deps1_dy,size(Aterm1a,1),1);
Aterm1=Aterm1a.*Aterm1b;
Aterm2=MB_eval((1-eta1)/2,deg);
Aterm_left=Aterm1.*Aterm2;


dVy1=[];
for ii=0:deg
    Aterm_rightL=OP_eval(eta1,deg-ii,2*ii+1,0);
    % note: Aterm_left(ii+1,:) is the component at degree "ii".
    Aterm_leftL=repmat(Aterm_left(ii+1,:),size(Aterm_rightL,1),1);
    dVy1L=Aterm_leftL.*Aterm_rightL;
    dVy1=[dVy1; dVy1L];
end


% Second term of the sums of dVy.
dVy2=[];
Bterm1=OP_eval(eps1,deg,0,0);
Bterm2Ahead=zeros(1,length(x));%Part.der.at degree 0 is equal to constant 1.
Bterm2Atail=MB_eval((1-eta1)/2,deg-1);
Bterm2Btail=repmat((-1)*(1:deg)',1,size(Bterm2Atail,2));
Bterm2=[Bterm2Ahead; Bterm2Atail.*Bterm2Btail];
Bterm_left=Bterm1.*Bterm2;
for ii=0:deg
    Bterm_rightL=OP_eval(eta1,deg-ii,2*ii+1,0);
    % note: term_left(ii+1,:) is the component at degree "ii".
    Bterm_leftL=repmat(Bterm_left(ii+1,:),size(Bterm_rightL,1),1);
    dVy2L=Bterm_leftL.*Bterm_rightL;
    dVy2=[dVy2; dVy2L];
end





% Third term of the sums of dVy.
Cterm1=Bterm1; % Bterm1=OP_eval(eta1,deg,0,0);
Cterm2=Aterm2; % Aterm2=MB_eval((1-eta1)/2,deg);
Cterm_left=Cterm1.*Cterm2;

dVy3=[];
for ii=0:deg
    % the term 2 is due to the chain rule.
    Cterm_rightL=2*diffOP_eval(eta1,deg-ii,2*ii+1,0);
    % note: term_left(ii+1,:) is the component at degree "ii".
    Cterm_leftL=repmat(Cterm_left(ii+1,:),size(Cterm_rightL,1),1);
    dVy3L=Cterm_leftL.*Cterm_rightL;
    dVy3=[dVy3; dVy3L];
end

dVy=dVy1+dVy2+dVy3;






function y=MB_eval(x,deg)

% evaluation at x of the monomial basis at any degree (y is a matrix,
% where the k-th row is the evaluation of the polynomial at degree k-1.


if size(x,2) == 1
    flag=0;
    x=x';
end

y=[];

for ii=0:deg
    y=[y; x.^ii];
end



function y=OP_eval(x,deg,a,b)

% evaluation at x of the Jacobi orthogonal polynomial with exponents "a", "b",
% at any degree (y is a matrix, where the k-th row is the evaluation of the
% polynomial at degree k-1).

if size(x,2) == 1
    flag=0;
    x=x';
end

if deg > 0
    [An,Bn,Cn]=recurr_coeffs(deg,a,b);
    
    L=length(An);
    y=[zeros(size(x)); ones(size(x))];
    for ii=1:L
        yL=(An(ii)*x+Bn(ii)).*y(ii+1,:)-Cn(ii)*y(ii,:);
        y=[y; yL];
    end
    
    y=y(2:end,:);
else
    y=ones(size(x));
end



function [An,Bn,Cn]=recurr_coeffs(deg,a,b)

% Rec. coeffs of Jacobi weight with exponents "a" and "b".
% p[n+1](x)=(A[n]*x+B[n])*p[n](x)-C[n]*p[n-1](x), p[1]=1, p[0]=0.


n=1:(deg-1); n=n';

An=((2*n+a+b+1).*(2*n+a+b+2))./(2*(n+1).*(n+a+b+1));
An=[1+(a+b)/2; An];

Bn=((a^2-b^2).*(2*n+a+b+1))./(2*(n+1).*(n+a+b+1).*(2*n+a+b));
Bn=[(a-b)/2; Bn];

Cn=((n+a).*(n+b).*(2*n+a+b+2))./((n+1).*(n+a+b+1).*(2*n+a+b));
Cn=[0; Cn];




function y=diffOP_eval(pts,deg,a,b)

L=length(pts);
if deg > 0
    Aterm1apre=OP_eval(pts,deg-1,a+1,b+1);
    Aterm1a=[zeros(1,L); Aterm1apre];
    Aterm1b_pre=(1/2)*((0:deg)+a+b+1)';
    Aterm1b=repmat(Aterm1b_pre,1,size(Aterm1a,2));
    y=Aterm1a.*Aterm1b;
else
    y=zeros(1,L);
end