
function rule=cub_square_mpx(ade)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Morrow-Patterson-Xu quasi minimal rule with degree of exactness "ade" on 
% the unit-square [-1,1] w.r.t. the weight function
%
%             w(x,y)=(1/pi^2)*(1-x^2)^(-1/2)*(1-y^2)^(-1/2).
%
% The procedure gives nodes and weights stored in the N x 3 matrix "rule",
% so that the nodes are "(rule(:,1),rule(:,2))" and the weights "rule(:,3)".
% Next
% * if "n" is odd then "N=(n+1)*(n+3)/2";
% * if "n" is even then "N=(n+2)^2/2".
%--------------------------------------------------------------------------
% REFERENCE
%--------------------------------------------------------------------------
% [1] Hyper2d: a numerical code for hyperinterpolation on rectangles
%     M. Caliari, S. De Marchi, R. Montagna, M. Vianello
%     Applied Mathematics and Computation
%     Volume 183, Issue 2, 15 December 2006, Pages 1138-1147
%     https://www.math.unipd.it/~marcov/pdf/hyper2d.pdf
%--------------------------------------------------------------------------
% NOTE
%--------------------------------------------------------------------------
% The code on the paper [1] is not correct (missing variable "z").
% The present routine has been tested and it works correctly. In particular
% it has been compared with a tensorial gaussian rule w.r.t. to the weight
% function "w".
%--------------------------------------------------------------------------
% EXAMPLES
%--------------------------------------------------------------------------
% >> rule=cub_square_mp(5)
% 
% rule =
% 
%     1.0000    1.0000    0.0556
%     1.0000   -0.5000    0.1111
%    -0.5000    1.0000    0.1111
%    -0.5000   -0.5000    0.2222
%     0.5000    0.5000    0.2222
%     0.5000   -1.0000    0.1111
%    -1.0000    0.5000    0.1111
%    -1.0000   -1.0000    0.0556
% 
% >> rule=cub_square_mp(7)
% 
% rule =
% 
%     1.0000    0.7071    0.0625
%     1.0000   -0.7071    0.0625
%     0.0000    0.7071    0.1250
%     0.0000   -0.7071    0.1250
%    -1.0000    0.7071    0.0625
%    -1.0000   -0.7071    0.0625
%     0.7071    1.0000    0.0625
%     0.7071    0.0000    0.1250
%     0.7071   -1.0000    0.0625
%    -0.7071    1.0000    0.0625
%    -0.7071    0.0000    0.1250
%    -0.7071   -1.0000    0.0625
%
% >>
%--------------------------------------------------------------------------
% AUTHORS OF THE PRESENT CODE
%--------------------------------------------------------------------------
% L. Rinaldi, A. Sommariva, M. Vianello.
%--------------------------------------------------------------------------
% TESTS
%--------------------------------------------------------------------------
% Tested on Matlab R2024B, on a PC running Intel(R) N150 (800 MHz) with 16
% GB of RAM.
%--------------------------------------------------------------------------
% LICENSE
%--------------------------------------------------------------------------
% Copyright (C) 2025 Laura Rinaldi, Alvise Sommariva, Marco Vianello.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%
% Authors:
%
% Laura Rinaldi    <laura.rinaldi@unipd.it>
% Alvise Sommariva <alvise@math.unipd.it>
% Marco Vianello   <marcov@math.unipd.it>
%
% Date: November 23, 2025
%--------------------------------------------------------------------------

n=floor(ade/2); % ade <= 2n+1 (see page 3 of [1])

k=0:n+1;
z=cos(k'*pi/(n+1));

if rem(n,2) == 0

    % Building A_even
    zE=z(1:2:end);
    [AEx,AEy]=meshgrid(zE);

    % Building B_even
    zO=z(2:2:end);
    [BEx,BEy]=meshgrid(zO);

    % Defining weights
    inner_val=2*(n+1)^(-2);
    boundary_val=(n+1)^(-2);
    corner_val=(n+1)^(-2)/2;

    WAE=inner_val*ones(size(AEx));
    WAE(:,1)=boundary_val*ones(size(AEx,1),1);
    WAE(1,:)=boundary_val*ones(1,size(AEx,2));
    WAE(1,1)=corner_val;

    WBE=inner_val*ones(size(BEx));
    WBE(:,end)=boundary_val*ones(size(BEx,1),1);
    WBE(end,:)=boundary_val*ones(1,size(BEx,2));
    WBE(end,end)=corner_val;

    X1=AEx(:); X2=BEx(:); X=[X1; X2];
    Y1=AEy(:); Y2=BEy(:); Y=[Y1; Y2];
    W1=WAE(:); W2=WBE(:); W=[W1; W2];

    % Assembling results.
    rule=[X Y W];

else

    % Building A_odd
    zE=z(1:2:end);
    zO=z(2:2:end);
    [AEx,AEy]=meshgrid(zE,zO);

    % Building B_odd
    [BEx,BEy]=meshgrid(zO,zE);

    % Defining weights
    inner_val=2*(n+1)^(-2);
    boundary_val=(n+1)^(-2);

    WAE=inner_val*ones(size(AEx));
    WAE(:,1)=boundary_val*ones(size(AEx,1),1);
    WAE(:,end)=boundary_val*ones(size(AEx,1),1);

    WBE=inner_val*ones(size(BEx));
    WBE(1,:)=boundary_val*ones(1,size(BEx,2));
    WBE(end,:)=boundary_val*ones(1,size(BEx,2));

    X1=AEx(:); X2=BEx(:); X=[X1; X2];
    Y1=AEy(:); Y2=BEy(:); Y=[Y1; Y2];
    W1=WAE(:); W2=WBE(:); W=[W1; W2];

    % Assembling results.
    rule=[X Y W];

end

