function V = chebvand2(deg,gmesh,rect)

% computes the Chebyshev-Vandermonde matrix on a 2d or 3d given mesh 
% in the total-degree product Chebyshev basis of a rectangle containing 
% the mesh points

% INPUT:

% deg = polynomial degree

% gmesh = 2- or 3-column array of mesh point coordinates

% rect = 4- or 6-component vector such that the rectangle 
% [rect(1),rect(2)] x [rect(3),rect(4)] in 2d
% or [rect(1),rect(2)] x [rect(3),rect(4)] x [rect(5),rect(6)] in 3d
% contains the mesh


% OUTPUT:

% V = Chebyshev-Vandermonde matrix 
% A = Change of basis for monomials

% FUNCTION BODY

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

if length(gmesh(1,:)) == 2 

% couples with length less or equal to deg
j=(0:1:deg);
[j1,j2]=meshgrid(j);
good=find(j1(:)+j2(:)<deg+1);
couples=[j1(good) j2(good)];
dimsp=size(couples,1);


% 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
V1=cos(couples(:,1)*acos(map(:,1)'));
V2=cos(couples(:,2)*acos(map(:,2)'));
V=(V1.*V2)';
end;



if length(gmesh(1,:)) == 3 
    
%triples with length less or equal to deg
j=(0:1:deg);
[j1,j2,j3]=meshgrid(j,j,j);
good=find(j1(:)+j2(:)+j3(:)<deg+1);
triples=[j1(good) j2(good) j3(good)];

%mapping the mesh in the cube [-1,1]^3
a=rect(1);b=rect(2);
c=rect(3);d=rect(4);
e=rect(5);f=rect(6);
map=[(2*gmesh(:,1)-b-a)/(b-a) (2*gmesh(:,2)-d-c)/(d-c) ... 
(2*gmesh(:,3)-f-e)/(f-e)];

%Chebyshev-Vandermonde matrix on the mesh
V1=cos(triples(:,1)*acos(map(:,1)'));
V2=cos(triples(:,2)*acos(map(:,2)'));
V3=cos(triples(:,3)*acos(map(:,3)'));
V=((V1.*V2).*V3)';

end; 
end