function leb = dLEBsph(deg,nodes,w,test_pts,jvec,dbox,dimpoly)

%--------------------------------------------------------------------------
% Object:
% This routine computes the Lebesgue constant on "nodes" of d-variate
% weighted least-squares polynomial fitting at "X".
%--------------------------------------------------------------------------
% Input:
% deg: polynomial degree;
% nodes: d-column array of points on which the Lebesgue constant is
%    required, at degree "deg";
% w: 1-column array of nonnegative weights or nonnegative scalar in
%    case of equal weights;
% test_pts: d-column array of control point coordinates, useful to estimate
%    the Lebesgue constant;
% * jvec: vector of column indexes, selects a polynomial basis;
% * dbox: variable that defines a hyperrectangle with sides parallel to the
%    axis, containing the domain.
%    If "dbox" is not provided, it is the smaller "hyperrectangle", with
%    sides parallel to the cartesian axes, containing the pointset "X".
%    It is a matrix with dimension "2 x d", where "d" is the dimension of
%    the space in which it is embedded the domain.
%    For example, for a 2-sphere, it is "d=3", for a 2 dimensional polygon
%    it is "d=2".
%    As example, the set "[-1,1] x [0,1]" is described as "[-1 0; 1 1]".
% * domain_structure: structure defining the domain,
%    domain_struct.domain: string with domain name
%    domain_struct.parms: domain parameters.
%
% Note:
% The variables with an asterisk "*" are not mandatory and can be also set
% as empty matrix.
%--------------------------------------------------------------------------
% Output:
% leb: Lebesgue constant estimate on "nodes", based on evaluations on "X".
%--------------------------------------------------------------------------
% Dates:
% Written on 26/07/2020 by M. Dessole, F. Marcuzzi, M. Vianello.
%
% Modified by:
% 29/10/2020: M. Vianello;
% 03/11/2020: M. Dessole, M. Vianello;
% 28/11/2020: A. Sommariva.
%--------------------------------------------------------------------------


% .........................  Function Body ................................

% ..... troubleshooting .....
if nargin < 8, dimpoly=[]; end
if nargin < 7, domain_structure.domain='generic'; end
if nargin < 6, dbox=[]; end
if nargin < 5, jvec=[]; end
if nargin < 4, test_pts=[]; end
if nargin < 3, w=[]; end
if isempty(test_pts), test_pts=nodes; end
if isempty(w), w=1; end
if isempty(domain_structure), domain_structure.domain='generic'; end

domain=domain_structure.domain;


% ........................ main code below ................................

% ..... Vandermonde matrix .....


if isempty(dbox), dbox = boxdef([test_pts; nodes]); end

[UX,~] = dCHEBVAND(deg,test_pts,dbox);
[VY,jvec0,~,R,~] = dORTHVANDsph(deg,nodes,w,jvec,[],dbox);
if isempty(jvec), jvec=jvec0; end
VX = UX(:,jvec)/R;

% ..... Lebesgue constant approximation .....

VYVXt = VY*VX';

if isscalar(w)
    V = w*VYVXt;
else
    V = zeros(size(VYVXt));
    for k=1:size(V,2), V(:,k)=VYVXt(:,k).*w; end
end

leb = norm(V,1);





