function [K,B,V_ONB,varargout] = ONB(V,z,varargin)
%--------------------------------------------------------------------------
% ONB compute orthonormal basis kernel and Bergman function of a linear 
% space with respect to L^2_w product.
%--------------------------------------------------------------------------
%
% INPUT PARAMETERS
% V         [npts nbasis] real Vandermonde matrix at sampling pts
% z         [npts 1] sq root of weights of discrete product
%
% OPTIONAL INPUT PARAMETER
% Veval     [npts_eval nbasis] real Vandermonde matrix at evaluation pts  
%
%--------------------------------------------------------------------------
%
% OUTPUT PARAMETERS
% K         [npts npts] reproducing kerlen at sampling pts
% B         [npts 1] Bergman function at sampling pts
% V_ONB     [npts nbasis] Vnd. matrix of orthonormal basis at sampling pts
%
% OPTIONAL OUTPUT PARAMETERS
% K_eval         [npts npts] reproducing kerlen at eval pts
% B_eval         [npts 1] Bergman function at eval pts
% V_ONB_eval     [npts nbasis] Vnd. matrix of orthonormal basis at eval pts
%
%--------------------------------------------------------------------------
W = sparse(diag(abs(z)));
[~,R] = qr(V,0);
V1 = V/R;
[~,R1] = qr(W*V1,0);
V_ONB = real(V1/R1);
K=real(V_ONB*V_ONB');
B=diag(K);
if nargin >2
    Veval=varargin{1};
    V_ONB_eval=real((Veval/R)/R1);
    K_eval=V_ONB_eval*V_ONB_eval';
    B_eval=diag(K_eval);
    if nargout >3
        varargout{1}=V_ONB_eval;
        varargout{2}=K_eval;
        varargout{3}=B_eval;
    end
end