function [V_OP,BK,varargout] = DiscreteOP2(deg,nodes,weights,varargin)
rect=[min(nodes(1,:)),max(nodes(1,:)),min(nodes(2,:)),max(nodes(2,:))];
V=chebvand(deg,nodes',rect);
    W = sparse(diag(sqrt(weights)));
% orthogonalization process
    [Q,R] = qr(V,0);
    V1 = V/R;
    [Q,R1] = qr(W*V1,0);
% Vandermonde in the OP basis on the cubature nodes
    V_OP = V1/R1;
    %norm(V1'*V1-eye(size(V1,2)))
    BK=sum(V_OP.^2,2);
if nargin <4
    1;
else
    eval_nodes=varargin{1};
    recteval=[min(eval_nodes(1,:)),max(eval_nodes(1,:)),...
        min(eval_nodes(2,:)),max(eval_nodes(2,:))];
    Veval=chebvand(deg,eval_nodes',recteval);
    V_OP_eval=(Veval/R)/R1;
    BK_eval=sum(V_OP_eval.^2,2);
    if nargout >2
        varargout{1}=R;
        varargout{2}=R1;
        varargout{3}=V_OP_eval;
        varargout{4}=BK_eval;
    end;
end
