function [Q,R1,R2] = wamdop(deg,wam,rect);

% by Alvise and Marco Vianello, University of Padova
% October 30, 2014

% computes discrete orthogonal polynomials up to degree deg 
% on the mesh wam: Q is the corresponding orthogonal matrix, 
% inv(R1)*inv(R2) the orthogonalization matrix 

if isempty(rect)
rect=[min(wam(:,1)) max(wam(:,1)) min(wam(:,2)) max(wam(:,2))];
end

V=chebvand(deg,wam,rect);
[Q1,R1]=qr(V,0);
[Q,R2]=qr(Q1,0);

end


