function [V,couples,N] = chebvandrp(deg,p,gmesh,rect)
%--------------------------------------------------------------------------
% computes the bivariate Chebyshev-Vandermonde matrix (with the Chebyshev
% basis defined on the rectangle rect) at the target points gmesh
% INPUT:
% deg = polynomial degree
% gmesh = 2-column array of mesh point coordinates
% rect = 4-component vector such that the rectangle
% [rect(1),rect(2)] x [rect(3),rect(4)] contains the mesh
% -------------------------------------------------------------------------
% OUTPUT:
% V = Chebyshev-Vandermonde matrix at gmesh
%--------------------------------------------------------------------------
% FUNCTION BODY
% rectangle containing the mesh
if isempty(rect)
    rect=[min(gmesh(:,1)) max(gmesh(:,1)) min(gmesh(:,2)) max(gmesh(:,2))];
end;
% couples with length less or equal to deg
% graded lexicographical order
j=(0:1:deg);
[j1,j2]=meshgrid(j);
if p<inf
    good=find((j1(:).^p+j2(:).^p).^(1/p)<=deg);
    couples=[j1(good) j2(good)];
elseif p==inf
    good=(find(max([j1(:)';j2(:)'])<=deg))';
    couples=[j1(good) j2(good)];
end
    
% mapping the mesh in the square [-1,1]^2
a=rect(1);b=rect(2);c=rect(3);d=rect(4);
map=[(2*gmesh(:,1)-b-a)/(b-a) (2*gmesh(:,2)-d-c)/(d-c)];
% Chebyshev-Vandermonde matrix on the mesh
for j=1:length(couples(:,1))
    V(:,j)=chebpoly(couples(j,1),map(:,1)).*chebpoly(couples(j,2),map(:,2));
end
N=size(couples,1);
end

function tn = chebpoly(deg,x)
% -------------------------------------------------------------------------
% computes the Chebyshev polynomials on the real line by recurrence
t0=ones(length(x),1);
t1=x;
if deg==0
    tn=t0;
end;
if deg>0
    for j=2:deg
        t2=2*x.*t1-t0;
        t0=t1;
        t1=t2;
    end;
    tn=t1;
end;
end

