function [V,dbox] = dLEGVAND_orthn(deg,X,dbox,basis_indices)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% This routine computes the Legendre-Vandermonde matrix for degree "deg"
% on a d-dimensional point cloud "X".
%
% Important: The Legendre basis is orthonormal when dbox is the
% unit-square, cube or more generally hypercube [-1,1]^d.
%
% The Legendre basis is the tensorial Legendre basis of total degree
% "deg", shifted on the hyperrectangle defined by "dbox".
%
% If "dbox" is not provided, the routine sets that variable to define the
% smaller "hyperrectangle" (box) with sides parallel to the cartesian
% axes and containing the pointset "X".
%
% Note: We checked the orthonormality of the rule in [-1,1]^2.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% deg: polynomial degree;
% X: d-column array of "m" points cloud (matrix "m x d");
% * dbox: variable that defines the smallest hyperectangle with sides
%     parallel to the axis, containing the domain.
%     If "dbox" is not provided, it defines 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 instance, for a 2-sphere, it is "d=3", while 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]".
%
% * basis_indices: nchoosek(n+d,n) x d matrix of indices describing the
%     tensorial basis. Example: if "d=2" then the i-th row is a couple
%     (h,k) and the basis element will be "T_h(x) T_k(y)" with "T_h" 
%     univariate orthonormal Legendre polynomial of the first kind.
%
% Note: the variables with an asterisk "*" are not mandatory and can be
% also set as empty matrix.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% V: shifted Legendre-Vandermonde matrix for degree "deg" on
%    the pointset "X", relatively to "dbox".
% dbox: variable that defines the hyperrectangle with sides parallel to the
%     axis, containing the domain.
%--------------------------------------------------------------------------
% Dates:
%--------------------------------------------------------------------------
% Written by A. Sommariva, M. Vianello
% Modified by A. Sommariva. December 5, 2025.
%--------------------------------------------------------------------------


% ........................... Function body ...............................

% ...... troubleshooting ......

% box containing the cloud
if nargin < 3, dbox=[]; end
if isempty(dbox)
    a=min(X); b=max(X); dbox=[a;b];
else
    a=dbox(1,:); b=dbox(2,:);
end

% d-uples of indices with sum less or equal to "deg" graded lexicographical
% order
if nargin < 4
    d=size(X,2); N = nchoosek(deg+d,d); basis_indices = zeros(N,d);
    for i=2:N
        basis_indices(i,:) = mono_next_grlex(d,basis_indices(i-1,:));
    end
end


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

d=size(X,2);

% mapping the mesh in the hypercube "[-1,1]^d"
map = zeros(size(X));
for i=1:d
    map(:,i)=(2*X(:,i)-b(i)-a(i))/(b(i)-a(i));
end

% orthonormal Legendre-Vandermonde matrix on the mesh
T=legpolys_orthn(deg,map(:,1));

V=T(:,basis_indices(:,1)+1);
for i=2:d
    T=legpolys_orthn(deg,map(:,i));
    V=V.*T(:,basis_indices(:,i)+1);
end










function Torthn=legpolys_orthn(deg,x)

%--------------------------------------------------------------------------
% Object:
% This routine computes the Legendre-Vandermonde matrix on the real line
% by recurrence.
% The polynomial is orthonormal w.r.t. Legendre measure.
%--------------------------------------------------------------------------
% Input:
% deg: maximum polynomial degree
% x: 1-column array of abscissas
%--------------------------------------------------------------------------
% Output:
% Torthn: orthonormal Legendre-Vandermonde matrix at "x".
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello
% University of Padova, December 05, 2025
%--------------------------------------------------------------------------

% Determining classical Legendre polynomials.

T=ones(length(x),deg+1); % Initialize matrix. 

% degree 1
T(:,2)=x; 

% higher degrees
for k=2:deg
    n=k-1;
    T(:,k+1)=((2*n+1)/(n+1)).*x.*T(:,k)-(n/(n+1))*T(:,k-1);
end



% Determining orthonormal Legendre polynomials.
Torthn=zeros(length(x),deg+1);
for k=1:deg+1
    n=k-1;
    term=sqrt(2/(2*n+1));
    Torthn(:,k)=T(:,k)/term;
end



