
function [pts_sq,w_sq]=tens_cub_sq(deg,a,b,gl)

%--------------------------------------------------------------------------
% TENS. GAUSS-JACOBI (LOBATTO) RULE ON [-1,1] x [-1,1].
%--------------------------------------------------------------------------
% INPUT.
%--------
% deg IS THE DEGREE OF THE RULE.
% a,b ARE THE GAUSS-JACOBI EXPONENTS.
% gl: 0: GAUSS POINTS. 1: GAUSS-LOBATTO POINTS.
%--------------------------------------------------------------------------
% OUTPUT.
%--------
% x_sq, w_sq ARE COLUMN VECTORS OF NODES AND WEIGHTS.
%--------------------------------------------------------------------------
% REQUIRED ROUTINES.
%---------------------
%
% 1. gauss_jacobi
%
%--------------------------------------------------------------------------

if nargin == 1
    a=0; b=0; gl=0;
end

if gl == 0
    N=floor(deg/2)+1; % GAUSS LEGENDRE.
else
    N=floor(deg/2); % GAUSS LEGENDRE LOBATTO.
end

if nargin < 2
    a=0; b=0;
end

if nargin < 4
    gl = 0;
end

[nodes,weights]=gauss_jacobi(N,a,b,gl);

[X,Y]=meshgrid(nodes);

pts_sq=[X(:) Y(:)];

[Wx,Wy]=meshgrid(weights);
WW=Wx.*Wy;
w_sq=WW(:);