

function moments=compute_moments(XYZW,ade,bbox,basis_type)

%--------------------------------------------------------------------------
% Object:
% Moments of the rule with nodes XYZW(:,1:3) and weights XYZW(:,4),
% relative to a certain polynomial basis of degree "ade" on hyperrectangle
% defined by "bbox".
% The basis is by default the shifted Chebyshev tensorial type on "bbox",
% of total degree. Setting "basis_type=1" one may use the monomial basis.
%--------------------------------------------------------------------------
% Input:
% ade: polynomial degree.
%
% bbox: 6-component vector defining the hyperrectangle:
% [bbox(1),bbox(2)] x [bbox(3),bbox(4)] x [bbox(5),bbox(6)]
%
% basis_type: 0: shifted Chebyshev tensorial type, 1: monomial bais
%--------------------------------------------------------------------------
% Output:
% moments: approximation of the moments of the Vandermonde matrix
% of total degree "ade" (relatively to the choosen basis).
%--------------------------------------------------------------------------
% External routines:
% 1. chebvand (and its subroutines)
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello
% University of Padova, September 7, 2020
%--------------------------------------------------------------------------

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

if nargin < 4
    basis_type=0;
end

X=XYZW(:,1); Y=XYZW(:,2); Z=XYZW(:,3); W=XYZW(:,4);

if (basis_type == 0) & nargin < 3
    bbox=[min(X) max(X) min(Y) max(Y) min(Z) max(Z)];
end

% .......................... moments computation ..........................

nodes=XYZW(:,1:3);

switch basis_type
    case 0
        
        V=chebvand(nodes,ade,bbox);
        moments=(XYZW(:,4))'*V;
        moments=moments'; % column vector
        
    case 1
        
        moments=[];
        for ii=0:ade
            for jj=0:(ade-ii)
                for kk=0: ade-jj-ii % computing contribution
                    FXYZC=(X.^ii).*(Y.^jj).*(Z.^kk);
                    moments_L=W'*FXYZC;
                    moments=[moments; moments_L];
                end
            end
        end
        
    otherwise
        
        V=chebvand(nodes,ade,bbox);
        moments=(XYZW(:,4))'*V;
        moments=moments'; % column vector
        
end









%--------------------------------------------------------------------------
% chebvand
%--------------------------------------------------------------------------

function V = chebvand(points,ade,bbox)

%--------------------------------------------------------------------------
% Object:
% Vandermonde matrix at "points" points, relative to Chebyshev basis of 
% degree "ade" on rectangle "bbox". The basis is of tensorial type, 
% of total degree.
%--------------------------------------------------------------------------
% Input:
% ade: polynomial degree.
% points: 3-column array of point coordinates.
% bbox: 6-component vector such that the cube:
% [bbox(1),bbox(2)] x [bbox(3),bbox(4)] x [bbox(5),bbox(6)] contains X.
%--------------------------------------------------------------------------
% Output:
% V: Chebyshev-Vandermonde matrix at "points", graded lexic. order.
%    size is M x N, where M is the number of nodes and N is
%    (ade+1)*(ade+2)*(ade+3)/6.
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello
% University of Padova, September 1, 2020
%--------------------------------------------------------------------------

if nargin < 3 
    bbox=[];
end

if isempty(bbox)
    bbox=[min(points(:,1)) max(points(:,1)) ...
        min(points(:,2)) max(points(:,2)) ...
        min(points(:,3)) max(points(:,3))];
end

degrees=[];
for ii=0:ade
    for jj=0:(ade-ii)
        for kk=0: ade-jj-ii % computing contribution
            degrees=[degrees; ii jj kk];
        end
    end
end

% mapping the mesh in the cube [-1,1]^3
map=[(2*points(:,1)-bbox(2)-bbox(1))/(bbox(2)-bbox(1)) ...
    (2*points(:,2)-bbox(4)-bbox(3))/(bbox(4)-bbox(3)) ...
    (2*points(:,3)-bbox(6)-bbox(5))/(bbox(6)-bbox(5))];

% Chebyshev-Vandermonde matrix on the mesh
T1=chebpolys(ade,map(:,1));
T2=chebpolys(ade,map(:,2));
T3=chebpolys(ade,map(:,3));
V=T1(:,degrees(:,1)+1).*T2(:,degrees(:,2)+1).*T3(:,degrees(:,3)+1);








%--------------------------------------------------------------------------
% chebpolys
%--------------------------------------------------------------------------

function T=chebpolys(deg,x)

%--------------------------------------------------------------------------
% Object:
% This routine computes the Chebyshev-Vandermonde matrix on the real line
% by recurrence.
%--------------------------------------------------------------------------
% Input:
% deg: maximum polynomial degree
% x: 1-column array of abscissas
%--------------------------------------------------------------------------
% Output:
% T: Chebyshev-Vandermonde matrix at x, T(i,j+1)=T_j(x_i), j=0,...,deg.
%    Its size is M x N where M is the number of nodes and N the dimension 
%    of the polynomial space (i.e. "deg+1").
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello
% University of Padova, September 1, 2020.
%--------------------------------------------------------------------------

T=zeros(length(x),deg+1);
t0=ones(length(x),1);
T(:,1)=t0;
t1=x;
T(:,2)=t1;

for j=2:deg
    t2=2*x.*t1-t0;
    T(:,j+1)=t2;
    t0=t1;
    t1=t2;
end
